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