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