diff options
Diffstat (limited to 'moon-abe/pbc-0.5.14/arith/field.c')
| -rw-r--r-- | moon-abe/pbc-0.5.14/arith/field.c | 889 |
1 files changed, 889 insertions, 0 deletions
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); +} |