moon-abe cleanup

Change-Id: Ie1259856db03f0b9e80de3e967ec6bd1f03191b3

1 files changed, 0 insertions, 1724 deletions

diff --git a/moon-abe/pbc-0.5.14/arith/poly.c b/moon-abe/pbc-0.5.14/arith/poly.c
deleted file mode 100644
index bd2dad33..00000000
--- a/moon-abe/pbc-0.5.14/arith/poly.c
+++ /dev/null
@@ -1,1724 +0,0 @@
-#include <ctype.h>
-#include <stdarg.h>
-#include <stdio.h>
-#include <stdint.h> // for intptr_t
-#include <stdlib.h>
-#include <gmp.h>
-#include "pbc_utils.h"
-#include "pbc_field.h"
-#include "pbc_multiz.h"
-#include "pbc_poly.h"
-#include "pbc_memory.h"
-#include "misc/darray.h"
-
-// == Polynomial rings ==
-//
-// Per-field data:
-typedef struct {
-  field_ptr field;  // Ring where coefficients live.
-  fieldmap mapbase; // Map element from underlying field to constant term.
-} *pfptr;
-
-// Per-element data:
-//TODO: Would we ever need any field besides coeff?
-typedef struct {
-  // The coefficients are held in a darray which is resized as needed.
-  // The last array entry represents the leading coefficient and should be
-  // nonzero. An empty darray represents 0.
-  darray_t coeff;
-} *peptr;
-
-// == Polynomial modulo rings ==
-//
-// Per-field data:
-typedef struct {
-  field_ptr field;   // Base field.
-  fieldmap mapbase;  // Similar to mapbase above.
-  int n;             // Degree of extension.
-  element_t poly;    // Polynomial of degree n.
-  element_t *xpwr;   // x^n,...,x^{2n-2} mod poly
-} *mfptr;
-// Per-element data: just a pointer to an array of element_t. This array always
-// has size n.
-
-// Add or remove coefficients until there are exactly n of them. Any new
-// coefficients are initialized to zero, which violates the invariant that the
-// leading coefficient must be nonzero. Thus routines calling this function
-// must check for this and fix the polynomial if necessary, e.g. by calling
-// poly_remove_leading_zeroes().
-static void poly_alloc(element_ptr e, int n) {
-  pfptr pdp = e->field->data;
-  peptr p = e->data;
-  element_ptr e0;
-  int k = p->coeff->count;
-  while (k < n) {
-    e0 = pbc_malloc(sizeof(element_t));
-    element_init(e0, pdp->field);
-    darray_append(p->coeff, e0);
-    k++;
-  }
-  while (k > n) {
-    k--;
-    e0 = darray_at(p->coeff, k);
-    element_clear(e0);
-    pbc_free(e0);
-    darray_remove_last(p->coeff);
-  }
-}
-
-static void poly_init(element_ptr e) {
-  peptr p = e->data = pbc_malloc(sizeof(*p));
-  darray_init(p->coeff);
-}
-
-static void poly_clear(element_ptr e) {
-  peptr p = e->data;
-
-  poly_alloc(e, 0);
-  darray_clear(p->coeff);
-  pbc_free(e->data);
-}
-
-// Some operations may zero a leading coefficient, which will cause other
-// routines to fail. After such an operation, this function should be called,
-// as it strips all leading zero coefficients and frees the memory they
-// occupied, reestablishing the guarantee that the last element of the array
-// is nonzero.
-static void poly_remove_leading_zeroes(element_ptr e) {
-  peptr p = e->data;
-  int n = p->coeff->count - 1;
-  while (n >= 0) {
-    element_ptr e0 = p->coeff->item[n];
-    if (!element_is0(e0)) return;
-    element_clear(e0);
-    pbc_free(e0);
-    darray_remove_last(p->coeff);
-    n--;
-  }
-}
-
-static void poly_set0(element_ptr e) {
-  poly_alloc(e, 0);
-}
-
-static void poly_set1(element_ptr e) {
-  peptr p = e->data;
-  element_ptr e0;
-
-  poly_alloc(e, 1);
-  e0 = p->coeff->item[0];
-  element_set1(e0);
-}
-
-static int poly_is0(element_ptr e) {
-  peptr p = e->data;
-  return !p->coeff->count;
-}
-
-static int poly_is1(element_ptr e) {
-  peptr p = e->data;
-  if (p->coeff->count == 1) {
-    return element_is1(p->coeff->item[0]);
-  }
-  return 0;
-}
-
-static void poly_set_si(element_ptr e, signed long int op) {
-  peptr p = e->data;
-  element_ptr e0;
-
-  poly_alloc(e, 1);
-  e0 = p->coeff->item[0];
-  element_set_si(e0, op);
-  poly_remove_leading_zeroes(e);
-}
-
-static void poly_set_mpz(element_ptr e, mpz_ptr op) {
-  peptr p = e->data;
-
-  poly_alloc(e, 1);
-  element_set_mpz(p->coeff->item[0], op);
-  poly_remove_leading_zeroes(e);
-}
-
-static void poly_set_multiz(element_ptr e, multiz op) {
-  if (multiz_is_z(op)) {
-    // TODO: Remove unnecessary copy.
-    mpz_t z;
-    mpz_init(z);
-    multiz_to_mpz(z, op);
-    poly_set_mpz(e, z);
-    mpz_clear(z);
-    return;
-  }
-  peptr p = e->data;
-  int n = multiz_count(op);
-  poly_alloc(e, n);
-  int i;
-  for(i = 0; i < n; i++) {
-    element_set_multiz(p->coeff->item[i], multiz_at(op, i));
-  }
-  poly_remove_leading_zeroes(e);
-}
-
-static void poly_set(element_ptr dst, element_ptr src) {
-  peptr psrc = src->data;
-  peptr pdst = dst->data;
-  int i;
-
-  poly_alloc(dst, psrc->coeff->count);
-  for (i=0; i<psrc->coeff->count; i++) {
-    element_set(pdst->coeff->item[i], psrc->coeff->item[i]);
-  }
-}
-
-static int poly_coeff_count(element_ptr e) {
-  return ((peptr) e->data)->coeff->count;
-}
-
-static element_ptr poly_coeff(element_ptr e, int n) {
-  peptr ep = e->data;
-  PBC_ASSERT(n < poly_coeff_count(e), "coefficient out of range");
-  return (element_ptr) ep->coeff->item[n];
-}
-
-static int poly_sgn(element_ptr f) {
-  int res = 0;
-  int i;
-  int n = poly_coeff_count(f);
-  for (i=0; i<n; i++) {
-    res = element_sgn(poly_coeff(f, i));
-    if (res) break;
-  }
-  return res;
-}
-
-static void poly_add(element_ptr sum, element_ptr f, element_ptr g) {
-  int i, n, n1;
-  element_ptr big;
-
-  n = poly_coeff_count(f);
-  n1 = poly_coeff_count(g);
-  if (n > n1) {
-    big = f;
-    n = n1;
-    n1 = poly_coeff_count(f);
-  } else {
-    big = g;
-  }
-
-  poly_alloc(sum, n1);
-  for (i=0; i<n; i++) {
-    element_add(poly_coeff(sum, i), poly_coeff(f, i), poly_coeff(g, i));
-  }
-  for (; i<n1; i++) {
-    element_set(poly_coeff(sum, i), poly_coeff(big, i));
-  }
-  poly_remove_leading_zeroes(sum);
-}
-
-static void poly_sub(element_ptr diff, element_ptr f, element_ptr g) {
-  int i, n, n1;
-  element_ptr big;
-
-  n = poly_coeff_count(f);
-  n1 = poly_coeff_count(g);
-  if (n > n1) {
-    big = f;
-    n = n1;
-    n1 = poly_coeff_count(f);
-  } else {
-    big = g;
-  }
-
-  poly_alloc(diff, n1);
-  for (i=0; i<n; i++) {
-    element_sub(poly_coeff(diff, i), poly_coeff(f, i), poly_coeff(g, i));
-  }
-  for (; i<n1; i++) {
-    if (big == f) {
-      element_set(poly_coeff(diff, i), poly_coeff(big, i));
-    } else {
-      element_neg(poly_coeff(diff, i), poly_coeff(big, i));
-    }
-  }
-  poly_remove_leading_zeroes(diff);
-}
-
-static void poly_neg(element_ptr f, element_ptr g) {
-  peptr pf = f->data;
-  peptr pg = g->data;
-  int i, n;
-
-  n = pg->coeff->count;
-  poly_alloc(f, n);
-  for (i=0; i<n; i++) {
-    element_neg(pf->coeff->item[i], pg->coeff->item[i]);
-  }
-}
-
-static void poly_double(element_ptr f, element_ptr g) {
-  peptr pf = f->data;
-  peptr pg = g->data;
-  int i, n;
-
-  n = pg->coeff->count;
-  poly_alloc(f, n);
-  for (i=0; i<n; i++) {
-    element_double(pf->coeff->item[i], pg->coeff->item[i]);
-  }
-}
-
-static void poly_mul_mpz(element_ptr f, element_ptr g, mpz_ptr z) {
-  peptr pf = f->data;
-  peptr pg = g->data;
-  int i, n;
-
-  n = pg->coeff->count;
-  poly_alloc(f, n);
-  for (i=0; i<n; i++) {
-    element_mul_mpz(pf->coeff->item[i], pg->coeff->item[i], z);
-  }
-}
-
-static void poly_mul_si(element_ptr f, element_ptr g, signed long int z) {
-  peptr pf = f->data;
-  peptr pg = g->data;
-  int i, n;
-
-  n = pg->coeff->count;
-  poly_alloc(f, n);
-  for (i=0; i<n; i++) {
-    element_mul_si(pf->coeff->item[i], pg->coeff->item[i], z);
-  }
-}
-
-static void poly_mul(element_ptr r, element_ptr f, element_ptr g) {
-  peptr pprod;
-  peptr pf = f->data;
-  peptr pg = g->data;
-  pfptr pdp = r->field->data;
-  int fcount = pf->coeff->count;
-  int gcount = pg->coeff->count;
-  int i, j, n;
-  element_t prod;
-  element_t e0;
-
-  if (!fcount || !gcount) {
-    element_set0(r);
-    return;
-  }
-  element_init(prod, r->field);
-  pprod = prod->data;
-  n = fcount + gcount - 1;
-  poly_alloc(prod, n);
-  element_init(e0, pdp->field);
-  for (i=0; i<n; i++) {
-    element_ptr x = pprod->coeff->item[i];
-    element_set0(x);
-    for (j=0; j<=i; j++) {
-      if (j < fcount && i - j < gcount) {
-        element_mul(e0, pf->coeff->item[j], pg->coeff->item[i - j]);
-        element_add(x, x, e0);
-      }
-    }
-  }
-  poly_remove_leading_zeroes(prod);
-  element_set(r, prod);
-  element_clear(e0);
-  element_clear(prod);
-}
-
-static void polymod_random(element_ptr e) {
-  element_t *coeff = e->data;
-  int i, n = polymod_field_degree(e->field);
-
-  for (i=0; i<n; i++) {
-    element_random(coeff[i]);
-  }
-}
-
-static void polymod_from_hash(element_ptr e, void *data, int len) {
-  // TODO: Improve this.
-  element_t *coeff = e->data;
-  int i, n = polymod_field_degree(e->field);
-  for (i=0; i<n; i++) {
-    element_from_hash(coeff[i], data, len);
-  }
-}
-
-static size_t poly_out_str(FILE *stream, int base, element_ptr e) {
-  int i;
-  int n = poly_coeff_count(e);
-  size_t result = 2, status;
-
-  /*
-  if (!n) {
-    if (EOF == fputs("[0]", stream)) return 0;
-    return 3;
-  }
-  */
-  if (EOF == fputc('[', stream)) return 0;
-  for (i=0; i<n; i++) {
-    if (i) {
-      if (EOF == fputs(", ", stream)) return 0;
-      result += 2;
-    }
-    status = element_out_str(stream, base, poly_coeff(e, i));
-    if (!status) return 0;
-    result += status;
-  }
-  if (EOF == fputc(']', stream)) return 0;
-  return result;
-}
-
-static int poly_snprint(char *s, size_t size, element_ptr e) {
-  int i;
-  int n = poly_coeff_count(e);
-  size_t result = 0, left;
-  int status;
-
-  #define clip_sub() {                         \
-    result += status;                          \
-    left = result >= size ? 0 : size - result; \
-  }
-
-  status = snprintf(s, size, "[");
-  if (status < 0) return status;
-  clip_sub();
-
-  for (i=0; i<n; i++) {
-    if (i) {
-      status = snprintf(s + result, left, ", ");
-      if (status < 0) return status;
-      clip_sub();
-    }
-    status = element_snprint(s + result, left, poly_coeff(e, i));
-    if (status < 0) return status;
-    clip_sub();
-  }
-  status = snprintf(s + result, left, "]");
-  if (status < 0) return status;
-  return result + status;
-  #undef clip_sub
-}
-
-static void poly_div(element_ptr quot, element_ptr rem,
-    element_ptr a, element_ptr b) {
-  peptr pq, pr;
-  pfptr pdp = a->field->data;
-  element_t q, r;
-  element_t binv, e0;
-  element_ptr qe;
-  int m, n;
-  int i, k;
-
-  if (element_is0(b)) pbc_die("division by zero");
-  n = poly_degree(b);
-  m = poly_degree(a);
-  if (n > m) {
-    element_set(rem, a);
-    element_set0(quot);
-    return;
-  }
-  element_init(r, a->field);
-  element_init(q, a->field);
-  element_init(binv, pdp->field);
-  element_init(e0, pdp->field);
-  pq = q->data;
-  pr = r->data;
-  element_set(r, a);
-  k = m - n;
-  poly_alloc(q, k + 1);
-  element_invert(binv, poly_coeff(b, n));
-  while (k >= 0) {
-    qe = pq->coeff->item[k];
-    element_mul(qe, binv, pr->coeff->item[m]);
-    for (i=0; i<=n; i++) {
-      element_mul(e0, qe, poly_coeff(b, i));
-      element_sub(pr->coeff->item[i + k], pr->coeff->item[i + k], e0);
-    }
-    k--;
-    m--;
-  }
-  poly_remove_leading_zeroes(r);
-  element_set(quot, q);
-  element_set(rem, r);
-
-  element_clear(q);
-  element_clear(r);
-  element_clear(e0);
-  element_clear(binv);
-}
-
-static void poly_invert(element_ptr res, element_ptr f, element_ptr m) {
-  element_t q, r0, r1, r2;
-  element_t b0, b1, b2;
-  element_t inv;
-
-  element_init(b0, res->field);
-  element_init(b1, res->field);
-  element_init(b2, res->field);
-  element_init(q, res->field);
-  element_init(r0, res->field);
-  element_init(r1, res->field);
-  element_init(r2, res->field);
-  element_init(inv, poly_base_field(res));
-  element_set0(b0);
-  element_set1(b1);
-  element_set(r0, m);
-  element_set(r1, f);
-
-  for (;;) {
-    poly_div(q, r2, r0, r1);
-    if (element_is0(r2)) break;
-    element_mul(b2, b1, q);
-    element_sub(b2, b0, b2);
-    element_set(b0, b1);
-    element_set(b1, b2);
-    element_set(r0, r1);
-    element_set(r1, r2);
-  }
-  element_invert(inv, poly_coeff(r1, 0));
-  poly_const_mul(res, inv, b1);
-  element_clear(inv);
-  element_clear(q);
-  element_clear(r0);
-  element_clear(r1);
-  element_clear(r2);
-  element_clear(b0);
-  element_clear(b1);
-  element_clear(b2);
-}
-
-static void poly_to_polymod_truncate(element_ptr e, element_ptr f) {
-  mfptr p = e->field->data;
-  element_t *coeff = e->data;
-  int i;
-  int n;
-  n = poly_coeff_count(f);
-  if (n > p->n) n = p->n;
-
-  for (i=0; i<n; i++) {
-    element_set(coeff[i], poly_coeff(f, i));
-  }
-  for (; i<p->n; i++) {
-    element_set0(coeff[i]);
-  }
-}
-
-static void polymod_to_poly(element_ptr f, element_ptr e) {
-  mfptr p = e->field->data;
-  element_t *coeff = e->data;
-  int i, n = p->n;
-  poly_alloc(f, n);
-  for (i=0; i<n; i++) {
-    element_set(poly_coeff(f, i), coeff[i]);
-  }
-  poly_remove_leading_zeroes(f);
-}
-
-static void polymod_invert(element_ptr r, element_ptr e) {
-  mfptr p = r->field->data;
-  element_ptr minpoly = p->poly;
-  element_t f, r1;
-
-  element_init(f, minpoly->field);
-  element_init(r1, minpoly->field);
-  polymod_to_poly(f, e);
-
-  poly_invert(r1, f, p->poly);
-
-  poly_to_polymod_truncate(r, r1);
-
-  element_clear(f);
-  element_clear(r1);
-}
-
-static int poly_cmp(element_ptr f, element_ptr g) {
-  int i;
-  int n = poly_coeff_count(f);
-  int n1 = poly_coeff_count(g);
-  if (n != n1) return 1;
-  for (i=0; i<n; i++) {
-    if (element_cmp(poly_coeff(f, i), poly_coeff(g, i))) return 1;
-  }
-  return 0;
-}
-
-static void field_clear_poly(field_ptr f) {
-  pfptr p = f->data;
-  pbc_free(p);
-}
-
-// 2 bytes hold the number of terms, then the terms follow.
-// Bad for sparse polynomials.
-static int poly_length_in_bytes(element_t p) {
-  int count = poly_coeff_count(p);
-  int result = 2;
-  int i;
-  for (i=0; i<count; i++) {
-    result += element_length_in_bytes(poly_coeff(p, i));
-  }
-  return result;
-}
-
-static int poly_to_bytes(unsigned char *buf, element_t p) {
-  int count = poly_coeff_count(p);
-  int result = 2;
-  int i;
-  buf[0] = (unsigned char) count;
-  buf[1] = (unsigned char) (count >> 8);
-  for (i=0; i<count; i++) {
-    result += element_to_bytes(&buf[result], poly_coeff(p, i));
-  }
-  return result;
-}
-
-static int poly_from_bytes(element_t p, unsigned char *buf) {
-  int result = 2;
-  int count = buf[0] + buf[1] * 256;
-  int i;
-  poly_alloc(p, count);
-  for (i=0; i<count; i++) {
-    result += element_from_bytes(poly_coeff(p, i), &buf[result]);
-  }
-  return result;
-}
-
-// Is this useful? This returns to_mpz(constant term).
-static void poly_to_mpz(mpz_t z, element_ptr e) {
-  if (!poly_coeff_count(e)) {
-    mpz_set_ui(z, 0);
-  } else {
-    element_to_mpz(z, poly_coeff(e, 0));
-  }
-}
-
-static void poly_out_info(FILE *str, field_ptr f) {
-  pfptr p = f->data;
-  fprintf(str, "Polynomial ring over ");
-  field_out_info(str, p->field);
-}
-
-static void field_clear_polymod(field_ptr f) {
-  mfptr p = f->data;
-  int i, n = p->n;
-
-  for (i=0; i<n; i++) {
-    element_clear(p->xpwr[i]);
-  }
-  pbc_free(p->xpwr);
-
-  element_clear(p->poly);
-  pbc_free(f->data);
-}
-
-static int polymod_is_sqr(element_ptr e) {
-  int res;
-  mpz_t z;
-  element_t e0;
-
-  element_init(e0, e->field);
-  mpz_init(z);
-  mpz_sub_ui(z, e->field->order, 1);
-  mpz_divexact_ui(z, z, 2);
-
-  element_pow_mpz(e0, e, z);
-  res = element_is1(e0);
-  element_clear(e0);
-  mpz_clear(z);
-  return res;
-}
-
-// Find a square root in a polynomial modulo ring using Cantor-Zassenhaus aka
-// Legendre's method.
-static void polymod_sqrt(element_ptr res, element_ptr a) {
-  // TODO: Use a faster method? See Bernstein.
-  field_t kx;
-  element_t f;
-  element_t r, s;
-  element_t e0;
-  mpz_t z;
-
-  field_init_poly(kx, a->field);
-  mpz_init(z);
-  element_init(f, kx);
-  element_init(r, kx);
-  element_init(s, kx);
-  element_init(e0, a->field);
-
-  poly_alloc(f, 3);
-  element_set1(poly_coeff(f, 2));
-  element_neg(poly_coeff(f, 0), a);
-
-  mpz_sub_ui(z, a->field->order, 1);
-  mpz_divexact_ui(z, z, 2);
-  for (;;) {
-    int i;
-    element_ptr x;
-    element_ptr e1, e2;
-
-    poly_alloc(r, 2);
-    element_set1(poly_coeff(r, 1));
-    x = poly_coeff(r, 0);
-    element_random(x);
-    element_mul(e0, x, x);
-    if (!element_cmp(e0, a)) {
-      element_set(res, x);
-      break;
-    }
-    element_set1(s);
-    //TODO: this can be optimized greatly
-    //since we know r has the form ax + b
-    for (i = mpz_sizeinbase(z, 2) - 1; i >=0; i--) {
-      element_mul(s, s, s);
-      if (poly_degree(s) == 2) {
-        e1 = poly_coeff(s, 0);
-        e2 = poly_coeff(s, 2);
-        element_mul(e0, e2, a);
-        element_add(e1, e1, e0);
-        poly_alloc(s, 2);
-        poly_remove_leading_zeroes(s);
-      }
-      if (mpz_tstbit(z, i)) {
-        element_mul(s, s, r);
-        if (poly_degree(s) == 2) {
-          e1 = poly_coeff(s, 0);
-          e2 = poly_coeff(s, 2);
-          element_mul(e0, e2, a);
-          element_add(e1, e1, e0);
-          poly_alloc(s, 2);
-          poly_remove_leading_zeroes(s);
-        }
-      }
-    }
-    if (poly_degree(s) < 1) continue;
-    element_set1(e0);
-    e1 = poly_coeff(s, 0);
-    e2 = poly_coeff(s, 1);
-    element_add(e1, e1, e0);
-    element_invert(e0, e2);
-    element_mul(e0, e0, e1);
-    element_mul(e2, e0, e0);
-    if (!element_cmp(e2, a)) {
-      element_set(res, e0);
-      break;
-    }
-  }
-
-  mpz_clear(z);
-  element_clear(f);
-  element_clear(r);
-  element_clear(s);
-  element_clear(e0);
-  field_clear(kx);
-}
-
-static int polymod_to_bytes(unsigned char *data, element_t f) {
-  mfptr p = f->field->data;
-  element_t *coeff = f->data;
-  int i, n = p->n;
-  int len = 0;
-  for (i=0; i<n; i++) {
-    len += element_to_bytes(data + len, coeff[i]);
-  }
-  return len;
-}
-
-static int polymod_length_in_bytes(element_t f) {
-  mfptr p = f->field->data;
-  element_t *coeff = f->data;
-  int i, n = p->n;
-  int res = 0;
-
-  for (i=0; i<n; i++) {
-    res += element_length_in_bytes(coeff[i]);
-  }
-
-  return res;
-}
-
-static int polymod_from_bytes(element_t f, unsigned char *data) {
-  mfptr p = f->field->data;
-  element_t *coeff = f->data;
-  int i, n = p->n;
-  int len = 0;
-
-  for (i=0; i<n; i++) {
-    len += element_from_bytes(coeff[i], data + len);
-  }
-  return len;
-}
-
-static void polymod_init(element_t e) {
-  int i;
-  mfptr p = e->field->data;
-  int n = p->n;
-  element_t *coeff;
-  coeff = e->data = pbc_malloc(sizeof(element_t) * n);
-
-  for (i=0; i<n; i++) {
-    element_init(coeff[i], p->field);
-  }
-}
-
-static void polymod_clear(element_t e) {
-  mfptr p = e->field->data;
-  element_t *coeff = e->data;
-  int i, n = p->n;
-  for (i=0; i<n; i++) {
-    element_clear(coeff[i]);
-  }
-  pbc_free(e->data);
-}
-
-static void polymod_set_si(element_t e, signed long int x) {
-  mfptr p = e->field->data;
-  element_t *coeff = e->data;
-  int i, n = p->n;
-  element_set_si(coeff[0], x);
-  for (i=1; i<n; i++) {
-    element_set0(coeff[i]);
-  }
-}
-
-static void polymod_set_mpz(element_t e, mpz_t z) {
-  mfptr p = e->field->data;
-  element_t *coeff = e->data;
-  int i, n = p->n;
-  element_set_mpz(coeff[0], z);
-  for (i=1; i<n; i++) {
-    element_set0(coeff[i]);
-  }
-}
-
-static void polymod_set(element_t e, element_t f) {
-  mfptr p = e->field->data;
-  element_t *dst = e->data, *src = f->data;
-  int i, n = p->n;
-  for (i=0; i<n; i++) {
-    element_set(dst[i], src[i]);
-  }
-}
-
-static void polymod_neg(element_t e, element_t f) {
-  mfptr p = e->field->data;
-  element_t *dst = e->data, *src = f->data;
-  int i, n = p->n;
-  for (i=0; i<n; i++) {
-    element_neg(dst[i], src[i]);
-  }
-}
-
-static int polymod_cmp(element_ptr f, element_ptr g) {
-  mfptr p = f->field->data;
-  element_t *c1 = f->data, *c2 = g->data;
-  int i, n = p->n;
-  for (i=0; i<n; i++) {
-    if (element_cmp(c1[i], c2[i])) return 1;
-  }
-  return 0;
-}
-
-static void polymod_add(element_t r, element_t e, element_t f) {
-  mfptr p = r->field->data;
-  element_t *dst = r->data, *s1 = e->data, *s2 = f->data;
-  int i, n = p->n;
-  for (i=0; i<n; i++) {
-    element_add(dst[i], s1[i], s2[i]);
-  }
-}
-
-static void polymod_double(element_t r, element_t f) {
-  mfptr p = r->field->data;
-  element_t *dst = r->data, *s1 = f->data;
-  int i, n = p->n;
-  for (i=0; i<n; i++) {
-    element_double(dst[i], s1[i]);
-  }
-}
-
-static void polymod_sub(element_t r, element_t e, element_t f) {
-  mfptr p = r->field->data;
-  element_t *dst = r->data, *s1 = e->data, *s2 = f->data;
-  int i, n = p->n;
-  for (i=0; i<n; i++) {
-    element_sub(dst[i], s1[i], s2[i]);
-  }
-}
-
-static void polymod_mul_mpz(element_t e, element_t f, mpz_ptr z) {
-  mfptr p = e->field->data;
-  element_t *dst = e->data, *src = f->data;
-  int i, n = p->n;
-  for (i=0; i<n; i++) {
-    element_mul_mpz(dst[i], src[i], z);
-  }
-}
-
-static void polymod_mul_si(element_t e, element_t f, signed long int z) {
-  mfptr p = e->field->data;
-  element_t *dst = e->data, *src = f->data;
-  int i, n = p->n;
-  for (i=0; i<n; i++) {
-    element_mul_si(dst[i], src[i], z);
-  }
-}
-
-// Karatsuba multiplication for degree 2 polynomials.
-static void kar_poly_2(element_t *dst, element_t c3, element_t c4, element_t *s1, element_t *s2, element_t *scratch) {
-  element_ptr c01, c02, c12;
-
-  c12 = scratch[0];
-  c02 = scratch[1];
-  c01 = scratch[2];
-
-  element_add(c3, s1[0], s1[1]);
-  element_add(c4, s2[0], s2[1]);
-  element_mul(c01, c3, c4);
-  element_add(c3, s1[0], s1[2]);
-  element_add(c4, s2[0], s2[2]);
-  element_mul(c02, c3, c4);
-  element_add(c3, s1[1], s1[2]);
-  element_add(c4, s2[1], s2[2]);
-  element_mul(c12, c3, c4);
-
-  element_mul(dst[1], s1[1], s2[1]);
-
-  // Constant term.
-  element_mul(dst[0], s1[0], s2[0]);
-
-  // Coefficient of x^4.
-  element_mul(c4, s1[2], s2[2]);
-
-  // Coefficient of x^3.
-  element_add(c3, dst[1], c4);
-  element_sub(c3, c12, c3);
-
-  // Coefficient of x^2.
-  element_add(dst[2], c4, dst[0]);
-  element_sub(c02, c02, dst[2]);
-  element_add(dst[2], dst[1], c02);
-
-  // Coefficient of x.
-  element_sub(c01, c01, dst[0]);
-  element_sub(dst[1], c01, dst[1]);
-}
-
-// Degree 3, 6 polynomial moduli have dedicated routines for multiplication.
-static void polymod_mul_degree3(element_ptr res, element_ptr e, element_ptr f) {
-  mfptr p = res->field->data;
-  element_t *dst = res->data, *s1 = e->data, *s2 = f->data;
-  element_t c3, c4;
-  element_t p0;
-
-  element_init(p0, res->field);
-  element_init(c3, p->field);
-  element_init(c4, p->field);
-
-  kar_poly_2(dst, c3, c4, s1, s2, p0->data);
-
-  polymod_const_mul(p0, c3, p->xpwr[0]);
-  element_add(res, res, p0);
-  polymod_const_mul(p0, c4, p->xpwr[1]);
-  element_add(res, res, p0);
-
-  element_clear(p0);
-  element_clear(c3);
-  element_clear(c4);
-}
-
-static void polymod_mul_degree6(element_ptr res, element_ptr e, element_ptr f) {
-  mfptr p = res->field->data;
-  element_t *dst = res->data, *s0, *s1 = e->data, *s2 = f->data;
-  element_t *a0, *a1, *b0, *b1;
-  element_t p0, p1, p2, p3;
-
-  a0 = s1;
-  a1 = &s1[3];
-  b0 = s2;
-  b1 = &s2[3];
-
-  element_init(p0, res->field);
-  element_init(p1, res->field);
-  element_init(p2, res->field);
-  element_init(p3, res->field);
-
-  s0 = p0->data;
-  s1 = p1->data;
-  s2 = p2->data;
-  element_add(s0[0], a0[0], a1[0]);
-  element_add(s0[1], a0[1], a1[1]);
-  element_add(s0[2], a0[2], a1[2]);
-
-  element_add(s1[0], b0[0], b1[0]);
-  element_add(s1[1], b0[1], b1[1]);
-  element_add(s1[2], b0[2], b1[2]);
-
-  kar_poly_2(s2, s2[3], s2[4], s0, s1, p3->data);
-  kar_poly_2(s0, s0[3], s0[4], a0, b0, p3->data);
-  kar_poly_2(s1, s1[3], s1[4], a1, b1, p3->data);
-
-  element_set(dst[0], s0[0]);
-  element_set(dst[1], s0[1]);
-  element_set(dst[2], s0[2]);
-
-  element_sub(dst[3], s0[3], s0[0]);
-  element_sub(dst[3], dst[3], s1[0]);
-  element_add(dst[3], dst[3], s2[0]);
-
-  element_sub(dst[4], s0[4], s0[1]);
-  element_sub(dst[4], dst[4], s1[1]);
-  element_add(dst[4], dst[4], s2[1]);
-
-  element_sub(dst[5], s2[2], s0[2]);
-  element_sub(dst[5], dst[5], s1[2]);
-
-  // Start reusing part of s0 as scratch space(!)
-  element_sub(s0[0], s2[3], s0[3]);
-  element_sub(s0[0], s0[0], s1[3]);
-  element_add(s0[0], s0[0], s1[0]);
-
-  element_sub(s0[1], s2[4], s0[4]);
-  element_sub(s0[1], s0[1], s1[4]);
-  element_add(s0[1], s0[1], s1[1]);
-
-  polymod_const_mul(p3, s0[0], p->xpwr[0]);
-  element_add(res, res, p3);
-  polymod_const_mul(p3, s0[1], p->xpwr[1]);
-  element_add(res, res, p3);
-  polymod_const_mul(p3, s1[2], p->xpwr[2]);
-  element_add(res, res, p3);
-  polymod_const_mul(p3, s1[3], p->xpwr[3]);
-  element_add(res, res, p3);
-  polymod_const_mul(p3, s1[4], p->xpwr[4]);
-  element_add(res, res, p3);
-
-  element_clear(p0);
-  element_clear(p1);
-  element_clear(p2);
-  element_clear(p3);
-}
-
-// General polynomial modulo ring multiplication.
-static void polymod_mul(element_ptr res, element_ptr e, element_ptr f) {
-  mfptr p = res->field->data;
-  int n = p->n;
-  element_t *dst;
-  element_t *s1 = e->data, *s2 = f->data;
-  element_t prod, p0, c0;
-  int i, j;
-  element_t *high;  // Coefficients of x^n, ..., x^{2n-2}.
-
-  high = pbc_malloc(sizeof(element_t) * (n - 1));
-  for (i=0; i<n-1; i++) {
-    element_init(high[i], p->field);
-    element_set0(high[i]);
-  }
-  element_init(prod, res->field);
-  dst = prod->data;
-  element_init(p0, res->field);
-  element_init(c0, p->field);
-
-  for (i=0; i<n; i++) {
-    int ni = n - i;
-    for (j=0; j<ni; j++) {
-      element_mul(c0, s1[i], s2[j]);
-      element_add(dst[i + j], dst[i + j], c0);
-    }
-    for (;j<n; j++) {
-      element_mul(c0, s1[i], s2[j]);
-      element_add(high[j - ni], high[j - ni], c0);
-    }
-  }
-
-  for (i=0; i<n-1; i++) {
-    polymod_const_mul(p0, high[i], p->xpwr[i]);
-    element_add(prod, prod, p0);
-    element_clear(high[i]);
-  }
-  pbc_free(high);
-
-  element_set(res, prod);
-  element_clear(prod);
-  element_clear(p0);
-  element_clear(c0);
-}
-
-static void polymod_square_degree3(element_ptr res, element_ptr e) {
-  // TODO: Investigate if squaring is significantly cheaper than
-  // multiplication. If so convert to Karatsuba.
-  element_t *dst = res->data;
-  element_t *src = e->data;
-  mfptr p = res->field->data;
-  element_t p0;
-  element_t c0, c2;
-  element_ptr c1, c3;
-
-  element_init(p0, res->field);
-  element_init(c0, p->field);
-  element_init(c2, p->field);
-
-  c3 = p0->data;
-  c1 = c3 + 1;
-
-  element_mul(c3, src[0], src[1]);
-  element_mul(c1, src[0], src[2]);
-  element_square(dst[0], src[0]);
-
-  element_mul(c2, src[1], src[2]);
-  element_square(c0, src[2]);
-  element_square(dst[2], src[1]);
-
-  element_add(dst[1], c3, c3);
-
-  element_add(c1, c1, c1);
-  element_add(dst[2], dst[2], c1);
-
-  polymod_const_mul(p0, c0, p->xpwr[1]);
-  element_add(res, res, p0);
-
-  element_add(c2, c2, c2);
-  polymod_const_mul(p0, c2, p->xpwr[0]);
-  element_add(res, res, p0);
-
-  element_clear(p0);
-  element_clear(c0);
-  element_clear(c2);
-}
-
-static void polymod_square(element_ptr res, element_ptr e) {
-  element_t *dst;
-  element_t *src = e->data;
-  mfptr p = res->field->data;
-  int n = p->n;
-  element_t prod, p0, c0;
-  int i, j;
-  element_t *high; // Coefficients of x^n,...,x^{2n-2}.
-
-  high = pbc_malloc(sizeof(element_t) * (n - 1));
-  for (i=0; i<n-1; i++) {
-    element_init(high[i], p->field);
-    element_set0(high[i]);
-  }
-
-  element_init(prod, res->field);
-  dst = prod->data;
-  element_init(p0, res->field);
-  element_init(c0, p->field);
-
-  for (i=0; i<n; i++) {
-    int twicei = 2 * i;
-    element_square(c0, src[i]);
-    if (twicei < n) {
-      element_add(dst[twicei], dst[twicei], c0);
-    } else {
-      element_add(high[twicei - n], high[twicei - n], c0);
-    }
-
-    for (j=i+1; j<n-i; j++) {
-      element_mul(c0, src[i], src[j]);
-      element_add(c0, c0, c0);
-      element_add(dst[i + j], dst[i + j], c0);
-    }
-    for (;j<n; j++) {
-      element_mul(c0, src[i], src[j]);
-      element_add(c0, c0, c0);
-      element_add(high[i + j - n], high[i + j - n], c0);
-    }
-  }
-
-  for (i=0; i<n-1; i++) {
-    polymod_const_mul(p0, high[i], p->xpwr[i]);
-    element_add(prod, prod, p0);
-    element_clear(high[i]);
-  }
-  pbc_free(high);
-
-  element_set(res, prod);
-  element_clear(prod);
-  element_clear(p0);
-  element_clear(c0);
-}
-
-static int polymod_is0(element_ptr e) {
-  mfptr p = e->field->data;
-  element_t *coeff = e->data;
-  int i, n = p->n;
-
-  for (i=0; i<n; i++) {
-    if (!element_is0(coeff[i])) return 0;
-  }
-  return 1;
-}
-
-static int polymod_is1(element_ptr e) {
-  mfptr p = e->field->data;
-  element_t *coeff = e->data;
-  int i, n = p->n;
-
-  if (!element_is1(coeff[0])) return 0;
-  for (i=1; i<n; i++) {
-    if (!element_is0(coeff[i])) return 0;
-  }
-  return 1;
-}
-
-static void polymod_set0(element_ptr e) {
-  mfptr p = e->field->data;
-  element_t *coeff = e->data;
-  int i, n = p->n;
-
-  for (i=0; i<n; i++) {
-    element_set0(coeff[i]);
-  }
-}
-
-static void polymod_set1(element_ptr e) {
-  mfptr p = e->field->data;
-  element_t *coeff = e->data;
-  int i, n = p->n;
-
-  element_set1(coeff[0]);
-  for (i=1; i<n; i++) {
-    element_set0(coeff[i]);
-  }
-}
-
-static int polymod_sgn(element_ptr e) {
-  mfptr p = e->field->data;
-  element_t *coeff = e->data;
-  int res = 0;
-  int i, n = p->n;
-  for (i=0; i<n; i++) {
-    res = element_sgn(coeff[i]);
-    if (res) break;
-  }
-  return res;
-}
-
-static size_t polymod_out_str(FILE *stream, int base, element_ptr e) {
-  size_t result = 2, status;
-  mfptr p = e->field->data;
-  element_t *coeff = e->data;
-  int i, n = p->n;
-
-  if (EOF == fputc('[', stream)) return 0;
-  for (i=0; i<n; i++) {
-    if (i) {
-      if (EOF == fputs(", ", stream)) return 0;
-      result += 2;
-    }
-    status = element_out_str(stream, base, coeff[i]);
-    if (!status) return 0;
-    result += status;
-  }
-  if (EOF == fputc(']', stream)) return 0;
-  return result;
-}
-
-static int polymod_snprint(char *s, size_t size, element_ptr e) {
-  mfptr p = e->field->data;
-  element_t *coeff = e->data;
-  int i, n = p->n;
-  size_t result = 0, left;
-  int status;
-
-  #define clip_sub(void) {                     \
-    result += status;                          \
-    left = result >= size ? 0 : size - result; \
-  }
-
-  status = snprintf(s, size, "[");
-  if (status < 0) return status;
-  clip_sub();
-
-  for (i=0; i<n; i++) {
-    if (i) {
-      status = snprintf(s + result, left, ", ");
-      if (status < 0) return status;
-      clip_sub();
-    }
-    status = element_snprint(s + result, left, coeff[i]);
-    if (status < 0) return status;
-    clip_sub();
-  }
-  status = snprintf(s + result, left, "]");
-  if (status < 0) return status;
-  return result + status;
-  #undef clip_sub
-}
-
-static void polymod_set_multiz(element_ptr e, multiz m) {
-  mfptr p = e->field->data;
-  element_t *coeff = e->data;
-  int i, n = p->n;
-  if (multiz_is_z(m)) {
-    element_set_multiz(coeff[0], m);
-    for (i = 1; i < n; i++) element_set0(coeff[i]);
-    return;
-  }
-  int max = multiz_count(m);
-  for (i = 0; i < n; i++) {
-    if (i >= max) element_set0(coeff[i]);
-    else element_set_multiz(coeff[i], multiz_at(m, i));
-  }
-}
-
-static int polymod_set_str(element_ptr e, const char *s, int base) {
-  mfptr p = e->field->data;
-  element_t *coeff = e->data;
-  int i, n = p->n;
-  const char *cp = s;
-  element_set0(e);
-  while (*cp && isspace(*cp)) cp++;
-  if (*cp++ != '[') return 0;
-  for (i=0; i<n; i++) {
-    cp += element_set_str(coeff[i], cp, base);
-    while (*cp && isspace(*cp)) cp++;
-    if (i<n-1 && *cp++ != ',') return 0;
-  }
-  if (*cp++ != ']') return 0;
-  return cp - s;
-}
-
-static int polymod_coeff_count(element_ptr e) {
-  UNUSED_VAR(e);
-  mfptr p = e->field->data;
-  return p->n;
-}
-
-static element_ptr polymod_coeff(element_ptr e, int i) {
-  element_t *coeff = e->data;
-  return coeff[i];
-}
-
-static void polymod_to_mpz(mpz_t z, element_ptr e) {
-  element_to_mpz(z, polymod_coeff(e, 0));
-}
-
-// Compute x^n,...,x^{2n-2} mod poly.
-static void compute_x_powers(field_ptr field, element_ptr poly) {
-  mfptr p = field->data;
-  element_t p0;
-  element_ptr pwrn;
-  element_t *coeff, *coeff1;
-  int i, j;
-  int n = p->n;
-  element_t *xpwr;
-
-  xpwr = p->xpwr;
-
-  element_init(p0, field);
-  for (i=0; i<n; i++) {
-    element_init(xpwr[i], field);
-  }
-  pwrn = xpwr[0];
-  poly_to_polymod_truncate(pwrn, poly);
-  element_neg(pwrn, pwrn);
-
-  for (i=1; i<n; i++) {
-    coeff = xpwr[i-1]->data;
-    coeff1 = xpwr[i]->data;
-
-    element_set0(coeff1[0]);
-    for (j=1; j<n; j++) {
-      element_set(coeff1[j], coeff[j - 1]);
-    }
-    polymod_const_mul(p0, coeff[n - 1], pwrn);
-    element_add(xpwr[i], xpwr[i], p0);
-  }
-  element_clear(p0);
-}
-
-static void polymod_out_info(FILE *str, field_ptr f) {
-  mfptr p = f->data;
-  element_fprintf(str, "Extension, poly = %B, base field = ", p->poly);
-  field_out_info(str, p->field);
-}
-
-// Sets d = gcd(f, g).
-static void poly_gcd(element_ptr d, element_ptr f, element_ptr g) {
-  element_t a, b, q, r;
-  element_init(a, d->field);
-  element_init(b, d->field);
-  element_init(q, d->field);
-  element_init(r, d->field);
-
-  element_set(a, f);
-  element_set(b, g);
-  for(;;) {
-    //TODO: don't care about q
-    poly_div(q, r, a, b);
-    if (element_is0(r)) break;
-    element_set(a, b);
-    element_set(b, r);
-  }
-  element_set(d, b);
-  element_clear(a);
-  element_clear(b);
-  element_clear(q);
-  element_clear(r);
-}
-
-// Sets f = c g where c is the inverse of the leading coefficient of g.
-static void poly_make_monic(element_t f, element_t g) {
-  int n = poly_coeff_count(g);
-  int i;
-  element_ptr e0;
-  poly_alloc(f, n);
-  if (!n) return;
-
-  e0 = poly_coeff(f, n - 1);
-  element_invert(e0, poly_coeff(g, n - 1));
-  for (i=0; i<n-1; i++) {
-    element_mul(poly_coeff(f, i), poly_coeff(g, i), e0);
-  }
-  element_set1(e0);
-}
-
-// The above should be static.
-
-void field_init_poly(field_ptr f, field_ptr base_field) {
-  field_init(f);
-  pfptr p = f->data = pbc_malloc(sizeof(*p));
-  p->field = base_field;
-  p->mapbase = element_field_to_poly;
-  f->field_clear = field_clear_poly;
-  f->init = poly_init;
-  f->clear = poly_clear;
-  f->set_si = poly_set_si;
-  f->set_multiz = poly_set_multiz;
-  f->set_mpz = poly_set_mpz;
-  f->to_mpz = poly_to_mpz;
-  f->out_str = poly_out_str;
-  f->snprint = poly_snprint;
-  f->set = poly_set;
-  f->sign = poly_sgn;
-  f->add = poly_add;
-  f->doub = poly_double;
-  f->is0 = poly_is0;
-  f->is1 = poly_is1;
-  f->set0 = poly_set0;
-  f->set1 = poly_set1;
-  f->sub = poly_sub;
-  f->neg = poly_neg;
-  f->mul = poly_mul;
-  f->mul_mpz = poly_mul_mpz;
-  f->mul_si = poly_mul_si;
-  f->cmp = poly_cmp;
-  f->out_info = poly_out_info;
-  f->item_count = poly_coeff_count;
-  f->item = poly_coeff;
-
-  f->to_bytes = poly_to_bytes;
-  f->from_bytes = poly_from_bytes;
-  f->fixed_length_in_bytes = -1;
-  f->length_in_bytes = poly_length_in_bytes;
-}
-
-void poly_set_coeff(element_ptr e, element_ptr a, int n) {
-  peptr p = e->data;
-  if (p->coeff->count < n + 1) {
-    poly_alloc(e, n + 1);
-  }
-  element_ptr e0 = p->coeff->item[n];
-  element_set(e0, a);
-  if (p->coeff->count == n + 1 && element_is0(a)) poly_remove_leading_zeroes(e);
-}
-
-void poly_set_coeff0(element_ptr e, int n) {
-  peptr p = e->data;
-  if (n < p->coeff->count) {
-    element_set0(p->coeff->item[n]);
-    if (n == p->coeff->count - 1) poly_remove_leading_zeroes(e);
-  }
-}
-
-void poly_set_coeff1(element_ptr e, int n) {
-  peptr p = e->data;
-  if (p->coeff->count < n + 1) {
-    poly_alloc(e, n + 1);
-  }
-  element_set1(p->coeff->item[n]);
-}
-
-void poly_setx(element_ptr f) {
-  poly_alloc(f, 2);
-  element_set1(poly_coeff(f, 1));
-  element_set0(poly_coeff(f, 0));
-}
-
-void poly_const_mul(element_ptr res, element_ptr a, element_ptr poly) {
-  int i, n = poly_coeff_count(poly);
-  poly_alloc(res, n);
-  for (i=0; i<n; i++) {
-    element_mul(poly_coeff(res, i), a, poly_coeff(poly, i));
-  }
-  poly_remove_leading_zeroes(res);
-}
-
-void poly_random_monic(element_ptr f, int deg) {
-  int i;
-  poly_alloc(f, deg + 1);
-  for (i=0; i<deg; i++) {
-    element_random(poly_coeff(f, i));
-  }
-  element_set1(poly_coeff(f, i));
-}
-
-int polymod_field_degree(field_t f) {
-  mfptr p = f->data;
-  return p->n;
-}
-
-void field_init_polymod(field_ptr f, element_ptr poly) {
-  pfptr pdp = poly->field->data;
-  field_init(f);
-  mfptr p = f->data = pbc_malloc(sizeof(*p));
-  p->field = pdp->field;
-  p->mapbase = element_field_to_poly;
-  element_init(p->poly, poly->field);
-  element_set(p->poly, poly);
-  int n = p->n = poly_degree(p->poly);
-  f->field_clear = field_clear_polymod;
-  f->init = polymod_init;
-  f->clear = polymod_clear;
-  f->set_si = polymod_set_si;
-  f->set_mpz = polymod_set_mpz;
-  f->out_str = polymod_out_str;
-  f->snprint = polymod_snprint;
-  f->set_multiz = polymod_set_multiz;
-  f->set_str = polymod_set_str;
-  f->set = polymod_set;
-  f->sign = polymod_sgn;
-  f->add = polymod_add;
-  f->doub = polymod_double;
-  f->sub = polymod_sub;
-  f->neg = polymod_neg;
-  f->is0 = polymod_is0;
-  f->is1 = polymod_is1;
-  f->set0 = polymod_set0;
-  f->set1 = polymod_set1;
-  f->cmp = polymod_cmp;
-  f->to_mpz = polymod_to_mpz;
-  f->item_count = polymod_coeff_count;
-  f->item = polymod_coeff;
-  switch(n) {
-    case 3:
-      f->mul = polymod_mul_degree3;
-      f->square = polymod_square_degree3;
-      break;
-    case 6:
-      f->mul = polymod_mul_degree6;
-      f->square = polymod_square;
-      break;
-    default:
-      f->mul = polymod_mul;
-      f->square = polymod_square;
-      break;
-  }
-
-  f->mul_mpz = polymod_mul_mpz;
-  f->mul_si = polymod_mul_si;
-  f->random = polymod_random;
-  f->from_hash = polymod_from_hash;
-  f->invert = polymod_invert;
-  f->is_sqr = polymod_is_sqr;
-  f->sqrt = polymod_sqrt;
-  f->to_bytes = polymod_to_bytes;
-  f->from_bytes = polymod_from_bytes;
-  f->out_info = polymod_out_info;
-
-  if (pdp->field->fixed_length_in_bytes < 0) {
-    f->fixed_length_in_bytes = -1;
-    f->length_in_bytes = polymod_length_in_bytes;
-  } else {
-    f->fixed_length_in_bytes = pdp->field->fixed_length_in_bytes * poly_degree(poly);
-  }
-  mpz_pow_ui(f->order, p->field->order, n);
-
-  p->xpwr = pbc_malloc(sizeof(element_t) * n);
-  compute_x_powers(f, poly);
-}
-
-field_ptr poly_base_field(element_t f) {
-  return ((pfptr) f->field->data)->field;
-}
-
-void polymod_const_mul(element_ptr res, element_ptr a, element_ptr e) {
-  // a lies in R, e in R[x].
-  element_t *coeff = e->data, *dst = res->data;
-  int i, n = polymod_field_degree(e->field);
-
-  for (i=0; i<n; i++) {
-    element_mul(dst[i], coeff[i], a);
-  }
-}
-
-struct checkgcd_scope_var {
-  mpz_ptr     z, deg;
-  field_ptr   basef;
-  element_ptr xpow, x, f, g;
-};
-
-// Returns 0 if gcd(x^q^{n/d} - x, f) = 1, 1 otherwise.
-static int checkgcd(mpz_ptr fac, unsigned int mul, struct checkgcd_scope_var *v) {
-  UNUSED_VAR(mul);
-  mpz_divexact(v->z, v->deg, fac);
-  mpz_pow_ui(v->z, v->basef->order, mpz_get_ui(v->z));
-  element_pow_mpz(v->xpow, v->x, v->z);
-  element_sub(v->xpow, v->xpow, v->x);
-  if (element_is0(v->xpow)) return 1;
-  polymod_to_poly(v->g, v->xpow);
-  poly_gcd(v->g, v->f, v->g);
-  return poly_degree(v->g) != 0;
-}
-
-// Returns 1 if polynomial is irreducible, 0 otherwise.
-// A polynomial f(x) is irreducible in F_q[x] if and only if:
-//  (1) f(x) | x^{q^n} - x, and
-//  (2) gcd(f(x), x^{q^{n/d}} - x) = 1 for all primes d | n.
-// (Recall GF(p) is the splitting field for x^p - x.)
-int poly_is_irred(element_ptr f) {
-  int res = 0;
-  element_t xpow, x, g;
-  field_ptr basef = poly_base_field(f);
-  field_t rxmod;
-
-  // 0, units are not irreducibles.
-  // Assume coefficients are from a field.
-  if (poly_degree(f) <= 0) return 0;
-  // Degree 1 polynomials are always irreducible.
-  if (poly_degree(f) == 1) return 1;
-
-  field_init_polymod(rxmod, f);
-  element_init(xpow, rxmod);
-  element_init(x, rxmod);
-  element_init(g, f->field);
-  element_set1(polymod_coeff(x, 1));
-
-  // The degree fits in an unsigned int but I'm lazy and want to use my
-  // mpz trial division code.
-  mpz_t deg, z;
-  mpz_init(deg);
-  mpz_init(z);
-  mpz_set_ui(deg, poly_degree(f));
-
-  struct checkgcd_scope_var v = {.z = z, .deg = deg, .basef = basef,
-                                 .xpow = xpow, .x = x, .f = f, .g = g};
-  if (!pbc_trial_divide((int(*)(mpz_t,unsigned,void*))checkgcd, &v, deg, NULL)) {
-    // By now condition (2) has been satisfied. Check (1).
-    mpz_pow_ui(z, basef->order, poly_degree(f));
-    element_pow_mpz(xpow, x, z);
-    element_sub(xpow, xpow, x);
-    if (element_is0(xpow)) res = 1;
-  }
-
-  mpz_clear(deg);
-  mpz_clear(z);
-  element_clear(g);
-  element_clear(xpow);
-  element_clear(x);
-  field_clear(rxmod);
-  return res;
-}
-
-void element_field_to_poly(element_ptr f, element_ptr g) {
-  poly_alloc(f, 1);
-  element_set(poly_coeff(f, 0), g);
-  poly_remove_leading_zeroes(f);
-}
-
-void element_field_to_polymod(element_ptr f, element_ptr g) {
-  mfptr p = f->field->data;
-  element_t *coeff = f->data;
-  int i, n = p->n;
-  element_set(coeff[0], g);
-  for (i=1; i<n; i++) {
-    element_set0(coeff[i]);
-  }
-}
-
-// Returns 0 when a root exists and sets root to one of the roots.
-int poly_findroot(element_ptr root, element_ptr poly) {
-  // Compute gcd(x^q - x, poly).
-  field_t fpxmod;
-  element_t p, x, r, fac, g;
-  mpz_t q;
-
-  mpz_init(q);
-  mpz_set(q, poly_base_field(poly)->order);
-
-  field_init_polymod(fpxmod, poly);
-  element_init(p, fpxmod);
-  element_init(x, fpxmod);
-  element_init(g, poly->field);
-  element_set1(((element_t *) x->data)[1]);
-pbc_info("findroot: degree %d...", poly_degree(poly));
-  element_pow_mpz(p, x, q);
-  element_sub(p, p, x);
-
-  polymod_to_poly(g, p);
-  element_clear(p);
-  poly_gcd(g, g, poly);
-  poly_make_monic(g, g);
-  element_clear(x);
-  field_clear(fpxmod);
-
-  if (!poly_degree(g)) {
-    printf("no roots!\n");
-    mpz_clear(q);
-    element_clear(g);
-    return -1;
-  }
-
-  // Cantor-Zassenhaus algorithm.
-  element_init(fac, g->field);
-  element_init(x, g->field);
-  element_set_si(x, 1);
-  mpz_sub_ui(q, q, 1);
-  mpz_divexact_ui(q, q, 2);
-  element_init(r, g->field);
-  for (;;) {
-    if (poly_degree(g) == 1) break;  // Found a root!
-step_random:
-    poly_random_monic(r, 1);
-    // TODO: evaluate at g instead of bothering with gcd
-    poly_gcd(fac, r, g);
-
-    if (poly_degree(fac) > 0) {
-      poly_make_monic(g, fac);
-    } else {
-      field_init_polymod(fpxmod, g);
-      int n;
-      element_init(p, fpxmod);
-
-      poly_to_polymod_truncate(p, r);
-pbc_info("findroot: degree %d...", poly_degree(g));
-      element_pow_mpz(p, p, q);
-
-      polymod_to_poly(r, p);
-      element_clear(p);
-      field_clear(fpxmod);
-
-      element_add(r, r, x);
-      poly_gcd(fac, r, g);
-      n = poly_degree(fac);
-      if (n > 0 && n < poly_degree(g)) {
-        poly_make_monic(g, fac);
-      } else {
-        goto step_random;
-      }
-    }
-  }
-pbc_info("findroot: found root");
-  element_neg(root, poly_coeff(g, 0));
-  element_clear(r);
-  mpz_clear(q);
-  element_clear(x);
-  element_clear(g);
-  element_clear(fac);
-  return 0;
-}