init attribute-based encryption

Change-Id: Iba1a3d722110abf747a0fba366f3ebc911d25b25

16 files changed, 7598 insertions, 0 deletions

diff --git a/moon-abe/pbc-0.5.14/arith/dlog.c b/moon-abe/pbc-0.5.14/arith/dlog.c
new file mode 100644
index 00000000..f77df1b7
--- /dev/null
+++ b/moon-abe/pbc-0.5.14/arith/dlog.c
@@ -0,0 +1,187 @@
+// Brute force and Pollard rho discrete log algorithms.
+
+#include <stdarg.h>
+#include <stdint.h> // for intptr_t
+#include <stdio.h>
+#include <gmp.h>
+#include "pbc_utils.h"
+#include "pbc_field.h"
+#include "pbc_memory.h"
+#include "misc/darray.h"
+
+struct snapshot_s {
+  element_t a;
+  element_t b;
+  element_t snark;
+};
+typedef struct snapshot_s *snapshot_ptr;
+
+static void record(element_t asum, element_t bsum, element_t snark,
+                   darray_t hole, mpz_t counter) {
+  snapshot_ptr ss = pbc_malloc(sizeof(struct snapshot_s));
+  element_init_same_as(ss->a, asum);
+  element_init_same_as(ss->b, bsum);
+  element_init_same_as(ss->snark, snark);
+  element_set(ss->a, asum);
+  element_set(ss->b, bsum);
+  element_set(ss->snark, snark);
+  darray_append(hole, ss);
+  element_printf("snark %Zd: %B\n", counter, snark);
+}
+
+// g, h in some group of order r
+// finds x such that g^x = h
+// will hang if no such x exists
+// x in some field_t that set_mpz makes sense for
+void element_dlog_brute_force(element_t x, element_t g, element_t h) {
+  element_t g0;
+  mpz_t count;
+
+  mpz_init(count);
+  element_init_same_as(g0, g);
+
+  element_set(g0, g);
+  mpz_set_ui(count, 1);
+  while (element_cmp(g0, h)) {
+    element_mul(g0, g0, g);
+//element_printf("g0^%Zd = %B\n", count, g0);
+    mpz_add_ui(count, count, 1);
+  }
+  element_set_mpz(x, count);
+  mpz_clear(count);
+  element_clear(g0);
+}
+
+// x in Z_r, g, h in some group of order r
+// finds x such that g^x = h
+void element_dlog_pollard_rho(element_t x, element_t g, element_t h) {
+// see Blake, Seroussi and Smart
+// only one snark for this implementation
+  int i, s = 20;
+  field_ptr Zr = x->field, G = g->field;
+  element_t asum;
+  element_t bsum;
+  element_t a[s];
+  element_t b[s];
+  element_t m[s];
+  element_t g0, snark;
+  darray_t hole;
+  int interval = 5;
+  mpz_t counter;
+  int found = 0;
+
+  mpz_init(counter);
+  element_init(g0, G);
+  element_init(snark, G);
+  element_init(asum, Zr);
+  element_init(bsum, Zr);
+  darray_init(hole);
+  //set up multipliers
+  for (i = 0; i < s; i++) {
+    element_init(a[i], Zr);
+    element_init(b[i], Zr);
+    element_init(m[i], G);
+    element_random(a[i]);
+    element_random(b[i]);
+    element_pow_zn(g0, g, a[i]);
+    element_pow_zn(m[i], h, b[i]);
+    element_mul(m[i], m[i], g0);
+  }
+
+  element_random(asum);
+  element_random(bsum);
+  element_pow_zn(g0, g, asum);
+  element_pow_zn(snark, h, bsum);
+  element_mul(snark, snark, g0);
+
+  record(asum, bsum, snark, hole, counter);
+  for (;;) {
+    int len = element_length_in_bytes(snark);
+    unsigned char *buf = pbc_malloc(len);
+    unsigned char hash = 0;
+
+    element_to_bytes(buf, snark);
+    for (i = 0; i < len; i++) {
+      hash += buf[i];
+    }
+    i = hash % s;
+    pbc_free(buf);
+
+    element_mul(snark, snark, m[i]);
+    element_add(asum, asum, a[i]);
+    element_add(bsum, bsum, b[i]);
+
+    for (i = 0; i < hole->count; i++) {
+      snapshot_ptr ss = hole->item[i];
+      if (!element_cmp(snark, ss->snark)) {
+        element_sub(bsum, bsum, ss->b);
+        element_sub(asum, ss->a, asum);
+        //answer is x such that x * bsum = asum
+        //complications arise if gcd(bsum, r) > 1
+        //which can happen if r is not prime
+        if (!mpz_probab_prime_p(Zr->order, 10)) {
+          mpz_t za, zb, zd, zm;
+
+          mpz_init(za);
+          mpz_init(zb);
+          mpz_init(zd);
+          mpz_init(zm);
+
+          element_to_mpz(za, asum);
+          element_to_mpz(zb, bsum);
+          mpz_gcd(zd, zb, Zr->order);
+          mpz_divexact(zm, Zr->order, zd);
+          mpz_divexact(zb, zb, zd);
+          //if zd does not divide za there is no solution
+          mpz_divexact(za, za, zd);
+          mpz_invert(zb, zb, zm);
+          mpz_mul(zb, za, zb);
+          mpz_mod(zb, zb, zm);
+          do {
+            element_pow_mpz(g0, g, zb);
+            if (!element_cmp(g0, h)) {
+              element_set_mpz(x, zb);
+              break;
+            }
+            mpz_add(zb, zb, zm);
+            mpz_sub_ui(zd, zd, 1);
+          } while (mpz_sgn(zd));
+          mpz_clear(zm);
+          mpz_clear(za);
+          mpz_clear(zb);
+          mpz_clear(zd);
+        } else {
+          element_div(x, asum, bsum);
+        }
+        found = 1;
+        break;
+      }
+    }
+    if (found) break;
+
+    mpz_add_ui(counter, counter, 1);
+    if (mpz_tstbit(counter, interval)) {
+      record(asum, bsum, snark, hole, counter);
+      interval++;
+    }
+  }
+
+  for (i = 0; i < s; i++) {
+    element_clear(a[i]);
+    element_clear(b[i]);
+    element_clear(m[i]);
+  }
+  element_clear(g0);
+  element_clear(snark);
+  for (i = 0; i < hole->count; i++) {
+    snapshot_ptr ss = hole->item[i];
+    element_clear(ss->a);
+    element_clear(ss->b);
+    element_clear(ss->snark);
+    pbc_free(ss);
+  }
+  darray_clear(hole);
+  element_clear(asum);
+  element_clear(bsum);
+  mpz_clear(counter);
+}
diff --git a/moon-abe/pbc-0.5.14/arith/fasterfp.c b/moon-abe/pbc-0.5.14/arith/fasterfp.c
new file mode 100644
index 00000000..5ce8243a
--- /dev/null
+++ b/moon-abe/pbc-0.5.14/arith/fasterfp.c
@@ -0,0 +1,546 @@
+// Naive implementation of F_p.
+// It uses lowlevel GMP routines (mpn_* functions) like fastfp.c, but also
+// has a flag for the value 0, avoiding many memsets.
+//
+// I'm thinking of using the flag to also represent 1, -1,
+// but that complicates the logic even more, and I believe I need more
+// control than GMP is willing to give in order to avoid expensive
+// checks for 1, -1 everywhere.
+//
+// NOTE: does not work for moduli of the form:
+//   2^(something * 8 * sizeof(mp_limb_t))
+// See comments in add, double code.
+// (This kind of integer mod ring deserves its own implementation anyway.)
+
+#include <stdarg.h>
+#include <stdio.h>
+#include <stdint.h> // for intptr_t
+#include <stdlib.h>
+#include <string.h>
+#include <gmp.h>
+#include "pbc_utils.h"
+#include "pbc_field.h"
+#include "pbc_random.h"
+#include "pbc_fp.h"
+#include "pbc_memory.h"
+
+struct fp_field_data_s {
+  size_t limbs;
+  size_t bytes;
+  mp_limb_t *primelimbs;
+};
+typedef struct fp_field_data_s fp_field_data_t[1];
+typedef struct fp_field_data_s *fp_field_data_ptr;
+
+struct data_s {
+  int flag;
+  mp_limb_t *d;
+};
+typedef struct data_s *dataptr;
+
+static void fp_init(element_ptr e) {
+  fp_field_data_ptr p = e->field->data;
+  dataptr dp = e->data = pbc_malloc(sizeof(struct data_s));
+  dp->flag = 0;
+  dp->d = pbc_malloc(p->bytes);
+}
+
+static void fp_clear(element_ptr e) {
+  dataptr dp = e->data;
+  pbc_free(dp->d);
+  pbc_free(e->data);
+}
+
+//assumes z is nonzero
+static inline void from_mpz(element_ptr e, mpz_ptr z) {
+  fp_field_data_ptr p = e->field->data;
+  size_t count;
+  dataptr dp = e->data;
+  mpz_export(dp->d, &count, -1, sizeof(mp_limb_t), 0, 0, z);
+  memset((void *) (((unsigned char *) dp->d) + count * sizeof(mp_limb_t)),
+         0, (p->limbs - count) * sizeof(mp_limb_t));
+}
+
+static void fp_set_mpz(element_ptr e, mpz_ptr z) {
+  dataptr dp = e->data;
+  if (!mpz_sgn(z)) {
+    dp->flag = 0;
+  } else {
+    mpz_t tmp;
+    mpz_init(tmp);
+    mpz_mod(tmp, z, e->field->order);
+    from_mpz(e, tmp);
+    mpz_clear(tmp);
+    dp->flag = 2;
+  }
+}
+
+static void fp_set_si(element_ptr e, signed long int op) {
+  dataptr dp = e->data;
+  if (!op) {
+    dp->flag = 0;
+  } else {
+    const fp_field_data_ptr p = e->field->data;
+    const size_t t = p->limbs;
+    if (op < 0) {
+      mpn_sub_1(dp->d, p->primelimbs, t, -op);
+    } else {
+      dp->d[0] = op;
+      memset(&dp->d[1], 0, sizeof(mp_limb_t) * (t - 1));
+    }
+    dp->flag = 2;
+  }
+}
+
+static void fp_to_mpz(mpz_ptr z, element_ptr e) {
+  dataptr dp = e->data;
+  if (!dp->flag) {
+    mpz_set_ui(z, 0);
+  } else {
+    fp_field_data_ptr p = e->field->data;
+    mpz_import(z, p->limbs, -1, sizeof(mp_limb_t), 0, 0, dp->d);
+  }
+}
+
+static void fp_set0(element_ptr e) {
+  dataptr dp = e->data;
+  dp->flag = 0;
+}
+
+static void fp_set1(element_ptr e) {
+  fp_field_data_ptr p = e->field->data;
+  dataptr dp = e->data;
+  dp->flag = 2;
+  memset(&dp->d[1], 0, p->bytes - sizeof(mp_limb_t));
+  dp->d[0] = 1;
+}
+
+static int fp_is1(element_ptr e) {
+  dataptr dp = e->data;
+  if (!dp->flag) return 0;
+  else {
+    fp_field_data_ptr p = e->field->data;
+    size_t i, t = p->limbs;
+    if (dp->d[0] != 1) return 0;
+    for (i = 1; i < t; i++) if (dp->d[i]) return 0;
+    return 1;
+  }
+}
+
+static int fp_is0(element_ptr e) {
+  dataptr dp = e->data;
+  return !dp->flag;
+}
+
+static size_t fp_out_str(FILE * stream, int base, element_ptr e) {
+  size_t result;
+  mpz_t z;
+  mpz_init(z);
+  fp_to_mpz(z, e);
+  result = mpz_out_str(stream, base, z);
+  mpz_clear(z);
+  return result;
+}
+
+static void fp_set(element_ptr c, element_ptr a) {
+  dataptr ad = a->data;
+  dataptr cd = c->data;
+  if (a == c) return;
+  if (!ad->flag) {
+    cd->flag = 0;
+  } else {
+    fp_field_data_ptr p = a->field->data;
+
+    //Assembly is faster here, but I don't want to stoop to that level.
+    //Instead of calling slower memcpy, wrap stuff so that GMP assembly
+    //gets called.
+    /*
+       memcpy(cd->d, ad->d, p->bytes);
+     */
+    mpz_t z1, z2;
+    z1->_mp_d = cd->d;
+    z2->_mp_d = ad->d;
+    z1->_mp_size = z1->_mp_alloc = z2->_mp_size = z2->_mp_alloc = p->limbs;
+    mpz_set(z1, z2);
+
+    cd->flag = 2;
+  }
+}
+
+static void fp_add(element_ptr c, element_ptr a, element_ptr b) {
+  dataptr ad = a->data, bd = b->data;
+
+  if (!ad->flag) {
+    fp_set(c, b);
+  } else if (!bd->flag) {
+    fp_set(c, a);
+  } else {
+    dataptr cd = c->data;
+    fp_field_data_ptr p = a->field->data;
+    const size_t t = p->limbs;
+    mp_limb_t carry;
+    carry = mpn_add_n(cd->d, ad->d, bd->d, t);
+
+    if (carry) {
+      //assumes result of following sub is not zero,
+      //i.e. modulus cannot be 2^(n * bits_per_limb)
+      mpn_sub_n(cd->d, cd->d, p->primelimbs, t);
+      cd->flag = 2;
+    } else {
+      int i = mpn_cmp(cd->d, p->primelimbs, t);
+      if (!i) {
+        cd->flag = 0;
+      } else {
+        cd->flag = 2;
+        if (i > 0) {
+          mpn_sub_n(cd->d, cd->d, p->primelimbs, t);
+        }
+      }
+    }
+  }
+}
+
+static void fp_double(element_ptr c, element_ptr a) {
+  dataptr ad = a->data, cd = c->data;
+  if (!ad->flag) {
+    cd->flag = 0;
+  } else {
+    fp_field_data_ptr p = c->field->data;
+    const size_t t = p->limbs;
+    if (mpn_lshift(cd->d, ad->d, t, 1)) {
+      cd->flag = 2;
+      //again, assumes result is not zero:
+      mpn_sub_n(cd->d, cd->d, p->primelimbs, t);
+    } else {
+      int i = mpn_cmp(cd->d, p->primelimbs, t);
+      if (!i) {
+        cd->flag = 0;
+      } else {
+        cd->flag = 2;
+        if (i > 0) {
+          mpn_sub_n(cd->d, cd->d, p->primelimbs, t);
+        }
+      }
+    }
+  }
+}
+
+static void fp_halve(element_ptr c, element_ptr a) {
+  dataptr ad = a->data, cd = c->data;
+  if (!ad->flag) {
+    cd->flag = 0;
+  } else {
+    fp_field_data_ptr p = c->field->data;
+    const size_t t = p->limbs;
+    int carry = 0;
+    mp_limb_t *alimb = ad->d;
+    mp_limb_t *climb = cd->d;
+    if (alimb[0] & 1) {
+      carry = mpn_add_n(climb, alimb, p->primelimbs, t);
+    } else fp_set(c, a);
+
+    mpn_rshift(climb, climb, t, 1);
+    if (carry) climb[t - 1] |= ((mp_limb_t) 1) << (sizeof(mp_limb_t) * 8 - 1);
+  }
+}
+
+static void fp_neg(element_ptr c, element_ptr a) {
+  dataptr ad = a->data, cd = c->data;
+  if (!ad->flag) cd->flag = 0;
+  else {
+    fp_field_data_ptr p = a->field->data;
+    mpn_sub_n(cd->d, p->primelimbs, ad->d, p->limbs);
+    cd->flag = 2;
+  }
+}
+
+static void fp_sub(element_ptr c, element_ptr a, element_ptr b) {
+  dataptr ad = a->data, bd = b->data;
+
+  if (!ad->flag) {
+    fp_neg(c, b);
+  } else if (!bd->flag) {
+    fp_set(c, a);
+  } else {
+    fp_field_data_ptr p = c->field->data;
+    size_t t = p->limbs;
+    dataptr cd = c->data;
+    int i = mpn_cmp(ad->d, bd->d, t);
+
+    if (i == 0) {
+      cd->flag = 0;
+    } else {
+      cd->flag = 2;
+      mpn_sub_n(cd->d, ad->d, bd->d, t);
+      if (i < 0) {
+        mpn_add_n(cd->d, cd->d, p->primelimbs, t);
+      }
+    }
+  }
+}
+
+static void fp_mul(element_ptr c, element_ptr a, element_ptr b) {
+  dataptr ad = a->data, bd = b->data;
+  dataptr cd = c->data;
+
+  if (!ad->flag || !bd->flag) {
+    cd->flag = 0;
+  } else {
+    fp_field_data_ptr p = c->field->data;
+    size_t t = p->limbs;
+    //mp_limb_t tmp[3 * t + 1];
+    //mp_limb_t *qp = &tmp[2 * t];
+    mp_limb_t tmp[2 * t];
+    mp_limb_t qp[t + 1];
+    //static mp_limb_t tmp[2 * 100];
+    //static mp_limb_t qp[100 + 1];
+
+    mpn_mul_n(tmp, ad->d, bd->d, t);
+
+    mpn_tdiv_qr(qp, cd->d, 0, tmp, 2 * t, p->primelimbs, t);
+    cd->flag = 2;
+  }
+}
+
+static void fp_square(element_ptr c, element_ptr a) {
+  const fp_field_data_ptr p = c->field->data;
+  mpz_t z1, z2;
+  size_t diff;
+  dataptr ad = a->data;
+  dataptr cd = c->data;
+
+  if (!ad->flag) {
+    cd->flag = 0;
+  } else {
+    cd->flag = 2;
+    z1->_mp_d = cd->d;
+    z1->_mp_size = z1->_mp_alloc = p->limbs;
+    if (c == a) {
+      mpz_powm_ui(z1, z1, 2, c->field->order);
+    } else {
+      z2->_mp_d = ad->d;
+      z2->_mp_size = z2->_mp_alloc = p->limbs;
+      mpz_powm_ui(z1, z2, 2, c->field->order);
+    }
+
+    diff = p->limbs - z1->_mp_size;
+    if (diff) memset(&z1->_mp_d[z1->_mp_size], 0, diff * sizeof(mp_limb_t));
+
+    //mpn_sqr_n() might make the code below faster than the code above
+    //but GMP doesn't expose this function
+    /*
+       const fp_field_data_ptr p = c->field->data;
+       const size_t t = p->limbs;
+       mp_limb_t tmp[2 * t];
+       mp_limb_t qp[t + 1];
+
+       mpn_mul_n(tmp, ad->d, ad->d, t);
+
+       mpn_tdiv_qr(qp, cd->d, 0, tmp, 2 * t, p->primelimbs, t);
+     */
+  }
+}
+
+static void fp_mul_si(element_ptr c, element_ptr a, signed long int op) {
+  dataptr ad = a->data;
+  dataptr cd = c->data;
+
+  if (!ad->flag || !op) {
+    cd->flag = 0;
+  } else {
+    cd->flag = 2;
+    fp_field_data_ptr p = a->field->data;
+    size_t t = p->limbs;
+    mp_limb_t tmp[t + 1];
+    mp_limb_t qp[2];
+
+    tmp[t] = mpn_mul_1(tmp, ad->d, t, labs(op));
+    mpn_tdiv_qr(qp, cd->d, 0, tmp, t + 1, p->primelimbs, t);
+    if (op < 0) {               //TODO: don't need to check c != 0 this time
+      fp_neg(c, c);
+    }
+  }
+}
+
+static void fp_pow_mpz(element_ptr c, element_ptr a, mpz_ptr op) {
+  dataptr ad = a->data;
+  dataptr cd = c->data;
+  if (!ad->flag) cd->flag = 0;
+  else {
+    mpz_t z;
+    mpz_init(z);
+    fp_to_mpz(z, a);
+    mpz_powm(z, z, op, a->field->order);
+    from_mpz(c, z);
+    mpz_clear(z);
+    cd->flag = 2;
+  }
+}
+
+static void fp_invert(element_ptr c, element_ptr a) {
+  //assumes a is invertible
+  dataptr cd = c->data;
+  mpz_t z;
+  mpz_init(z);
+  fp_to_mpz(z, a);
+  mpz_invert(z, z, a->field->order);
+  from_mpz(c, z);
+  mpz_clear(z);
+  cd->flag = 2;
+}
+
+static void fp_random(element_ptr a) {
+  dataptr ad = a->data;
+  mpz_t z;
+  mpz_init(z);
+  pbc_mpz_random(z, a->field->order);
+  if (mpz_sgn(z)) {
+    from_mpz(a, z);
+    ad->flag = 2;
+  } else {
+    ad->flag = 0;
+  }
+  mpz_clear(z);
+}
+
+static void fp_from_hash(element_ptr a, void *data, int len) {
+  mpz_t z;
+
+  mpz_init(z);
+  pbc_mpz_from_hash(z, a->field->order, data, len);
+  fp_set_mpz(a, z);
+  mpz_clear(z);
+}
+
+static int fp_cmp(element_ptr a, element_ptr b) {
+  dataptr ad = a->data, bd = b->data;
+  if (!ad->flag) {
+    return bd->flag;
+  } else {
+    fp_field_data_ptr p = a->field->data;
+    return mpn_cmp(ad->d, bd->d, p->limbs);
+    //return memcmp(ad->d, bd->d, p->limbs);
+  }
+}
+
+static int fp_sgn_odd(element_ptr a) {
+  dataptr ad = a->data;
+  if (!ad->flag) return 0;
+  return ad->d[0] & 1 ? 1 : -1;
+}
+
+static int fp_sgn_even(element_ptr a) {
+  fp_field_data_ptr p = a->field->data;
+  dataptr ad = a->data;
+  if (!ad->flag) return 0;
+  mp_limb_t sum[p->limbs];
+
+  int carry = mpn_add_n(sum, ad->d, ad->d, p->limbs);
+  if (carry) return 1;
+  return mpn_cmp(sum, p->primelimbs, p->limbs);
+}
+
+
+static int fp_is_sqr(element_ptr a) {
+  dataptr ad = a->data;
+  int res;
+  mpz_t z;
+  mpz_init(z);
+  //0 is a square
+  if (!ad->flag) return 1;
+  fp_to_mpz(z, a);
+  res = mpz_legendre(z, a->field->order) == 1;
+  mpz_clear(z);
+  return res;
+}
+
+static int fp_to_bytes(unsigned char *data, element_t a) {
+  dataptr ad = a->data;
+  int n = a->field->fixed_length_in_bytes;
+  if (!ad->flag) {
+    memset(data, 0, n);
+  } else {
+    mpz_t z;
+
+    mpz_init(z);
+    fp_to_mpz(z, a);
+    pbc_mpz_out_raw_n(data, n, z);
+    mpz_clear(z);
+  }
+  return n;
+}
+
+static int fp_from_bytes(element_t a, unsigned char *data) {
+  dataptr ad = a->data;
+  int n;
+  mpz_t z;
+
+  mpz_init(z);
+
+  n = a->field->fixed_length_in_bytes;
+  mpz_import(z, n, 1, 1, 1, 0, data);
+  if (!mpz_sgn(z)) ad->flag = 0;
+  else {
+    ad->flag = 2;
+    from_mpz(a, z);
+  }
+  mpz_clear(z);
+  return n;
+}
+
+static void fp_out_info(FILE* str, field_ptr f) {
+  element_fprintf(str, "GF(%Zd): zero flag + mpn", f->order);
+}
+
+static void fp_field_clear(field_t f) {
+  fp_field_data_ptr p = f->data;
+  pbc_free(p->primelimbs);
+  pbc_free(p);
+}
+
+void field_init_faster_fp(field_ptr f, mpz_t prime) {
+  PBC_ASSERT(!mpz_fits_ulong_p(prime), "modulus too small");
+  fp_field_data_ptr p;
+  field_init(f);
+  f->init = fp_init;
+  f->clear = fp_clear;
+  f->set_si = fp_set_si;
+  f->set_mpz = fp_set_mpz;
+  f->out_str = fp_out_str;
+  f->add = fp_add;
+  f->sub = fp_sub;
+  f->set = fp_set;
+  f->mul = fp_mul;
+  f->mul_si = fp_mul_si;
+  f->square = fp_square;
+  f->doub = fp_double;
+  f->halve = fp_halve;
+  f->pow_mpz = fp_pow_mpz;
+  f->neg = fp_neg;
+  f->cmp = fp_cmp;
+  f->sign = mpz_odd_p(prime) ? fp_sgn_odd : fp_sgn_even;
+  f->invert = fp_invert;
+  f->random = fp_random;
+  f->from_hash = fp_from_hash;
+  f->is1 = fp_is1;
+  f->is0 = fp_is0;
+  f->set0 = fp_set0;
+  f->set1 = fp_set1;
+  f->is_sqr = fp_is_sqr;
+  f->sqrt = element_tonelli;
+  f->field_clear = fp_field_clear;
+  f->to_bytes = fp_to_bytes;
+  f->from_bytes = fp_from_bytes;
+  f->to_mpz = fp_to_mpz;
+
+  f->out_info = fp_out_info;
+
+  p = f->data = pbc_malloc(sizeof(fp_field_data_t));
+  p->limbs = mpz_size(prime);
+  p->bytes = p->limbs * sizeof(mp_limb_t);
+  p->primelimbs = pbc_malloc(p->bytes);
+  mpz_export(p->primelimbs, &p->limbs, -1, sizeof(mp_limb_t), 0, 0, prime);
+
+  mpz_set(f->order, prime);
+  f->fixed_length_in_bytes = (mpz_sizeinbase(prime, 2) + 7) / 8;
+}
diff --git a/moon-abe/pbc-0.5.14/arith/fastfp.c b/moon-abe/pbc-0.5.14/arith/fastfp.c
new file mode 100644
index 00000000..13c6fb87
--- /dev/null
+++ b/moon-abe/pbc-0.5.14/arith/fastfp.c
@@ -0,0 +1,382 @@
+// Naive implementation of F_p.
+// Uses lowlevel GMP routines (mpn_* functions).
+//
+// Within an element_t, ''data'' field of element holds pointer to array of
+// mp_limb_t, which is allocated on init and freed on clear.
+// Its size is fixed and determined by the number of limbs in the modulus.
+// This simplifies code but is inefficient for storing values like 0 and 1.
+
+#include <stdarg.h>
+#include <stdio.h>
+#include <stdint.h> // for intptr_t
+#include <stdlib.h>
+#include <string.h>
+#include <gmp.h>
+#include "pbc_utils.h"
+#include "pbc_field.h"
+#include "pbc_random.h"
+#include "pbc_fp.h"
+#include "pbc_memory.h"
+
+struct fp_field_data_s {
+  size_t limbs;
+  size_t bytes;
+  mp_limb_t *primelimbs;
+};
+typedef struct fp_field_data_s fp_field_data_t[1];
+typedef struct fp_field_data_s *fp_field_data_ptr;
+
+static void fp_init(element_ptr e) {
+  fp_field_data_ptr p = e->field->data;
+  e->data = pbc_malloc(p->bytes);
+  memset(e->data, 0, p->bytes);
+  // e->data = pbc_calloc(sizeof(mp_limb_t), p->limbs);
+}
+
+static void fp_clear(element_ptr e) {
+  pbc_free(e->data);
+}
+
+static inline void from_mpz(element_ptr e, mpz_ptr z) {
+  fp_field_data_ptr p = e->field->data;
+  size_t count;
+  mpz_export(e->data, &count, -1, sizeof(mp_limb_t), 0, 0, z);
+  memset((void *) (((unsigned char *) e->data) + count * sizeof(mp_limb_t)), 0,
+         (p->limbs - count) * sizeof(mp_limb_t));
+}
+
+static void fp_set_mpz(element_ptr e, mpz_ptr z) {
+  mpz_t tmp;
+  mpz_init(tmp);
+  mpz_mod(tmp, z, e->field->order);
+  from_mpz(e, tmp);
+  mpz_clear(tmp);
+}
+
+static void fp_set_si(element_ptr e, signed long int op) {
+  const fp_field_data_ptr p = e->field->data;
+  const size_t t = p->limbs;
+  mp_limb_t *d = e->data;
+  if (op < 0) {
+    mpn_sub_1(d, p->primelimbs, t, -op);
+  } else {
+    d[0] = op;
+    memset(&d[1], 0, sizeof(mp_limb_t) * (t - 1));
+  }
+}
+
+static void fp_to_mpz(mpz_ptr z, element_ptr a) {
+  fp_field_data_ptr p = a->field->data;
+  mpz_import(z, p->limbs, -1, sizeof(mp_limb_t), 0, 0, a->data);
+}
+
+static void fp_set0(element_ptr e) {
+  fp_field_data_ptr p = e->field->data;
+  memset(e->data, 0, p->bytes);
+}
+
+static void fp_set1(element_ptr e) {
+  fp_field_data_ptr p = e->field->data;
+  mp_limb_t *d = e->data;
+  memset(&d[1], 0, p->bytes - sizeof(mp_limb_t));
+  d[0] = 1;
+}
+
+static int fp_is1(element_ptr e) {
+  fp_field_data_ptr p = e->field->data;
+  size_t i, t = p->limbs;
+  mp_limb_t *d = e->data;
+  if (d[0] != 1) return 0;
+  for (i = 1; i < t; i++) if (d[i]) return 0;
+  return 1;
+}
+
+static int fp_is0(element_ptr e) {
+  fp_field_data_ptr p = e->field->data;
+  size_t i, t = p->limbs;
+  mp_limb_t *d = e->data;
+  for (i = 0; i < t; i++) if (d[i]) return 0;
+  return 1;
+}
+
+static size_t fp_out_str(FILE * stream, int base, element_ptr e) {
+  size_t result;
+  mpz_t z;
+  mpz_init(z);
+  fp_to_mpz(z, e);
+  result = mpz_out_str(stream, base, z);
+  mpz_clear(z);
+  return result;
+}
+
+static void fp_add(element_ptr r, element_ptr a, element_ptr b) {
+  fp_field_data_ptr p = r->field->data;
+  const size_t t = p->limbs;
+  mp_limb_t carry;
+  carry = mpn_add_n(r->data, a->data, b->data, t);
+
+  if (carry || mpn_cmp(r->data, p->primelimbs, t) >= 0) {
+    mpn_sub_n(r->data, r->data, p->primelimbs, t);
+  }
+}
+
+static void fp_double(element_ptr r, element_ptr a) {
+  fp_field_data_ptr p = r->field->data;
+  const size_t t = p->limbs;
+  if (mpn_lshift(r->data, a->data, t, 1)
+      || mpn_cmp(r->data, p->primelimbs, t) >= 0) {
+    mpn_sub_n(r->data, r->data, p->primelimbs, t);
+  }
+}
+
+static void fp_set(element_ptr c, element_ptr a) {
+  fp_field_data_ptr p = a->field->data;
+  if (c == a) return;
+
+  // Assembly is faster here, but I don't want to stoop to that level.
+  // Instead of calling slower memcpy, wrap stuff so that GMP assembly
+  // gets called.
+  /*
+     memcpy(c->data, a->data, p->bytes);
+   */
+  mpz_t z1, z2;
+  z1->_mp_d = c->data;
+  z2->_mp_d = a->data;
+  z1->_mp_size = z1->_mp_alloc = z2->_mp_size = z2->_mp_alloc = p->limbs;
+  mpz_set(z1, z2);
+}
+
+static void fp_halve(element_ptr r, element_ptr a) {
+  fp_field_data_ptr p = r->field->data;
+  const size_t t = p->limbs;
+  int carry = 0;
+  mp_limb_t *alimb = a->data;
+  mp_limb_t *rlimb = r->data;
+  if (alimb[0] & 1) carry = mpn_add_n(rlimb, alimb, p->primelimbs, t);
+  else fp_set(r, a);
+
+  mpn_rshift(rlimb, rlimb, t, 1);
+  if (carry) rlimb[t - 1] |= ((mp_limb_t) 1) << (sizeof(mp_limb_t) * 8 - 1);
+}
+
+static void fp_sub(element_ptr r, element_ptr a, element_ptr b) {
+  fp_field_data_ptr p = r->field->data;
+  size_t t = p->limbs;
+  if (mpn_sub_n(r->data, a->data, b->data, t)) {
+    mpn_add_n(r->data, r->data, p->primelimbs, t);
+  }
+}
+
+static void fp_mul(element_ptr c, element_ptr a, element_ptr b) {
+  fp_field_data_ptr p = c->field->data;
+  size_t t = p->limbs;
+  //mp_limb_t tmp[3 * t + 1];
+  //mp_limb_t *qp = &tmp[2 * t];
+  mp_limb_t tmp[2 * t];
+  mp_limb_t qp[t + 1];
+  //static mp_limb_t tmp[2 * 100];
+  //static mp_limb_t qp[100 + 1];
+
+  mpn_mul_n(tmp, a->data, b->data, t);
+
+  mpn_tdiv_qr(qp, c->data, 0, tmp, 2 * t, p->primelimbs, t);
+}
+
+static void fp_square(element_ptr c, element_ptr a) {
+  const fp_field_data_ptr r = c->field->data;
+  mpz_t z1, z2;
+  size_t diff;
+
+  z1->_mp_d = c->data;
+  z1->_mp_size = z1->_mp_alloc = r->limbs;
+  if (c == a) {
+    mpz_powm_ui(z1, z1, 2, c->field->order);
+  } else {
+    z2->_mp_d = a->data;
+    z2->_mp_size = z2->_mp_alloc = r->limbs;
+    mpz_powm_ui(z1, z2, 2, c->field->order);
+  }
+
+  diff = r->limbs - z1->_mp_size;
+  if (diff) memset(&z1->_mp_d[z1->_mp_size], 0, diff * sizeof(mp_limb_t));
+
+  //mpn_sqr_n() might make the code below faster than the code above
+  //but GMP doesn't expose this function
+  /*
+     const fp_field_data_ptr r = c->field->data;
+     const size_t t = r->limbs;
+     mp_limb_t tmp[2 * t];
+     mp_limb_t qp[t + 1];
+
+     mpn_mul_n(tmp, a->data, a->data, t);
+
+     mpn_tdiv_qr(qp, c->data, 0, tmp, 2 * t, r->primelimbs, t);
+   */
+}
+
+static void fp_neg(element_ptr n, element_ptr a) {
+  if (fp_is0(a)) {
+    fp_set0(n);
+  } else {
+    fp_field_data_ptr p = a->field->data;
+    mpn_sub_n(n->data, p->primelimbs, a->data, p->limbs);
+  }
+}
+
+static void fp_mul_si(element_ptr e, element_ptr a, signed long int op) {
+  fp_field_data_ptr p = e->field->data;
+  size_t t = p->limbs;
+  mp_limb_t tmp[t + 1];
+  mp_limb_t qp[2];
+
+  tmp[t] = mpn_mul_1(tmp, a->data, t, labs(op));
+  mpn_tdiv_qr(qp, e->data, 0, tmp, t + 1, p->primelimbs, t);
+  if (op < 0) {
+    fp_neg(e, e);
+  }
+}
+
+static void fp_pow_mpz(element_ptr c, element_ptr a, mpz_ptr op) {
+  mpz_t z;
+  mpz_init(z);
+  fp_to_mpz(z, a);
+  mpz_powm(z, z, op, c->field->order);
+  from_mpz(c, z);
+  mpz_clear(z);
+}
+
+static void fp_invert(element_ptr e, element_ptr a) {
+  mpz_t z;
+  mpz_init(z);
+  fp_to_mpz(z, a);
+  mpz_invert(z, z, e->field->order);
+  from_mpz(e, z);
+  mpz_clear(z);
+}
+
+static void fp_random(element_ptr a) {
+  mpz_t z;
+  mpz_init(z);
+  pbc_mpz_random(z, a->field->order);
+  from_mpz(a, z);
+  mpz_clear(z);
+}
+
+static void fp_from_hash(element_ptr a, void *data, int len) {
+  mpz_t z;
+
+  mpz_init(z);
+  pbc_mpz_from_hash(z, a->field->order, data, len);
+  fp_set_mpz(a, z);
+  mpz_clear(z);
+}
+
+static int fp_cmp(element_ptr a, element_ptr b) {
+  fp_field_data_ptr p = a->field->data;
+  return mpn_cmp(a->data, b->data, p->limbs);
+  //return memcmp(a->data, b->data, p->limbs);
+}
+
+static int fp_sgn_odd(element_ptr a) {
+  if (fp_is0(a)) return 0;
+  mp_limb_t *lp = a->data;
+  return lp[0] & 1 ? 1 : -1;
+}
+
+static int fp_sgn_even(element_ptr a) {
+  fp_field_data_ptr p = a->field->data;
+  if (fp_is0(a)) return 0;
+  mp_limb_t sum[p->limbs];
+
+  int carry = mpn_add_n(sum, a->data, a->data, p->limbs);
+  if (carry) return 1;
+  return mpn_cmp(sum, p->primelimbs, p->limbs);
+}
+
+static int fp_is_sqr(element_ptr a) {
+  int res;
+  mpz_t z;
+  mpz_init(z);
+  //0 is a square
+  if (fp_is0(a)) return 1;
+  fp_to_mpz(z, a);
+  res = mpz_legendre(z, a->field->order) == 1;
+  mpz_clear(z);
+  return res;
+}
+
+static int fp_to_bytes(unsigned char *data, element_t e) {
+  mpz_t z;
+  int n;
+
+  mpz_init(z);
+  fp_to_mpz(z, e);
+  n = e->field->fixed_length_in_bytes;
+  pbc_mpz_out_raw_n(data, n, z);
+  mpz_clear(z);
+  return n;
+}
+
+static int fp_from_bytes(element_t e, unsigned char *data) {
+  int n;
+  mpz_t z;
+
+  mpz_init(z);
+
+  n = e->field->fixed_length_in_bytes;
+  mpz_import(z, n, 1, 1, 1, 0, data);
+  fp_set_mpz(e, z);
+  mpz_clear(z);
+  return n;
+}
+
+static void fp_field_clear(field_t f) {
+  fp_field_data_ptr p = f->data;
+  pbc_free(p->primelimbs);
+  pbc_free(p);
+}
+
+void field_init_fast_fp(field_ptr f, mpz_t prime) {
+  PBC_ASSERT(!mpz_fits_ulong_p(prime), "modulus too small");
+  fp_field_data_ptr p;
+  field_init(f);
+  f->init = fp_init;
+  f->clear = fp_clear;
+  f->set_si = fp_set_si;
+  f->set_mpz = fp_set_mpz;
+  f->out_str = fp_out_str;
+  f->add = fp_add;
+  f->sub = fp_sub;
+  f->set = fp_set;
+  f->mul = fp_mul;
+  f->mul_si = fp_mul_si;
+  f->square = fp_square;
+  f->doub = fp_double;
+  f->halve = fp_halve;
+  f->pow_mpz = fp_pow_mpz;
+  f->neg = fp_neg;
+  f->cmp = fp_cmp;
+  f->sign = mpz_odd_p(prime) ? fp_sgn_odd : fp_sgn_even;
+  f->invert = fp_invert;
+  f->random = fp_random;
+  f->from_hash = fp_from_hash;
+  f->is1 = fp_is1;
+  f->is0 = fp_is0;
+  f->set0 = fp_set0;
+  f->set1 = fp_set1;
+  f->is_sqr = fp_is_sqr;
+  f->sqrt = element_tonelli;
+  f->field_clear = fp_field_clear;
+  f->to_bytes = fp_to_bytes;
+  f->from_bytes = fp_from_bytes;
+  f->to_mpz = fp_to_mpz;
+
+  p = f->data = pbc_malloc(sizeof(fp_field_data_t));
+  p->limbs = mpz_size(prime);
+  p->bytes = p->limbs * sizeof(mp_limb_t);
+  p->primelimbs = pbc_malloc(p->bytes);
+  mpz_export(p->primelimbs, &p->limbs, -1, sizeof(mp_limb_t), 0, 0, prime);
+
+  mpz_set(f->order, prime);
+  f->fixed_length_in_bytes = (mpz_sizeinbase(prime, 2) + 7) / 8;
+}
diff --git a/moon-abe/pbc-0.5.14/arith/field.c b/moon-abe/pbc-0.5.14/arith/field.c
new file mode 100644
index 00000000..af94e37f
--- /dev/null
+++ b/moon-abe/pbc-0.5.14/arith/field.c
@@ -0,0 +1,889 @@
+#include <ctype.h>
+#include <stdarg.h>
+#include <stdint.h> // for intptr_t
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h> // for memcmp()
+#include <gmp.h>
+#include "pbc_utils.h"
+#include "pbc_field.h"
+#include "pbc_multiz.h"
+#include "pbc_memory.h"
+
+// returns recommended window size.  n is exponent.
+static int optimal_pow_window_size(mpz_ptr n) {
+  int exp_bits;
+
+  exp_bits = mpz_sizeinbase(n, 2);
+
+  // try to minimize 2^k + n/(k+1).
+  return exp_bits > 9065 ? 8 :
+         exp_bits > 3529 ? 7 :
+         exp_bits > 1324 ? 6 :
+         exp_bits > 474 ? 5 :
+         exp_bits > 157 ? 4 :
+         exp_bits > 47 ? 3 :
+         2;
+}
+
+/* builds k-bit lookup window for base a */
+static element_t *build_pow_window(element_ptr a, int k) {
+  int s;
+  int lookup_size;
+  element_t *lookup;
+
+  if (k < 1) return NULL;  // no window
+
+  /* build 2^k lookup table.  lookup[i] = x^i. */
+  /* TODO: a more careful word-finding algorithm would allow
+   *       us to avoid calculating even lookup entries > 2
+   */
+  lookup_size = 1 << k;
+  lookup = pbc_malloc(lookup_size * sizeof(element_t));
+  element_init(lookup[0], a->field);
+  element_set1(lookup[0]);
+  for (s = 1; s < lookup_size; s++) {
+    element_init(lookup[s], a->field);
+    element_mul(lookup[s], lookup[s - 1], a);
+  }
+
+  return lookup;
+}
+
+static void clear_pow_window(int k, element_t * lookup) {
+  int s;
+  int lookup_size = 1 << k;
+
+  for (s = 0; s < lookup_size; s++) {
+    element_clear(lookup[s]);
+  }
+  pbc_free(lookup);
+}
+
+/*
+ * left-to-right exponentiation with k-bit window.
+ * NB. must have k >= 1.
+ */
+static void element_pow_wind(element_ptr x, mpz_ptr n,
+                             int k, element_t * a_lookup) {
+  int s;
+  int bit;
+
+  int inword;                   // boolean: currently reading word?
+  int word = 0;                 // the word to look up. 0<word<base
+  int wbits = 0;                // # of bits so far in word. wbits<=k.
+
+  element_t result;
+
+  // early abort if raising to power 0
+  if (!mpz_sgn(n)) {
+    element_set1(x);
+    return;
+  }
+
+  element_init(result, x->field);
+  element_set1(result);
+
+  for (inword = 0, s = mpz_sizeinbase(n, 2) - 1; s >= 0; s--) {
+    element_square(result, result);
+    bit = mpz_tstbit(n, s);
+
+    if (!inword && !bit)
+      continue;                 // keep scanning.  note continue.
+
+    if (!inword) {              // was scanning, just found word
+      inword = 1;               // so, start new word
+      word = 1;
+      wbits = 1;
+    } else {
+      word = (word << 1) + bit;
+      wbits++;                  // continue word
+    }
+
+    if (wbits == k || s == 0) {
+      element_mul(result, result, a_lookup[word]);
+      inword = 0;
+    }
+  }
+
+  element_set(x, result);
+  element_clear(result);
+}
+
+static void generic_pow_mpz(element_ptr x, element_ptr a, mpz_ptr n) {
+  int k;
+  element_t *a_lookup;
+
+  if (mpz_is0(n)) {
+    element_set1(x);
+    return;
+  }
+
+  k = optimal_pow_window_size(n);
+  a_lookup = build_pow_window(a, k);
+  element_pow_wind(x, n, k, a_lookup);
+  clear_pow_window(k, a_lookup);
+}
+
+/* TODO: Allow fields to choose this exponentiation routine so we can compare.
+static void naive_generic_pow_mpz(element_ptr x, element_ptr a, mpz_ptr n) {
+  int s;
+
+  element_t result;
+
+  if (mpz_is0(n)) {
+    element_set1(x);
+    return;
+  }
+
+  element_init(result, x->field);
+  element_set1(result);
+
+  for (s = mpz_sizeinbase(n, 2) - 1; s >= 0; s--) {
+    element_square(result, result);
+    if (mpz_tstbit(n, s)) {
+      element_mul(result, result, a);
+    }
+  }
+  element_set(x, result);
+  element_clear(result);
+}
+*/
+
+void element_pow2_mpz(element_ptr x, element_ptr a1, mpz_ptr n1,
+                      element_ptr a2, mpz_ptr n2) {
+  int s, s1, s2;
+  int b1, b2;
+
+  element_t result, a1a2;
+
+  if (mpz_is0(n1) && mpz_is0(n2)) {
+    element_set1(x);
+    return;
+  }
+
+  element_init(result, x->field);
+  element_set1(result);
+
+  element_init(a1a2, x->field);
+  element_mul(a1a2, a1, a2);
+
+  s1 = mpz_sizeinbase(n1, 2) - 1;
+  s2 = mpz_sizeinbase(n2, 2) - 1;
+  for (s = (s1 > s2) ? s1 : s2; s >= 0; s--) {
+    element_mul(result, result, result);
+    b1 = mpz_tstbit(n1, s);
+    b2 = mpz_tstbit(n2, s);
+    if (b1 && b2) {
+      element_mul(result, result, a1a2);
+    } else if (b1) {
+      element_mul(result, result, a1);
+    } else if (b2) {
+      element_mul(result, result, a2);
+    }
+  }
+
+  element_set(x, result);
+  element_clear(result);
+  element_clear(a1a2);
+}
+
+void element_pow3_mpz(element_ptr x, element_ptr a1, mpz_ptr n1,
+                      element_ptr a2, mpz_ptr n2,
+                      element_ptr a3, mpz_ptr n3) {
+  int s, s1, s2, s3;
+  int b;
+  int i;
+
+  element_t result;
+  element_t lookup[8];
+
+  if (mpz_is0(n1) && mpz_is0(n2) && mpz_is0(n3)) {
+    element_set1(x);
+    return;
+  }
+
+  element_init(result, x->field);
+  element_set1(result);
+
+  for (i = 0; i < 8; i++)
+    element_init(lookup[i], x->field);
+
+  // build lookup table.
+  element_set1(lookup[0]);
+  element_set(lookup[1], a1);
+  element_set(lookup[2], a2);
+  element_set(lookup[4], a3);
+  element_mul(lookup[3], a1, a2);
+  element_mul(lookup[5], a1, a3);
+  element_mul(lookup[6], a2, a3);
+  element_mul(lookup[7], lookup[6], a1);
+
+  // calculate largest exponent bitsize
+  s1 = mpz_sizeinbase(n1, 2) - 1;
+  s2 = mpz_sizeinbase(n2, 2) - 1;
+  s3 = mpz_sizeinbase(n3, 2) - 1;
+  s = (s1 > s2) ? ((s1 > s3) ? s1 : s3)
+      : ((s2 > s3) ? s2 : s3);
+
+  for (; s >= 0; s--) {
+    element_mul(result, result, result);
+    b = (mpz_tstbit(n1, s))
+        + (mpz_tstbit(n2, s) << 1)
+        + (mpz_tstbit(n3, s) << 2);
+    element_mul(result, result, lookup[b]);
+  }
+
+  element_set(x, result);
+  element_clear(result);
+  for (i = 0; i < 8; i++)
+    element_clear(lookup[i]);
+}
+
+struct element_base_table {
+  int k;
+  int bits;
+  int num_lookups;
+  element_t **table;
+};
+
+/* build k-bit base table for n-bit exponentiation w/ base a */
+static void *element_build_base_table(element_ptr a, int bits, int k) {
+  struct element_base_table *base_table;
+  element_t multiplier;
+  int i, j;
+  int lookup_size;
+
+  element_t *lookup;
+
+  // pbc_info("building %d bits %d k", bits, k);
+  lookup_size = 1 << k;
+
+  base_table = pbc_malloc(sizeof(struct element_base_table));
+  base_table->num_lookups = bits / k + 1;
+  base_table->k = k;
+  base_table->bits = bits;
+  base_table->table =
+      pbc_malloc(base_table->num_lookups * sizeof(element_t *));
+
+  element_init(multiplier, a->field);
+  element_set(multiplier, a);
+
+  for (i = 0; i < base_table->num_lookups; i++) {
+    lookup = pbc_malloc(lookup_size * sizeof(element_t));
+    element_init(lookup[0], a->field);
+    element_set1(lookup[0]);
+    for (j = 1; j < lookup_size; j++) {
+      element_init(lookup[j], a->field);
+      element_mul(lookup[j], multiplier, lookup[j - 1]);
+    }
+    element_mul(multiplier, multiplier, lookup[lookup_size - 1]);
+    base_table->table[i] = lookup;
+  }
+
+  element_clear(multiplier);
+  return base_table;
+}
+
+/*
+ * exponentiation using aggressive base lookup table
+ * must have k >= 1.
+ */
+static void element_pow_base_table(element_ptr x, mpz_ptr power,
+                                   struct element_base_table *base_table) {
+  int word;                     /* the word to look up. 0<word<base */
+  int row, s;                   /* row and col in base table */
+  int num_lookups;
+
+  element_t result;
+  mpz_t n;
+  mpz_init_set(n, power);
+
+  // Early abort if raising to power 0.
+  if (!mpz_sgn(n)) {
+    element_set1(x);
+    return;
+  }
+  if (mpz_cmp(n, x->field->order) > 0) {
+    mpz_mod(n, n, x->field->order);
+  }
+
+  element_init(result, x->field);
+  element_set1(result);
+
+  num_lookups = mpz_sizeinbase(n, 2) / base_table->k + 1;
+
+  for (row = 0; row < num_lookups; row++) {
+    word = 0;
+    for (s = 0; s < base_table->k; s++) {
+      word |= mpz_tstbit(n, base_table->k * row + s) << s;
+    }
+    if (word > 0) {
+      element_mul(result, result, base_table->table[row][word]);
+    }
+  }
+
+  element_set(x, result);
+  element_clear(result);
+  mpz_clear(n);
+}
+
+static void default_element_pp_init(element_pp_t p, element_t in) {
+  p->data =
+      element_build_base_table(in, mpz_sizeinbase(in->field->order, 2), 5);
+}
+
+static void default_element_pp_pow(element_t out, mpz_ptr power, element_pp_t p) {
+  element_pow_base_table(out, power, p->data);
+}
+
+static void default_element_pp_clear(element_pp_t p) {
+  struct element_base_table *base_table = p->data;
+  int lookup_size = 1 << base_table->k;
+  element_t *lookup;
+  int i, j;
+
+  element_t **epp = base_table->table;
+
+  for (i = 0; i < base_table->num_lookups; i++) {
+    lookup = epp[i];
+    for (j = 0; j < lookup_size; j++) {
+      element_clear(lookup[j]);
+    }
+    pbc_free(lookup);
+  }
+  pbc_free(epp);
+
+  pbc_free(base_table);
+}
+
+void field_set_nqr(field_ptr f, element_t nqr) {
+  if (!f->nqr) {
+    f->nqr = pbc_malloc(sizeof(element_t));
+    element_init(f->nqr, f);
+  }
+  element_set(f->nqr, nqr);
+}
+
+void field_gen_nqr(field_ptr f) {
+  f->nqr = pbc_malloc(sizeof(element_t));
+  element_init(f->nqr, f);
+  do {
+    element_random(f->nqr);
+  } while (element_is_sqr(f->nqr));
+}
+
+element_ptr field_get_nqr(field_ptr f) {
+  if (!f->nqr) field_gen_nqr(f);
+  return f->nqr;
+}
+
+static void generic_square(element_ptr r, element_ptr a) {
+  element_mul(r, a, a);
+}
+static void generic_mul_mpz(element_ptr r, element_ptr a, mpz_ptr z) {
+  element_t e0;
+  element_init(e0, r->field);
+  element_set_mpz(e0, z);
+  element_mul(r, a, e0);
+  element_clear(e0);
+}
+
+static void generic_mul_si(element_ptr r, element_ptr a, signed long int n) {
+  element_t e0;
+  element_init(e0, r->field);
+  element_set_si(e0, n);
+  element_mul(r, a, e0);
+  element_clear(e0);
+}
+
+static void generic_double(element_ptr r, element_ptr a) {
+  element_add(r, a, a);
+}
+
+static void generic_halve(element_ptr r, element_ptr a) {
+  element_t e0;
+  element_init(e0, r->field);
+  element_set_si(e0, 2);
+  element_invert(e0, e0);
+  element_mul(r, a, e0);
+  element_clear(e0);
+}
+
+static void zero_to_mpz(mpz_t z, element_ptr a) {
+  UNUSED_VAR(a);
+  mpz_set_ui(z, 0);
+}
+
+static void zero_set_mpz(element_ptr a, mpz_t z) {
+  UNUSED_VAR(z);
+  element_set0(a);
+}
+
+static void zero_random(element_ptr a) {
+  element_set0(a);
+}
+
+static void generic_set_si(element_ptr a, long int si) {
+  mpz_t z;
+  mpz_init(z);
+  mpz_set_si(z, si);
+  element_set_mpz(a, z);
+  mpz_clear(z);
+}
+
+static void generic_set_multiz(element_ptr a, multiz m) {
+  mpz_t z;
+  mpz_init(z);
+  multiz_to_mpz(z, m);
+  element_set_mpz(a, z);
+  mpz_clear(z);
+}
+
+static void generic_sub(element_ptr c, element_ptr a, element_ptr b) {
+  if (c != a) {
+    element_neg(c, b);
+    element_add(c, c, a);
+  } else {
+    element_t tmp;
+    element_init(tmp, a->field);
+    element_neg(tmp, b);
+    element_add(c, tmp, a);
+    element_clear(tmp);
+  }
+}
+
+static void generic_div(element_ptr c, element_ptr a, element_ptr b) {
+  if (c != a) {
+    element_invert(c, b);
+    element_mul(c, c, a);
+  } else {
+    element_t tmp;
+    element_init(tmp, a->field);
+    element_invert(tmp, b);
+    element_mul(c, tmp, a);
+    element_clear(tmp);
+  }
+}
+
+static void generic_add_ui(element_ptr c, element_ptr a,
+                           unsigned long int b) {
+  element_t e;
+  mpz_t z;
+  element_init(e, c->field);
+  mpz_init(z);
+  mpz_set_ui(z, b);
+  element_set_mpz(e, z);
+  element_add(c, a, e);
+  mpz_clear(z);
+  element_clear(e);
+}
+
+static int generic_cmp(element_ptr a, element_ptr b) {
+  int result;
+  unsigned char *buf1, *buf2;
+  int len;
+  if (a == b) return 0;
+  len = element_length_in_bytes(a);
+  if (len != element_length_in_bytes(b)) return 1;
+  buf1 = pbc_malloc(len);
+  buf2 = pbc_malloc(len);
+  element_to_bytes(buf1, a);
+  element_to_bytes(buf2, b);
+  result = memcmp(buf1, buf2, len);
+  pbc_free(buf1);
+  pbc_free(buf2);
+  return result;
+}
+
+static int generic_is0(element_ptr a) {
+  int result;
+  element_t b;
+  element_init(b, a->field);
+  result = !element_cmp(a, b); // element_cmp returns 0 if 'a' and 'b' are the same, nonzero otherwise. generic_is0 returns true if 'a' is 0.
+  element_clear(b);
+  return result;
+}
+
+static int generic_is1(element_ptr a) {
+  int result;
+  element_t b;
+  element_init(b, a->field);
+  element_set1(b);
+  result = !element_cmp(a, b); // element_cmp returns 0 if 'a' and 'b' are the same, nonzero otherwise. generic_is1 returns true if 'a' is 1.
+  element_clear(b);
+  return result;
+}
+
+static void generic_out_info(FILE * out, field_ptr f) {
+  element_fprintf(out, "unknown field %p, order = %Zd", f, f->order);
+}
+
+static int generic_item_count(element_ptr e) {
+  UNUSED_VAR(e);
+  return 0;
+}
+
+static element_ptr generic_item(element_ptr e, int i) {
+  UNUSED_VAR(e);
+  UNUSED_VAR(i);
+  return NULL;
+}
+
+static element_ptr generic_get_x(element_ptr e) {
+  return element_item(e, 0);
+}
+
+static element_ptr generic_get_y(element_ptr e) {
+  return element_item(e, 1);
+}
+
+static int default_element_snprint(char *s, size_t n, element_t e) {
+  UNUSED_VAR(e);
+  if (n == 1) {
+    s[0] = '0';
+  } else if (n >= 2) {
+    s[0] = '?';
+    s[1] = '\0';
+  }
+  return 1;
+}
+
+static int default_element_set_str(element_t e, const char *s, int base) {
+  UNUSED_VAR(s);
+  UNUSED_VAR(base);
+  element_set0(e);
+  return 0;
+}
+
+static void warn_field_clear(field_ptr f) {
+  pbc_warn("field %p has no clear function", f);
+}
+
+void field_out_info(FILE* out, field_ptr f) {
+  f->out_info(out, f);
+}
+
+void field_init(field_ptr f) {
+  // should be called by each field_init_*
+  f->nqr = NULL;
+  mpz_init(f->order);
+
+  // this should later be set
+  f->field_clear = warn_field_clear;
+
+  // and this to something more helpful
+  f->out_info = generic_out_info;
+
+  // many of these can usually be optimized for particular fields
+  // provided for developer's convenience
+  f->halve = generic_halve;
+  f->doub = generic_double;
+  f->square = generic_square;
+  f->mul_mpz = generic_mul_mpz;
+  f->mul_si = generic_mul_si;
+  f->cmp = generic_cmp;
+  f->sub = generic_sub;
+  f->div = generic_div;
+  f->add_ui = generic_add_ui;
+
+  // default: converts all elements to integer 0
+  // reads all integers as 0
+  // random always outputs 0
+  f->to_mpz = zero_to_mpz;
+  f->set_mpz = zero_set_mpz;
+  f->set_multiz = generic_set_multiz;
+  f->random = zero_random;
+  f->set_si = generic_set_si;
+  f->is1 = generic_is1;
+  f->is0 = generic_is0;
+
+  // By default, an element has no components.
+  f->item_count = generic_item_count;
+  f->item = generic_item;
+  f->get_x = generic_get_x;
+  f->get_y = generic_get_y;
+
+  // these are fast, thanks to Hovav
+  f->pow_mpz = generic_pow_mpz;
+  f->pp_init = default_element_pp_init;
+  f->pp_clear = default_element_pp_clear;
+  f->pp_pow = default_element_pp_pow;
+
+  f->snprint = default_element_snprint;
+  f->set_str = default_element_set_str;
+  f->pairing = NULL;
+}
+
+void field_clear(field_ptr f) {
+  if (f->nqr) {
+    element_clear(f->nqr);
+    pbc_free(f->nqr);
+  }
+  mpz_clear(f->order);
+  f->field_clear(f);
+}
+
+void pbc_mpz_out_raw_n(unsigned char *data, int n, mpz_t z) {
+  size_t count;
+  if (mpz_sgn(z)) {
+    count = (mpz_sizeinbase(z, 2) + 7) / 8;
+    mpz_export(&data[n - count], NULL, 1, 1, 1, 0, z);
+    memset(data, 0, n - count);
+  } else {
+    memset(data, 0, n);
+  }
+}
+
+//for short hashes H, do
+//  buf = H || 0 || H || 1 || H || ...
+//before calling mpz_import
+void pbc_mpz_from_hash(mpz_t z, mpz_t limit,
+                       unsigned char *data, unsigned int len) {
+  size_t i = 0, n, count = (mpz_sizeinbase(limit, 2) + 7) / 8;
+  unsigned char buf[count];
+  unsigned char counter = 0;
+  int done = 0;
+  for (;;) {
+    if (len >= count - i) {
+      n = count - i;
+      done = 1;
+    } else n = len;
+    memcpy(buf + i, data, n);
+    i += n;
+    if (done) break;
+    buf[i] = counter;
+    counter++;
+    i++;
+    if (i == count) break;
+  }
+  PBC_ASSERT(i == count, "did not read whole buffer");
+  mpz_import(z, count, 1, 1, 1, 0, buf);
+  while (mpz_cmp(z, limit) > 0) {
+    mpz_tdiv_q_2exp(z, z, 1);
+  }
+}
+
+// Square root algorithm for Fp.
+// TODO: What happens if this is run on other kinds of fields?
+void element_tonelli(element_ptr x, element_ptr a) {
+  int s;
+  int i;
+  mpz_t e;
+  mpz_t t, t0;
+  element_t ginv, e0;
+  element_ptr nqr;
+
+  mpz_init(t);
+  mpz_init(e);
+  mpz_init(t0);
+  element_init(ginv, a->field);
+  element_init(e0, a->field);
+  nqr = field_get_nqr(a->field);
+
+  element_invert(ginv, nqr);
+
+  //let q be the order of the field
+  //q - 1 = 2^s t, t odd
+  mpz_sub_ui(t, a->field->order, 1);
+  s = mpz_scan1(t, 0);
+  mpz_tdiv_q_2exp(t, t, s);
+  mpz_set_ui(e, 0);
+  for (i = 2; i <= s; i++) {
+    mpz_sub_ui(t0, a->field->order, 1);
+    mpz_tdiv_q_2exp(t0, t0, i);
+    element_pow_mpz(e0, ginv, e);
+    element_mul(e0, e0, a);
+    element_pow_mpz(e0, e0, t0);
+    if (!element_is1(e0)) mpz_setbit(e, i - 1);
+  }
+  element_pow_mpz(e0, ginv, e);
+  element_mul(e0, e0, a);
+  mpz_add_ui(t, t, 1);
+  mpz_tdiv_q_2exp(t, t, 1);
+  mpz_tdiv_q_2exp(e, e, 1);
+
+  // (suggested by Hovav Shacham) replace next three lines with
+  //   element_pow2_mpz(x, e0, t, nqr, e);
+  // once sliding windows are implemented for pow2.
+  element_pow_mpz(e0, e0, t);
+  element_pow_mpz(x, nqr, e);
+  element_mul(x, x, e0);
+
+  mpz_clear(t);
+  mpz_clear(e);
+  mpz_clear(t0);
+  element_clear(ginv);
+  element_clear(e0);
+}
+
+// Like mpz_set_str except returns number of bytes read and allows trailing
+// junk. This simplifies code for parsing elements like "[123, 456]".
+// TODO: Handle 0x, 0X and 0 conventions for hexadecimal and octal.
+int pbc_mpz_set_str(mpz_t z, const char *s, int base) {
+  int b, i = 0;
+  mpz_set_ui(z, 0);
+  if (!base) b = 10;
+  else if (base < 2 || base > 36) return 0;
+  else b = base;
+
+  for (;;) {
+    int j;
+    char c = s[i];
+    if (!c) break;
+    if (isspace(c)) {
+      i++;
+      continue;
+    }
+    if (isdigit(c)) {
+      j = c - '0';
+    } else if (c >= 'A' && c <= 'Z') {
+      j = c - 'A';
+    } else if (c >= 'a' && c <= 'z') {
+      j = c - 'a';
+    } else break;
+
+    if (j >= b) break;
+
+    mpz_mul_ui(z, z, b);
+    mpz_add_ui(z, z, j);
+    i++;
+  }
+  return i;
+}
+
+// Divides `n` with primes up to `limit`. For each factor found,
+// call `fun`. If the callback returns nonzero, then aborts and returns 1.
+// Otherwise returns 0.
+int pbc_trial_divide(int (*fun)(mpz_t factor,
+                                unsigned int multiplicity,
+                                void *scope_ptr),
+                     void *scope_ptr,
+                     mpz_t n,
+                     mpz_ptr limit) {
+  mpz_t p, m;
+  mpz_t fac;
+  unsigned int mul;
+
+  mpz_init(fac);
+  mpz_init(p);
+  mpz_init(m);
+  mpz_set(m ,n);
+  mpz_set_ui(p, 2);
+
+  while (mpz_cmp_ui(m, 1)) {
+    if (mpz_probab_prime_p(m, 10)) {
+      mpz_set(p, m);
+    }
+    if (limit && mpz_cmp(p, limit) > 0) {
+      mpz_set(p, m);
+    }
+    if (mpz_divisible_p(m, p)) {
+      mul = 0;
+      mpz_set(fac, p);
+      do {
+        mpz_divexact(m, m, p);
+        mul++;
+      } while (mpz_divisible_p(m, p));
+      if (fun(fac, mul, scope_ptr)) {
+        mpz_clear(fac);
+        mpz_clear(m);
+        mpz_clear(p);
+        return 1;
+      }
+    }
+    mpz_nextprime(p, p);
+  }
+
+  mpz_clear(fac);
+  mpz_clear(m);
+  mpz_clear(p);
+  return 0;
+}
+
+// For each digit of 'n', call fun(). If it returns 1, then return 1 and
+// abort. Otherwise return 0.
+int pbc_mpz_trickle(int (*fun)(char), int base, mpz_t n) {
+  // TODO: Support different bases.
+  if (!base) base = 10;
+  if (base < 2 || base > 10) {
+    pbc_warn("only bases 2 to 10 supported");
+    return 1;
+  }
+  mpz_t d, z, q;
+  mpz_init(d);
+  mpz_init(z);
+  mpz_init(q);
+  mpz_set(z, n);
+  int res;
+  int len;
+  mpz_ui_pow_ui(d, base, len = mpz_sizeinbase(z, base));
+  if (mpz_cmp(d, z) > 0) {
+    len--;
+    mpz_divexact_ui(d, d, base);
+  }
+  while (mpz_cmp_ui(z, base) >= 0) {
+    mpz_fdiv_qr(q, z, z, d);
+    res = fun('0' + mpz_get_ui(q));
+    if (res) goto clean;
+    mpz_divexact_ui(d, d, base);
+    len--;
+  }
+  while (len) {
+    res = fun('0');
+    if (res) goto clean;
+    len--;
+  }
+  res = fun('0' + mpz_get_ui(z));
+clean:
+  mpz_clear(q);
+  mpz_clear(z);
+  mpz_clear(d);
+  return res;
+}
+
+void element_multi_double(element_t n[], element_t a[], int m) {
+  element_ptr *temp1 = pbc_malloc(sizeof(*temp1)*m);
+  element_ptr *temp2 = pbc_malloc(sizeof(*temp2)*m);
+  int i;
+
+  for(i=0; i<m; i++) {
+    PBC_ASSERT_MATCH2(n[i], a[i]);
+    temp1[i] = n[i];
+    temp2[i] = a[i];
+  }
+  n[0]->field->multi_doub(temp1, temp2, m);
+  pbc_free(temp1);
+  pbc_free(temp2);
+}
+
+void element_multi_add(element_t n[], element_t a[],element_t b[], int m) {
+  size_t size = sizeof(element_ptr)*m;
+  element_ptr *temp1 = pbc_malloc(size);
+  element_ptr *temp2 = pbc_malloc(size);
+  element_ptr *temp3 = pbc_malloc(size);
+
+  int i;
+  for(i=0; i<m; i++){
+    PBC_ASSERT_MATCH3(n[i], a[i], b[i]);
+    temp1[i] = n[i];
+    temp2[i] = a[i];
+    temp3[i] = b[i];
+  }
+
+  n[0]->field->multi_add(temp1, temp2, temp3, m);
+  pbc_free(temp1);
+  pbc_free(temp2);
+  pbc_free(temp3);
+}
+
+element_ptr element_new(field_ptr f) {
+  element_ptr e = pbc_malloc(sizeof(*e));
+  element_init(e, f);
+  return e;
+}
+
+void element_free(element_ptr e) {
+  element_clear(e);
+  pbc_free(e);
+}
diff --git a/moon-abe/pbc-0.5.14/arith/fieldquadratic.c b/moon-abe/pbc-0.5.14/arith/fieldquadratic.c
new file mode 100644
index 00000000..bfb46027
--- /dev/null
+++ b/moon-abe/pbc-0.5.14/arith/fieldquadratic.c
@@ -0,0 +1,692 @@
+// Quadratic extension fields.
+//
+// The fq_ functions are for general quadratic extensions.
+// The fi_ functions are faster versions of some of these functions specialized
+// for fields extended by sqrt(-1).
+// TODO: Instead of lazily generating a quadratic nonresidue, in this case
+// we can use sqrt(base field nqr) as the nqr of the extension.
+
+#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_fieldquadratic.h"
+#include "pbc_memory.h"
+
+// Per-element data.
+typedef struct {
+  // Elements have the form x + ya, where a is the square root of a quadratic
+  // nonresidue in the base field.
+  element_t x;
+  element_t y;
+} *eptr;
+
+// Per-field data: we use ''data'' as a field_ptr to the base field.
+
+// Return the quadratic nonresidue used to build this field.
+// Should only be called from routines used exclusively by the generic quadratic
+// extension code.
+static inline element_ptr fq_nqr(field_ptr f) {
+  return field_get_nqr((field_ptr) f->data);
+}
+
+static void fq_init(element_ptr e) {
+  eptr p = e->data = pbc_malloc(sizeof(*p));
+  field_ptr f = e->field->data;
+  element_init(p->x, f);
+  element_init(p->y, f);
+}
+
+static void fq_clear(element_ptr e) {
+  eptr p = e->data;
+  element_clear(p->x);
+  element_clear(p->y);
+  pbc_free(e->data);
+}
+
+static void fq_set_si(element_ptr e, signed long int i) {
+  eptr p = e->data;
+  element_set_si(p->x, i);
+  element_set0(p->y);
+}
+
+static void fq_set_mpz(element_ptr e, mpz_t z) {
+  eptr p = e->data;
+  element_set_mpz(p->x, z);
+  element_set0(p->y);
+}
+
+// Projection: attempts to convert Re(e) to mpz.
+static void fq_to_mpz(mpz_t z, element_ptr e) {
+  eptr p = e->data;
+  element_to_mpz(z, p->x);
+}
+
+static void fq_set0(element_ptr e) {
+  eptr p = e->data;
+  element_set0(p->x);
+  element_set0(p->y);
+}
+
+static void fq_set1(element_ptr e) {
+  eptr p = e->data;
+  element_set1(p->x);
+  element_set0(p->y);
+}
+
+static int fq_is0(element_ptr e) {
+  eptr p = e->data;
+  return element_is0(p->x) && element_is0(p->y);
+}
+
+static int fq_is1(element_ptr e) {
+  eptr p = e->data;
+  return element_is1(p->x) && element_is0(p->y);
+}
+
+static size_t fq_out_str(FILE *stream, int base, element_ptr e) {
+  size_t result = 4, status;
+  eptr p = e->data;
+  if (EOF == fputc('[', stream)) return 0;
+  result = element_out_str(stream, base, p->x);
+  if (!result) return 0;
+  if (EOF == fputs(", ", stream)) return 0;
+  status = element_out_str(stream, base, p->y);
+  if (!status) return 0;
+  if (EOF == fputc(']', stream)) return 0;
+  return result + status;
+}
+
+static int fq_snprint(char *s, size_t n, element_ptr e) {
+  eptr p = e->data;
+  size_t result = 0, left;
+  int status;
+
+  #define clip_sub() {                   \
+    result += status;                    \
+    left = result >= n ? 0 : n - result; \
+  }
+
+  status = snprintf(s, n, "[");
+  if (status < 0) return status;
+  clip_sub();
+  status = element_snprint(s + result, left, p->x);
+  if (status < 0) return status;
+  clip_sub();
+  status = snprintf(s + result, left, ", ");
+  if (status < 0) return status;
+  clip_sub();
+  status = element_snprint(s + result, left, p->y);
+  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 fq_set_multiz(element_ptr e, multiz m) {
+  eptr p = e->data;
+  if (multiz_is_z(m)) {
+    element_set_multiz(p->x, m);
+    element_set0(p->y);
+    return;
+  }
+  element_set_multiz(p->x, multiz_at(m, 0));
+  if (2 > multiz_count(m)) element_set0(p->y);
+  else element_set_multiz(p->y, multiz_at(m, 1));
+}
+
+static int fq_set_str(element_ptr e, const char *s, int base) {
+  const char *cp = s;
+  element_set0(e);
+  while (*cp && isspace(*cp)) cp++;
+  if (*cp++ != '[') return 0;
+  eptr p = e->data;
+  cp += element_set_str(p->x, cp, base);
+  while (*cp && isspace(*cp)) cp++;
+  if (*cp++ != ',') return 0;
+  cp += element_set_str(p->y, cp, base);
+  if (*cp++ != ']') return 0;
+  return cp - s;
+}
+
+static int fq_sign(element_ptr n) {
+  int res;
+  eptr r = n->data;
+  res = element_sign(r->x);
+  if (!res) return element_sign(r->y);
+  return res;
+}
+
+static void fq_add(element_ptr n, element_ptr a, element_ptr b) {
+  eptr p = a->data;
+  eptr q = b->data;
+  eptr r = n->data;
+  element_add(r->x, p->x, q->x);
+  element_add(r->y, p->y, q->y);
+}
+
+static void fq_double(element_ptr n, element_ptr a) {
+  eptr p = a->data;
+  eptr r = n->data;
+  element_double(r->x, p->x);
+  element_double(r->y, p->y);
+}
+
+static void fq_sub(element_ptr n, element_ptr a, element_ptr b) {
+  eptr p = a->data;
+  eptr q = b->data;
+  eptr r = n->data;
+  element_sub(r->x, p->x, q->x);
+  element_sub(r->y, p->y, q->y);
+}
+
+static void fq_set(element_ptr n, element_ptr a) {
+  eptr p = a->data;
+  eptr r = n->data;
+  element_set(r->x, p->x);
+  element_set(r->y, p->y);
+}
+
+static void fq_mul(element_ptr n, element_ptr a, element_ptr b) {
+  eptr p = a->data;
+  eptr q = b->data;
+  eptr r = n->data;
+
+  element_ptr nqr = fq_nqr(n->field);
+  element_t e0, e1, e2;
+
+  element_init(e0, p->x->field);
+  element_init(e1, e0->field);
+  element_init(e2, e0->field);
+  /* naive:
+  element_mul(e0, p->x, q->x);
+  element_mul(e1, p->y, q->y);
+  element_mul(e1, e1, nqr);
+  element_add(e0, e0, e1);
+  element_mul(e1, p->x, q->y);
+  element_mul(e2, p->y, q->x);
+  element_add(e1, e1, e2);
+  element_set(r->x, e0);
+  element_set(r->y, e1);
+  */
+  // Karatsuba:
+  element_add(e0, p->x, p->y);
+  element_add(e1, q->x, q->y);
+  element_mul(e2, e0, e1);
+  element_mul(e0, p->x, q->x);
+  element_mul(e1, p->y, q->y);
+  element_mul(r->x, e1, nqr);
+  element_add(r->x, r->x, e0);
+  element_sub(e2, e2, e0);
+  element_sub(r->y, e2, e1);
+
+  element_clear(e0);
+  element_clear(e1);
+  element_clear(e2);
+}
+
+static void fq_mul_mpz(element_ptr n, element_ptr a, mpz_ptr z) {
+  eptr p = a->data;
+  eptr r = n->data;
+  element_mul_mpz(r->x, p->x, z);
+  element_mul_mpz(r->y, p->y, z);
+}
+
+static void fq_mul_si(element_ptr n, element_ptr a, signed long int z) {
+  eptr p = a->data;
+  eptr r = n->data;
+  element_mul_si(r->x, p->x, z);
+  element_mul_si(r->y, p->y, z);
+}
+
+static void fq_square(element_ptr n, element_ptr a) {
+  eptr p = a->data;
+  eptr r = n->data;
+  element_ptr nqr = fq_nqr(n->field);
+  element_t e0, e1;
+
+  element_init(e0, p->x->field);
+  element_init(e1, e0->field);
+  element_square(e0, p->x);
+  element_square(e1, p->y);
+  element_mul(e1, e1, nqr);
+  element_add(e0, e0, e1);
+  element_mul(e1, p->x, p->y);
+  //TODO: which is faster?
+  //element_add(e1, e1, e1);
+  element_double(e1, e1);
+  element_set(r->x, e0);
+  element_set(r->y, e1);
+  element_clear(e0);
+  element_clear(e1);
+}
+
+static void fq_neg(element_ptr n, element_ptr a) {
+  eptr p = a->data;
+  eptr r = n->data;
+  element_neg(r->x, p->x);
+  element_neg(r->y, p->y);
+}
+
+static void fq_random(element_ptr e) {
+  eptr p = e->data;
+  element_random(p->x);
+  element_random(p->y);
+}
+
+static int fq_cmp(element_ptr a, element_ptr b) {
+  eptr p = a->data;
+  eptr q = b->data;
+  return element_cmp(p->x, q->x) || element_cmp(p->y, q->y);
+}
+
+static void fq_invert(element_ptr n, element_ptr a) {
+  eptr p = a->data;
+  eptr r = n->data;
+  element_ptr nqr = fq_nqr(n->field);
+  element_t e0, e1;
+
+  element_init(e0, p->x->field);
+  element_init(e1, e0->field);
+  element_square(e0, p->x);
+  element_square(e1, p->y);
+  element_mul(e1, e1, nqr);
+  element_sub(e0, e0, e1);
+  element_invert(e0, e0);
+  element_mul(r->x, p->x, e0);
+  element_neg(e0, e0);
+  element_mul(r->y, p->y, e0);
+
+  element_clear(e0);
+  element_clear(e1);
+}
+
+static void fq_from_hash(element_ptr n, void *data, int len) {
+  eptr r = n->data;
+  int k = len / 2;
+  element_from_hash(r->x, data, k);
+  element_from_hash(r->y, (char *)data + k, len - k);
+}
+
+static int fq_length_in_bytes(element_ptr e) {
+  eptr p = e->data;
+  return element_length_in_bytes(p->x) + element_length_in_bytes(p->y);
+}
+
+static int fq_to_bytes(unsigned char *data, element_t e) {
+  eptr p = e->data;
+  int len;
+  len = element_to_bytes(data, p->x);
+  len += element_to_bytes(data + len, p->y);
+  return len;
+}
+
+static int fq_from_bytes(element_t e, unsigned char *data) {
+  eptr p = e->data;
+  int len;
+  len = element_from_bytes(p->x, data);
+  len += element_from_bytes(p->y, data + len);
+  return len;
+}
+
+static int fq_is_sqr(element_ptr e) {
+  //x + y sqrt(nqr) is a square iff x^2 - nqr y^2 is (in the base field)
+  eptr p = e->data;
+  element_t e0, e1;
+  element_ptr nqr = fq_nqr(e->field);
+  int result;
+  element_init(e0, p->x->field);
+  element_init(e1, e0->field);
+  element_square(e0, p->x);
+  element_square(e1, p->y);
+  element_mul(e1, e1, nqr);
+  element_sub(e0, e0, e1);
+  result = element_is_sqr(e0);
+  element_clear(e0);
+  element_clear(e1);
+  return result;
+}
+
+static void fq_sqrt(element_ptr n, element_ptr e) {
+  eptr p = e->data;
+  eptr r = n->data;
+  element_ptr nqr = fq_nqr(n->field);
+  element_t e0, e1, e2;
+
+  //if (a+b sqrt(nqr))^2 = x+y sqrt(nqr) then
+  //2a^2 = x +- sqrt(x^2 - nqr y^2)
+  //(take the sign which allows a to exist)
+  //and 2ab = y
+  element_init(e0, p->x->field);
+  element_init(e1, e0->field);
+  element_init(e2, e0->field);
+  element_square(e0, p->x);
+  element_square(e1, p->y);
+  element_mul(e1, e1, nqr);
+  element_sub(e0, e0, e1);
+  element_sqrt(e0, e0);
+  //e0 = sqrt(x^2 - nqr y^2)
+  element_add(e1, p->x, e0);
+  element_set_si(e2, 2);
+  element_invert(e2, e2);
+  element_mul(e1, e1, e2);
+  //e1 = (x + sqrt(x^2 - nqr y^2))/2
+  if (!element_is_sqr(e1)) {
+    element_sub(e1, e1, e0);
+    //e1 should be a square
+  }
+  element_sqrt(e0, e1);
+  element_add(e1, e0, e0);
+  element_invert(e1, e1);
+  element_mul(r->y, p->y, e1);
+  element_set(r->x, e0);
+  element_clear(e0);
+  element_clear(e1);
+  element_clear(e2);
+}
+
+static int fq_item_count(element_ptr e) {
+  UNUSED_VAR(e);
+  return 2;
+}
+
+static element_ptr fq_item(element_ptr e, int i) {
+  eptr p = e->data;
+  switch(i) {
+    case 0:
+      return p->x;
+    case 1:
+      return p->y;
+    default:
+      return NULL;
+  }
+}
+
+static void field_clear_fq(field_ptr f) {
+  UNUSED_VAR(f);
+  //f->order gets cleared automatically
+}
+
+static void fq_out_info(FILE *out, field_ptr f) {
+  field_ptr fbase = f->data;
+  element_fprintf(out, "extension x^2 + %B, base field: ", fq_nqr(f));
+  field_out_info(out, fbase);
+}
+
+// Specialized versions of some of the above for the case K[i].
+
+static void fi_mul(element_ptr n, element_ptr a, element_ptr b) {
+  eptr p = a->data;
+  eptr q = b->data;
+  eptr r = n->data;
+  element_t e0, e1, e2;
+
+  element_init(e0, p->x->field);
+  element_init(e1, e0->field);
+  element_init(e2, e0->field);
+  /* Naive method:
+  element_mul(e0, p->x, q->x);
+  element_mul(e1, p->y, q->y);
+  element_sub(e0, e0, e1);
+  element_mul(e1, p->x, q->y);
+  element_mul(e2, p->y, q->x);
+  element_add(e1, e1, e2);
+  element_set(r->x, e0);
+  element_set(r->y, e1);
+  */
+  // Karatsuba multiplicaiton:
+  element_add(e0, p->x, p->y);
+  element_add(e1, q->x, q->y);
+  element_mul(e2, e0, e1);
+  element_mul(e0, p->x, q->x);
+  element_sub(e2, e2, e0);
+  element_mul(e1, p->y, q->y);
+  element_sub(r->x, e0, e1);
+  element_sub(r->y, e2, e1);
+
+  element_clear(e0);
+  element_clear(e1);
+  element_clear(e2);
+}
+
+static void fi_square(element_ptr n, element_ptr a) {
+  eptr p = a->data;
+  eptr r = n->data;
+  element_t e0, e1;
+
+  element_init(e0, p->x->field);
+  element_init(e1, e0->field);
+  // Re(n) = x^2 - y^2 = (x+y)(x-y)
+  element_add(e0, p->x, p->y);
+  element_sub(e1, p->x, p->y);
+  element_mul(e0, e0, e1);
+  // Im(n) = 2xy
+  element_mul(e1, p->x, p->y);
+  element_add(e1, e1, e1);
+  element_set(r->x, e0);
+  element_set(r->y, e1);
+  element_clear(e0);
+  element_clear(e1);
+}
+
+static void fi_invert(element_ptr n, element_ptr a) {
+  eptr p = a->data;
+  eptr r = n->data;
+  element_t e0, e1;
+
+  element_init(e0, p->x->field);
+  element_init(e1, e0->field);
+  element_square(e0, p->x);
+  element_square(e1, p->y);
+  element_add(e0, e0, e1);
+  element_invert(e0, e0);
+  element_mul(r->x, p->x, e0);
+  element_neg(e0, e0);
+  element_mul(r->y, p->y, e0);
+
+  element_clear(e0);
+  element_clear(e1);
+}
+
+static int fi_is_sqr(element_ptr e) {
+  // x + yi is a square <=> x^2 + y^2 is (in the base field).
+
+  // Proof: (=>) if x+yi = (a+bi)^2, then a^2 - b^2 = x, 2ab = y,
+  // thus (a^2 + b^2)^2 = (a^2 - b^2)^2 + (2ab)^2 =  x^2 + y^2
+
+  // (<=) Suppose A^2 = x^2 + y^2. If there exist a, b satisfying:
+  //   a^2 = (+-A + x)/2, b^2 = (+-A - x)/2
+  // then (a + bi)^2 = x + yi.
+  //
+  // We show that exactly one of (A + x)/2, (-A + x)/2 is a quadratic residue
+  // (thus a, b do exist). Suppose not. Then the product (x^2 - A^2) / 4 is
+  // some quadratic residue, a contradiction since this would imply x^2 - A^2 =
+  // -y^2 is also a quadratic residue, but we know -1 is not a quadratic
+  // residue. QED.
+  eptr p = e->data;
+  element_t e0, e1;
+  int result;
+  element_init(e0, p->x->field);
+  element_init(e1, e0->field);
+  element_square(e0, p->x);
+  element_square(e1, p->y);
+  element_add(e0, e0, e1);
+  result = element_is_sqr(e0);
+  element_clear(e0);
+  element_clear(e1);
+  return result;
+}
+
+static void fi_sqrt(element_ptr n, element_ptr e) {
+  eptr p = e->data;
+  eptr r = n->data;
+  element_t e0, e1, e2;
+
+  // If (a+bi)^2 = x+yi then 2a^2 = x +- sqrt(x^2 + y^2)
+  // where we choose the sign so that a exists, and 2ab = y.
+  // Thus 2b^2 = - (x -+ sqrt(x^2 + y^2)).
+  element_init(e0, p->x->field);
+  element_init(e1, e0->field);
+  element_init(e2, e0->field);
+  element_square(e0, p->x);
+  element_square(e1, p->y);
+  element_add(e0, e0, e1);
+  element_sqrt(e0, e0);
+  // e0 = sqrt(x^2 + y^2)
+  element_add(e1, p->x, e0);
+  element_set_si(e2, 2);
+  element_invert(e2, e2);
+  element_mul(e1, e1, e2);
+  // e1 = (x + sqrt(x^2 + y^2))/2
+  if (!element_is_sqr(e1)) {
+    element_sub(e1, e1, e0);
+    // e1 should be a square.
+  }
+  element_sqrt(e0, e1);
+  element_add(e1, e0, e0);
+  element_invert(e1, e1);
+  element_mul(r->y, p->y, e1);
+  element_set(r->x, e0);
+  element_clear(e0);
+  element_clear(e1);
+  element_clear(e2);
+}
+
+static void fi_out_info(FILE *out, field_ptr f) {
+  field_ptr fbase = f->data;
+  fprintf(out, "extension x^2 + 1, base field: ");
+  field_out_info(out, fbase);
+}
+
+static void field_clear_fi(field_ptr f) {
+  UNUSED_VAR(f);
+}
+
+// All the above should be static.
+
+void element_field_to_quadratic(element_ptr r, element_ptr a) {
+  eptr p = r->data;
+  element_set(p->x, a);
+  element_set0(p->y);
+}
+
+void element_field_to_fi(element_ptr a, element_ptr b) {
+  element_field_to_quadratic(a, b);
+}
+
+static element_ptr fq_get_x(element_ptr a) {
+  return ((eptr) a->data)->x;
+}
+
+static element_ptr fq_get_y(element_ptr a) {
+  return ((eptr) a->data)->y;
+}
+
+void field_init_quadratic(field_ptr f, field_ptr fbase) {
+  field_init(f);
+
+  f->field_clear = field_clear_fq;
+  f->data = fbase;
+
+  f->init = fq_init;
+  f->clear = fq_clear;
+  f->set_si = fq_set_si;
+  f->set_mpz = fq_set_mpz;
+  f->to_mpz = fq_to_mpz;
+  f->out_str = fq_out_str;
+  f->snprint = fq_snprint;
+  f->set_multiz = fq_set_multiz;
+  f->set_str = fq_set_str;
+  f->sign = fq_sign;
+  f->add = fq_add;
+  f->sub = fq_sub;
+  f->set = fq_set;
+  f->mul = fq_mul;
+  f->mul_mpz = fq_mul_mpz;
+  f->mul_si = fq_mul_si;
+  f->square = fq_square;
+  f->doub = fq_double;
+  f->neg = fq_neg;
+  f->cmp = fq_cmp;
+  f->invert = fq_invert;
+  f->random = fq_random;
+  f->from_hash = fq_from_hash;
+  f->is1 = fq_is1;
+  f->is0 = fq_is0;
+  f->set0 = fq_set0;
+  f->set1 = fq_set1;
+  f->is_sqr = fq_is_sqr;
+  f->sqrt = fq_sqrt;
+  f->to_bytes = fq_to_bytes;
+  f->from_bytes = fq_from_bytes;
+  f->out_info = fq_out_info;
+  f->item_count = fq_item_count;
+  f->item = fq_item;
+  f->get_x = fq_get_x;
+  f->get_y = fq_get_y;
+
+  mpz_mul(f->order, fbase->order, fbase->order);
+  if (fbase->fixed_length_in_bytes < 0) {
+    f->length_in_bytes = fq_length_in_bytes;
+    f->fixed_length_in_bytes = -1;
+  } else {
+    f->fixed_length_in_bytes = 2 * fbase->fixed_length_in_bytes;
+  }
+}
+
+void field_init_fi(field_ptr f, field_ptr fbase) {
+  field_init(f);
+  f->field_clear = field_clear_fi;
+  f->data = fbase;
+  f->init = fq_init;
+  f->clear = fq_clear;
+  f->set_si = fq_set_si;
+  f->set_mpz = fq_set_mpz;
+  f->to_mpz = fq_to_mpz;
+  f->out_str = fq_out_str;
+  f->snprint = fq_snprint;
+  f->set_multiz = fq_set_multiz;
+  f->set_str = fq_set_str;
+  f->sign = fq_sign;
+  f->add = fq_add;
+  f->sub = fq_sub;
+  f->set = fq_set;
+  f->mul = fi_mul;
+  f->mul_mpz = fq_mul_mpz;
+  f->mul_si = fq_mul_si;
+  f->square = fi_square;
+  f->doub = fq_double;
+  f->neg = fq_neg;
+  f->cmp = fq_cmp;
+  f->invert = fi_invert;
+  f->random = fq_random;
+  f->from_hash = fq_from_hash;
+  f->is1 = fq_is1;
+  f->is0 = fq_is0;
+  f->set0 = fq_set0;
+  f->set1 = fq_set1;
+  f->is_sqr = fi_is_sqr;
+  f->sqrt = fi_sqrt;
+  f->to_bytes = fq_to_bytes;
+  f->from_bytes = fq_from_bytes;
+  f->out_info = fi_out_info;
+  f->item_count = fq_item_count;
+  f->item = fq_item;
+  f->get_x = fq_get_x;
+  f->get_y = fq_get_y;
+
+  mpz_mul(f->order, fbase->order, fbase->order);
+  if (fbase->fixed_length_in_bytes < 0) {
+    f->length_in_bytes = fq_length_in_bytes;
+    f->fixed_length_in_bytes = -1;
+  } else {
+    f->fixed_length_in_bytes = 2 * fbase->fixed_length_in_bytes;
+  }
+}
diff --git a/moon-abe/pbc-0.5.14/arith/fp.c b/moon-abe/pbc-0.5.14/arith/fp.c
new file mode 100644
index 00000000..e0127a8e
--- /dev/null
+++ b/moon-abe/pbc-0.5.14/arith/fp.c
@@ -0,0 +1,49 @@
+// F_p initialization.
+//
+// Specific implementations of F_p are found in naivefp.c, fastfp.c, fasterfp.c
+// and montfp.c. For pairing-based cryptosystems, montfp.c is the fastest.
+// I keep all versions around for testing, and also to show off the modularity
+// of the code.
+
+#include <stdarg.h>
+#include <stdio.h>
+#include <stdint.h> // for intptr_t
+#include <gmp.h>
+#include <string.h>
+#include "pbc_utils.h"
+#include "pbc_field.h"
+#include "pbc_fp.h"
+
+// By default, use the montfp.c implementation of F_p. After
+// pbc_tweak_use_fp(), future field_init_fp calls will use the specified
+// implementation. This is useful for benchmarking and testing.
+static void (*option_fpinit) (field_ptr f, mpz_t prime) = field_init_mont_fp;
+
+void pbc_tweak_use_fp(char *s) {
+  if (!strcmp(s, "naive")) {
+    option_fpinit = field_init_naive_fp;
+  } else if (!strcmp(s, "fast")) {
+    option_fpinit = field_init_fast_fp;
+  } else if (!strcmp(s, "faster")) {
+    option_fpinit = field_init_faster_fp;
+  } else if (!strcmp(s, "mont")) {
+    option_fpinit = field_init_mont_fp;
+  } else {
+    pbc_error("no such Fp implementation: %s", s);
+  }
+}
+
+void field_init_fp(field_ptr f, mpz_t modulus) {
+  if (mpz_fits_ulong_p(modulus)) {
+    // If this case mattered, I'd have written a F_p implementation specialized
+    // for moduli that fits into machine words.
+    field_init_naive_fp(f, modulus);
+  } else {
+    if (mpz_odd_p(modulus)) {
+      option_fpinit(f, modulus);
+    } else {
+      // montfp.c only supports odd moduli.
+      field_init_faster_fp(f, modulus);
+    }
+  }
+}
diff --git a/moon-abe/pbc-0.5.14/arith/init_random.c b/moon-abe/pbc-0.5.14/arith/init_random.c
new file mode 100644
index 00000000..bd040a38
--- /dev/null
+++ b/moon-abe/pbc-0.5.14/arith/init_random.c
@@ -0,0 +1,18 @@
+#include <stdio.h>
+#include <stdint.h> // for intptr_t
+#include <stdlib.h>
+#include <gmp.h>
+#include "pbc_utils.h"
+#include "pbc_random.h"
+
+void pbc_init_random(void) {
+  FILE *fp;
+  fp = fopen("/dev/urandom", "rb");
+  if (!fp) {
+    pbc_warn("could not open /dev/urandom, using deterministic random number generator");
+    pbc_random_set_deterministic(0);
+  } else {
+    pbc_random_set_file("/dev/urandom");
+    fclose(fp);
+  }
+}
diff --git a/moon-abe/pbc-0.5.14/arith/init_random.win32.c b/moon-abe/pbc-0.5.14/arith/init_random.win32.c
new file mode 100644
index 00000000..ec7f8732
--- /dev/null
+++ b/moon-abe/pbc-0.5.14/arith/init_random.win32.c
@@ -0,0 +1,52 @@
+// Win32 Compatibility Code added by Yulian Kalev and Stefan Georg Weber.
+#include <stdio.h>
+#include <stdint.h> // for intptr_t
+#include <stdlib.h>
+#include <windows.h>
+#include <wincrypt.h>
+#include <gmp.h>
+#include "pbc_random.h"
+#include "pbc_utils.h"
+#include "pbc_memory.h"
+
+static void win32_mpz_random(mpz_t r, mpz_t limit, void *data) {
+  UNUSED_VAR (data);
+  HCRYPTPROV phProv;
+  unsigned int error;
+  if (!CryptAcquireContext(&phProv,NULL,NULL,PROV_RSA_FULL,0)) {
+    error = GetLastError();
+    if (error == 0x80090016) { //need to create a new keyset
+      if (!CryptAcquireContext(&phProv,NULL,NULL,PROV_RSA_FULL,CRYPT_NEWKEYSET)) {
+        pbc_error("Couldn't create CryptContext: %x", (int)GetLastError());
+        return;
+      }
+    } else {
+      pbc_error("Couldn't create CryptContext: %x", error);
+      return;
+    }
+  }
+  int n, bytecount, leftover;
+  unsigned char *bytes;
+  mpz_t z;
+  mpz_init(z);
+  n = mpz_sizeinbase(limit, 2);
+  bytecount = (n + 7) / 8;
+  leftover = n % 8;
+  bytes = (unsigned char *) pbc_malloc(bytecount);
+  for (;;) {
+    CryptGenRandom(phProv,bytecount,(byte *)bytes);
+    if (leftover) {
+      *bytes = *bytes % (1 << leftover);
+    }
+    mpz_import(z, bytecount, 1, 1, 0, 0, bytes);
+    if (mpz_cmp(z, limit) < 0) break;
+  }
+  CryptReleaseContext(phProv,0);
+  mpz_set(r, z);
+  mpz_clear(z);
+  pbc_free(bytes);
+}
+
+void pbc_init_random(void) {
+  pbc_random_set_function(win32_mpz_random, NULL);
+}
diff --git a/moon-abe/pbc-0.5.14/arith/montfp.c b/moon-abe/pbc-0.5.14/arith/montfp.c
new file mode 100644
index 00000000..c79bb72b
--- /dev/null
+++ b/moon-abe/pbc-0.5.14/arith/montfp.c
@@ -0,0 +1,596 @@
+// F_p using Montgomery representation.
+//
+// Let b = 256^sizeof(mp_limb_t).
+// Let R = b^t be the smallest power of b greater than the modulus p.
+// Then x is stored as xR (mod p).
+// Addition: same as naive implementation.
+// Multipication: Montgomery reduction.
+// Code assumes the modulus p is odd.
+//
+// TODO: mul_2exp(x, p->bytes * 8) could be replaced with
+// faster code that messes with GMP internals
+
+#include <stdarg.h>
+#include <stdio.h>
+#include <stdint.h> // for intptr_t
+#include <stdlib.h>
+#include <string.h>
+#include <gmp.h>
+#include "pbc_utils.h"
+#include "pbc_field.h"
+#include "pbc_random.h"
+#include "pbc_fp.h"
+#include "pbc_memory.h"
+
+// Per-field data.
+typedef struct {
+  size_t limbs;           // Number of limbs per element.
+  size_t bytes;           // Number of bytes per element.
+  mp_limb_t *primelimbs;  // Points to an array of limbs holding the modulus.
+  mp_limb_t negpinv;      // -p^-1 mod b
+  mp_limb_t *R;           // R mod p
+  mp_limb_t *R3;          // R^3 mod p
+} *fptr;
+
+// Per-element data.
+typedef struct {
+  char flag;     // flag == 0 means the element is zero.
+  mp_limb_t *d;  // Otherwise d points to an array holding the element.
+} *eptr;
+
+// Copies limbs of z into dst and zeroes any leading limbs, where n is the
+// total number of limbs.
+// Requires z to have at most n limbs.
+static inline void set_limbs(mp_limb_t *dst, mpz_t z, size_t n) {
+  size_t count;
+  mpz_export(dst, &count, -1, sizeof(mp_limb_t), 0, 0, z);
+  memset((void *) (((unsigned char *) dst) + count * sizeof(mp_limb_t)),
+         0, (n - count) * sizeof(mp_limb_t));
+}
+
+static void fp_init(element_ptr e) {
+  fptr p = e->field->data;
+  eptr ep = e->data = pbc_malloc(sizeof(*ep));
+  ep->flag = 0;
+  ep->d = pbc_malloc(p->bytes);
+}
+
+static void fp_clear(element_ptr e) {
+  eptr ep = e->data;
+  pbc_free(ep->d);
+  pbc_free(e->data);
+}
+
+static void fp_set_mpz(element_ptr e, mpz_ptr z) {
+  fptr p = e->field->data;
+  eptr ep = e->data;
+  if (!mpz_sgn(z)) ep->flag = 0;
+  else {
+    mpz_t tmp;
+    mpz_init(tmp);
+    mpz_mul_2exp(tmp, z, p->bytes * 8);
+    mpz_mod(tmp, tmp, e->field->order);
+    if (!mpz_sgn(tmp)) ep->flag = 0;
+    else {
+      set_limbs(ep->d, tmp, p->limbs);
+      ep->flag = 2;
+    }
+    mpz_clear(tmp);
+  }
+}
+
+static void fp_set_si(element_ptr e, signed long int op) {
+  fptr p = e->field->data;
+  eptr ep = e->data;
+  if (!op) ep->flag = 0;
+  else {
+    mpz_t tmp;
+    mpz_init(tmp);
+    // TODO: Could be optimized.
+    mpz_set_si(tmp, op);
+    mpz_mul_2exp(tmp, tmp, p->bytes * 8);
+    mpz_mod(tmp, tmp, e->field->order);
+    if (!mpz_sgn(tmp)) ep->flag = 0;
+    else {
+      set_limbs(ep->d, tmp, p->limbs);
+      ep->flag = 2;
+    }
+    mpz_clear(tmp);
+  }
+}
+
+// Montgomery reduction.
+// Algorithm II.4 from Blake, Seroussi and Smart.
+static void mont_reduce(mp_limb_t *x, mp_limb_t *y, fptr p) {
+  size_t t = p->limbs;
+  size_t i;
+  mp_limb_t flag = 0;
+  for (i = 0; i < t; i++) {
+    mp_limb_t u = y[i] * p->negpinv;
+    mp_limb_t carry = mpn_addmul_1(&y[i], p->primelimbs, t, u);
+    //mpn_add_1(&y[i+t], &y[i+t], t - i + 1, carry);
+    flag += mpn_add_1(&y[i + t], &y[i + t], t - i, carry);
+  }
+  if (flag || mpn_cmp(&y[t], p->primelimbs, t) >= 0) {
+    mpn_sub_n(x, &y[t], p->primelimbs, t);
+  } else {
+    // TODO: GMP set might be faster.
+    memcpy(x, &y[t], t * sizeof(mp_limb_t));
+  }
+}
+
+static void fp_to_mpz(mpz_ptr z, element_ptr e) {
+  eptr ep = e->data;
+  if (!ep->flag) mpz_set_ui(z, 0);
+  else {
+    // x is stored as xR.
+    // We must divide out R to convert to standard representation.
+    fptr p = e->field->data;
+    mp_limb_t tmp[2 * p->limbs];
+
+    memcpy(tmp, ep->d, p->limbs * sizeof(mp_limb_t));
+    memset(&tmp[p->limbs], 0, p->limbs * sizeof(mp_limb_t));
+    _mpz_realloc(z, p->limbs);
+    mont_reduce(z->_mp_d, tmp, p);
+    // Remove leading zero limbs.
+    for (z->_mp_size = p->limbs; !z->_mp_d[z->_mp_size - 1]; z->_mp_size--);
+  }
+}
+
+static void fp_set0(element_ptr e) {
+  eptr ep = e->data;
+  ep->flag = 0;
+}
+
+static void fp_set1(element_ptr e) {
+  fptr p = e->field->data;
+  eptr ep = e->data;
+  ep->flag = 2;
+  memcpy(ep->d, p->R, p->bytes);
+}
+
+static int fp_is1(element_ptr e) {
+  eptr ep = e->data;
+  if (!ep->flag) return 0;
+  else {
+    fptr p = e->field->data;
+    return !mpn_cmp(ep->d, p->R, p->limbs);
+  }
+}
+
+static int fp_is0(element_ptr e) {
+  eptr ep = e->data;
+  return !ep->flag;
+}
+
+static size_t fp_out_str(FILE * stream, int base, element_ptr e) {
+  size_t result;
+  mpz_t z;
+  mpz_init(z);
+  fp_to_mpz(z, e);
+  result = mpz_out_str(stream, base, z);
+  mpz_clear(z);
+  return result;
+}
+
+static int fp_snprint(char *s, size_t n, element_ptr e) {
+  int result;
+  mpz_t z;
+  mpz_init(z);
+  fp_to_mpz(z, e);
+  result = gmp_snprintf(s, n, "%Zd", z);
+  mpz_clear(z);
+  return result;
+}
+
+static int fp_set_str(element_ptr e, const char *s, int base) {
+  mpz_t z;
+  mpz_init(z);
+  int result = pbc_mpz_set_str(z, s, base);
+  mpz_mod(z, z, e->field->order);
+  fp_set_mpz(e, z);
+  mpz_clear(z);
+  return result;
+}
+
+static void fp_set(element_ptr c, element_ptr a) {
+  eptr ad = a->data;
+  eptr cd = c->data;
+  if (a == c) return;
+  if (!ad->flag) cd->flag = 0;
+  else {
+    fptr p = a->field->data;
+
+    // Assembly is faster, but I don't want to stoop to that level.
+    // Instead of memcpy(), we rewrite so GMP assembly ends up being invoked.
+    /*
+       memcpy(cd->d, ad->d, p->bytes);
+     */
+    mpz_t z1, z2;
+    z1->_mp_d = cd->d;
+    z2->_mp_d = ad->d;
+    z1->_mp_size = z1->_mp_alloc = z2->_mp_size = z2->_mp_alloc = p->limbs;
+    mpz_set(z1, z2);
+
+    cd->flag = 2;
+  }
+}
+
+static void fp_add(element_ptr c, element_ptr a, element_ptr b) {
+  eptr ad = a->data, bd = b->data;
+
+  if (!ad->flag) {
+    fp_set(c, b);
+  } else if (!bd->flag) {
+    fp_set(c, a);
+  } else {
+    eptr cd = c->data;
+    fptr p = a->field->data;
+    const size_t t = p->limbs;
+    mp_limb_t carry;
+    carry = mpn_add_n(cd->d, ad->d, bd->d, t);
+
+    if (carry) {
+      // Assumes result of following sub is not zero,
+      // i.e. modulus cannot be 2^(n * bits_per_limb).
+      mpn_sub_n(cd->d, cd->d, p->primelimbs, t);
+      cd->flag = 2;
+    } else {
+      int i = mpn_cmp(cd->d, p->primelimbs, t);
+      if (!i) {
+        cd->flag = 0;
+      } else {
+        cd->flag = 2;
+        if (i > 0) {
+          mpn_sub_n(cd->d, cd->d, p->primelimbs, t);
+        }
+      }
+    }
+  }
+}
+
+static void fp_double(element_ptr c, element_ptr a) {
+  eptr ad = a->data, cd = c->data;
+  if (!ad->flag) {
+    cd->flag = 0;
+  } else {
+    fptr p = c->field->data;
+    const size_t t = p->limbs;
+    if (mpn_lshift(cd->d, ad->d, t, 1)) {
+      cd->flag = 2;
+      // Again, assumes result is not zero.
+      mpn_sub_n(cd->d, cd->d, p->primelimbs, t);
+    } else {
+      int i = mpn_cmp(cd->d, p->primelimbs, t);
+      if (!i) {
+        cd->flag = 0;
+      } else {
+        cd->flag = 2;
+        if (i > 0) {
+          mpn_sub_n(cd->d, cd->d, p->primelimbs, t);
+        }
+      }
+    }
+  }
+}
+
+static void fp_halve(element_ptr c, element_ptr a) {
+  eptr ad = a->data, cd = c->data;
+  if (!ad->flag) {
+    cd->flag = 0;
+  } else {
+    fptr p = c->field->data;
+    const size_t t = p->limbs;
+    int carry = 0;
+    mp_limb_t *alimb = ad->d;
+    mp_limb_t *climb = cd->d;
+    if (alimb[0] & 1) {
+      carry = mpn_add_n(climb, alimb, p->primelimbs, t);
+    } else fp_set(c, a);
+
+    mpn_rshift(climb, climb, t, 1);
+    if (carry) climb[t - 1] |= ((mp_limb_t) 1) << (sizeof(mp_limb_t) * 8 - 1);
+  }
+}
+
+static void fp_neg(element_ptr c, element_ptr a) {
+  eptr ad = a->data, cd = c->data;
+  if (!ad->flag) cd->flag = 0;
+  else {
+    fptr p = a->field->data;
+    mpn_sub_n(cd->d, p->primelimbs, ad->d, p->limbs);
+    cd->flag = 2;
+  }
+}
+
+static void fp_sub(element_ptr c, element_ptr a, element_ptr b) {
+  eptr ad = a->data, bd = b->data;
+
+  if (!ad->flag) {
+    fp_neg(c, b);
+  } else if (!bd->flag) {
+    fp_set(c, a);
+  } else {
+    fptr p = c->field->data;
+    size_t t = p->limbs;
+    eptr cd = c->data;
+    int i = mpn_cmp(ad->d, bd->d, t);
+
+    if (i == 0) {
+      cd->flag = 0;
+    } else {
+      cd->flag = 2;
+      mpn_sub_n(cd->d, ad->d, bd->d, t);
+      if (i < 0) {
+        mpn_add_n(cd->d, cd->d, p->primelimbs, t);
+      }
+    }
+  }
+}
+
+// Montgomery multiplication.
+// See Blake, Seroussi and Smart.
+static inline void mont_mul(mp_limb_t *c, mp_limb_t *a, mp_limb_t *b,
+                            fptr p) {
+  // Instead of right shifting every iteration
+  // I allocate more room for the z array.
+  size_t i, t = p->limbs;
+  mp_limb_t z[2 * t + 1];
+  mp_limb_t u = (a[0] * b[0]) * p->negpinv;
+  mp_limb_t v = z[t] = mpn_mul_1(z, b, t, a[0]);
+  z[t] += mpn_addmul_1(z, p->primelimbs, t, u);
+  z[t + 1] = z[t] < v;  // Handle overflow.
+  for (i = 1; i < t; i++) {
+    u = (z[i] + a[i] * b[0]) * p->negpinv;
+    v = z[t + i] += mpn_addmul_1(z + i, b, t, a[i]);
+    z[t + i] += mpn_addmul_1(z + i, p->primelimbs, t, u);
+    z[t + i + 1] = z[t + i] < v;
+  }
+  if (z[t * 2] || mpn_cmp(z + t, p->primelimbs, t) >= 0) {
+    mpn_sub_n(c, z + t, p->primelimbs, t);
+  } else {
+    memcpy(c, z + t, t * sizeof(mp_limb_t));
+    // Doesn't seem to make a difference:
+    /*
+       mpz_t z1, z2;
+       z1->_mp_d = c;
+       z2->_mp_d = z + t;
+       z1->_mp_size = z1->_mp_alloc = z2->_mp_size = z2->_mp_alloc = t;
+       mpz_set(z1, z2);
+     */
+  }
+}
+
+static void fp_mul(element_ptr c, element_ptr a, element_ptr b) {
+  eptr ad = a->data, bd = b->data;
+  eptr cd = c->data;
+
+  if (!ad->flag || !bd->flag) {
+    cd->flag = 0;
+  } else {
+    fptr p = c->field->data;
+    mont_mul(cd->d, ad->d, bd->d, p);
+    cd->flag = 2;
+  }
+}
+
+static void fp_pow_mpz(element_ptr c, element_ptr a, mpz_ptr op) {
+  // Alternative: rewrite GMP mpz_powm().
+  fptr p = a->field->data;
+  eptr ad = a->data;
+  eptr cd = c->data;
+  if (!ad->flag) cd->flag = 0;
+  else {
+    mpz_t z;
+    mpz_init(z);
+    fp_to_mpz(z, a);
+    mpz_powm(z, z, op, a->field->order);
+    mpz_mul_2exp(z, z, p->bytes * 8);
+    mpz_mod(z, z, a->field->order);
+    set_limbs(cd->d, z, p->limbs);
+    mpz_clear(z);
+    cd->flag = 2;
+  }
+}
+
+// Inversion is slower than in a naive Fp implementation because of an extra
+// multiplication.
+// Requires nonzero a.
+static void fp_invert(element_ptr c, element_ptr a) {
+  eptr ad = a->data;
+  eptr cd = c->data;
+  fptr p = a->field->data;
+  mp_limb_t tmp[p->limbs];
+  mpz_t z;
+
+  mpz_init(z);
+
+  // Copy the limbs into a regular mpz_t so we can invert using the standard
+  // mpz_invert().
+  mpz_import(z, p->limbs, -1, sizeof(mp_limb_t), 0, 0, ad->d);
+  mpz_invert(z, z, a->field->order);
+  set_limbs(tmp, z, p->limbs);
+
+  // Normalize.
+  mont_mul(cd->d, tmp, p->R3, p);
+  cd->flag = 2;
+  mpz_clear(z);
+}
+
+static void fp_random(element_ptr a) {
+  fptr p = a->field->data;
+  eptr ad = a->data;
+  mpz_t z;
+  mpz_init(z);
+  pbc_mpz_random(z, a->field->order);
+  if (mpz_sgn(z)) {
+    mpz_mul_2exp(z, z, p->bytes * 8);
+    mpz_mod(z, z, a->field->order);
+    set_limbs(ad->d, z, p->limbs);
+    ad->flag = 2;
+  } else {
+    ad->flag = 0;
+  }
+  mpz_clear(z);
+}
+
+static void fp_from_hash(element_ptr a, void *data, int len) {
+  mpz_t z;
+
+  mpz_init(z);
+  pbc_mpz_from_hash(z, a->field->order, data, len);
+  fp_set_mpz(a, z);
+  mpz_clear(z);
+}
+
+static int fp_cmp(element_ptr a, element_ptr b) {
+  eptr ad = a->data, bd = b->data;
+  if (!ad->flag) return bd->flag;
+  else {
+    fptr p = a->field->data;
+    return mpn_cmp(ad->d, bd->d, p->limbs);
+    //return memcmp(ad->d, bd->d, p->limbs);
+  }
+}
+
+static int fp_sgn_odd(element_ptr a) {
+  eptr ad = a->data;
+  if (!ad->flag) return 0;
+  else {
+    mpz_t z;
+    mpz_init(z);
+    int res;
+    fp_to_mpz(z, a);
+    res = mpz_odd_p(z) ? 1 : -1;
+    mpz_clear(z);
+    return res;
+  }
+}
+
+static int fp_is_sqr(element_ptr a) {
+  eptr ad = a->data;
+  int res;
+  mpz_t z;
+  mpz_init(z);
+  // 0 is a square.
+  if (!ad->flag) return 1;
+  fp_to_mpz(z, a);
+  res = mpz_legendre(z, a->field->order) == 1;
+  mpz_clear(z);
+  return res;
+}
+
+static int fp_to_bytes(unsigned char *data, element_t a) {
+  mpz_t z;
+  int n = a->field->fixed_length_in_bytes;
+
+  mpz_init(z);
+  fp_to_mpz(z, a);
+  pbc_mpz_out_raw_n(data, n, z);
+  mpz_clear(z);
+  return n;
+}
+
+static int fp_from_bytes(element_t a, unsigned char *data) {
+  fptr p = a->field->data;
+  eptr ad = a->data;
+  int n;
+  mpz_t z;
+
+  mpz_init(z);
+
+  n = a->field->fixed_length_in_bytes;
+  mpz_import(z, n, 1, 1, 1, 0, data);
+  if (!mpz_sgn(z)) ad->flag = 0;
+  else {
+    ad->flag = 2;
+    mpz_mul_2exp(z, z, p->bytes * 8);
+    mpz_mod(z, z, a->field->order);
+    set_limbs(ad->d, z, p->limbs);
+  }
+  mpz_clear(z);
+  return n;
+}
+
+static void fp_field_clear(field_t f) {
+  fptr p = f->data;
+  pbc_free(p->primelimbs);
+  pbc_free(p->R);
+  pbc_free(p->R3);
+  pbc_free(p);
+}
+
+// The only public functions. All the above should be static.
+
+static void fp_out_info(FILE * out, field_ptr f) {
+  element_fprintf(out, "GF(%Zd): Montgomery representation", f->order);
+}
+
+void field_init_mont_fp(field_ptr f, mpz_t prime) {
+  PBC_ASSERT(!mpz_fits_ulong_p(prime), "modulus too small");
+  fptr p;
+  field_init(f);
+  f->init = fp_init;
+  f->clear = fp_clear;
+  f->set_si = fp_set_si;
+  f->set_mpz = fp_set_mpz;
+  f->out_str = fp_out_str;
+  f->snprint = fp_snprint;
+  f->set_str = fp_set_str;
+  f->add = fp_add;
+  f->sub = fp_sub;
+  f->set = fp_set;
+  f->mul = fp_mul;
+  f->doub = fp_double;
+  f->halve = fp_halve;
+  f->pow_mpz = fp_pow_mpz;
+  f->neg = fp_neg;
+  f->sign = fp_sgn_odd;
+  f->cmp = fp_cmp;
+  f->invert = fp_invert;
+  f->random = fp_random;
+  f->from_hash = fp_from_hash;
+  f->is1 = fp_is1;
+  f->is0 = fp_is0;
+  f->set0 = fp_set0;
+  f->set1 = fp_set1;
+  f->is_sqr = fp_is_sqr;
+  f->sqrt = element_tonelli;
+  f->field_clear = fp_field_clear;
+  f->to_bytes = fp_to_bytes;
+  f->from_bytes = fp_from_bytes;
+  f->to_mpz = fp_to_mpz;
+  f->out_info = fp_out_info;
+
+  // Initialize per-field data specific to this implementation.
+  p = f->data = pbc_malloc(sizeof(*p));
+  p->limbs = mpz_size(prime);
+  p->bytes = p->limbs * sizeof(mp_limb_t);
+  p->primelimbs = pbc_malloc(p->bytes);
+  mpz_export(p->primelimbs, &p->limbs, -1, sizeof(mp_limb_t), 0, 0, prime);
+
+  mpz_set(f->order, prime);
+  f->fixed_length_in_bytes = (mpz_sizeinbase(prime, 2) + 7) / 8;
+
+  // Compute R, R3 and negpinv.
+  mpz_t z;
+  mpz_init(z);
+
+  p->R = pbc_malloc(p->bytes);
+  p->R3 = pbc_malloc(p->bytes);
+  mpz_setbit(z, p->bytes * 8);
+  mpz_mod(z, z, prime);
+  set_limbs(p->R, z, p->limbs);
+
+  mpz_powm_ui(z, z, 3, prime);
+  set_limbs(p->R3, z, p->limbs);
+
+  mpz_set_ui(z, 0);
+
+  // Algorithm II.5 in Blake, Seroussi and Smart is better but this suffices
+  // since we're only doing it once.
+  mpz_setbit(z, p->bytes * 8);
+  mpz_invert(z, prime, z);
+  p->negpinv = -mpz_get_ui(z);
+  mpz_clear(z);
+}
diff --git a/moon-abe/pbc-0.5.14/arith/multiz.c b/moon-abe/pbc-0.5.14/arith/multiz.c
new file mode 100644
index 00000000..6c8b43cc
--- /dev/null
+++ b/moon-abe/pbc-0.5.14/arith/multiz.c
@@ -0,0 +1,589 @@
+// Multinomials over Z.
+// e.g. [[1, 2], 3, [4, [5, 6]]] means
+// (1 + 2y) + 3 x + (4 + (5 + 6z)y)x^2
+// Convenient interchange format for different groups, rings, and fields.
+
+// TODO: Canonicalize, e.g. [[1]], 0, 0] --> 1.
+
+#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_random.h"
+#include "pbc_fp.h"
+#include "pbc_memory.h"
+#include "misc/darray.h"
+
+// Per-element data.
+struct multiz_s {
+  // Either it's an mpz, or a list of mpzs.
+  char type;
+  union {
+    mpz_t z;
+    darray_t a;
+  };
+};
+
+enum {
+  T_MPZ,
+  T_ARR,
+};
+
+static multiz multiz_new_empty_list(void) {
+  multiz ep = pbc_malloc(sizeof(*ep));
+  ep->type = T_ARR;
+  darray_init(ep->a);
+  return ep;
+}
+
+void multiz_append(element_ptr x, element_ptr e) {
+  multiz l = x->data;
+  darray_append(l->a, e->data);
+}
+
+static multiz multiz_new(void) {
+  multiz ep = pbc_malloc(sizeof(*ep));
+  ep->type = T_MPZ;
+  mpz_init(ep->z);
+  return ep;
+}
+
+static void f_init(element_ptr e) {
+  e->data = multiz_new();
+}
+
+static void multiz_free(multiz ep) {
+  switch(ep->type) {
+    case T_MPZ:
+      mpz_clear(ep->z);
+      break;
+    default:
+      PBC_ASSERT(T_ARR == ep->type, "no such type");
+      darray_forall(ep->a, (void(*)(void*))multiz_free);
+      darray_clear(ep->a);
+      break;
+  }
+  pbc_free(ep);
+}
+
+static void f_clear(element_ptr e) {
+  multiz_free(e->data);
+}
+
+element_ptr multiz_new_list(element_ptr e) {
+  element_ptr x = pbc_malloc(sizeof(*x));
+  element_init_same_as(x, e);
+  multiz_free(x->data);
+  x->data = multiz_new_empty_list();
+  multiz_append(x, e);
+  return x;
+}
+
+static void f_set_si(element_ptr e, signed long int op) {
+  multiz_free(e->data);
+  f_init(e);
+  multiz ep = e->data;
+  mpz_set_si(ep->z, op);
+}
+
+static void f_set_mpz(element_ptr e, mpz_ptr z) {
+  multiz_free(e->data);
+  f_init(e);
+  multiz ep = e->data;
+  mpz_set(ep->z, z);
+}
+
+static void f_set0(element_ptr e) {
+  multiz_free(e->data);
+  f_init(e);
+}
+
+static void f_set1(element_ptr e) {
+  multiz_free(e->data);
+  f_init(e);
+  multiz ep = e->data;
+  mpz_set_ui(ep->z, 1);
+}
+
+static size_t multiz_out_str(FILE *stream, int base, multiz ep) {
+  switch(ep->type) {
+    case T_MPZ:
+      return mpz_out_str(stream, base, ep->z);
+    default:
+      PBC_ASSERT(T_ARR == ep->type, "no such type");
+      fputc('[', stream);
+      size_t res = 1;
+      int n = darray_count(ep->a);
+      int i;
+      for(i = 0; i < n; i++) {
+        if (i) res += 2, fputs(", ", stream);
+        res += multiz_out_str(stream, base, darray_at(ep->a, i));
+      }
+      fputc(']', stream);
+      res++;
+      return res;
+  }
+}
+
+static size_t f_out_str(FILE *stream, int base, element_ptr e) {
+  return multiz_out_str(stream, base, e->data);
+}
+
+void multiz_to_mpz(mpz_ptr z, multiz ep) {
+  while(ep->type == T_ARR) ep = darray_at(ep->a, 0);
+  PBC_ASSERT(T_MPZ == ep->type, "no such type");
+  mpz_set(z, ep->z);
+}
+
+static void f_to_mpz(mpz_ptr z, element_ptr a) {
+  multiz_to_mpz(z, a->data);
+}
+
+static int multiz_sgn(multiz ep) {
+  while(ep->type == T_ARR) ep = darray_at(ep->a, 0);
+  PBC_ASSERT(T_MPZ == ep->type, "no such type");
+  return mpz_sgn(ep->z);
+}
+
+static int f_sgn(element_ptr a) {
+  return multiz_sgn(a->data);
+}
+
+static void add_to_x(void *data,
+                     multiz x,
+                     void (*fun)(mpz_t, const mpz_t, void *scope_ptr),
+                     void *scope_ptr);
+
+static multiz multiz_new_unary(const multiz y,
+    void (*fun)(mpz_t, const mpz_t, void *scope_ptr), void *scope_ptr) {
+  multiz x = pbc_malloc(sizeof(*x));
+  switch(y->type) {
+    case T_MPZ:
+      x->type = T_MPZ;
+      mpz_init(x->z);
+      fun(x->z, y->z, scope_ptr);
+      break;
+    default:
+      PBC_ASSERT(T_ARR == ep->type, "no such type");
+      x->type = T_ARR;
+      darray_init(x->a);
+      darray_forall4(y->a,
+                     (void(*)(void*,void*,void*,void*))add_to_x,
+                     x,
+                     fun,
+                     scope_ptr);
+      break;
+  }
+  return x;
+}
+
+static void add_to_x(void *data,
+                     multiz x,
+                     void (*fun)(mpz_t, const mpz_t, void *scope_ptr),
+                     void *scope_ptr) {
+  darray_append(x->a, multiz_new_unary(data, fun, scope_ptr));
+}
+
+static void mpzset(mpz_t dst, const mpz_t src, void *scope_ptr) {
+  UNUSED_VAR(scope_ptr);
+  mpz_set(dst, src);
+}
+
+static multiz multiz_clone(multiz y) {
+  return multiz_new_unary(y, (void(*)(mpz_t, const mpz_t, void *))mpzset, NULL);
+}
+
+static multiz multiz_new_bin(const multiz a, const multiz b,
+    void (*fun)(mpz_t, const mpz_t, const mpz_t)) {
+  if (T_MPZ == a->type) {
+    if (T_MPZ == b->type) {
+      multiz x = multiz_new();
+      fun(x->z, a->z, b->z);
+      return x;
+    } else {
+      multiz x = multiz_clone(b);
+      multiz z = x;
+      PBC_ASSERT(T_ARR == z->type, "no such type");
+      while(z->type == T_ARR) z = darray_at(z->a, 0);
+      fun(z->z, a->z, z->z);
+      return x;
+    }
+  } else {
+    PBC_ASSERT(T_ARR == a->type, "no such type");
+    if (T_MPZ == b->type) {
+      multiz x = multiz_clone(a);
+      multiz z = x;
+      PBC_ASSERT(T_ARR == z->type, "no such type");
+      while(z->type == T_ARR) z = darray_at(z->a, 0);
+      fun(z->z, b->z, z->z);
+      return x;
+    } else {
+      PBC_ASSERT(T_ARR == b->type, "no such type");
+      int m = darray_count(a->a);
+      int n = darray_count(b->a);
+      int min = m < n ? m : n;
+      int max = m > n ? m : n;
+      multiz x = multiz_new_empty_list();
+      int i;
+      for(i = 0; i < min; i++) {
+        multiz z = multiz_new_bin(darray_at(a->a, i), darray_at(b->a, i), fun);
+        darray_append(x->a, z);
+      }
+      multiz zero = multiz_new();
+      for(; i < max; i++) {
+        multiz z = multiz_new_bin(m > n ? darray_at(a->a, i) : zero,
+                                  n > m ? darray_at(b->a, i) : zero,
+                                  fun);
+        darray_append(x->a, z);
+      }
+      multiz_free(zero);
+      return x;
+    }
+  }
+}
+static multiz multiz_new_add(const multiz a, const multiz b) {
+  return multiz_new_bin(a, b, mpz_add);
+}
+
+static void f_add(element_ptr n, element_ptr a, element_ptr b) {
+  multiz delme = n->data;
+  n->data = multiz_new_add(a->data, b->data);
+  multiz_free(delme);
+}
+
+static multiz multiz_new_sub(const multiz a, const multiz b) {
+  return multiz_new_bin(a, b, mpz_sub);
+}
+static void f_sub(element_ptr n, element_ptr a, element_ptr b) {
+  multiz delme = n->data;
+  n->data = multiz_new_sub(a->data, b->data);
+  multiz_free(delme);
+}
+
+static void mpzmul(mpz_t x, const mpz_t y, const mpz_t z) {
+  mpz_mul(x, y, z);
+}
+
+static multiz multiz_new_mul(const multiz a, const multiz b) {
+  if (T_MPZ == a->type) {
+    // Multiply each coefficient of b by a->z.
+    return multiz_new_unary(b, (void(*)(mpz_t, const mpz_t, void *))mpzmul, a->z);
+  } else {
+    PBC_ASSERT(T_ARR == a->type, "no such type");
+    if (T_MPZ == b->type) {
+      // Multiply each coefficient of a by b->z.
+      return multiz_new_unary(a, (void(*)(mpz_t, const mpz_t, void *))mpzmul, b->z);
+    } else {
+      PBC_ASSERT(T_ARR == b->type, "no such type");
+      int m = darray_count(a->a);
+      int n = darray_count(b->a);
+      int max = m + n - 1;
+      multiz x = multiz_new_empty_list();
+      int i;
+      multiz zero = multiz_new();
+      for(i = 0; i < max; i++) {
+        multiz z = multiz_new();
+        int j;
+        for (j = 0; j <= i; j++) {
+          multiz y = multiz_new_mul(j < m ? darray_at(a->a, j) : zero,
+                                    i - j < n ? darray_at(b->a, i - j) : zero);
+          multiz t = multiz_new_add(z, y);
+          multiz_free(y);
+          multiz_free(z);
+          z = t;
+        }
+        darray_append(x->a, z);
+      }
+      multiz_free(zero);
+      return x;
+    }
+  }
+}
+static void f_mul(element_ptr n, element_ptr a, element_ptr b) {
+  multiz delme = n->data;
+  n->data = multiz_new_mul(a->data, b->data);
+  multiz_free(delme);
+}
+
+static void f_mul_mpz(element_ptr n, element_ptr a, mpz_ptr z) {
+  multiz delme = n->data;
+  n->data = multiz_new_unary(a->data, (void(*)(mpz_t, const mpz_t, void *))mpzmul, z);
+  multiz_free(delme);
+}
+
+static void mulsi(mpz_t x, const mpz_t y, signed long *i) {
+  mpz_mul_si(x, y, *i);
+}
+
+static void f_mul_si(element_ptr n, element_ptr a, signed long int z) {
+  multiz delme = n->data;
+  n->data = multiz_new_unary(a->data, (void(*)(mpz_t, const mpz_t, void *))mulsi, &z);
+  multiz_free(delme);
+}
+
+static void mpzneg(mpz_t dst, const mpz_t src, void *scope_ptr) {
+  UNUSED_VAR(scope_ptr);
+  mpz_neg(dst, src);
+}
+
+static multiz multiz_new_neg(multiz z) {
+  return multiz_new_unary(z, (void(*)(mpz_t, const mpz_t, void *))mpzneg, NULL);
+}
+
+static void f_set(element_ptr n, element_ptr a) {
+  multiz delme = n->data;
+  n->data = multiz_clone(a->data);
+  multiz_free(delme);
+}
+
+static void f_neg(element_ptr n, element_ptr a) {
+  multiz delme = n->data;
+  n->data = multiz_new_neg(a->data);
+  multiz_free(delme);
+}
+
+static void f_div(element_ptr c, element_ptr a, element_ptr b) {
+  mpz_t d;
+  mpz_init(d);
+  element_to_mpz(d, b);
+  multiz delme = c->data;
+  c->data = multiz_new_unary(a->data, (void(*)(mpz_t, const mpz_t, void *))mpz_tdiv_q, d);
+  mpz_clear(d);
+  multiz_free(delme);
+}
+
+// Doesn't make sense if order is infinite.
+static void f_random(element_ptr n) {
+  multiz delme = n->data;
+  f_init(n);
+  multiz_free(delme);
+}
+
+static void f_from_hash(element_ptr n, void *data, int len) {
+  mpz_t z;
+  mpz_init(z);
+  mpz_import(z, len, -1, 1, -1, 0, data);
+  f_set_mpz(n, z);
+  mpz_clear(z);
+}
+
+static int f_is1(element_ptr n) {
+  multiz ep = n->data;
+  return ep->type == T_MPZ && !mpz_cmp_ui(ep->z, 1);
+}
+
+int multiz_is0(multiz m) {
+  return m->type == T_MPZ && mpz_is0(m->z);
+}
+
+static int f_is0(element_ptr n) {
+  return multiz_is0(n->data);
+}
+
+static int f_item_count(element_ptr e) {
+  multiz z = e->data;
+  if (T_MPZ == z->type) return 0;
+  return darray_count(z->a);
+}
+
+// TODO: Redesign multiz so this doesn't leak.
+static element_ptr f_item(element_ptr e, int i) {
+  multiz z = e->data;
+  if (T_MPZ == z->type) return NULL;
+  element_ptr r = malloc(sizeof(*r));
+  r->field = e->field;
+  r->data = darray_at(z->a, i);
+  return r;
+}
+
+// Usual meaning when both are integers.
+// Otherwise, compare coefficients.
+static int multiz_cmp(multiz a, multiz b) {
+  if (T_MPZ == a->type) {
+    if (T_MPZ == b->type) {
+      // Simplest case: both are integers.
+      return mpz_cmp(a->z, b->z);
+    }
+    // Leading coefficient of b.
+    while(T_ARR == b->type) b = darray_last(b->a);
+    PBC_ASSERT(T_MPZ == b->type, "no such type");
+    return -mpz_sgn(b->z);
+  }
+  PBC_ASSERT(T_ARR == a->type, "no such type");
+  if (T_MPZ == b->type) {
+    // Leading coefficient of a.
+    while(T_ARR == a->type) a = darray_last(a->a);
+    PBC_ASSERT(T_MPZ == a->type, "no such type");
+    return mpz_sgn(a->z);
+  }
+  PBC_ASSERT(T_ARR == b->type, "no such type");
+  int m = darray_count(a->a);
+  int n = darray_count(b->a);
+  if (m > n) {
+    // Leading coefficient of a.
+    while(T_ARR == a->type) a = darray_last(a->a);
+    PBC_ASSERT(T_MPZ == a->type, "no such type");
+    return mpz_sgn(a->z);
+  }
+  if (n > m) {
+    // Leading coefficient of b.
+    while(T_ARR == b->type) b = darray_last(b->a);
+    PBC_ASSERT(T_MPZ == b->type, "no such type");
+    return -mpz_sgn(b->z);
+  }
+  for(n--; n >= 0; n--) {
+    int i = multiz_cmp(darray_at(a->a, n), darray_at(b->a, n));
+    if (i) return i;
+  }
+  return 0;
+}
+static int f_cmp(element_ptr x, element_ptr y) {
+  return multiz_cmp(x->data, y->data);
+}
+
+static void f_field_clear(field_t f) { UNUSED_VAR (f); }
+
+// OpenSSL convention:
+//   4 bytes containing length
+//   followed by number in big-endian, most-significant bit set if negative
+//   (prepending null byte if necessary)
+// Positive numbers also the same as mpz_out_raw.
+static int z_to_bytes(unsigned char *data, element_t e) {
+  mpz_ptr z = e->data;
+  size_t msb = mpz_sizeinbase(z, 2);
+  size_t n = 4;
+  size_t i;
+
+  if (!(msb % 8)) {
+    data[4] = 0;
+    n++;
+  }
+  if (mpz_sgn(z) < 0) {
+    mpz_export(data + n, NULL, 1, 1, 1, 0, z);
+    data[4] |= 128;
+  } else {
+    mpz_export(data + n, NULL, 1, 1, 1, 0, z);
+  }
+  n += (msb + 7) / 8 - 4;
+  for (i=0; i<4; i++) {
+    data[i] = (n >> 8 * (3 - i));
+  }
+  n += 4;
+
+  return n;
+}
+
+static int z_from_bytes(element_t e, unsigned char *data) {
+  unsigned char *ptr;
+  size_t i, n;
+  mpz_ptr z = e->data;
+  mpz_t z1;
+  int neg = 0;
+
+  mpz_init(z1);
+  mpz_set_ui(z, 0);
+
+  ptr = data;
+  n = 0;
+  for (i=0; i<4; i++) {
+    n += ((unsigned int) *ptr) << 8 * (3 - i);
+    ptr++;
+  }
+  if (data[4] & 128) {
+    neg = 1;
+    data[4] &= 127;
+  }
+  for (i=0; i<n; i++) {
+    mpz_set_ui(z1, *ptr);
+    mpz_mul_2exp(z1, z1, 8 * (n - 1 - i));
+    ptr++;
+    mpz_add(z, z, z1);
+  }
+  mpz_clear(z1);
+  if (neg) mpz_neg(z, z);
+  return n;
+}
+
+static int z_length_in_bytes(element_ptr a) {
+  return (mpz_sizeinbase(a->data, 2) + 7) / 8 + 4;
+}
+
+static void f_out_info(FILE *out, field_ptr f) {
+  UNUSED_VAR(f);
+  fprintf(out, "Z multinomials");
+}
+
+static int f_set_str(element_ptr e, const char *s, int base) {
+  // TODO: Square brackets.
+  mpz_t z;
+  mpz_init(z);
+  int result = pbc_mpz_set_str(z, s, base);
+  f_set_mpz(e, z);
+  mpz_clear(z);
+  return result;
+}
+
+static void f_set_multiz(element_ptr e, multiz m) {
+  multiz delme = e->data;
+  e->data = multiz_clone(m);
+  multiz_free(delme);
+}
+
+void field_init_multiz(field_ptr f) {
+  field_init(f);
+  f->init = f_init;
+  f->clear = f_clear;
+  f->set_si = f_set_si;
+  f->set_mpz = f_set_mpz;
+  f->set_multiz = f_set_multiz;
+  f->set_str = f_set_str;
+  f->out_str = f_out_str;
+  f->sign = f_sgn;
+  f->add = f_add;
+  f->sub = f_sub;
+  f->set = f_set;
+  f->mul = f_mul;
+  f->mul_mpz = f_mul_mpz;
+  f->mul_si = f_mul_si;
+  f->neg = f_neg;
+  f->cmp = f_cmp;
+  f->div = f_div;
+  f->random = f_random;
+  f->from_hash = f_from_hash;
+  f->is1 = f_is1;
+  f->is0 = f_is0;
+  f->set0 = f_set0;
+  f->set1 = f_set1;
+  f->field_clear = f_field_clear;
+  f->to_bytes = z_to_bytes;
+  f->from_bytes = z_from_bytes;
+  f->to_mpz = f_to_mpz;
+  f->length_in_bytes = z_length_in_bytes;
+  f->item = f_item;
+  f->item_count = f_item_count;
+
+  f->out_info = f_out_info;
+
+  mpz_set_ui(f->order, 0);
+  f->data = NULL;
+  f->fixed_length_in_bytes = -1;
+}
+
+int multiz_is_z(multiz m) {
+  return T_MPZ == m->type;
+}
+
+int multiz_count(multiz m) {
+  if (T_ARR != m->type) return -1;
+  return darray_count(m->a);
+}
+
+multiz multiz_at(multiz m, int i) {
+  PBC_ASSERT(T_ARR == m->type, "wrong type");
+  PBC_ASSERT(darray_count(m->a) > i, "out of bounds");
+  return darray_at(m->a, i);
+}
diff --git a/moon-abe/pbc-0.5.14/arith/naivefp.c b/moon-abe/pbc-0.5.14/arith/naivefp.c
new file mode 100644
index 00000000..ceb1b7fb
--- /dev/null
+++ b/moon-abe/pbc-0.5.14/arith/naivefp.c
@@ -0,0 +1,270 @@
+// Naive implementation of F_p.
+// Little more than wrappers around GMP mpz functions.
+
+#include <stdarg.h>
+#include <stdio.h>
+#include <stdint.h> // for intptr_t
+#include <stdlib.h>
+#include <string.h>
+#include <gmp.h>
+#include "pbc_utils.h"
+#include "pbc_field.h"
+#include "pbc_random.h"
+#include "pbc_fp.h"
+#include "pbc_memory.h"
+
+static void zp_init(element_ptr e) {
+  e->data = pbc_malloc(sizeof(mpz_t));
+  mpz_init(e->data);
+}
+
+static void zp_clear(element_ptr e) {
+  mpz_clear(e->data);
+  pbc_free(e->data);
+}
+
+static void zp_set_si(element_ptr e, signed long int op) {
+  mpz_set_si(e->data, op);
+  mpz_mod(e->data, e->data, e->field->order);
+}
+
+static void zp_set_mpz(element_ptr e, mpz_ptr z) {
+  mpz_set(e->data, z);
+  mpz_mod(e->data, e->data, e->field->order);
+}
+
+static void zp_set0(element_ptr e) {
+  mpz_set_si(e->data, 0);
+}
+
+static void zp_set1(element_ptr e) {
+  mpz_set_si(e->data, 1);
+}
+
+static size_t zp_out_str(FILE * stream, int base, element_ptr e) {
+  return mpz_out_str(stream, base, e->data);
+}
+
+static int zp_snprint(char *s, size_t n, element_ptr e) {
+  return gmp_snprintf(s, n, "%Zd", e->data);
+}
+
+static int zp_set_str(element_ptr e, const char *s, int base) {
+  int result = pbc_mpz_set_str(e->data, s, base);
+  mpz_mod(e->data, e->data, e->field->order);
+  return result;
+}
+
+static int zp_sgn_odd(element_ptr a) {
+  mpz_ptr z = a->data;
+
+  return mpz_is0(z) ? 0 : (mpz_odd_p(z) ? 1 : -1);
+}
+
+static int zp_sgn_even(element_ptr a) {
+  mpz_t z;
+  mpz_init(z);
+  int res;
+
+  if (mpz_is0(a->data)) {
+    res = 0;
+  } else {
+    mpz_add(z, a->data, a->data);
+    res = mpz_cmp(z, a->field->order);
+  }
+  mpz_clear(z);
+  return res;
+}
+
+static void zp_add(element_ptr n, element_ptr a, element_ptr b) {
+  /*
+     mpz_add(n->data, a->data, b->data);
+     mpz_mod(n->data, n->data, n->field->order);
+   */
+  //This seems faster:
+  mpz_add(n->data, a->data, b->data);
+  if (mpz_cmp(n->data, n->field->order) >= 0) {
+    mpz_sub(n->data, n->data, n->field->order);
+  }
+}
+
+static void zp_sub(element_ptr n, element_ptr a, element_ptr b) {
+  //mpz_sub(n->data, a->data, b->data);
+  //mpz_mod(n->data, n->data, n->field->order);
+  mpz_sub(n->data, a->data, b->data);
+  if (mpz_sgn((mpz_ptr) n->data) < 0) {
+    mpz_add(n->data, n->data, n->field->order);
+  }
+}
+
+static void zp_square(element_ptr c, element_ptr a) {
+  /*
+     mpz_mul(c->data, a->data, a->data);
+     mpz_mod(c->data, c->data, c->field->order);
+   */
+  mpz_powm_ui(c->data, a->data, 2, c->field->order);
+
+  /*
+     const mpz_ptr prime = c->field->order;
+     const size_t t = prime->_mp_size;
+     const mpz_ptr p = a->data;
+     const mpz_ptr r = c->data;
+     mp_limb_t tmp[2 * t];
+     mp_limb_t qp[t + 1];
+
+     mpn_mul_n(tmp, p->_mp_d, p->_mp_d, t);
+
+     mpn_tdiv_qr(qp, r->_mp_d, 0, tmp, 2 * t, prime->_mp_d, t);
+   */
+}
+
+static void zp_double(element_ptr n, element_ptr a) {
+  //mpz_add(n->data, a->data, a->data);
+  mpz_mul_2exp(n->data, a->data, 1);
+  if (mpz_cmp(n->data, n->field->order) >= 0) {
+    mpz_sub(n->data, n->data, n->field->order);
+  }
+}
+
+static void zp_halve(element_ptr n, element_ptr a) {
+  mpz_ptr z = a->data;
+  if (mpz_odd_p(z)) {
+    mpz_add(n->data, z, a->field->order);
+    mpz_tdiv_q_2exp(n->data, n->data, 1);
+  } else {
+    mpz_tdiv_q_2exp(n->data, a->data, 1);
+  }
+}
+
+static void zp_mul(element_ptr n, element_ptr a, element_ptr b) {
+  mpz_mul(n->data, a->data, b->data);
+  mpz_mod(n->data, n->data, n->field->order);
+}
+
+static void zp_mul_mpz(element_ptr n, element_ptr a, mpz_ptr z) {
+  mpz_mul(n->data, a->data, z);
+  mpz_mod(n->data, n->data, n->field->order);
+}
+
+static void zp_mul_si(element_ptr n, element_ptr a, signed long int z) {
+  mpz_mul_si(n->data, a->data, z);
+  mpz_mod(n->data, n->data, n->field->order);
+}
+
+static void zp_pow_mpz(element_ptr n, element_ptr a, mpz_ptr z) {
+  mpz_powm(n->data, a->data, z, n->field->order);
+}
+
+static void zp_set(element_ptr n, element_ptr a) {
+  mpz_set(n->data, a->data);
+}
+
+static void zp_neg(element_ptr n, element_ptr a) {
+  if (mpz_is0(a->data)) {
+    mpz_set_ui(n->data, 0);
+  } else {
+    mpz_sub(n->data, n->field->order, a->data);
+  }
+}
+
+static void zp_invert(element_ptr n, element_ptr a) {
+  mpz_invert(n->data, a->data, n->field->order);
+}
+
+static void zp_random(element_ptr n) {
+  pbc_mpz_random(n->data, n->field->order);
+}
+
+static void zp_from_hash(element_ptr n, void *data, int len) {
+  pbc_mpz_from_hash(n->data, n->field->order, data, len);
+}
+
+static int zp_is1(element_ptr n) {
+  return !mpz_cmp_ui((mpz_ptr) n->data, 1);
+}
+
+static int zp_is0(element_ptr n) {
+  return mpz_is0(n->data);
+}
+
+static int zp_cmp(element_ptr a, element_ptr b) {
+  return mpz_cmp((mpz_ptr) a->data, (mpz_ptr) b->data);
+}
+
+static int zp_is_sqr(element_ptr a) {
+  //0 is a square
+  if (mpz_is0(a->data)) return 1;
+  return mpz_legendre(a->data, a->field->order) == 1;
+}
+
+static void zp_field_clear(field_t f) {
+  UNUSED_VAR(f);
+}
+
+static int zp_to_bytes(unsigned char *data, element_t e) {
+  int n;
+
+  n = e->field->fixed_length_in_bytes;
+
+  pbc_mpz_out_raw_n(data, n, e->data);
+  return n;
+}
+
+static int zp_from_bytes(element_t e, unsigned char *data) {
+  mpz_ptr z = e->data;
+  int n;
+  n = e->field->fixed_length_in_bytes;
+  mpz_import(z, n, 1, 1, 1, 0, data);
+  return n;
+}
+
+static void zp_to_mpz(mpz_ptr z, element_ptr a) {
+  mpz_set(z, a->data);
+}
+
+static void zp_out_info(FILE * out, field_ptr f) {
+  element_fprintf(out, "GF(%Zd), GMP wrapped", f->order);
+}
+
+void field_init_naive_fp(field_ptr f, mpz_t prime) {
+  field_init(f);
+  f->init = zp_init;
+  f->clear = zp_clear;
+  f->set_si = zp_set_si;
+  f->set_mpz = zp_set_mpz;
+  f->out_str = zp_out_str;
+  f->snprint = zp_snprint;
+  f->set_str = zp_set_str;
+  f->sign = mpz_odd_p(prime) ? zp_sgn_odd : zp_sgn_even;
+  f->add = zp_add;
+  f->sub = zp_sub;
+  f->set = zp_set;
+  f->square = zp_square;
+  f->doub = zp_double;
+  f->halve = zp_halve;
+  f->mul = zp_mul;
+  f->mul_mpz = zp_mul_mpz;
+  f->mul_si = zp_mul_si;
+  f->pow_mpz = zp_pow_mpz;
+  f->neg = zp_neg;
+  f->cmp = zp_cmp;
+  f->invert = zp_invert;
+  f->random = zp_random;
+  f->from_hash = zp_from_hash;
+  f->is1 = zp_is1;
+  f->is0 = zp_is0;
+  f->set0 = zp_set0;
+  f->set1 = zp_set1;
+  f->is_sqr = zp_is_sqr;
+  f->sqrt = element_tonelli;
+  f->field_clear = zp_field_clear;
+  f->to_bytes = zp_to_bytes;
+  f->from_bytes = zp_from_bytes;
+  f->to_mpz = zp_to_mpz;
+
+  f->out_info = zp_out_info;
+
+  mpz_set(f->order, prime);
+  f->data = NULL;
+  f->fixed_length_in_bytes = (mpz_sizeinbase(prime, 2) + 7) / 8;
+}
diff --git a/moon-abe/pbc-0.5.14/arith/poly.c b/moon-abe/pbc-0.5.14/arith/poly.c
new file mode 100644
index 00000000..bd2dad33
--- /dev/null
+++ b/moon-abe/pbc-0.5.14/arith/poly.c
@@ -0,0 +1,1724 @@
+#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;
+}
diff --git a/moon-abe/pbc-0.5.14/arith/random.c b/moon-abe/pbc-0.5.14/arith/random.c
new file mode 100644
index 00000000..68228b3f
--- /dev/null
+++ b/moon-abe/pbc-0.5.14/arith/random.c
@@ -0,0 +1,87 @@
+#include <stdio.h>
+#include <stdint.h> // for intptr_t
+#include <stdlib.h>
+#include <gmp.h>
+#include "pbc_random.h"
+#include "pbc_utils.h"
+#include "pbc_memory.h"
+
+void pbc_init_random(void);
+
+// Must use pointer due to lack of gmp_randstate_ptr.
+static gmp_randstate_t *get_rs(void) {
+  static int rs_is_ready;
+  static gmp_randstate_t rs;
+  if (!rs_is_ready) {
+    gmp_randinit_default(rs);
+    rs_is_ready = 1;
+  }
+  return &rs;
+}
+
+static void deterministic_mpz_random(mpz_t z, mpz_t limit, void *data) {
+  UNUSED_VAR (data);
+  mpz_urandomm(z, *get_rs(), limit);
+}
+
+static void file_mpz_random(mpz_t r, mpz_t limit, void *data) {
+  char *filename = (char *) data;
+  FILE *fp;
+  int n, bytecount, leftover;
+  unsigned char *bytes;
+  mpz_t z;
+  mpz_init(z);
+  fp = fopen(filename, "rb");
+  if (!fp) return;
+  n = mpz_sizeinbase(limit, 2);
+  bytecount = (n + 7) / 8;
+  leftover = n % 8;
+  bytes = (unsigned char *) pbc_malloc(bytecount);
+  for (;;) {
+    if (!fread(bytes, 1, bytecount, fp)) {
+      pbc_warn("error reading source of random bits");
+      return;
+    }
+    if (leftover) {
+      *bytes = *bytes % (1 << leftover);
+    }
+    mpz_import(z, bytecount, 1, 1, 0, 0, bytes);
+    if (mpz_cmp(z, limit) < 0) break;
+  }
+  fclose(fp);
+  mpz_set(r, z);
+  mpz_clear(z);
+  pbc_free(bytes);
+}
+
+static void (*current_mpz_random)(mpz_t, mpz_t, void *);
+static void *current_random_data;
+static int random_function_ready = 0;
+
+void pbc_random_set_function(void (*fun)(mpz_t, mpz_t, void *), void *data) {
+  current_mpz_random = fun;
+  current_random_data = data;
+  random_function_ready = 1;
+}
+
+void pbc_mpz_random(mpz_t z, mpz_t limit) {
+  if (!random_function_ready) pbc_init_random();
+  current_mpz_random(z, limit, current_random_data);
+}
+
+void pbc_mpz_randomb(mpz_t z, unsigned int bits) {
+  mpz_t limit;
+  mpz_init(limit);
+  mpz_setbit(limit, bits);
+  pbc_mpz_random(z, limit);
+  mpz_clear(limit);
+}
+
+void pbc_random_set_deterministic(unsigned int seed) {
+  gmp_randseed_ui(*get_rs(), seed);
+  pbc_random_set_function(deterministic_mpz_random, NULL);
+}
+
+void pbc_random_set_file(char *filename) {
+  pbc_random_set_function(file_mpz_random, filename);
+}
diff --git a/moon-abe/pbc-0.5.14/arith/ternary_extension_field.c b/moon-abe/pbc-0.5.14/arith/ternary_extension_field.c
new file mode 100644
index 00000000..3c79e3bd
--- /dev/null
+++ b/moon-abe/pbc-0.5.14/arith/ternary_extension_field.c
@@ -0,0 +1,950 @@
+/* $GF(3^m)     = GF(3)[x]/(x^m + x^t + 2)$
+   $GF(3^{2*m}) = GF(3^m)[x]/(x^2 + 1)$
+   $GF(3^{3*m}) = GF(3^m)[x]/(x^3 - x -1)$
+   $GF(3^{6*m}) = GF(3^{2*m})[x]/(x^3 - x -1)$
+
+   The "gf3_*" functions are for $GF(3)$.
+   The "gf3m_*" functions are for $GF(3^m)$.
+   The "gf32m_*" functions are for $GF(3^{2*m})$.
+   The "gf33m_*" functions are for $GF(3^{3*m})$ and $GF(3^{6*m})$.
+
+   (gf3m field_t).data is a pointer of struct params
+   (gf3m element_t).data is a pointer of unsigned long
+   (gf32m element_t).data is gf32m_ptr
+   (gf33m element_t).data is gf33m_ptr */
+
+#include <stdlib.h>
+#include <string.h>
+#include <stdarg.h>
+#include <stdio.h>
+#include <stdint.h>
+#include <gmp.h>
+#include "pbc_utils.h"
+#include "pbc_memory.h"
+#include "pbc_field.h"
+
+typedef unsigned long gf3;
+
+typedef struct { /* private data of $GF(3^m)$ */
+    unsigned int len; /* the number of native machine integers required to represent one GF(3^m) element */
+    unsigned int m; /* the irreducible polynomial is $x^m + x^t + 2$ */
+    unsigned int t; /* the irreducible polynomial is $x^m + x^t + 2$ */
+    element_ptr p; /* $p$ is the irreducible polynomial. */
+} params;
+
+typedef struct {
+    element_t _0, _1;
+} gf32m_s;
+
+typedef gf32m_s *gf32m_ptr;
+
+typedef struct {
+    element_t _0, _1, _2;
+} gf33m_s;
+
+typedef gf33m_s *gf33m_ptr;
+
+#define W (sizeof(unsigned long)*8) /* number of GF(3) elements in one processor integer */
+#define PARAM(e) ((params *)e->field->data)
+#define LEN(e) (PARAM(e)->len)
+#define SIZE(e) (LEN(e) * 2 * sizeof(unsigned long))
+#define DATA1(e) ((unsigned long*)e->data)
+#define DATA2(e) ((unsigned long*)e->data + LEN(e))
+#define GF32M(e) ((gf32m_s *)e->data)
+#define GF33M(e) ((gf33m_s *)e->data)
+#define BASE(e) ((field_ptr)e->field->data)
+#define print(e) {printf(#e": "); element_out_str(stdout, 10, e); printf("\n");}
+
+static size_t gf3m_out_str(FILE *stream, int base, element_t e) {
+    if (base != 10 && base != 16)
+        pbc_die("only support base 10 and base 16");
+    size_t size = 0;
+    unsigned i;
+    unsigned long *d = DATA1(e);
+    for (i = 0; i < LEN(e) * 2; i++) {
+        if (base == 16)
+            size += fprintf(stream, "0x%lx,", d[i]);
+        else
+            size += fprintf(stream, "%lu,", d[i]);
+    }
+    return size;
+}
+
+/* $a <- 0$ */
+static void gf3m_zero(element_t a) {
+    memset(a->data, 0, SIZE(a));
+}
+
+static void gf3m_init(element_t e) {
+    e->data = pbc_malloc(SIZE(e));
+    gf3m_zero(e);
+}
+
+static void gf3m_clear(element_t e) {
+    pbc_free(e->data);
+}
+
+/* $e <- a$ */
+static void gf3m_assign(element_t e, element_t a) {
+    memcpy(e->data, a->data, SIZE(a));
+}
+
+/* $a <- a/x$. $len$ is the number of elements in $a$ */
+static void shift_down(unsigned int len, unsigned long a[]) {
+    unsigned long h = 0;
+    const unsigned long x = 1ul << (W - 1);
+    int i;
+    for (i = len - 1; i >= 0; i--) {
+        unsigned long l = a[i] & 1;
+        a[i] >>= 1;
+        if (h)
+            a[i] |= x;
+        h = l;
+    }
+}
+
+/* $e <- e/x$ */
+static void gf3m_shift_down(element_t e) {
+    shift_down(LEN(e), DATA1(e));
+    shift_down(LEN(e), DATA2(e));
+}
+
+/* $a <- a*x$. $len$ is the number of elements in $a$ */
+static void shift_up(unsigned int len, unsigned long a[]) {
+    unsigned long l = 0;
+    const unsigned long x = 1ul << (W - 1), y = x - 1;
+    unsigned i;
+    for (i = 0; i < len; i++) {
+        unsigned long h = a[i] & x;
+        a[i] = ((a[i] & y) << 1) | l;
+        l = h ? 1 : 0;
+    }
+}
+
+/* $e <- e*x$ */
+static void gf3m_shift_up(element_t e) {
+    shift_up(LEN(e), DATA1(e));
+    shift_up(LEN(e), DATA2(e));
+}
+
+/* return the coefficient of $x^pos$ in $e$ */
+static unsigned gf3m_get(element_t e, unsigned pos) {
+    unsigned long *a1 = DATA1(e), *a2 = DATA2(e);
+    unsigned x = pos / W;
+    unsigned long y = 1ul << (pos % W), v1 = a1[x] & y, v2 = a2[x] & y;
+    return v1 ? 1 : (v2 ? 2 : 0);
+}
+
+/* set the coefficient of $x^pos$ as 1 */
+static void gf3m_set(element_t e, unsigned pos, unsigned value) {
+    unsigned long *a = DATA1(e);
+    /* assert value == 0, 1 or 2 */
+    if (value == 2)
+        a = DATA2(e);
+    if (value)
+        a[pos / W] |= 1ul << (pos % W);
+}
+
+/* $e <- a+b$ */
+static void gf3m_add(element_t e, element_t a, element_t b) {
+    unsigned long *e1 = DATA1(e), *e2 = DATA2(e), *a1 = DATA1(a),
+            *a2 = DATA2(a), *b1 = DATA1(b), *b2 = DATA2(b);
+    unsigned i;
+    for (i = 0; i < LEN(e); i++, e1++, e2++, a1++, a2++, b1++, b2++) {
+        unsigned long t = (*a1 | *a2) & (*b1 | *b2), c1 = t ^ (*a1 | *b1), c2 =
+                t ^ (*a2 | *b2);
+        *e1 = c1;
+        *e2 = c2;
+    }
+}
+
+/* $e <- x-y$ */
+static void gf3m_sub(element_t e, element_t a, element_t b) {
+    unsigned long *e1 = DATA1(e), *e2 = DATA2(e), *a1 = DATA1(a),
+            *a2 = DATA2(a), *b1 = DATA2(b), *b2 = DATA1(b);
+    unsigned i;
+    for (i = 0; i < LEN(e); i++, e1++, e2++, a1++, a2++, b1++, b2++) {
+        unsigned long t = (*a1 | *a2) & (*b1 | *b2), c1 = t ^ (*a1 | *b1), c2 =
+                t ^ (*a2 | *b2);
+        *e1 = c1;
+        *e2 = c2;
+    }
+}
+
+/* return 0 if $a == b$ in $GF(3^m)$, 1 otherwise. */
+static int gf3m_cmp(element_t a, element_t b) {
+    unsigned long *pa = DATA1(a), *pb = DATA1(b);
+    unsigned i;
+    for (i = 0; i < LEN(a) * 2; i++, pa++, pb++)
+        if (*pa != *pb)
+            return 1;
+    return 0;
+}
+
+/* $a <- 1$ */
+static void gf3m_one(element_t a) {
+    gf3m_zero(a);
+    *DATA1(a) = 1;
+}
+
+static int gf3m_is0(element_t e) {
+    unsigned i;
+    for (i = 0; i < LEN(e) * 2; i++)
+        if (DATA1(e)[i])
+            return 0;
+    return 1;
+}
+
+static int gf3m_is1(element_t e) {
+    unsigned i;
+    if (DATA1(e)[0] != 1)
+        return 0;
+    for (i = 1; i < LEN(e) * 2; i++)
+        if (DATA1(e)[i])
+            return 0;
+    return 1;
+}
+
+/* set $a$ to be a random element in $GF(3^m)$ */
+static void gf3m_random(element_t a) {
+    /* TODO: use uniform distribution? */
+    params *c = PARAM(a);
+    unsigned rm = c->m % W;
+    const unsigned long i1 = ~0ul;
+    unsigned long i2 = (1ul << rm) - 1;
+    unsigned long *a1 = DATA1(a), *a2 = DATA2(a);
+    unsigned i;
+    for (i = 0; i < c->len - 1; i++, a1++, a2++) { /* TODO: if $RAND_MAX < i1$ ? */
+        *a1 = rand() & i1;
+        *a2 = rand() & i1 & ~(*a1); /* assuring there is no bit that a1[x] & a2[x] == 1 */
+    }
+    unsigned long x = rm ? i2 : i1;
+    *a1 = rand() & x;
+    *a2 = rand() & x & ~(*a1);
+}
+
+static void swap(unsigned long *a, unsigned long *b) {
+    *a ^= *b;
+    *b ^= *a;
+    *a ^= *b;
+}
+
+/* $y <- (-x)$ */
+static void gf3m_neg(element_t y, element_t x) {
+    unsigned long *a1 = DATA1(x), *a2 = DATA2(x), *c1 = DATA1(y),
+            *c2 = DATA2(y);
+    if (a1 == c1) {
+        unsigned i;
+        for (i = 0; i < LEN(y); i++, a1++, a2++)
+            swap(a1, a2);
+    } else {
+        memcpy(c1, a2, SIZE(y) / 2);
+        memcpy(c2, a1, SIZE(y) / 2);
+    }
+}
+
+/* doing reduction
+ * The function returns the value of $a$ modulo $the irreducible trinomial$.
+ * $degree$ equals the degree of $a$.
+ * $2*len$ is the number of elements in $a$ */
+static void gf3m_reduct(element_t e, unsigned len, unsigned degree) {
+    // the $len$ argument exists because sometimes $len != p->len$
+    params *p = PARAM(e);
+    unsigned old = p->len;
+    p->len = len;
+    element_t px;
+    element_init(px, e->field);
+    gf3m_set(px, degree, 1);
+    gf3m_set(px, degree - p->m + p->t, 1);
+    gf3m_set(px, degree - p->m, 2);
+    while (degree >= p->m) {
+        unsigned v = gf3m_get(e, degree);
+        if (v == 1)
+            gf3m_sub(e, e, px);
+        else if (v == 2)
+            gf3m_add(e, e, px);
+        degree--;
+        gf3m_shift_down(px);
+    }
+    element_clear(px);
+    p->len = old;
+}
+
+/* doing multiplication of $n \in \{0,1,2\}$ and $a$ in $GF(3^m)$
+ * The function sets $e <- n * a$. */
+static void gf3m_f1(element_t e, unsigned n, element_t a) {
+    /* assert $e$ is not $a$ */
+    if (n == 0)
+        memset(DATA1(e), 0, SIZE(e));
+    else if (n == 1)
+        memcpy(DATA1(e), DATA1(a), SIZE(e));
+    else {
+        memcpy(DATA1(e), DATA2(a), SIZE(e) / 2);
+        memcpy(DATA2(e), DATA1(a), SIZE(e) / 2);
+    }
+}
+
+/* $e <- e*x mod p(x)$ */
+static void gf3m_f2(element_t e) {
+    params *p = PARAM(e);
+    gf3m_shift_up(e);
+    unsigned v = gf3m_get(e, p->m);
+    if (v == 1)
+        gf3m_sub(e, e, p->p);
+    else if (v == 2)
+        gf3m_add(e, e, p->p);
+}
+
+/* doing multiplication in GF(3^m)
+ * The function sets $e == a*b \in GF(3^m)$ */
+static void gf3m_mult(element_t e, element_ptr a, element_t b) {
+    params *p = PARAM(a);
+    element_t aa, t, c;
+    element_init(aa, a->field);
+    element_set(aa, a);
+    a = aa; // clone $a$
+    element_init(t, a->field);
+    element_init(c, a->field);
+    unsigned i;
+    for (i = 0; i < p->m; i++) {
+        unsigned v = gf3m_get(b, i);
+        gf3m_f1(t, v, a); /* t == b[i]*a in GF(3^m) */
+        gf3m_add(c, c, t); /* c += b[i]*a in GF(3^m) */
+        gf3m_f2(a); /* a == a*x in GF(3^m) */
+    }
+    element_set(e, c);
+    element_clear(t);
+    element_clear(c);
+    element_clear(aa);
+}
+
+/* $e <- x^3$ */
+static void gf3m_cubic(element_t e, element_t x) {
+    /* TODO: faster algorithm */
+    params *p = PARAM(x);
+    unsigned old = p->len;
+    unsigned len = (3 * p->m - 2 + W - 1) / W; /* length of $b1 */
+    p->len = len;
+    element_t a;
+    element_init(a, x->field);
+    unsigned i;
+    for (i = 0; i < p->m; i++) {
+        p->len = old;
+        unsigned v = gf3m_get(x, i);
+        p->len = len;
+        gf3m_set(a, 3 * i, v);
+    }
+    gf3m_reduct(a, len, 3 * p->m - 3);
+    p->len = old;
+    memcpy(DATA1(e), DATA1(a), SIZE(e) / 2);
+    memcpy(DATA2(e), DATA1(a) + len, SIZE(e) / 2);
+    element_clear(a);
+}
+
+/* multiplication modulo 3 of two elements in GF(3)
+ * for example, $mult(2,2) == 1$, and $mult(1,2) == 2$ */
+static unsigned gf3_mult(unsigned a, unsigned b) {
+    static const unsigned l[] = { 0, 1, 2, 0, 1 };
+    return l[a * b];
+}
+
+static void gf3m_swap(element_t a, element_t b) {
+    unsigned long *p = DATA1(a);
+    a->data = b->data;
+    b->data = p;
+}
+
+/* computing the inversion of an element $a$ in GF(3^m), i.e., $e <- a^{-1}$
+ The algorithm is by Tim Kerins, Emanuel Popovici and William Marnane
+ in the paper of "Algorithms and Architectures for use in FPGA",
+ Lecture Notes in Computer Science, 2004, Volume 3203/2004, 74-83.
+ Note that $U$ must have an extra bit, i.e, (_m + W - 1) // W == (_m + W) // W */
+static void gf3m_invert(element_t e, element_t a) {
+    struct field_s *f = a->field;
+    params *p = PARAM(a);
+    unsigned lenA = p->len;
+    unsigned lenS = (3 * p->m + W - 1) / W;
+    p->len = lenS;
+    element_t S, R, t, U, V, t2;
+    element_init(S, f);
+    element_init(R, f);
+    element_init(t, f);
+    memcpy(DATA1(S), DATA1(p->p), lenA * sizeof(unsigned long)); /* S = p(x) */
+    memcpy(DATA1(S) + lenS, DATA1(p->p) + lenA, lenA * sizeof(unsigned long));
+    memcpy(DATA1(R), DATA1(a), lenA * sizeof(unsigned long)); /* R = _clone(a) */
+    memcpy(DATA1(R) + lenS, DATA1(a) + lenA, lenA * sizeof(unsigned long));
+    p->len = lenA;
+    element_init(U, f);
+    gf3m_one(U);
+    element_init(V, f);
+    element_init(t2, f);
+    unsigned d = 0, i, r_m, s_m, q, x;
+    for (i = 0; i < p->m * 2; i++) {
+        p->len = lenS;
+        r_m = gf3m_get(R, p->m), s_m = gf3m_get(S, p->m);
+        if (r_m == 0) {
+            gf3m_shift_up(R); /* R = xR */
+            p->len = lenA;
+            gf3m_f2(U); /* U = xU mod p */
+            d++;
+        } else {
+            q = gf3_mult(r_m, s_m);
+            gf3m_f1(t, q, R);
+            gf3m_sub(S, S, t); /* S = S-qR */
+            gf3m_shift_up(S); /* S = xS */
+            p->len = lenA;
+            gf3m_f1(t2, q, U);
+            gf3m_sub(V, V, t2); /* V = V-qU */
+            if (d == 0) {
+                gf3m_swap(S, R);
+                gf3m_swap(U, V);
+                gf3m_f2(U); /* U = xU mod p*/
+                d++;
+            } else {
+                x = gf3m_get(U, 0);
+                if (x == 1) /* assuring x|U */
+                    gf3m_add(U, U, p->p);
+                else if (x == 2)
+                    gf3m_sub(U, U, p->p);
+                gf3m_shift_down(U); /* divide U by $x$ */
+                d--;
+            }
+        }
+    }
+    p->len = lenS;
+    r_m = gf3m_get(R, p->m); /* assume r_m is not zero */
+    p->len = lenA;
+    if (r_m == 2)
+        gf3m_neg(U, U);
+    memcpy(e->data, U->data, lenA * 2 * sizeof(unsigned long));
+    element_clear(S);
+    element_clear(R);
+    element_clear(U);
+    element_clear(V);
+    element_clear(t);
+    element_clear(t2);
+}
+
+static void gf3m_sqrt(element_t e, element_t a) {
+    field_ptr f = e->field;
+    mpz_t t;
+    mpz_init(t); // t == (field_order  + 1) / 4
+    mpz_set(t, f->order);
+    mpz_add_ui(t, t, 1);
+    mpz_tdiv_q_2exp(t, t, 2);
+    element_pow_mpz(e, a, t);
+    mpz_clear(t);
+}
+
+int gf3m_to_bytes(unsigned char *d, element_ptr e) {
+    unsigned long *a = DATA1(e), *b = DATA2(e);
+    unsigned long i, j;
+    for (i = 0; i < LEN(e); i++, a++, b++) {
+        for (j = 0; j < sizeof(unsigned long) * 8; j += 8) {
+            *(d++) = (unsigned char) ((*a) >> j);
+            *(d++) = (unsigned char) ((*b) >> j);
+        }
+    }
+    return SIZE(e);
+}
+
+int gf3m_from_bytes(element_ptr e, unsigned char *d) {
+    unsigned long *a = DATA1(e), *b = DATA2(e);
+    unsigned i;
+    int j;
+    for (i = 0; i < LEN(e); i++, a++, b++, d += sizeof(unsigned long) * 2) {
+        *a = 0, *b = 0;
+        j = 2 * sizeof(unsigned long) - 2;
+        while (j >= 0) {
+            *a <<= 8, *b <<= 8;
+            *a += d[j];
+            *b += d[j + 1];
+            j -= 2;
+        }
+    }
+    return SIZE(e);
+}
+
+static void field_clear_gf3m(field_t f) {
+    params *p = f->data;
+    gf3m_clear(p->p);
+    pbc_free(p->p);
+    pbc_free(p);
+}
+
+/* initialize the finite field as $GF(3^m)$, whose irreducible polynomial is with the degree of $m$ */
+void field_init_gf3m(field_t f, unsigned m, unsigned t) {
+    params *p = pbc_malloc(sizeof(*p));
+    p->len = (m + (W - 1) + 1) / W; /* extra one bit for $_p$ */
+    p->m = m;
+    p->t = t;
+    p->p = pbc_malloc(sizeof(*(p->p)));
+    p->p->field = f;
+    p->p->data = pbc_malloc(2 * sizeof(unsigned long) * p->len);
+    memset(p->p->data, 0, 2 * sizeof(unsigned long) * p->len);
+    unsigned long *p1 = p->p->data, *p2 = p1 + p->len;
+    p2[0] = 1; /* _p == x^m+x^t+2 */
+    unsigned int p_t = p->t;
+    p1[p_t / W] |= 1ul << (p_t % W);
+    p1[m / W] |= 1ul << (m % W);
+
+    field_init(f);
+    f->field_clear = field_clear_gf3m;
+    f->init = gf3m_init;
+    f->clear = gf3m_clear;
+    f->set = gf3m_assign;
+    f->set0 = gf3m_zero;
+    f->set1 = gf3m_one;
+    f->is0 = gf3m_is0;
+    f->is1 = gf3m_is1;
+    f->add = gf3m_add;
+    f->sub = gf3m_sub;
+    f->mul = gf3m_mult;
+    f->cubic = gf3m_cubic;
+    f->invert = gf3m_invert;
+    f->neg = gf3m_neg;
+    f->random = gf3m_random;
+    f->cmp = gf3m_cmp;
+    f->sqrt = gf3m_sqrt;
+    f->from_bytes = gf3m_from_bytes;
+    f->to_bytes = gf3m_to_bytes;
+    f->out_str = gf3m_out_str;
+    f->fixed_length_in_bytes = 2 * sizeof(unsigned long) * p->len;
+    f->data = p;
+    f->name = "GF(3^m)";
+
+    mpz_set_ui(f->order, 3);
+    mpz_pow_ui(f->order, f->order, p->m);
+}
+
+static size_t gf32m_out_str(FILE *stream, int base, element_t e) {
+    UNUSED_VAR(base);
+    element_ptr e0 = GF32M(e)->_0, e1 = GF32M(e)->_1;
+    size_t size = 0;
+    size += element_out_str(stream, base, e0);
+    size += element_out_str(stream, base, e1);
+    return size;
+}
+
+static void gf32m_init(element_t e) {
+    e->data = pbc_malloc(sizeof(gf32m_s));
+    gf32m_ptr p = (gf32m_ptr) e->data;
+    field_ptr base = BASE(e);
+    element_init(p->_0, base);
+    element_init(p->_1, base);
+}
+
+static void gf32m_clear(element_t e) {
+    gf32m_ptr p = (gf32m_ptr) e->data;
+    element_clear(p->_0);
+    element_clear(p->_1);
+    pbc_free(e->data);
+}
+
+static void gf32m_set0(element_t e) {
+    element_ptr e0 = GF32M(e)->_0, e1 = GF32M(e)->_1;
+    element_set0(e0);
+    element_set0(e1);
+}
+
+static void gf32m_set1(element_t e) {
+    element_ptr e0 = GF32M(e)->_0, e1 = GF32M(e)->_1;
+    element_set1(e0);
+    element_set0(e1);
+}
+
+static int gf32m_item_count(element_t e) {
+    UNUSED_VAR(e);
+    return 2;
+}
+
+static element_ptr gf32m_item(element_t a, int i) {
+    element_ptr a0 = GF32M(a)->_0, a1 = GF32M(a)->_1;
+    return i == 0 ? a0 : a1;
+}
+
+static void gf32m_assign(element_t e, element_t a) {
+    element_ptr a0 = GF32M(a)->_0, a1 = GF32M(a)->_1, e0 = GF32M(e)->_0, e1 =
+            GF32M(e)->_1;
+    element_set(e0, a0);
+    element_set(e1, a1);
+}
+
+static void gf32m_random(element_t e) {
+    element_ptr e0 = GF32M(e)->_0, e1 = GF32M(e)->_1;
+    element_random(e0);
+    element_random(e1);
+}
+
+/* return 0 if $a == b$, 1 otherwise */
+static int gf32m_cmp(element_t a, element_t b) {
+    element_ptr a0 = GF32M(a)->_0, a1 = GF32M(a)->_1, b0 = GF32M(b)->_0, b1 =
+            GF32M(b)->_1;
+    return element_cmp(a0, b0) || element_cmp(a1, b1);
+}
+
+/* $c <- a+b$ */
+static void gf32m_add(element_t c, element_t a, element_t b) {
+    element_ptr a0 = GF32M(a)->_0, a1 = GF32M(a)->_1, b0 = GF32M(b)->_0, b1 =
+            GF32M(b)->_1, c0 = GF32M(c)->_0, c1 = GF32M(c)->_1;
+    element_add(c0, a0, b0);
+    element_add(c1, a1, b1);
+}
+
+/* $c <- a-b$ */
+static void gf32m_sub(element_t c, element_t a, element_t b) {
+    element_ptr a0 = GF32M(a)->_0, a1 = GF32M(a)->_1, b0 = GF32M(b)->_0, b1 =
+            GF32M(b)->_1, c0 = GF32M(c)->_0, c1 = GF32M(c)->_1;
+    element_sub(c0, a0, b0);
+    element_sub(c1, a1, b1);
+}
+
+/* $c <- (-a)$ */
+static void gf32m_neg(element_t c, element_t a) {
+    element_ptr a0 = GF32M(a)->_0, a1 = GF32M(a)->_1, c0 = GF32M(c)->_0, c1 =
+            GF32M(c)->_1;
+    element_neg(c0, a0);
+    element_neg(c1, a1);
+}
+
+/* $e<- a*b$ */
+static void gf32m_mult(element_t e, element_t a, element_t b) {
+    element_ptr a0 = GF32M(a)->_0, a1 = GF32M(a)->_1, b0 = GF32M(b)->_0, b1 =
+            GF32M(b)->_1, e0 = GF32M(e)->_0, e1 = GF32M(e)->_1;
+    field_ptr base = BASE(a);
+    element_t a0b0, a1b1, t0, t1, c1;
+    element_init(a0b0, base);
+    element_init(a1b1, base);
+    element_init(t0, base);
+    element_init(t1, base);
+    element_init(c1, base);
+    element_mul(a0b0, a0, b0);
+    element_mul(a1b1, a1, b1);
+    element_add(t0, a1, a0);
+    element_add(t1, b1, b0);
+    element_mul(c1, t0, t1); // c1 == (a1+a0)*(b1+b0)
+    element_sub(c1, c1, a1b1);
+    element_sub(c1, c1, a0b0);
+    element_ptr c0 = a0b0;
+    element_sub(c0, c0, a1b1); // c0 == a0*b0 - a1*b1
+    element_set(e0, c0);
+    element_set(e1, c1);
+    element_clear(a0b0);
+    element_clear(a1b1);
+    element_clear(t0);
+    element_clear(t1);
+    element_clear(c1);
+}
+
+/* $e <- a^3$ */
+static void gf32m_cubic(element_t e, element_t a) {
+    element_ptr a0 = GF32M(a)->_0, a1 = GF32M(a)->_1, e0 = GF32M(e)->_0, e1 =
+            GF32M(e)->_1;
+    field_ptr base = BASE(a);
+    element_t c0, c1;
+    element_init(c0, base);
+    element_init(c1, base);
+    element_cubic(c0, a0);
+    element_cubic(c1, a1);
+    element_neg(c1, c1); // c1 == -(a1^3)
+    element_set(e0, c0);
+    element_set(e1, c1);
+    element_clear(c0);
+    element_clear(c1);
+}
+
+void field_clear_gf32m(field_t f) {
+    UNUSED_VAR(f);
+}
+
+/* initialize the finite field as $base_field[x]/(x^2 + 1)$, whose base field is $b$ */
+void field_init_gf32m(field_t f, field_t b) {
+    field_init(f);
+    f->data = b;
+    f->field_clear = field_clear_gf32m;
+    f->init = gf32m_init;
+    f->clear = gf32m_clear;
+    f->set = gf32m_assign;
+    f->set0 = gf32m_set0;
+    f->set1 = gf32m_set1;
+    f->random = gf32m_random;
+    f->cmp = gf32m_cmp;
+    f->add = gf32m_add;
+    f->sub = gf32m_sub;
+    f->neg = gf32m_neg;
+    f->mul = gf32m_mult;
+    f->cubic = gf32m_cubic;
+    f->item_count = gf32m_item_count;
+    f->item = gf32m_item;
+    f->out_str = gf32m_out_str;
+    mpz_pow_ui(f->order, b->order, 2);
+    f->name = "GF(3^{2*m})";
+}
+
+static size_t gf33m_out_str(FILE *stream, int base, element_t e) {
+    UNUSED_VAR(base);
+    element_ptr e0 = GF33M(e)->_0, e1 = GF33M(e)->_1, e2 = GF33M(e)->_2;
+    size_t size = 0;
+    size += element_out_str(stream, base, e0);
+    size += element_out_str(stream, base, e1);
+    size += element_out_str(stream, base, e2);
+    return size;
+}
+
+static void gf33m_init(element_t e) {
+    e->data = pbc_malloc(sizeof(gf33m_s));
+    gf33m_ptr p = (gf33m_ptr) e->data;
+    field_ptr base = BASE(e);
+    element_init(p->_0, base);
+    element_init(p->_1, base);
+    element_init(p->_2, base);
+}
+
+static void gf33m_clear(element_t e) {
+    gf33m_ptr p = (gf33m_ptr) e->data;
+    element_clear(p->_0);
+    element_clear(p->_1);
+    element_clear(p->_2);
+    pbc_free(e->data);
+}
+
+static void gf33m_set0(element_t e) {
+    element_ptr e0 = GF33M(e)->_0, e1 = GF33M(e)->_1, e2 = GF33M(e)->_2;
+    element_set0(e0);
+    element_set0(e1);
+    element_set0(e2);
+}
+
+static void gf33m_set1(element_t e) {
+    element_ptr e0 = GF33M(e)->_0, e1 = GF33M(e)->_1, e2 = GF33M(e)->_2;
+    element_set1(e0);
+    element_set0(e1);
+    element_set0(e2);
+}
+
+static int gf33m_item_count(element_t e) {
+    UNUSED_VAR(e);
+    return 3;
+}
+
+static element_ptr gf33m_item(element_t a, int i) {
+    element_ptr a0 = GF33M(a)->_0, a1 = GF33M(a)->_1, a2 = GF33M(a)->_2;
+    return i == 0 ? a0 : (i == 1 ? a1 : a2);
+}
+
+static void gf33m_assign(element_t e, element_t a) {
+    element_ptr a0 = GF33M(a)->_0, a1 = GF33M(a)->_1, a2 = GF33M(a)->_2, e0 =
+            GF33M(e)->_0, e1 = GF33M(e)->_1, e2 = GF33M(e)->_2;
+    element_set(e0, a0);
+    element_set(e1, a1);
+    element_set(e2, a2);
+}
+
+static void gf33m_random(element_t e) {
+    element_ptr e0 = GF33M(e)->_0, e1 = GF33M(e)->_1, e2 = GF33M(e)->_2;
+    element_random(e0);
+    element_random(e1);
+    element_random(e2);
+}
+
+/* return 0 if $a == b$, 1 otherwise */
+static int gf33m_cmp(element_t a, element_t b) {
+    element_ptr a0 = GF33M(a)->_0, a1 = GF33M(a)->_1, a2 = GF33M(a)->_2, b0 =
+            GF33M(b)->_0, b1 = GF33M(b)->_1, b2 = GF33M(b)->_2;
+    return element_cmp(a0, b0) || element_cmp(a1, b1) || element_cmp(a2, b2);
+}
+
+/* $c <- a+b$ */
+static void gf33m_add(element_t c, element_t a, element_t b) {
+    element_ptr a0 = GF33M(a)->_0, a1 = GF33M(a)->_1, a2 = GF33M(a)->_2, b0 =
+            GF33M(b)->_0, b1 = GF33M(b)->_1, b2 = GF33M(b)->_2, c0 =
+            GF33M(c)->_0, c1 = GF33M(c)->_1, c2 = GF33M(c)->_2;
+    element_add(c0, a0, b0);
+    element_add(c1, a1, b1);
+    element_add(c2, a2, b2);
+}
+
+/* $c <- a-b$ */
+static void gf33m_sub(element_t c, element_t a, element_t b) {
+    element_ptr a0 = GF33M(a)->_0, a1 = GF33M(a)->_1, a2 = GF33M(a)->_2, b0 =
+            GF33M(b)->_0, b1 = GF33M(b)->_1, b2 = GF33M(b)->_2, c0 =
+            GF33M(c)->_0, c1 = GF33M(c)->_1, c2 = GF33M(c)->_2;
+    element_sub(c0, a0, b0);
+    element_sub(c1, a1, b1);
+    element_sub(c2, a2, b2);
+}
+
+/* $c <- a*b$ */
+static void gf33m_mult(element_t e, element_t a, element_t b) {
+    element_ptr a0 = GF33M(a)->_0, a1 = GF33M(a)->_1, a2 = GF33M(a)->_2, b0 =
+            GF33M(b)->_0, b1 = GF33M(b)->_1, b2 = GF33M(b)->_2, e0 =
+            GF33M(e)->_0, e1 = GF33M(e)->_1, e2 = GF33M(e)->_2;
+    field_ptr base = BASE(e);
+    element_t t0, t1, c1, a0b0, a1b1, a2b2;
+    element_init(t0, base);
+    element_init(t1, base);
+    element_init(c1, base);
+    element_init(a0b0, base);
+    element_init(a1b1, base);
+    element_init(a2b2, base);
+    element_mul(a0b0, a0, b0);
+    element_mul(a1b1, a1, b1);
+    element_mul(a2b2, a2, b2);
+    element_ptr d0 = a0b0;
+    element_add(t0, a1, a0);
+    element_add(t1, b1, b0);
+    element_t d1;
+    element_init(d1, base);
+    element_mul(d1, t0, t1);
+    element_sub(d1, d1, a1b1);
+    element_sub(d1, d1, a0b0);
+    element_add(t0, a2, a0);
+    element_add(t1, b2, b0);
+    element_t d2;
+    element_init(d2, base);
+    element_mul(d2, t0, t1);
+    element_add(d2, d2, a1b1);
+    element_sub(d2, d2, a2b2);
+    element_sub(d2, d2, a0b0);
+    element_add(t0, a2, a1);
+    element_add(t1, b2, b1);
+    element_t d3;
+    element_init(d3, base);
+    element_mul(d3, t0, t1);
+    element_sub(d3, d3, a2b2);
+    element_sub(d3, d3, a1b1);
+    element_ptr d4 = a2b2;
+    element_add(t0, d0, d3);
+    element_ptr c0 = t0;
+    element_add(c1, d1, d3);
+    element_add(c1, c1, d4);
+    element_add(t1, d2, d4);
+    element_ptr c2 = t1;
+    element_set(e0, c0);
+    element_set(e1, c1);
+    element_set(e2, c2);
+    element_clear(t0);
+    element_clear(t1);
+    element_clear(c1);
+    element_clear(a0b0);
+    element_clear(a1b1);
+    element_clear(a2b2);
+    element_clear(d1);
+    element_clear(d2);
+    element_clear(d3);
+}
+
+/* $e <- a^3$ */
+static void gf33m_cubic(element_t e, element_t a) {
+    field_ptr base = BASE(a);
+    element_ptr a0 = GF33M(a)->_0, a1 = GF33M(a)->_1, a2 = GF33M(a)->_2, e0 =
+            GF33M(e)->_0, e1 = GF33M(e)->_1, e2 = GF33M(e)->_2;
+    element_t a03, a13, a23;
+    element_init(a03, base);
+    element_init(a13, base);
+    element_init(a23, base);
+    element_cubic(a03, a0);
+    element_cubic(a13, a1);
+    element_cubic(a23, a2);
+    element_add(a03, a03, a13);
+    element_add(a03, a03, a23);
+    element_ptr c0 = a03;
+    element_sub(a13, a13, a23);
+    element_ptr c1 = a13;
+    element_ptr c2 = a23;
+    element_set(e0, c0);
+    element_set(e1, c1);
+    element_set(e2, c2);
+    element_clear(a03);
+    element_clear(a13);
+    element_clear(a23);
+}
+
+/* $e <- a^{-1}$ */
+static void gf33m_invert(element_t e, element_t a) {
+    element_ptr a0 = GF33M(a)->_0, a1 = GF33M(a)->_1, a2 = GF33M(a)->_2, e0 =
+            GF33M(e)->_0, e1 = GF33M(e)->_1, e2 = GF33M(e)->_2;
+    field_ptr base = BASE(e);
+    element_t a02, a12, a22;
+    element_init(a02, base);
+    element_init(a12, base);
+    element_init(a22, base);
+    element_mul(a02, a0, a0);
+    element_mul(a12, a1, a1);
+    element_mul(a22, a2, a2);
+    element_t v0;
+    element_init(v0, base);
+    element_sub(v0, a0, a2); // v0 == a0-a2
+    element_t delta;
+    element_init(delta, base);
+    element_mul(delta, v0, a02); // delta = (a0-a2)*(a0^2), free
+    element_sub(v0, a1, a0); // v0 == a1-a0
+    element_t c0;
+    element_init(c0, base);
+    element_mul(c0, v0, a12); // c0 == (a1-a0)*(a1^2)
+    element_add(delta, delta, c0); // delta = (a0-a2)*(a0^2) + (a1-a0)*(a1^2)
+    element_sub(v0, a2, v0); // v0 == a2-(a1-a0) = a0-a1+a2
+    element_t c1;
+    element_init(c1, base);
+    element_mul(c1, v0, a22); // c1 == (a0-a1+a2)*(a2^2)
+    element_add(delta, delta, c1); // delta = (a0-a2)*(a0^2) + (a1-a0)*(a1^2) + (a0-a1+a2)*(a2^2)
+    element_invert(delta, delta); // delta = [(a0-a2)*(a0^2) + (a1-a0)*(a1^2) + (a0-a1+a2)*(a2^2)] ^ {-1}
+    element_add(v0, a02, a22); // v0 == a0^2+a2^2
+    element_t c2;
+    element_init(c2, base);
+    element_mul(c2, a0, a2); // c2 == a0*a2
+    element_sub(c0, v0, c2); // c0 == a0^2+a2^2-a0*a2
+    element_add(v0, a1, a2); // v0 == a1+a2
+    element_t c3;
+    element_init(c3, base);
+    element_mul(c3, a1, v0); // c3 == a1*(a1+a2)
+    element_sub(c0, c0, c3); // c0 == a0^2+a2^2-a0*a2-a1*(a1+a2)
+    element_mul(c0, c0, delta); // c0 *= delta
+    element_mul(c1, a0, a1); // c1 == a0*a1
+    element_sub(c1, a22, c1); // c1 == a2^2-a0*a1
+    element_mul(c1, c1, delta); // c1 *= delta
+    element_sub(c2, a12, c2); // c2 == a1^2-a0*a2
+    element_sub(c2, c2, a22); // c2 == a1^2-a0*a2-a2^2
+    element_mul(c2, c2, delta); // c2 *= delta
+    element_set(e0, c0);
+    element_set(e1, c1);
+    element_set(e2, c2);
+    element_clear(a02);
+    element_clear(a12);
+    element_clear(a22);
+    element_clear(v0);
+    element_clear(delta);
+    element_clear(c0);
+    element_clear(c1);
+    element_clear(c2);
+    element_clear(c3);
+}
+
+void field_clear_gf33m(field_t f) {
+    UNUSED_VAR(f);
+}
+
+/* initialize the finite field as $base_field[x]/(x^3 - x - 1)$, whose base field is $b$ */
+void field_init_gf33m(field_t f, field_t b) {
+    field_init(f);
+    f->data = b;
+    f->field_clear = field_clear_gf33m;
+    f->init = gf33m_init;
+    f->clear = gf33m_clear;
+    f->set = gf33m_assign;
+    f->set0 = gf33m_set0;
+    f->set1 = gf33m_set1;
+    f->random = gf33m_random;
+    f->cmp = gf33m_cmp;
+    f->add = gf33m_add;
+    f->sub = gf33m_sub;
+    f->mul = gf33m_mult;
+    f->cubic = gf33m_cubic;
+    f->invert = gf33m_invert;
+    f->item_count = gf33m_item_count;
+    f->item = gf33m_item;
+    f->out_str = gf33m_out_str;
+    mpz_pow_ui(f->order, b->order, 3);
+    f->name = "GF(3^{3*m})";
+}
+
diff --git a/moon-abe/pbc-0.5.14/arith/tinyfp.c b/moon-abe/pbc-0.5.14/arith/tinyfp.c
new file mode 100644
index 00000000..50e883e1
--- /dev/null
+++ b/moon-abe/pbc-0.5.14/arith/tinyfp.c
@@ -0,0 +1,304 @@
+// F_p for small p, i.e. at most sizeof(long) bytes long.
+// Assumes long long is at least twice long.
+
+// TODO: Fix outstanding bugs and use in PBC.
+
+#include <stdarg.h>
+#include <stdio.h>
+#include <stdint.h> // for intptr_t
+#include <stdlib.h>
+#include <string.h>
+#include <gmp.h>
+#include "pbc_utils.h"
+#include "pbc_field.h"
+#include "pbc_random.h"
+#include "pbc_fp.h"
+#include "pbc_memory.h"
+
+// Mostly wrappers. We use GMP routines for pow_mpz and invert.
+
+static void fp_init(element_ptr e) {
+  unsigned long *p = e->data = pbc_malloc(sizeof(unsigned long));
+  *p = 0;
+}
+
+static void fp_clear(element_ptr e) {
+  pbc_free(e->data);
+}
+
+static void fp_set_mpz(element_ptr e, mpz_ptr z) {
+  mpz_t r;
+  mpz_init(r);
+  unsigned long *p = e->field->data;
+  unsigned long *l = e->data;
+  mpz_fdiv_r_ui(r, z, *p);
+  *l = mpz_get_ui(r);
+  mpz_clear(r);
+}
+
+static void fp_set_si(element_ptr e, signed long int op) {
+  unsigned long int *d = e->data;
+  unsigned long *p = e->field->data;
+  if (op < 0) {
+    *d = (-op) % *p;
+    *d = *p - *d;
+  } else {
+    *d = op % *p;
+  }
+}
+
+static void fp_to_mpz(mpz_ptr z, element_ptr e) {
+  unsigned long int *l = e->data;
+  mpz_set_ui(z, *l);
+}
+
+static void fp_set0(element_ptr e) {
+  unsigned long int *l = e->data;
+  *l = 0;
+}
+
+static void fp_set1(element_ptr e) {
+  unsigned long int *l = e->data;
+  *l = 1;
+}
+
+static int fp_is1(element_ptr e) {
+  unsigned long int *l = e->data;
+  return *l == 1;
+}
+
+static int fp_is0(element_ptr e) {
+  unsigned long int *l = e->data;
+  return *l == 0;
+}
+
+static size_t fp_out_str(FILE *stream, int base, element_ptr e) {
+  size_t result;
+  mpz_t z;
+  mpz_init(z);
+  fp_to_mpz(z, e);
+  result = mpz_out_str(stream, base, z);
+  mpz_clear(z);
+  return result;
+}
+
+static void fp_add(element_ptr c, element_ptr a, element_ptr b) {
+  unsigned long *prime = a->field->data;
+  unsigned long *p = a->data;
+  unsigned long *q = b->data;
+  unsigned long *r = c->data;
+  unsigned long l0;
+  l0 = *p + *q;
+  if (l0 < *p) {
+    //overflow
+    l0 -= *prime;
+  }
+  *r = l0 % *prime;
+}
+
+static void fp_double(element_ptr c, element_ptr a) {
+  unsigned long *prime = a->field->data;
+  unsigned long *p = a->data;
+  unsigned long *r = c->data;
+  *r = 2 * *p;
+  if (*r < *p) {
+    //overflow
+    *r -= *prime;
+  }
+  *r = *r % *prime;
+}
+
+static void fp_sub(element_ptr c, element_ptr a, element_ptr b) {
+  unsigned long *prime = a->field->data;
+  unsigned long *p = a->data;
+  unsigned long *q = b->data;
+  unsigned long *r = c->data;
+
+  if (*p >= *q) {
+    *r = *p - *q;
+  } else {
+    *r = *prime - *q + *p;
+  }
+}
+
+static void fp_mul(element_ptr c, element_ptr a, element_ptr b) {
+  unsigned long *prime = a->field->data;
+  unsigned long *p = a->data;
+  unsigned long *q = b->data;
+  unsigned long long ll;
+  unsigned long *r = c->data;
+
+  ll = *p * *q;
+  *r = ll % *prime;
+}
+
+static void fp_square(element_ptr c, element_ptr a) {
+  fp_mul(c, a, a);
+}
+
+static void fp_neg(element_ptr c, element_ptr a) {
+  unsigned long *prime = a->field->data;
+  unsigned long *r = c->data;
+  unsigned long *p = a->data;
+  if (*p) {
+    *r = *prime - *p;
+  } else {
+    *r = 0;
+  }
+}
+
+static void fp_mul_si(element_ptr c, element_ptr a, signed long int op) {
+  unsigned long *prime = a->field->data;
+  unsigned long *p = a->data;
+  unsigned long long ll;
+  unsigned long *r = c->data;
+
+  ll = *p * op;
+  *r = ll % *prime;
+}
+
+static void fp_pow_mpz(element_ptr c, element_ptr a, mpz_ptr op) {
+  unsigned long *r = c->data;
+  mpz_t z;
+  mpz_init(z);
+  fp_to_mpz(z, a);
+  mpz_powm(z, z, op, a->field->order);
+  *r = mpz_get_ui(z);
+  mpz_clear(z);
+}
+
+static void fp_set(element_ptr c, element_ptr a) {
+  unsigned long *p = a->data;
+  unsigned long *r = c->data;
+  *r = *p;
+}
+
+static void fp_invert(element_ptr c, element_ptr a) {
+  unsigned long *r = c->data;
+  mpz_t z;
+  mpz_init(z);
+  fp_to_mpz(z, a);
+  mpz_invert(z, z, a->field->order);
+  *r = mpz_get_ui(z);
+  mpz_clear(z);
+}
+
+static void fp_random(element_ptr c) {
+  unsigned long *r = c->data;
+  mpz_t z;
+  mpz_init(z);
+  pbc_mpz_random(z, c->field->order);
+  *r = mpz_get_ui(z);
+  mpz_clear(z);
+}
+
+static void fp_from_hash(element_ptr n, void *data, int len) {
+  mpz_t z;
+
+  mpz_init(z);
+  mpz_import(z, len, -1, 1, -1, 0, data);
+  fp_set_mpz(n, z);
+  mpz_clear(z);
+}
+
+static int fp_cmp(element_ptr a, element_ptr b) {
+  unsigned long *p = a->data;
+  unsigned long *q = b->data;
+  return *p != *q;
+}
+
+static int fp_sgn_odd(element_ptr a) {
+  unsigned long *p = a->data;
+  if (!*p) return 0;
+  return *p & 1 ? 1 : -1;
+}
+
+static int fp_is_sqr(element_ptr a) {
+  int res;
+  mpz_t z;
+  mpz_init(z);
+  //0 is a square
+  if (fp_is0(a)) return 1;
+  fp_to_mpz(z, a);
+  res = mpz_legendre(z, a->field->order) == 1;
+  mpz_clear(z);
+  return res;
+}
+
+static int fp_to_bytes(unsigned char *data, element_t e) {
+  unsigned long *p = e->data;
+  unsigned long l = *p;
+  int i, n = e->field->fixed_length_in_bytes;
+  for (i = 0; i < n; i++) {
+    data[n - i - 1] = (unsigned char) l;
+    l >>= 8;
+  }
+  return n;
+}
+
+static int fp_from_bytes(element_t e, unsigned char *data) {
+  unsigned char *ptr = data;
+  unsigned long *p = e->data;
+  int i, n = e->field->fixed_length_in_bytes;
+  *p = 0;
+  for (i=0; i<n; i++) {
+    *p <<= 8;
+    *p += *ptr;
+    ptr++;
+  }
+  return n;
+}
+
+static void fp_field_clear(field_t f) {
+  pbc_free(f->data);
+}
+
+void field_init_tiny_fp(field_ptr f, mpz_t prime) {
+  unsigned long *p;
+
+  PBC_ASSERT(mpz_fits_ulong_p(prime), "modulus too big");
+
+  field_init(f);
+  f->init = fp_init;
+  f->clear = fp_clear;
+  f->set_si = fp_set_si;
+  f->set_mpz = fp_set_mpz;
+  f->out_str = fp_out_str;
+  f->add = fp_add;
+  f->sub = fp_sub;
+  f->set = fp_set;
+  f->mul = fp_mul;
+  f->mul_si = fp_mul_si;
+  f->square = fp_square;
+  f->doub = fp_double;
+  f->pow_mpz = fp_pow_mpz;
+  f->neg = fp_neg;
+  f->cmp = fp_cmp;
+  f->sign = fp_sgn_odd;
+  f->invert = fp_invert;
+  f->random = fp_random;
+  f->from_hash = fp_from_hash;
+  f->is1 = fp_is1;
+  f->is0 = fp_is0;
+  f->set0 = fp_set0;
+  f->set1 = fp_set1;
+  f->is_sqr = fp_is_sqr;
+  f->sqrt = element_tonelli;
+  f->field_clear = fp_field_clear;
+  f->to_bytes = fp_to_bytes;
+  f->from_bytes = fp_from_bytes;
+  f->to_mpz = fp_to_mpz;
+
+  p = f->data = pbc_malloc(sizeof(long));
+  *p = mpz_get_ui(prime);
+  {
+    unsigned long int l = 255;
+    f->fixed_length_in_bytes = 1;
+    while (l < *p) {
+      f->fixed_length_in_bytes++;
+      l <<= 8;
+      l += 255;
+    }
+  }
+  mpz_set(f->order, prime);
+}
diff --git a/moon-abe/pbc-0.5.14/arith/z.c b/moon-abe/pbc-0.5.14/arith/z.c
new file mode 100644
index 00000000..ff5a4a97
--- /dev/null
+++ b/moon-abe/pbc-0.5.14/arith/z.c
@@ -0,0 +1,263 @@
+// The ring Z.
+//
+// Wrappers around GMP mpz functions.
+#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_z.h"
+#include "pbc_random.h"
+#include "pbc_fp.h"
+#include "pbc_memory.h"
+
+static void z_init(element_ptr e) {
+  e->data = pbc_malloc(sizeof(mpz_t));
+  mpz_init(e->data);
+}
+
+static void z_clear(element_ptr e) {
+  mpz_clear(e->data);
+  pbc_free(e->data);
+}
+
+static void z_set_si(element_ptr e, signed long int op) {
+  mpz_set_si(e->data, op);
+}
+
+static void z_set_mpz(element_ptr e, mpz_ptr z) {
+  mpz_set(e->data, z);
+}
+
+static void z_set0(element_ptr e) {
+  mpz_set_ui(e->data, 0);
+}
+
+static void z_set1(element_ptr e) {
+  mpz_set_ui(e->data, 1);
+}
+
+static size_t z_out_str(FILE *stream, int base, element_ptr e) {
+  return mpz_out_str(stream, base, e->data);
+}
+
+static int z_sgn(element_ptr a) {
+  mpz_ptr z = a->data;
+  return mpz_sgn(z);
+}
+
+static void z_add(element_ptr n, element_ptr a, element_ptr b) {
+  mpz_add(n->data, a->data, b->data);
+}
+
+static void z_sub(element_ptr n, element_ptr a, element_ptr b) {
+  mpz_sub(n->data, a->data, b->data);
+}
+
+static void z_square(element_ptr c, element_ptr a) {
+  mpz_mul(c->data, a->data, a->data);
+}
+
+static void z_double(element_ptr n, element_ptr a) {
+  mpz_mul_2exp(n->data, a->data, 1);
+}
+
+static void z_halve(element_ptr n, element_ptr a) {
+  mpz_tdiv_q_2exp(n->data, a->data, -1);
+}
+
+static void z_mul(element_ptr n, element_ptr a, element_ptr b) {
+  mpz_mul(n->data, a->data, b->data);
+}
+
+static void z_mul_mpz(element_ptr n, element_ptr a, mpz_ptr z) {
+  mpz_mul(n->data, a->data, z);
+}
+
+static void z_mul_si(element_ptr n, element_ptr a, signed long int z) {
+  mpz_mul_si(n->data, a->data, z);
+}
+
+static void z_pow_mpz(element_ptr n, element_ptr a, mpz_ptr z) {
+  mpz_pow_ui(n->data, a->data, mpz_get_ui(z));
+}
+
+static void z_set(element_ptr n, element_ptr a) {
+  mpz_set(n->data, a->data);
+}
+
+static void z_neg(element_ptr n, element_ptr a) {
+  mpz_neg(n->data, a->data);
+}
+
+static void z_invert(element_ptr n, element_ptr a) {
+  if (!mpz_cmpabs_ui(a->data, 1)) {
+    mpz_set(n->data, a->data);
+  } else mpz_set_ui(n->data, 0);
+}
+
+static void z_div(element_ptr c, element_ptr a, element_ptr b) {
+  mpz_tdiv_q(c->data, a->data, b->data);
+}
+
+//(doesn't make sense if order is infinite)
+static void z_random(element_ptr n) {
+  mpz_set_ui(n->data, 0);
+}
+
+static void z_from_hash(element_ptr n, void *data, int len) {
+  mpz_import(n->data, len, -1, 1, -1, 0, data);
+}
+
+static int z_is1(element_ptr n) {
+  return !mpz_cmp_ui((mpz_ptr) n->data, 1);
+}
+
+static int z_is0(element_ptr n) {
+  return mpz_is0(n->data);
+}
+
+static int z_cmp(element_ptr a, element_ptr b) {
+  return mpz_cmp((mpz_ptr) a->data, (mpz_ptr) b->data);
+}
+
+static int z_is_sqr(element_ptr a) {
+  return mpz_perfect_power_p(a->data);
+}
+
+static void z_sqrt(element_ptr c, element_ptr a) {
+  mpz_sqrt(c->data, a->data);
+}
+
+static void z_field_clear(field_t f) {
+  UNUSED_VAR (f);
+}
+
+// OpenSSL convention:
+//   4 bytes containing length
+//   followed by number in big-endian, most-significant bit set if negative
+//   (prepending null byte if necessary)
+// Positive numbers also the same as mpz_out_raw.
+static int z_to_bytes(unsigned char *data, element_t e) {
+  mpz_ptr z = e->data;
+  size_t msb = mpz_sizeinbase(z, 2);
+  size_t n = 4;
+  size_t i;
+
+  if (!(msb % 8)) {
+    data[4] = 0;
+    n++;
+  }
+  if (mpz_sgn(z) < 0) {
+    mpz_export(data + n, NULL, 1, 1, 1, 0, z);
+    data[4] |= 128;
+  } else {
+    mpz_export(data + n, NULL, 1, 1, 1, 0, z);
+  }
+  n += (msb + 7) / 8 - 4;
+  for (i=0; i<4; i++) {
+    data[i] = (n >> 8 * (3 - i));
+  }
+  n += 4;
+
+  return n;
+}
+
+static int z_from_bytes(element_t e, unsigned char *data) {
+  unsigned char *ptr;
+  size_t i, n;
+  mpz_ptr z = e->data;
+  mpz_t z1;
+  int neg = 0;
+
+  mpz_init(z1);
+  mpz_set_ui(z, 0);
+
+  ptr = data;
+  n = 0;
+  for (i=0; i<4; i++) {
+    n += ((unsigned int) *ptr) << 8 * (3 - i);
+    ptr++;
+  }
+  if (data[4] & 128) {
+    neg = 1;
+    data[4] &= 127;
+  }
+  for (i=0; i<n; i++) {
+    mpz_set_ui(z1, *ptr);
+    mpz_mul_2exp(z1, z1, 8 * (n - 1 - i));
+    ptr++;
+    mpz_add(z, z, z1);
+  }
+  mpz_clear(z1);
+  if (neg) mpz_neg(z, z);
+  return n;
+}
+
+static void z_to_mpz(mpz_ptr z, element_ptr a) {
+  mpz_set(z, a->data);
+}
+
+static int z_length_in_bytes(element_ptr a) {
+  return (mpz_sizeinbase(a->data, 2) + 7) / 8 + 4;
+}
+
+static void z_out_info(FILE *out, field_ptr f) {
+  UNUSED_VAR(f);
+  fprintf(out, "Z: wrapped GMP");
+}
+
+static int z_set_str(element_ptr e, const char *s, int base) {
+  mpz_t z;
+  mpz_init(z);
+  int result = pbc_mpz_set_str(z, s, base);
+  z_set_mpz(e, z);
+  mpz_clear(z);
+  return result;
+}
+
+void field_init_z(field_ptr f) {
+  field_init(f);
+  f->init = z_init;
+  f->clear = z_clear;
+  f->set_si = z_set_si;
+  f->set_mpz = z_set_mpz;
+  f->set_str = z_set_str;
+  f->out_str = z_out_str;
+  f->sign = z_sgn;
+  f->add = z_add;
+  f->sub = z_sub;
+  f->set = z_set;
+  f->square = z_square;
+  f->doub = z_double;
+  f->halve = z_halve;
+  f->mul = z_mul;
+  f->mul_mpz = z_mul_mpz;
+  f->mul_si = z_mul_si;
+  f->pow_mpz = z_pow_mpz;
+  f->neg = z_neg;
+  f->cmp = z_cmp;
+  f->invert = z_invert;
+  f->div = z_div;
+  f->random = z_random;
+  f->from_hash = z_from_hash;
+  f->is1 = z_is1;
+  f->is0 = z_is0;
+  f->set0 = z_set0;
+  f->set1 = z_set1;
+  f->is_sqr = z_is_sqr;
+  f->sqrt = z_sqrt;
+  f->field_clear = z_field_clear;
+  f->to_bytes = z_to_bytes;
+  f->from_bytes = z_from_bytes;
+  f->to_mpz = z_to_mpz;
+  f->length_in_bytes = z_length_in_bytes;
+
+  f->out_info = z_out_info;
+
+  mpz_set_ui(f->order, 0);
+  f->data = NULL;
+  f->fixed_length_in_bytes = -1;
+}