diff options
Diffstat (limited to 'moon-abe/pbc-0.5.14/arith')
| -rw-r--r-- | moon-abe/pbc-0.5.14/arith/dlog.c | 187 | ||||
| -rw-r--r-- | moon-abe/pbc-0.5.14/arith/fasterfp.c | 546 | ||||
| -rw-r--r-- | moon-abe/pbc-0.5.14/arith/fastfp.c | 382 | ||||
| -rw-r--r-- | moon-abe/pbc-0.5.14/arith/field.c | 889 | ||||
| -rw-r--r-- | moon-abe/pbc-0.5.14/arith/fieldquadratic.c | 692 | ||||
| -rw-r--r-- | moon-abe/pbc-0.5.14/arith/fp.c | 49 | ||||
| -rw-r--r-- | moon-abe/pbc-0.5.14/arith/init_random.c | 18 | ||||
| -rw-r--r-- | moon-abe/pbc-0.5.14/arith/init_random.win32.c | 52 | ||||
| -rw-r--r-- | moon-abe/pbc-0.5.14/arith/montfp.c | 596 | ||||
| -rw-r--r-- | moon-abe/pbc-0.5.14/arith/multiz.c | 589 | ||||
| -rw-r--r-- | moon-abe/pbc-0.5.14/arith/naivefp.c | 270 | ||||
| -rw-r--r-- | moon-abe/pbc-0.5.14/arith/poly.c | 1724 | ||||
| -rw-r--r-- | moon-abe/pbc-0.5.14/arith/random.c | 87 | ||||
| -rw-r--r-- | moon-abe/pbc-0.5.14/arith/ternary_extension_field.c | 950 | ||||
| -rw-r--r-- | moon-abe/pbc-0.5.14/arith/tinyfp.c | 304 | ||||
| -rw-r--r-- | moon-abe/pbc-0.5.14/arith/z.c | 263 |
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; +} |