aboutsummaryrefslogtreecommitdiffstats
path: root/moon-abe/pbc-0.5.14/arith/poly.c
diff options
context:
space:
mode:
authorwukong <rebirthmonkey@gmail.com>2015-11-23 17:48:48 +0100
committerwukong <rebirthmonkey@gmail.com>2015-11-23 17:48:48 +0100
commitfca74d4bc3569506a6659880a89aa009dc11f552 (patch)
tree4cefd06af989608ea8ebd3bc6306889e2a1ad175 /moon-abe/pbc-0.5.14/arith/poly.c
parent840ac3ebca7af381132bf7e93c1e4c0430d6b16a (diff)
moon-abe cleanup
Change-Id: Ie1259856db03f0b9e80de3e967ec6bd1f03191b3
Diffstat (limited to 'moon-abe/pbc-0.5.14/arith/poly.c')
-rw-r--r--moon-abe/pbc-0.5.14/arith/poly.c1724
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;
-}