1 #include <ctype.h>
2 #include <stdarg.h>
3 #include <stdio.h>
4 #include <stdint.h> // for intptr_t
5 #include <stdlib.h>
6 #include <gmp.h>
7 #include "pbc_utils.h"
8 #include "pbc_field.h"
9 #include "pbc_multiz.h"
10 #include "pbc_poly.h"
11 #include "pbc_memory.h"
12 #include "misc/darray.h"
13
14 // == Polynomial rings ==
15 //
16 // Per-field data:
17 typedef struct {
18 field_ptr field; // Ring where coefficients live.
19 fieldmap mapbase; // Map element from underlying field to constant term.
20 } *pfptr;
21
22 // Per-element data:
23 //TODO: Would we ever need any field besides coeff?
24 typedef struct {
25 // The coefficients are held in a darray which is resized as needed.
26 // The last array entry represents the leading coefficient and should be
27 // nonzero. An empty darray represents 0.
28 darray_t coeff;
29 } *peptr;
30
31 // == Polynomial modulo rings ==
32 //
33 // Per-field data:
34 typedef struct {
35 field_ptr field; // Base field.
36 fieldmap mapbase; // Similar to mapbase above.
37 int n; // Degree of extension.
38 element_t poly; // Polynomial of degree n.
39 element_t *xpwr; // x^n,...,x^{2n-2} mod poly
40 } *mfptr;
41 // Per-element data: just a pointer to an array of element_t. This array always
42 // has size n.
43
44 // Add or remove coefficients until there are exactly n of them. Any new
45 // coefficients are initialized to zero, which violates the invariant that the
46 // leading coefficient must be nonzero. Thus routines calling this function
47 // must check for this and fix the polynomial if necessary, e.g. by calling
48 // poly_remove_leading_zeroes().
poly_alloc(element_ptr e,int n)49 static void poly_alloc(element_ptr e, int n) {
50 pfptr pdp = e->field->data;
51 peptr p = e->data;
52 element_ptr e0;
53 int k = p->coeff->count;
54 while (k < n) {
55 e0 = pbc_malloc(sizeof(element_t));
56 element_init(e0, pdp->field);
57 darray_append(p->coeff, e0);
58 k++;
59 }
60 while (k > n) {
61 k--;
62 e0 = darray_at(p->coeff, k);
63 element_clear(e0);
64 pbc_free(e0);
65 darray_remove_last(p->coeff);
66 }
67 }
68
poly_init(element_ptr e)69 static void poly_init(element_ptr e) {
70 peptr p = e->data = pbc_malloc(sizeof(*p));
71 darray_init(p->coeff);
72 }
73
poly_clear(element_ptr e)74 static void poly_clear(element_ptr e) {
75 peptr p = e->data;
76
77 poly_alloc(e, 0);
78 darray_clear(p->coeff);
79 pbc_free(e->data);
80 }
81
82 // Some operations may zero a leading coefficient, which will cause other
83 // routines to fail. After such an operation, this function should be called,
84 // as it strips all leading zero coefficients and frees the memory they
85 // occupied, reestablishing the guarantee that the last element of the array
86 // is nonzero.
poly_remove_leading_zeroes(element_ptr e)87 static void poly_remove_leading_zeroes(element_ptr e) {
88 peptr p = e->data;
89 int n = p->coeff->count - 1;
90 while (n >= 0) {
91 element_ptr e0 = p->coeff->item[n];
92 if (!element_is0(e0)) return;
93 element_clear(e0);
94 pbc_free(e0);
95 darray_remove_last(p->coeff);
96 n--;
97 }
98 }
99
poly_set0(element_ptr e)100 static void poly_set0(element_ptr e) {
101 poly_alloc(e, 0);
102 }
103
poly_set1(element_ptr e)104 static void poly_set1(element_ptr e) {
105 peptr p = e->data;
106 element_ptr e0;
107
108 poly_alloc(e, 1);
109 e0 = p->coeff->item[0];
110 element_set1(e0);
111 }
112
poly_is0(element_ptr e)113 static int poly_is0(element_ptr e) {
114 peptr p = e->data;
115 return !p->coeff->count;
116 }
117
poly_is1(element_ptr e)118 static int poly_is1(element_ptr e) {
119 peptr p = e->data;
120 if (p->coeff->count == 1) {
121 return element_is1(p->coeff->item[0]);
122 }
123 return 0;
124 }
125
poly_set_si(element_ptr e,signed long int op)126 static void poly_set_si(element_ptr e, signed long int op) {
127 peptr p = e->data;
128 element_ptr e0;
129
130 poly_alloc(e, 1);
131 e0 = p->coeff->item[0];
132 element_set_si(e0, op);
133 poly_remove_leading_zeroes(e);
134 }
135
poly_set_mpz(element_ptr e,mpz_ptr op)136 static void poly_set_mpz(element_ptr e, mpz_ptr op) {
137 peptr p = e->data;
138
139 poly_alloc(e, 1);
140 element_set_mpz(p->coeff->item[0], op);
141 poly_remove_leading_zeroes(e);
142 }
143
poly_set_multiz(element_ptr e,multiz op)144 static void poly_set_multiz(element_ptr e, multiz op) {
145 if (multiz_is_z(op)) {
146 // TODO: Remove unnecessary copy.
147 mpz_t z;
148 mpz_init(z);
149 multiz_to_mpz(z, op);
150 poly_set_mpz(e, z);
151 mpz_clear(z);
152 return;
153 }
154 peptr p = e->data;
155 int n = multiz_count(op);
156 poly_alloc(e, n);
157 int i;
158 for(i = 0; i < n; i++) {
159 element_set_multiz(p->coeff->item[i], multiz_at(op, i));
160 }
161 poly_remove_leading_zeroes(e);
162 }
163
poly_set(element_ptr dst,element_ptr src)164 static void poly_set(element_ptr dst, element_ptr src) {
165 peptr psrc = src->data;
166 peptr pdst = dst->data;
167 int i;
168
169 poly_alloc(dst, psrc->coeff->count);
170 for (i=0; i<psrc->coeff->count; i++) {
171 element_set(pdst->coeff->item[i], psrc->coeff->item[i]);
172 }
173 }
174
poly_coeff_count(element_ptr e)175 static int poly_coeff_count(element_ptr e) {
176 return ((peptr) e->data)->coeff->count;
177 }
178
poly_coeff(element_ptr e,int n)179 static element_ptr poly_coeff(element_ptr e, int n) {
180 peptr ep = e->data;
181 PBC_ASSERT(n < poly_coeff_count(e), "coefficient out of range");
182 return (element_ptr) ep->coeff->item[n];
183 }
184
poly_sgn(element_ptr f)185 static int poly_sgn(element_ptr f) {
186 int res = 0;
187 int i;
188 int n = poly_coeff_count(f);
189 for (i=0; i<n; i++) {
190 res = element_sgn(poly_coeff(f, i));
191 if (res) break;
192 }
193 return res;
194 }
195
poly_add(element_ptr sum,element_ptr f,element_ptr g)196 static void poly_add(element_ptr sum, element_ptr f, element_ptr g) {
197 int i, n, n1;
198 element_ptr big;
199
200 n = poly_coeff_count(f);
201 n1 = poly_coeff_count(g);
202 if (n > n1) {
203 big = f;
204 n = n1;
205 n1 = poly_coeff_count(f);
206 } else {
207 big = g;
208 }
209
210 poly_alloc(sum, n1);
211 for (i=0; i<n; i++) {
212 element_add(poly_coeff(sum, i), poly_coeff(f, i), poly_coeff(g, i));
213 }
214 for (; i<n1; i++) {
215 element_set(poly_coeff(sum, i), poly_coeff(big, i));
216 }
217 poly_remove_leading_zeroes(sum);
218 }
219
poly_sub(element_ptr diff,element_ptr f,element_ptr g)220 static void poly_sub(element_ptr diff, element_ptr f, element_ptr g) {
221 int i, n, n1;
222 element_ptr big;
223
224 n = poly_coeff_count(f);
225 n1 = poly_coeff_count(g);
226 if (n > n1) {
227 big = f;
228 n = n1;
229 n1 = poly_coeff_count(f);
230 } else {
231 big = g;
232 }
233
234 poly_alloc(diff, n1);
235 for (i=0; i<n; i++) {
236 element_sub(poly_coeff(diff, i), poly_coeff(f, i), poly_coeff(g, i));
237 }
238 for (; i<n1; i++) {
239 if (big == f) {
240 element_set(poly_coeff(diff, i), poly_coeff(big, i));
241 } else {
242 element_neg(poly_coeff(diff, i), poly_coeff(big, i));
243 }
244 }
245 poly_remove_leading_zeroes(diff);
246 }
247
poly_neg(element_ptr f,element_ptr g)248 static void poly_neg(element_ptr f, element_ptr g) {
249 peptr pf = f->data;
250 peptr pg = g->data;
251 int i, n;
252
253 n = pg->coeff->count;
254 poly_alloc(f, n);
255 for (i=0; i<n; i++) {
256 element_neg(pf->coeff->item[i], pg->coeff->item[i]);
257 }
258 }
259
poly_double(element_ptr f,element_ptr g)260 static void poly_double(element_ptr f, element_ptr g) {
261 peptr pf = f->data;
262 peptr pg = g->data;
263 int i, n;
264
265 n = pg->coeff->count;
266 poly_alloc(f, n);
267 for (i=0; i<n; i++) {
268 element_double(pf->coeff->item[i], pg->coeff->item[i]);
269 }
270 }
271
poly_mul_mpz(element_ptr f,element_ptr g,mpz_ptr z)272 static void poly_mul_mpz(element_ptr f, element_ptr g, mpz_ptr z) {
273 peptr pf = f->data;
274 peptr pg = g->data;
275 int i, n;
276
277 n = pg->coeff->count;
278 poly_alloc(f, n);
279 for (i=0; i<n; i++) {
280 element_mul_mpz(pf->coeff->item[i], pg->coeff->item[i], z);
281 }
282 }
283
poly_mul_si(element_ptr f,element_ptr g,signed long int z)284 static void poly_mul_si(element_ptr f, element_ptr g, signed long int z) {
285 peptr pf = f->data;
286 peptr pg = g->data;
287 int i, n;
288
289 n = pg->coeff->count;
290 poly_alloc(f, n);
291 for (i=0; i<n; i++) {
292 element_mul_si(pf->coeff->item[i], pg->coeff->item[i], z);
293 }
294 }
295
poly_mul(element_ptr r,element_ptr f,element_ptr g)296 static void poly_mul(element_ptr r, element_ptr f, element_ptr g) {
297 peptr pprod;
298 peptr pf = f->data;
299 peptr pg = g->data;
300 pfptr pdp = r->field->data;
301 int fcount = pf->coeff->count;
302 int gcount = pg->coeff->count;
303 int i, j, n;
304 element_t prod;
305 element_t e0;
306
307 if (!fcount || !gcount) {
308 element_set0(r);
309 return;
310 }
311 element_init(prod, r->field);
312 pprod = prod->data;
313 n = fcount + gcount - 1;
314 poly_alloc(prod, n);
315 element_init(e0, pdp->field);
316 for (i=0; i<n; i++) {
317 element_ptr x = pprod->coeff->item[i];
318 element_set0(x);
319 for (j=0; j<=i; j++) {
320 if (j < fcount && i - j < gcount) {
321 element_mul(e0, pf->coeff->item[j], pg->coeff->item[i - j]);
322 element_add(x, x, e0);
323 }
324 }
325 }
326 poly_remove_leading_zeroes(prod);
327 element_set(r, prod);
328 element_clear(e0);
329 element_clear(prod);
330 }
331
polymod_random(element_ptr e)332 static void polymod_random(element_ptr e) {
333 element_t *coeff = e->data;
334 int i, n = polymod_field_degree(e->field);
335
336 for (i=0; i<n; i++) {
337 element_random(coeff[i]);
338 }
339 }
340
polymod_from_hash(element_ptr e,void * data,int len)341 static void polymod_from_hash(element_ptr e, void *data, int len) {
342 // TODO: Improve this.
343 element_t *coeff = e->data;
344 int i, n = polymod_field_degree(e->field);
345 for (i=0; i<n; i++) {
346 element_from_hash(coeff[i], data, len);
347 }
348 }
349
poly_out_str(FILE * stream,int base,element_ptr e)350 static size_t poly_out_str(FILE *stream, int base, element_ptr e) {
351 int i;
352 int n = poly_coeff_count(e);
353 size_t result = 2, status;
354
355 /*
356 if (!n) {
357 if (EOF == fputs("[0]", stream)) return 0;
358 return 3;
359 }
360 */
361 if (EOF == fputc('[', stream)) return 0;
362 for (i=0; i<n; i++) {
363 if (i) {
364 if (EOF == fputs(", ", stream)) return 0;
365 result += 2;
366 }
367 status = element_out_str(stream, base, poly_coeff(e, i));
368 if (!status) return 0;
369 result += status;
370 }
371 if (EOF == fputc(']', stream)) return 0;
372 return result;
373 }
374
poly_snprint(char * s,size_t size,element_ptr e)375 static int poly_snprint(char *s, size_t size, element_ptr e) {
376 int i;
377 int n = poly_coeff_count(e);
378 size_t result = 0, left;
379 int status;
380
381 #define clip_sub() { \
382 result += status; \
383 left = result >= size ? 0 : size - result; \
384 }
385
386 status = snprintf(s, size, "[");
387 if (status < 0) return status;
388 clip_sub();
389
390 for (i=0; i<n; i++) {
391 if (i) {
392 status = snprintf(s + result, left, ", ");
393 if (status < 0) return status;
394 clip_sub();
395 }
396 status = element_snprint(s + result, left, poly_coeff(e, i));
397 if (status < 0) return status;
398 clip_sub();
399 }
400 status = snprintf(s + result, left, "]");
401 if (status < 0) return status;
402 return result + status;
403 #undef clip_sub
404 }
405
poly_div(element_ptr quot,element_ptr rem,element_ptr a,element_ptr b)406 static void poly_div(element_ptr quot, element_ptr rem,
407 element_ptr a, element_ptr b) {
408 peptr pq, pr;
409 pfptr pdp = a->field->data;
410 element_t q, r;
411 element_t binv, e0;
412 element_ptr qe;
413 int m, n;
414 int i, k;
415
416 if (element_is0(b)) pbc_die("division by zero");
417 n = poly_degree(b);
418 m = poly_degree(a);
419 if (n > m) {
420 element_set(rem, a);
421 element_set0(quot);
422 return;
423 }
424 element_init(r, a->field);
425 element_init(q, a->field);
426 element_init(binv, pdp->field);
427 element_init(e0, pdp->field);
428 pq = q->data;
429 pr = r->data;
430 element_set(r, a);
431 k = m - n;
432 poly_alloc(q, k + 1);
433 element_invert(binv, poly_coeff(b, n));
434 while (k >= 0) {
435 qe = pq->coeff->item[k];
436 element_mul(qe, binv, pr->coeff->item[m]);
437 for (i=0; i<=n; i++) {
438 element_mul(e0, qe, poly_coeff(b, i));
439 element_sub(pr->coeff->item[i + k], pr->coeff->item[i + k], e0);
440 }
441 k--;
442 m--;
443 }
444 poly_remove_leading_zeroes(r);
445 element_set(quot, q);
446 element_set(rem, r);
447
448 element_clear(q);
449 element_clear(r);
450 element_clear(e0);
451 element_clear(binv);
452 }
453
poly_invert(element_ptr res,element_ptr f,element_ptr m)454 static void poly_invert(element_ptr res, element_ptr f, element_ptr m) {
455 element_t q, r0, r1, r2;
456 element_t b0, b1, b2;
457 element_t inv;
458
459 element_init(b0, res->field);
460 element_init(b1, res->field);
461 element_init(b2, res->field);
462 element_init(q, res->field);
463 element_init(r0, res->field);
464 element_init(r1, res->field);
465 element_init(r2, res->field);
466 element_init(inv, poly_base_field(res));
467 element_set0(b0);
468 element_set1(b1);
469 element_set(r0, m);
470 element_set(r1, f);
471
472 for (;;) {
473 poly_div(q, r2, r0, r1);
474 if (element_is0(r2)) break;
475 element_mul(b2, b1, q);
476 element_sub(b2, b0, b2);
477 element_set(b0, b1);
478 element_set(b1, b2);
479 element_set(r0, r1);
480 element_set(r1, r2);
481 }
482 element_invert(inv, poly_coeff(r1, 0));
483 poly_const_mul(res, inv, b1);
484 element_clear(inv);
485 element_clear(q);
486 element_clear(r0);
487 element_clear(r1);
488 element_clear(r2);
489 element_clear(b0);
490 element_clear(b1);
491 element_clear(b2);
492 }
493
poly_to_polymod_truncate(element_ptr e,element_ptr f)494 static void poly_to_polymod_truncate(element_ptr e, element_ptr f) {
495 mfptr p = e->field->data;
496 element_t *coeff = e->data;
497 int i;
498 int n;
499 n = poly_coeff_count(f);
500 if (n > p->n) n = p->n;
501
502 for (i=0; i<n; i++) {
503 element_set(coeff[i], poly_coeff(f, i));
504 }
505 for (; i<p->n; i++) {
506 element_set0(coeff[i]);
507 }
508 }
509
polymod_to_poly(element_ptr f,element_ptr e)510 static void polymod_to_poly(element_ptr f, element_ptr e) {
511 mfptr p = e->field->data;
512 element_t *coeff = e->data;
513 int i, n = p->n;
514 poly_alloc(f, n);
515 for (i=0; i<n; i++) {
516 element_set(poly_coeff(f, i), coeff[i]);
517 }
518 poly_remove_leading_zeroes(f);
519 }
520
polymod_invert(element_ptr r,element_ptr e)521 static void polymod_invert(element_ptr r, element_ptr e) {
522 mfptr p = r->field->data;
523 element_ptr minpoly = p->poly;
524 element_t f, r1;
525
526 element_init(f, minpoly->field);
527 element_init(r1, minpoly->field);
528 polymod_to_poly(f, e);
529
530 poly_invert(r1, f, p->poly);
531
532 poly_to_polymod_truncate(r, r1);
533
534 element_clear(f);
535 element_clear(r1);
536 }
537
poly_cmp(element_ptr f,element_ptr g)538 static int poly_cmp(element_ptr f, element_ptr g) {
539 int i;
540 int n = poly_coeff_count(f);
541 int n1 = poly_coeff_count(g);
542 if (n != n1) return 1;
543 for (i=0; i<n; i++) {
544 if (element_cmp(poly_coeff(f, i), poly_coeff(g, i))) return 1;
545 }
546 return 0;
547 }
548
field_clear_poly(field_ptr f)549 static void field_clear_poly(field_ptr f) {
550 pfptr p = f->data;
551 pbc_free(p);
552 }
553
554 // 2 bytes hold the number of terms, then the terms follow.
555 // Bad for sparse polynomials.
poly_length_in_bytes(element_t p)556 static int poly_length_in_bytes(element_t p) {
557 int count = poly_coeff_count(p);
558 int result = 2;
559 int i;
560 for (i=0; i<count; i++) {
561 result += element_length_in_bytes(poly_coeff(p, i));
562 }
563 return result;
564 }
565
poly_to_bytes(unsigned char * buf,element_t p)566 static int poly_to_bytes(unsigned char *buf, element_t p) {
567 int count = poly_coeff_count(p);
568 int result = 2;
569 int i;
570 buf[0] = (unsigned char) count;
571 buf[1] = (unsigned char) (count >> 8);
572 for (i=0; i<count; i++) {
573 result += element_to_bytes(&buf[result], poly_coeff(p, i));
574 }
575 return result;
576 }
577
poly_from_bytes(element_t p,unsigned char * buf)578 static int poly_from_bytes(element_t p, unsigned char *buf) {
579 int result = 2;
580 int count = buf[0] + buf[1] * 256;
581 int i;
582 poly_alloc(p, count);
583 for (i=0; i<count; i++) {
584 result += element_from_bytes(poly_coeff(p, i), &buf[result]);
585 }
586 return result;
587 }
588
589 // Is this useful? This returns to_mpz(constant term).
poly_to_mpz(mpz_t z,element_ptr e)590 static void poly_to_mpz(mpz_t z, element_ptr e) {
591 if (!poly_coeff_count(e)) {
592 mpz_set_ui(z, 0);
593 } else {
594 element_to_mpz(z, poly_coeff(e, 0));
595 }
596 }
597
poly_out_info(FILE * str,field_ptr f)598 static void poly_out_info(FILE *str, field_ptr f) {
599 pfptr p = f->data;
600 fprintf(str, "Polynomial ring over ");
601 field_out_info(str, p->field);
602 }
603
field_clear_polymod(field_ptr f)604 static void field_clear_polymod(field_ptr f) {
605 mfptr p = f->data;
606 int i, n = p->n;
607
608 for (i=0; i<n; i++) {
609 element_clear(p->xpwr[i]);
610 }
611 pbc_free(p->xpwr);
612
613 element_clear(p->poly);
614 pbc_free(f->data);
615 }
616
polymod_is_sqr(element_ptr e)617 static int polymod_is_sqr(element_ptr e) {
618 int res;
619 mpz_t z;
620 element_t e0;
621
622 element_init(e0, e->field);
623 mpz_init(z);
624 mpz_sub_ui(z, e->field->order, 1);
625 mpz_divexact_ui(z, z, 2);
626
627 element_pow_mpz(e0, e, z);
628 res = element_is1(e0);
629 element_clear(e0);
630 mpz_clear(z);
631 return res;
632 }
633
634 // Find a square root in a polynomial modulo ring using Cantor-Zassenhaus aka
635 // Legendre's method.
polymod_sqrt(element_ptr res,element_ptr a)636 static void polymod_sqrt(element_ptr res, element_ptr a) {
637 // TODO: Use a faster method? See Bernstein.
638 field_t kx;
639 element_t f;
640 element_t r, s;
641 element_t e0;
642 mpz_t z;
643
644 field_init_poly(kx, a->field);
645 mpz_init(z);
646 element_init(f, kx);
647 element_init(r, kx);
648 element_init(s, kx);
649 element_init(e0, a->field);
650
651 poly_alloc(f, 3);
652 element_set1(poly_coeff(f, 2));
653 element_neg(poly_coeff(f, 0), a);
654
655 mpz_sub_ui(z, a->field->order, 1);
656 mpz_divexact_ui(z, z, 2);
657 for (;;) {
658 int i;
659 element_ptr x;
660 element_ptr e1, e2;
661
662 poly_alloc(r, 2);
663 element_set1(poly_coeff(r, 1));
664 x = poly_coeff(r, 0);
665 element_random(x);
666 element_mul(e0, x, x);
667 if (!element_cmp(e0, a)) {
668 element_set(res, x);
669 break;
670 }
671 element_set1(s);
672 //TODO: this can be optimized greatly
673 //since we know r has the form ax + b
674 for (i = mpz_sizeinbase(z, 2) - 1; i >=0; i--) {
675 element_mul(s, s, s);
676 if (poly_degree(s) == 2) {
677 e1 = poly_coeff(s, 0);
678 e2 = poly_coeff(s, 2);
679 element_mul(e0, e2, a);
680 element_add(e1, e1, e0);
681 poly_alloc(s, 2);
682 poly_remove_leading_zeroes(s);
683 }
684 if (mpz_tstbit(z, i)) {
685 element_mul(s, s, r);
686 if (poly_degree(s) == 2) {
687 e1 = poly_coeff(s, 0);
688 e2 = poly_coeff(s, 2);
689 element_mul(e0, e2, a);
690 element_add(e1, e1, e0);
691 poly_alloc(s, 2);
692 poly_remove_leading_zeroes(s);
693 }
694 }
695 }
696 if (poly_degree(s) < 1) continue;
697 element_set1(e0);
698 e1 = poly_coeff(s, 0);
699 e2 = poly_coeff(s, 1);
700 element_add(e1, e1, e0);
701 element_invert(e0, e2);
702 element_mul(e0, e0, e1);
703 element_mul(e2, e0, e0);
704 if (!element_cmp(e2, a)) {
705 element_set(res, e0);
706 break;
707 }
708 }
709
710 mpz_clear(z);
711 element_clear(f);
712 element_clear(r);
713 element_clear(s);
714 element_clear(e0);
715 field_clear(kx);
716 }
717
polymod_to_bytes(unsigned char * data,element_t f)718 static int polymod_to_bytes(unsigned char *data, element_t f) {
719 mfptr p = f->field->data;
720 element_t *coeff = f->data;
721 int i, n = p->n;
722 int len = 0;
723 for (i=0; i<n; i++) {
724 len += element_to_bytes(data + len, coeff[i]);
725 }
726 return len;
727 }
728
polymod_length_in_bytes(element_t f)729 static int polymod_length_in_bytes(element_t f) {
730 mfptr p = f->field->data;
731 element_t *coeff = f->data;
732 int i, n = p->n;
733 int res = 0;
734
735 for (i=0; i<n; i++) {
736 res += element_length_in_bytes(coeff[i]);
737 }
738
739 return res;
740 }
741
polymod_from_bytes(element_t f,unsigned char * data)742 static int polymod_from_bytes(element_t f, unsigned char *data) {
743 mfptr p = f->field->data;
744 element_t *coeff = f->data;
745 int i, n = p->n;
746 int len = 0;
747
748 for (i=0; i<n; i++) {
749 len += element_from_bytes(coeff[i], data + len);
750 }
751 return len;
752 }
753
polymod_init(element_t e)754 static void polymod_init(element_t e) {
755 int i;
756 mfptr p = e->field->data;
757 int n = p->n;
758 element_t *coeff;
759 coeff = e->data = pbc_malloc(sizeof(element_t) * n);
760
761 for (i=0; i<n; i++) {
762 element_init(coeff[i], p->field);
763 }
764 }
765
polymod_clear(element_t e)766 static void polymod_clear(element_t e) {
767 mfptr p = e->field->data;
768 element_t *coeff = e->data;
769 int i, n = p->n;
770 for (i=0; i<n; i++) {
771 element_clear(coeff[i]);
772 }
773 pbc_free(e->data);
774 }
775
polymod_set_si(element_t e,signed long int x)776 static void polymod_set_si(element_t e, signed long int x) {
777 mfptr p = e->field->data;
778 element_t *coeff = e->data;
779 int i, n = p->n;
780 element_set_si(coeff[0], x);
781 for (i=1; i<n; i++) {
782 element_set0(coeff[i]);
783 }
784 }
785
polymod_set_mpz(element_t e,mpz_t z)786 static void polymod_set_mpz(element_t e, mpz_t z) {
787 mfptr p = e->field->data;
788 element_t *coeff = e->data;
789 int i, n = p->n;
790 element_set_mpz(coeff[0], z);
791 for (i=1; i<n; i++) {
792 element_set0(coeff[i]);
793 }
794 }
795
polymod_set(element_t e,element_t f)796 static void polymod_set(element_t e, element_t f) {
797 mfptr p = e->field->data;
798 element_t *dst = e->data, *src = f->data;
799 int i, n = p->n;
800 for (i=0; i<n; i++) {
801 element_set(dst[i], src[i]);
802 }
803 }
804
polymod_neg(element_t e,element_t f)805 static void polymod_neg(element_t e, element_t f) {
806 mfptr p = e->field->data;
807 element_t *dst = e->data, *src = f->data;
808 int i, n = p->n;
809 for (i=0; i<n; i++) {
810 element_neg(dst[i], src[i]);
811 }
812 }
813
polymod_cmp(element_ptr f,element_ptr g)814 static int polymod_cmp(element_ptr f, element_ptr g) {
815 mfptr p = f->field->data;
816 element_t *c1 = f->data, *c2 = g->data;
817 int i, n = p->n;
818 for (i=0; i<n; i++) {
819 if (element_cmp(c1[i], c2[i])) return 1;
820 }
821 return 0;
822 }
823
polymod_add(element_t r,element_t e,element_t f)824 static void polymod_add(element_t r, element_t e, element_t f) {
825 mfptr p = r->field->data;
826 element_t *dst = r->data, *s1 = e->data, *s2 = f->data;
827 int i, n = p->n;
828 for (i=0; i<n; i++) {
829 element_add(dst[i], s1[i], s2[i]);
830 }
831 }
832
polymod_double(element_t r,element_t f)833 static void polymod_double(element_t r, element_t f) {
834 mfptr p = r->field->data;
835 element_t *dst = r->data, *s1 = f->data;
836 int i, n = p->n;
837 for (i=0; i<n; i++) {
838 element_double(dst[i], s1[i]);
839 }
840 }
841
polymod_sub(element_t r,element_t e,element_t f)842 static void polymod_sub(element_t r, element_t e, element_t f) {
843 mfptr p = r->field->data;
844 element_t *dst = r->data, *s1 = e->data, *s2 = f->data;
845 int i, n = p->n;
846 for (i=0; i<n; i++) {
847 element_sub(dst[i], s1[i], s2[i]);
848 }
849 }
850
polymod_mul_mpz(element_t e,element_t f,mpz_ptr z)851 static void polymod_mul_mpz(element_t e, element_t f, mpz_ptr z) {
852 mfptr p = e->field->data;
853 element_t *dst = e->data, *src = f->data;
854 int i, n = p->n;
855 for (i=0; i<n; i++) {
856 element_mul_mpz(dst[i], src[i], z);
857 }
858 }
859
polymod_mul_si(element_t e,element_t f,signed long int z)860 static void polymod_mul_si(element_t e, element_t f, signed long int z) {
861 mfptr p = e->field->data;
862 element_t *dst = e->data, *src = f->data;
863 int i, n = p->n;
864 for (i=0; i<n; i++) {
865 element_mul_si(dst[i], src[i], z);
866 }
867 }
868
869 // Karatsuba multiplication for degree 2 polynomials.
kar_poly_2(element_t * dst,element_t c3,element_t c4,element_t * s1,element_t * s2,element_t * scratch)870 static void kar_poly_2(element_t *dst, element_t c3, element_t c4, element_t *s1, element_t *s2, element_t *scratch) {
871 element_ptr c01, c02, c12;
872
873 c12 = scratch[0];
874 c02 = scratch[1];
875 c01 = scratch[2];
876
877 element_add(c3, s1[0], s1[1]);
878 element_add(c4, s2[0], s2[1]);
879 element_mul(c01, c3, c4);
880 element_add(c3, s1[0], s1[2]);
881 element_add(c4, s2[0], s2[2]);
882 element_mul(c02, c3, c4);
883 element_add(c3, s1[1], s1[2]);
884 element_add(c4, s2[1], s2[2]);
885 element_mul(c12, c3, c4);
886
887 element_mul(dst[1], s1[1], s2[1]);
888
889 // Constant term.
890 element_mul(dst[0], s1[0], s2[0]);
891
892 // Coefficient of x^4.
893 element_mul(c4, s1[2], s2[2]);
894
895 // Coefficient of x^3.
896 element_add(c3, dst[1], c4);
897 element_sub(c3, c12, c3);
898
899 // Coefficient of x^2.
900 element_add(dst[2], c4, dst[0]);
901 element_sub(c02, c02, dst[2]);
902 element_add(dst[2], dst[1], c02);
903
904 // Coefficient of x.
905 element_sub(c01, c01, dst[0]);
906 element_sub(dst[1], c01, dst[1]);
907 }
908
909 // Degree 3, 6 polynomial moduli have dedicated routines for multiplication.
polymod_mul_degree3(element_ptr res,element_ptr e,element_ptr f)910 static void polymod_mul_degree3(element_ptr res, element_ptr e, element_ptr f) {
911 mfptr p = res->field->data;
912 element_t *dst = res->data, *s1 = e->data, *s2 = f->data;
913 element_t c3, c4;
914 element_t p0;
915
916 element_init(p0, res->field);
917 element_init(c3, p->field);
918 element_init(c4, p->field);
919
920 kar_poly_2(dst, c3, c4, s1, s2, p0->data);
921
922 polymod_const_mul(p0, c3, p->xpwr[0]);
923 element_add(res, res, p0);
924 polymod_const_mul(p0, c4, p->xpwr[1]);
925 element_add(res, res, p0);
926
927 element_clear(p0);
928 element_clear(c3);
929 element_clear(c4);
930 }
931
polymod_mul_degree6(element_ptr res,element_ptr e,element_ptr f)932 static void polymod_mul_degree6(element_ptr res, element_ptr e, element_ptr f) {
933 mfptr p = res->field->data;
934 element_t *dst = res->data, *s0, *s1 = e->data, *s2 = f->data;
935 element_t *a0, *a1, *b0, *b1;
936 element_t p0, p1, p2, p3;
937
938 a0 = s1;
939 a1 = &s1[3];
940 b0 = s2;
941 b1 = &s2[3];
942
943 element_init(p0, res->field);
944 element_init(p1, res->field);
945 element_init(p2, res->field);
946 element_init(p3, res->field);
947
948 s0 = p0->data;
949 s1 = p1->data;
950 s2 = p2->data;
951 element_add(s0[0], a0[0], a1[0]);
952 element_add(s0[1], a0[1], a1[1]);
953 element_add(s0[2], a0[2], a1[2]);
954
955 element_add(s1[0], b0[0], b1[0]);
956 element_add(s1[1], b0[1], b1[1]);
957 element_add(s1[2], b0[2], b1[2]);
958
959 kar_poly_2(s2, s2[3], s2[4], s0, s1, p3->data);
960 kar_poly_2(s0, s0[3], s0[4], a0, b0, p3->data);
961 kar_poly_2(s1, s1[3], s1[4], a1, b1, p3->data);
962
963 element_set(dst[0], s0[0]);
964 element_set(dst[1], s0[1]);
965 element_set(dst[2], s0[2]);
966
967 element_sub(dst[3], s0[3], s0[0]);
968 element_sub(dst[3], dst[3], s1[0]);
969 element_add(dst[3], dst[3], s2[0]);
970
971 element_sub(dst[4], s0[4], s0[1]);
972 element_sub(dst[4], dst[4], s1[1]);
973 element_add(dst[4], dst[4], s2[1]);
974
975 element_sub(dst[5], s2[2], s0[2]);
976 element_sub(dst[5], dst[5], s1[2]);
977
978 // Start reusing part of s0 as scratch space(!)
979 element_sub(s0[0], s2[3], s0[3]);
980 element_sub(s0[0], s0[0], s1[3]);
981 element_add(s0[0], s0[0], s1[0]);
982
983 element_sub(s0[1], s2[4], s0[4]);
984 element_sub(s0[1], s0[1], s1[4]);
985 element_add(s0[1], s0[1], s1[1]);
986
987 polymod_const_mul(p3, s0[0], p->xpwr[0]);
988 element_add(res, res, p3);
989 polymod_const_mul(p3, s0[1], p->xpwr[1]);
990 element_add(res, res, p3);
991 polymod_const_mul(p3, s1[2], p->xpwr[2]);
992 element_add(res, res, p3);
993 polymod_const_mul(p3, s1[3], p->xpwr[3]);
994 element_add(res, res, p3);
995 polymod_const_mul(p3, s1[4], p->xpwr[4]);
996 element_add(res, res, p3);
997
998 element_clear(p0);
999 element_clear(p1);
1000 element_clear(p2);
1001 element_clear(p3);
1002 }
1003
1004 // General polynomial modulo ring multiplication.
polymod_mul(element_ptr res,element_ptr e,element_ptr f)1005 static void polymod_mul(element_ptr res, element_ptr e, element_ptr f) {
1006 mfptr p = res->field->data;
1007 int n = p->n;
1008 element_t *dst;
1009 element_t *s1 = e->data, *s2 = f->data;
1010 element_t prod, p0, c0;
1011 int i, j;
1012 element_t *high; // Coefficients of x^n, ..., x^{2n-2}.
1013
1014 high = pbc_malloc(sizeof(element_t) * (n - 1));
1015 for (i=0; i<n-1; i++) {
1016 element_init(high[i], p->field);
1017 element_set0(high[i]);
1018 }
1019 element_init(prod, res->field);
1020 dst = prod->data;
1021 element_init(p0, res->field);
1022 element_init(c0, p->field);
1023
1024 for (i=0; i<n; i++) {
1025 int ni = n - i;
1026 for (j=0; j<ni; j++) {
1027 element_mul(c0, s1[i], s2[j]);
1028 element_add(dst[i + j], dst[i + j], c0);
1029 }
1030 for (;j<n; j++) {
1031 element_mul(c0, s1[i], s2[j]);
1032 element_add(high[j - ni], high[j - ni], c0);
1033 }
1034 }
1035
1036 for (i=0; i<n-1; i++) {
1037 polymod_const_mul(p0, high[i], p->xpwr[i]);
1038 element_add(prod, prod, p0);
1039 element_clear(high[i]);
1040 }
1041 pbc_free(high);
1042
1043 element_set(res, prod);
1044 element_clear(prod);
1045 element_clear(p0);
1046 element_clear(c0);
1047 }
1048
polymod_square_degree3(element_ptr res,element_ptr e)1049 static void polymod_square_degree3(element_ptr res, element_ptr e) {
1050 // TODO: Investigate if squaring is significantly cheaper than
1051 // multiplication. If so convert to Karatsuba.
1052 element_t *dst = res->data;
1053 element_t *src = e->data;
1054 mfptr p = res->field->data;
1055 element_t p0;
1056 element_t c0, c2;
1057 element_ptr c1, c3;
1058
1059 element_init(p0, res->field);
1060 element_init(c0, p->field);
1061 element_init(c2, p->field);
1062
1063 c3 = p0->data;
1064 c1 = c3 + 1;
1065
1066 element_mul(c3, src[0], src[1]);
1067 element_mul(c1, src[0], src[2]);
1068 element_square(dst[0], src[0]);
1069
1070 element_mul(c2, src[1], src[2]);
1071 element_square(c0, src[2]);
1072 element_square(dst[2], src[1]);
1073
1074 element_add(dst[1], c3, c3);
1075
1076 element_add(c1, c1, c1);
1077 element_add(dst[2], dst[2], c1);
1078
1079 polymod_const_mul(p0, c0, p->xpwr[1]);
1080 element_add(res, res, p0);
1081
1082 element_add(c2, c2, c2);
1083 polymod_const_mul(p0, c2, p->xpwr[0]);
1084 element_add(res, res, p0);
1085
1086 element_clear(p0);
1087 element_clear(c0);
1088 element_clear(c2);
1089 }
1090
polymod_square(element_ptr res,element_ptr e)1091 static void polymod_square(element_ptr res, element_ptr e) {
1092 element_t *dst;
1093 element_t *src = e->data;
1094 mfptr p = res->field->data;
1095 int n = p->n;
1096 element_t prod, p0, c0;
1097 int i, j;
1098 element_t *high; // Coefficients of x^n,...,x^{2n-2}.
1099
1100 high = pbc_malloc(sizeof(element_t) * (n - 1));
1101 for (i=0; i<n-1; i++) {
1102 element_init(high[i], p->field);
1103 element_set0(high[i]);
1104 }
1105
1106 element_init(prod, res->field);
1107 dst = prod->data;
1108 element_init(p0, res->field);
1109 element_init(c0, p->field);
1110
1111 for (i=0; i<n; i++) {
1112 int twicei = 2 * i;
1113 element_square(c0, src[i]);
1114 if (twicei < n) {
1115 element_add(dst[twicei], dst[twicei], c0);
1116 } else {
1117 element_add(high[twicei - n], high[twicei - n], c0);
1118 }
1119
1120 for (j=i+1; j<n-i; j++) {
1121 element_mul(c0, src[i], src[j]);
1122 element_add(c0, c0, c0);
1123 element_add(dst[i + j], dst[i + j], c0);
1124 }
1125 for (;j<n; j++) {
1126 element_mul(c0, src[i], src[j]);
1127 element_add(c0, c0, c0);
1128 element_add(high[i + j - n], high[i + j - n], c0);
1129 }
1130 }
1131
1132 for (i=0; i<n-1; i++) {
1133 polymod_const_mul(p0, high[i], p->xpwr[i]);
1134 element_add(prod, prod, p0);
1135 element_clear(high[i]);
1136 }
1137 pbc_free(high);
1138
1139 element_set(res, prod);
1140 element_clear(prod);
1141 element_clear(p0);
1142 element_clear(c0);
1143 }
1144
polymod_is0(element_ptr e)1145 static int polymod_is0(element_ptr e) {
1146 mfptr p = e->field->data;
1147 element_t *coeff = e->data;
1148 int i, n = p->n;
1149
1150 for (i=0; i<n; i++) {
1151 if (!element_is0(coeff[i])) return 0;
1152 }
1153 return 1;
1154 }
1155
polymod_is1(element_ptr e)1156 static int polymod_is1(element_ptr e) {
1157 mfptr p = e->field->data;
1158 element_t *coeff = e->data;
1159 int i, n = p->n;
1160
1161 if (!element_is1(coeff[0])) return 0;
1162 for (i=1; i<n; i++) {
1163 if (!element_is0(coeff[i])) return 0;
1164 }
1165 return 1;
1166 }
1167
polymod_set0(element_ptr e)1168 static void polymod_set0(element_ptr e) {
1169 mfptr p = e->field->data;
1170 element_t *coeff = e->data;
1171 int i, n = p->n;
1172
1173 for (i=0; i<n; i++) {
1174 element_set0(coeff[i]);
1175 }
1176 }
1177
polymod_set1(element_ptr e)1178 static void polymod_set1(element_ptr e) {
1179 mfptr p = e->field->data;
1180 element_t *coeff = e->data;
1181 int i, n = p->n;
1182
1183 element_set1(coeff[0]);
1184 for (i=1; i<n; i++) {
1185 element_set0(coeff[i]);
1186 }
1187 }
1188
polymod_sgn(element_ptr e)1189 static int polymod_sgn(element_ptr e) {
1190 mfptr p = e->field->data;
1191 element_t *coeff = e->data;
1192 int res = 0;
1193 int i, n = p->n;
1194 for (i=0; i<n; i++) {
1195 res = element_sgn(coeff[i]);
1196 if (res) break;
1197 }
1198 return res;
1199 }
1200
polymod_out_str(FILE * stream,int base,element_ptr e)1201 static size_t polymod_out_str(FILE *stream, int base, element_ptr e) {
1202 size_t result = 2, status;
1203 mfptr p = e->field->data;
1204 element_t *coeff = e->data;
1205 int i, n = p->n;
1206
1207 if (EOF == fputc('[', stream)) return 0;
1208 for (i=0; i<n; i++) {
1209 if (i) {
1210 if (EOF == fputs(", ", stream)) return 0;
1211 result += 2;
1212 }
1213 status = element_out_str(stream, base, coeff[i]);
1214 if (!status) return 0;
1215 result += status;
1216 }
1217 if (EOF == fputc(']', stream)) return 0;
1218 return result;
1219 }
1220
polymod_snprint(char * s,size_t size,element_ptr e)1221 static int polymod_snprint(char *s, size_t size, element_ptr e) {
1222 mfptr p = e->field->data;
1223 element_t *coeff = e->data;
1224 int i, n = p->n;
1225 size_t result = 0, left;
1226 int status;
1227
1228 #define clip_sub(void) { \
1229 result += status; \
1230 left = result >= size ? 0 : size - result; \
1231 }
1232
1233 status = snprintf(s, size, "[");
1234 if (status < 0) return status;
1235 clip_sub();
1236
1237 for (i=0; i<n; i++) {
1238 if (i) {
1239 status = snprintf(s + result, left, ", ");
1240 if (status < 0) return status;
1241 clip_sub();
1242 }
1243 status = element_snprint(s + result, left, coeff[i]);
1244 if (status < 0) return status;
1245 clip_sub();
1246 }
1247 status = snprintf(s + result, left, "]");
1248 if (status < 0) return status;
1249 return result + status;
1250 #undef clip_sub
1251 }
1252
polymod_set_multiz(element_ptr e,multiz m)1253 static void polymod_set_multiz(element_ptr e, multiz m) {
1254 mfptr p = e->field->data;
1255 element_t *coeff = e->data;
1256 int i, n = p->n;
1257 if (multiz_is_z(m)) {
1258 element_set_multiz(coeff[0], m);
1259 for (i = 1; i < n; i++) element_set0(coeff[i]);
1260 return;
1261 }
1262 int max = multiz_count(m);
1263 for (i = 0; i < n; i++) {
1264 if (i >= max) element_set0(coeff[i]);
1265 else element_set_multiz(coeff[i], multiz_at(m, i));
1266 }
1267 }
1268
polymod_set_str(element_ptr e,const char * s,int base)1269 static int polymod_set_str(element_ptr e, const char *s, int base) {
1270 mfptr p = e->field->data;
1271 element_t *coeff = e->data;
1272 int i, n = p->n;
1273 const char *cp = s;
1274 element_set0(e);
1275 while (*cp && isspace(*cp)) cp++;
1276 if (*cp++ != '[') return 0;
1277 for (i=0; i<n; i++) {
1278 cp += element_set_str(coeff[i], cp, base);
1279 while (*cp && isspace(*cp)) cp++;
1280 if (i<n-1 && *cp++ != ',') return 0;
1281 }
1282 if (*cp++ != ']') return 0;
1283 return cp - s;
1284 }
1285
polymod_coeff_count(element_ptr e)1286 static int polymod_coeff_count(element_ptr e) {
1287 UNUSED_VAR(e);
1288 mfptr p = e->field->data;
1289 return p->n;
1290 }
1291
polymod_coeff(element_ptr e,int i)1292 static element_ptr polymod_coeff(element_ptr e, int i) {
1293 element_t *coeff = e->data;
1294 return coeff[i];
1295 }
1296
polymod_to_mpz(mpz_t z,element_ptr e)1297 static void polymod_to_mpz(mpz_t z, element_ptr e) {
1298 element_to_mpz(z, polymod_coeff(e, 0));
1299 }
1300
1301 // Compute x^n,...,x^{2n-2} mod poly.
compute_x_powers(field_ptr field,element_ptr poly)1302 static void compute_x_powers(field_ptr field, element_ptr poly) {
1303 mfptr p = field->data;
1304 element_t p0;
1305 element_ptr pwrn;
1306 element_t *coeff, *coeff1;
1307 int i, j;
1308 int n = p->n;
1309 element_t *xpwr;
1310
1311 xpwr = p->xpwr;
1312
1313 element_init(p0, field);
1314 for (i=0; i<n; i++) {
1315 element_init(xpwr[i], field);
1316 }
1317 pwrn = xpwr[0];
1318 poly_to_polymod_truncate(pwrn, poly);
1319 element_neg(pwrn, pwrn);
1320
1321 for (i=1; i<n; i++) {
1322 coeff = xpwr[i-1]->data;
1323 coeff1 = xpwr[i]->data;
1324
1325 element_set0(coeff1[0]);
1326 for (j=1; j<n; j++) {
1327 element_set(coeff1[j], coeff[j - 1]);
1328 }
1329 polymod_const_mul(p0, coeff[n - 1], pwrn);
1330 element_add(xpwr[i], xpwr[i], p0);
1331 }
1332 element_clear(p0);
1333 }
1334
polymod_out_info(FILE * str,field_ptr f)1335 static void polymod_out_info(FILE *str, field_ptr f) {
1336 mfptr p = f->data;
1337 element_fprintf(str, "Extension, poly = %B, base field = ", p->poly);
1338 field_out_info(str, p->field);
1339 }
1340
1341 // Sets d = gcd(f, g).
poly_gcd(element_ptr d,element_ptr f,element_ptr g)1342 static void poly_gcd(element_ptr d, element_ptr f, element_ptr g) {
1343 element_t a, b, q, r;
1344 element_init(a, d->field);
1345 element_init(b, d->field);
1346 element_init(q, d->field);
1347 element_init(r, d->field);
1348
1349 element_set(a, f);
1350 element_set(b, g);
1351 for(;;) {
1352 //TODO: don't care about q
1353 poly_div(q, r, a, b);
1354 if (element_is0(r)) break;
1355 element_set(a, b);
1356 element_set(b, r);
1357 }
1358 element_set(d, b);
1359 element_clear(a);
1360 element_clear(b);
1361 element_clear(q);
1362 element_clear(r);
1363 }
1364
1365 // Sets f = c g where c is the inverse of the leading coefficient of g.
poly_make_monic(element_t f,element_t g)1366 static void poly_make_monic(element_t f, element_t g) {
1367 int n = poly_coeff_count(g);
1368 int i;
1369 element_ptr e0;
1370 poly_alloc(f, n);
1371 if (!n) return;
1372
1373 e0 = poly_coeff(f, n - 1);
1374 element_invert(e0, poly_coeff(g, n - 1));
1375 for (i=0; i<n-1; i++) {
1376 element_mul(poly_coeff(f, i), poly_coeff(g, i), e0);
1377 }
1378 element_set1(e0);
1379 }
1380
1381 // The above should be static.
1382
field_init_poly(field_ptr f,field_ptr base_field)1383 void field_init_poly(field_ptr f, field_ptr base_field) {
1384 field_init(f);
1385 pfptr p = f->data = pbc_malloc(sizeof(*p));
1386 p->field = base_field;
1387 p->mapbase = element_field_to_poly;
1388 f->field_clear = field_clear_poly;
1389 f->init = poly_init;
1390 f->clear = poly_clear;
1391 f->set_si = poly_set_si;
1392 f->set_multiz = poly_set_multiz;
1393 f->set_mpz = poly_set_mpz;
1394 f->to_mpz = poly_to_mpz;
1395 f->out_str = poly_out_str;
1396 f->snprint = poly_snprint;
1397 f->set = poly_set;
1398 f->sign = poly_sgn;
1399 f->add = poly_add;
1400 f->doub = poly_double;
1401 f->is0 = poly_is0;
1402 f->is1 = poly_is1;
1403 f->set0 = poly_set0;
1404 f->set1 = poly_set1;
1405 f->sub = poly_sub;
1406 f->neg = poly_neg;
1407 f->mul = poly_mul;
1408 f->mul_mpz = poly_mul_mpz;
1409 f->mul_si = poly_mul_si;
1410 f->cmp = poly_cmp;
1411 f->out_info = poly_out_info;
1412 f->item_count = poly_coeff_count;
1413 f->item = poly_coeff;
1414
1415 f->to_bytes = poly_to_bytes;
1416 f->from_bytes = poly_from_bytes;
1417 f->fixed_length_in_bytes = -1;
1418 f->length_in_bytes = poly_length_in_bytes;
1419 }
1420
poly_set_coeff(element_ptr e,element_ptr a,int n)1421 void poly_set_coeff(element_ptr e, element_ptr a, int n) {
1422 peptr p = e->data;
1423 if (p->coeff->count < n + 1) {
1424 poly_alloc(e, n + 1);
1425 }
1426 element_ptr e0 = p->coeff->item[n];
1427 element_set(e0, a);
1428 if (p->coeff->count == n + 1 && element_is0(a)) poly_remove_leading_zeroes(e);
1429 }
1430
poly_set_coeff0(element_ptr e,int n)1431 void poly_set_coeff0(element_ptr e, int n) {
1432 peptr p = e->data;
1433 if (n < p->coeff->count) {
1434 element_set0(p->coeff->item[n]);
1435 if (n == p->coeff->count - 1) poly_remove_leading_zeroes(e);
1436 }
1437 }
1438
poly_set_coeff1(element_ptr e,int n)1439 void poly_set_coeff1(element_ptr e, int n) {
1440 peptr p = e->data;
1441 if (p->coeff->count < n + 1) {
1442 poly_alloc(e, n + 1);
1443 }
1444 element_set1(p->coeff->item[n]);
1445 }
1446
poly_setx(element_ptr f)1447 void poly_setx(element_ptr f) {
1448 poly_alloc(f, 2);
1449 element_set1(poly_coeff(f, 1));
1450 element_set0(poly_coeff(f, 0));
1451 }
1452
poly_const_mul(element_ptr res,element_ptr a,element_ptr poly)1453 void poly_const_mul(element_ptr res, element_ptr a, element_ptr poly) {
1454 int i, n = poly_coeff_count(poly);
1455 poly_alloc(res, n);
1456 for (i=0; i<n; i++) {
1457 element_mul(poly_coeff(res, i), a, poly_coeff(poly, i));
1458 }
1459 poly_remove_leading_zeroes(res);
1460 }
1461
poly_random_monic(element_ptr f,int deg)1462 void poly_random_monic(element_ptr f, int deg) {
1463 int i;
1464 poly_alloc(f, deg + 1);
1465 for (i=0; i<deg; i++) {
1466 element_random(poly_coeff(f, i));
1467 }
1468 element_set1(poly_coeff(f, i));
1469 }
1470
polymod_field_degree(field_t f)1471 int polymod_field_degree(field_t f) {
1472 mfptr p = f->data;
1473 return p->n;
1474 }
1475
field_init_polymod(field_ptr f,element_ptr poly)1476 void field_init_polymod(field_ptr f, element_ptr poly) {
1477 pfptr pdp = poly->field->data;
1478 field_init(f);
1479 mfptr p = f->data = pbc_malloc(sizeof(*p));
1480 p->field = pdp->field;
1481 p->mapbase = element_field_to_poly;
1482 element_init(p->poly, poly->field);
1483 element_set(p->poly, poly);
1484 int n = p->n = poly_degree(p->poly);
1485 f->field_clear = field_clear_polymod;
1486 f->init = polymod_init;
1487 f->clear = polymod_clear;
1488 f->set_si = polymod_set_si;
1489 f->set_mpz = polymod_set_mpz;
1490 f->out_str = polymod_out_str;
1491 f->snprint = polymod_snprint;
1492 f->set_multiz = polymod_set_multiz;
1493 f->set_str = polymod_set_str;
1494 f->set = polymod_set;
1495 f->sign = polymod_sgn;
1496 f->add = polymod_add;
1497 f->doub = polymod_double;
1498 f->sub = polymod_sub;
1499 f->neg = polymod_neg;
1500 f->is0 = polymod_is0;
1501 f->is1 = polymod_is1;
1502 f->set0 = polymod_set0;
1503 f->set1 = polymod_set1;
1504 f->cmp = polymod_cmp;
1505 f->to_mpz = polymod_to_mpz;
1506 f->item_count = polymod_coeff_count;
1507 f->item = polymod_coeff;
1508 switch(n) {
1509 case 3:
1510 f->mul = polymod_mul_degree3;
1511 f->square = polymod_square_degree3;
1512 break;
1513 case 6:
1514 f->mul = polymod_mul_degree6;
1515 f->square = polymod_square;
1516 break;
1517 default:
1518 f->mul = polymod_mul;
1519 f->square = polymod_square;
1520 break;
1521 }
1522
1523 f->mul_mpz = polymod_mul_mpz;
1524 f->mul_si = polymod_mul_si;
1525 f->random = polymod_random;
1526 f->from_hash = polymod_from_hash;
1527 f->invert = polymod_invert;
1528 f->is_sqr = polymod_is_sqr;
1529 f->sqrt = polymod_sqrt;
1530 f->to_bytes = polymod_to_bytes;
1531 f->from_bytes = polymod_from_bytes;
1532 f->out_info = polymod_out_info;
1533
1534 if (pdp->field->fixed_length_in_bytes < 0) {
1535 f->fixed_length_in_bytes = -1;
1536 f->length_in_bytes = polymod_length_in_bytes;
1537 } else {
1538 f->fixed_length_in_bytes = pdp->field->fixed_length_in_bytes * poly_degree(poly);
1539 }
1540 mpz_pow_ui(f->order, p->field->order, n);
1541
1542 p->xpwr = pbc_malloc(sizeof(element_t) * n);
1543 compute_x_powers(f, poly);
1544 }
1545
poly_base_field(element_t f)1546 field_ptr poly_base_field(element_t f) {
1547 return ((pfptr) f->field->data)->field;
1548 }
1549
polymod_const_mul(element_ptr res,element_ptr a,element_ptr e)1550 void polymod_const_mul(element_ptr res, element_ptr a, element_ptr e) {
1551 // a lies in R, e in R[x].
1552 element_t *coeff = e->data, *dst = res->data;
1553 int i, n = polymod_field_degree(e->field);
1554
1555 for (i=0; i<n; i++) {
1556 element_mul(dst[i], coeff[i], a);
1557 }
1558 }
1559
1560 struct checkgcd_scope_var {
1561 mpz_ptr z, deg;
1562 field_ptr basef;
1563 element_ptr xpow, x, f, g;
1564 };
1565
1566 // Returns 0 if gcd(x^q^{n/d} - x, f) = 1, 1 otherwise.
checkgcd(mpz_ptr fac,unsigned int mul,struct checkgcd_scope_var * v)1567 static int checkgcd(mpz_ptr fac, unsigned int mul, struct checkgcd_scope_var *v) {
1568 UNUSED_VAR(mul);
1569 mpz_divexact(v->z, v->deg, fac);
1570 mpz_pow_ui(v->z, v->basef->order, mpz_get_ui(v->z));
1571 element_pow_mpz(v->xpow, v->x, v->z);
1572 element_sub(v->xpow, v->xpow, v->x);
1573 if (element_is0(v->xpow)) return 1;
1574 polymod_to_poly(v->g, v->xpow);
1575 poly_gcd(v->g, v->f, v->g);
1576 return poly_degree(v->g) != 0;
1577 }
1578
1579 // Returns 1 if polynomial is irreducible, 0 otherwise.
1580 // A polynomial f(x) is irreducible in F_q[x] if and only if:
1581 // (1) f(x) | x^{q^n} - x, and
1582 // (2) gcd(f(x), x^{q^{n/d}} - x) = 1 for all primes d | n.
1583 // (Recall GF(p) is the splitting field for x^p - x.)
poly_is_irred(element_ptr f)1584 int poly_is_irred(element_ptr f) {
1585 int res = 0;
1586 element_t xpow, x, g;
1587 field_ptr basef = poly_base_field(f);
1588 field_t rxmod;
1589
1590 // 0, units are not irreducibles.
1591 // Assume coefficients are from a field.
1592 if (poly_degree(f) <= 0) return 0;
1593 // Degree 1 polynomials are always irreducible.
1594 if (poly_degree(f) == 1) return 1;
1595
1596 field_init_polymod(rxmod, f);
1597 element_init(xpow, rxmod);
1598 element_init(x, rxmod);
1599 element_init(g, f->field);
1600 element_set1(polymod_coeff(x, 1));
1601
1602 // The degree fits in an unsigned int but I'm lazy and want to use my
1603 // mpz trial division code.
1604 mpz_t deg, z;
1605 mpz_init(deg);
1606 mpz_init(z);
1607 mpz_set_ui(deg, poly_degree(f));
1608
1609 struct checkgcd_scope_var v = {.z = z, .deg = deg, .basef = basef,
1610 .xpow = xpow, .x = x, .f = f, .g = g};
1611 if (!pbc_trial_divide((int(*)(mpz_t,unsigned,void*))checkgcd, &v, deg, NULL)) {
1612 // By now condition (2) has been satisfied. Check (1).
1613 mpz_pow_ui(z, basef->order, poly_degree(f));
1614 element_pow_mpz(xpow, x, z);
1615 element_sub(xpow, xpow, x);
1616 if (element_is0(xpow)) res = 1;
1617 }
1618
1619 mpz_clear(deg);
1620 mpz_clear(z);
1621 element_clear(g);
1622 element_clear(xpow);
1623 element_clear(x);
1624 field_clear(rxmod);
1625 return res;
1626 }
1627
element_field_to_poly(element_ptr f,element_ptr g)1628 void element_field_to_poly(element_ptr f, element_ptr g) {
1629 poly_alloc(f, 1);
1630 element_set(poly_coeff(f, 0), g);
1631 poly_remove_leading_zeroes(f);
1632 }
1633
element_field_to_polymod(element_ptr f,element_ptr g)1634 void element_field_to_polymod(element_ptr f, element_ptr g) {
1635 mfptr p = f->field->data;
1636 element_t *coeff = f->data;
1637 int i, n = p->n;
1638 element_set(coeff[0], g);
1639 for (i=1; i<n; i++) {
1640 element_set0(coeff[i]);
1641 }
1642 }
1643
1644 // Returns 0 when a root exists and sets root to one of the roots.
poly_findroot(element_ptr root,element_ptr poly)1645 int poly_findroot(element_ptr root, element_ptr poly) {
1646 // Compute gcd(x^q - x, poly).
1647 field_t fpxmod;
1648 element_t p, x, r, fac, g;
1649 mpz_t q;
1650
1651 mpz_init(q);
1652 mpz_set(q, poly_base_field(poly)->order);
1653
1654 field_init_polymod(fpxmod, poly);
1655 element_init(p, fpxmod);
1656 element_init(x, fpxmod);
1657 element_init(g, poly->field);
1658 element_set1(((element_t *) x->data)[1]);
1659 pbc_info("findroot: degree %d...", poly_degree(poly));
1660 element_pow_mpz(p, x, q);
1661 element_sub(p, p, x);
1662
1663 polymod_to_poly(g, p);
1664 element_clear(p);
1665 poly_gcd(g, g, poly);
1666 poly_make_monic(g, g);
1667 element_clear(x);
1668 field_clear(fpxmod);
1669
1670 if (!poly_degree(g)) {
1671 printf("no roots!\n");
1672 mpz_clear(q);
1673 element_clear(g);
1674 return -1;
1675 }
1676
1677 // Cantor-Zassenhaus algorithm.
1678 element_init(fac, g->field);
1679 element_init(x, g->field);
1680 element_set_si(x, 1);
1681 mpz_sub_ui(q, q, 1);
1682 mpz_divexact_ui(q, q, 2);
1683 element_init(r, g->field);
1684 for (;;) {
1685 if (poly_degree(g) == 1) break; // Found a root!
1686 step_random:
1687 poly_random_monic(r, 1);
1688 // TODO: evaluate at g instead of bothering with gcd
1689 poly_gcd(fac, r, g);
1690
1691 if (poly_degree(fac) > 0) {
1692 poly_make_monic(g, fac);
1693 } else {
1694 field_init_polymod(fpxmod, g);
1695 int n;
1696 element_init(p, fpxmod);
1697
1698 poly_to_polymod_truncate(p, r);
1699 pbc_info("findroot: degree %d...", poly_degree(g));
1700 element_pow_mpz(p, p, q);
1701
1702 polymod_to_poly(r, p);
1703 element_clear(p);
1704 field_clear(fpxmod);
1705
1706 element_add(r, r, x);
1707 poly_gcd(fac, r, g);
1708 n = poly_degree(fac);
1709 if (n > 0 && n < poly_degree(g)) {
1710 poly_make_monic(g, fac);
1711 } else {
1712 goto step_random;
1713 }
1714 }
1715 }
1716 pbc_info("findroot: found root");
1717 element_neg(root, poly_coeff(g, 0));
1718 element_clear(r);
1719 mpz_clear(q);
1720 element_clear(x);
1721 element_clear(g);
1722 element_clear(fac);
1723 return 0;
1724 }
1725