5 #include <mLib/alloc.h>
14 /*----- Polynomials over finite fields ------------------------------------*/
21 #define POLY_INIT { 0, 0, 0 }
23 #define POLY_ZEROP(p) (!(p)->len)
24 #define POLY_CONSTANTP(p) ((p)->len == 1)
25 #define POLY_DEGREE(p) ((p)->len - 1)
27 void poly_free(poly *p)
31 for (i = 0; i < p->sz; i++)
33 xfree(p->v); p->v = 0;
37 void poly_ensure(field *f, poly *p, unsigned sz)
45 n = p->sz; if (!n) n = 2; do n <<= 1; while (n < sz);
46 x = xmalloc(n * sizeof(mp *));
47 for (i = 0; i < p->sz; i++)
50 x[i] = MP_COPY(f->zero);
56 void poly_fix(field *f, poly *p, unsigned n)
58 while (n && F_ZEROP(f, p->v[n - 1]))
63 #define MAX(a, b) ((a) > (b) ? (a) : (b))
65 static void putmphex(field *f, const char *name, mp *x)
69 printf("%s = 0x", name);
70 z = F_OUT(f, MP_NEW, x);
71 mp_writefile(z, stdout, 16);
76 void poly_dump(field *f, const char *name, poly *p)
81 printf("%s = (", name);
82 for (i = p->len; i--; ) {
84 z = F_OUT(f, z, p->v[i]);
85 mp_writefile(z, stdout, 16);
92 void poly_add(field *f, poly *d, poly *x, poly *y)
95 unsigned n = MAX(x->len, y->len);
97 /* poly_dump(f, "add x", x); */
98 /* poly_dump(f, "add y", y); */
100 for (i = 0; i < n; i++) {
101 mp *a = (i < x->len) ? x->v[i] : f->zero;
102 mp *b = (i < y->len) ? y->v[i] : f->zero;
103 d->v[i] = F_ADD(f, d->v[i], a, b);
106 /* poly_dump(f, "add d", d); */
109 void poly_sub(field *f, poly *d, poly *x, poly *y)
112 unsigned n = MAX(x->len, y->len);
114 /* poly_dump(f, "sub x", x); */
115 /* poly_dump(f, "sub y", y); */
116 poly_ensure(f, d, n);
117 for (i = 0; i < n; i++) {
118 mp *a = (i < x->len) ? x->v[i] : f->zero;
119 mp *b = (i < y->len) ? y->v[i] : f->zero;
120 d->v[i] = F_SUB(f, d->v[i], a, b);
123 /* poly_dump(f, "sub d", d); */
126 static void omuladd(field *f, poly *x, poly *y, mp *z, unsigned o)
129 unsigned n = MAX(x->len, y->len + o);
132 /* poly_dump(f, "omuladd x", x); */
133 /* poly_dump(f, "omuladd y", y); */
134 /* putmphex(f, "omuladd z", z); */
135 /* printf("omuladd o = %u\n", o); */
136 poly_ensure(f, x, n);
137 for (i = o; i < n; i++) {
138 if (i < y->len + o) {
139 mp *a = (i < x->len) ? x->v[i] : f->zero;
140 t = F_MUL(f, t, y->v[i - o], z);
141 x->v[i] = F_ADD(f, x->v[i], a, t);
146 /* poly_dump(f, "omuladd d", x); */
149 static void omulsub(field *f, poly *x, poly *y, mp *z, unsigned o)
152 unsigned n = MAX(x->len, y->len + o);
155 /* poly_dump(f, "omulsub x", x); */
156 /* poly_dump(f, "omulsub y", y); */
157 /* putmphex(f, "omulsub z", z); */
158 /* printf("omulsub o = %u\n", o); */
159 poly_ensure(f, x, n);
160 for (i = o; i < n; i++) {
161 if (i < y->len + o) {
162 mp *a = (i < x->len) ? x->v[i] : f->zero;
163 t = F_MUL(f, t, y->v[i - o], z);
164 x->v[i] = F_SUB(f, x->v[i], a, t);
169 /* poly_dump(f, "omulsub d", x); */
172 void poly_mul(field *f, poly *d, poly *x, poly *y)
177 /* poly_dump(f, "mul x", x); */
178 /* poly_dump(f, "mul y", y); */
179 poly_ensure(f, &dd, x->len + y->len);
180 for (i = 0; i < y->len; i++) {
181 if (!F_ZEROP(f, y->v[i]))
182 omuladd(f, &dd, x, y->v[i], i);
186 /* poly_dump(f, "mul d", d); */
189 void poly_copy(field *f, poly *d, poly *p)
193 poly_ensure(f, d, p->len);
194 for (i = 0; i < p->len; i++) {
196 d->v[i] = MP_COPY(p->v[i]);
201 void poly_div(field *f, poly *q, poly *r, poly *x, poly *y)
203 poly qq = POLY_INIT, rr = POLY_INIT;
209 /* poly_dump(f, "div x", x); */
210 /* poly_dump(f, "div y", y); */
212 i = F_INV(f, MP_NEW, y->v[y->len - 1]);
213 poly_copy(f, &rr, x);
214 if (rr.len >= y->len) {
215 o = oo = rr.len - y->len + 1;
217 if (q) poly_ensure(f, &qq, o);
220 if (!F_ZEROP(f, rr.v[j])) {
221 z = F_MUL(f, z, rr.v[j], i);
222 if (q) qq.v[o] = MP_COPY(z);
223 omulsub(f, &rr, y, z, o);
228 poly_fix(f, &qq, oo);
231 /* poly_dump(f, "div q", q); */
233 /* poly_dump(f, "div r", &rr); */
244 void poly_monic(field *f, poly *d, poly *p)
251 /* poly_dump(f, "monic p", p); */
252 poly_ensure(f, d, n);
253 z = F_INV(f, MP_NEW, p->v[n - 1]);
254 for (i = 0; i < n; i++)
255 d->v[i] = F_MUL(f, d->v[i], p->v[i], z);
257 /* poly_dump(f, "monic d", d); */
261 void poly_gcd(field *f, poly *g, poly *x, poly *y)
263 poly uu = POLY_INIT, vv = POLY_INIT, rr = POLY_INIT;
264 poly *u = &uu, *v = &vv, *r = &rr;
267 /* poly_dump(f, "gcd x", x); */
268 /* poly_dump(f, "gcd y", y); */
269 poly_copy(f, u, x); poly_copy(f, v, y);
270 if (u->len < v->len) { t = u; u = v; v = t; }
271 while (!POLY_ZEROP(v)) {
272 poly_div(f, 0, r, u, v);
273 t = u; u = v; v = r; r = t;
279 /* poly_dump(f, "gcd g", g); */
282 mp *poly_solve(field *f, mp *d, mp *p, grand *r)
289 unsigned m = f->nbits;
292 poly_ensure(f, &g, m + 1);
293 for (i = 0; i <= m; i++)
294 g.v[i] = mp_testbit(p, i) ? MP_COPY(f->one) : MP_COPY(f->zero);
295 poly_fix(f, &g, m + 1);
296 poly_ensure(f, &ut, 2);
297 while (POLY_DEGREE(&g) > 1) {
298 ut.v[1] = mprand(ut.v[1], m, r, 0);
300 poly_copy(f, &c, &ut);
301 /* poly_dump(f, "c-in", &c); */
303 for (i = 1; i < m; i++) {
304 poly_mul(f, &c, &c, &c);
305 poly_add(f, &c, &c, &ut);
306 poly_div(f, 0, &c, &c, &g);
307 /* putchar('.'); fflush(stdout); */
309 /* poly_dump(f, "c-out", &c); */
310 poly_gcd(f, &h, &c, &g);
311 /* poly_dump(f, "h", &h); */
312 if (POLY_CONSTANTP(&h) || POLY_DEGREE(&h) == POLY_DEGREE(&g))
314 if (2 * POLY_DEGREE(&h) > POLY_DEGREE(&g))
315 poly_div(f, &g, 0, &g, &h);
317 poly_copy(f, &g, &h);
328 /*----- Other stuff -------------------------------------------------------*/
330 static int primep(unsigned x)
334 for (i = 2; i*i < x; i++) {
341 static unsigned gcd(unsigned u, unsigned v)
345 if (u < v) { r = u; u = v; v = r; }
353 static unsigned order(unsigned x, unsigned p)
357 if (!x || x == 1) return (0);
358 for (y = x, k = 1; y != 1; y = (y*x)%p, k++);
362 static int onbtype(unsigned m)
383 #define PI 3.1415926535897932384626433832795028841971693993751
385 static mp *fieldpoly(unsigned m, int t, grand *rr)
393 for (i = 0; i <= m; i++)
394 p = mp_setbit(p, p, i);
400 for (i = 1; i < m; i++) {
403 z = r; r = q; q = p; p = z;
410 unsigned pp = t*m + 1;
413 struct cplx { double r, i; };
417 do uu = GR_RANGE(rr, pp); while (order(uu, pp) != t);
418 f = xmalloc(sizeof(struct cplx) * (m + 1));
419 for (i = 0; i <= m; i++) f[i].r = f[i].i = 0;
421 printf("poly init; type = %u\n", t);
422 for (i = m + 1; i--; )
423 printf("%3u: %g + %g i\n", i, f[i].r, f[i].i);
425 for (i = 1; i <= m; i++) {
427 for (j = 0; j < t; j++) {
428 double z = (pow(2, i) * pow(uu, j) * PI)/pp;
429 e.r -= cos(z); e.i -= sin(z);
431 printf("new factor: %g + %g i\n", e.r, e.i);
433 f[j + 1].r += f[j].r;
434 f[j + 1].i += f[j].i;
435 n.r = f[j].r * e.r - f[j].i * e.i;
436 n.i = f[j].r * e.i + f[j].i * e.r;
439 printf("poly after %u iters\n", i);
440 for (j = m + 1; j--; )
441 printf("%3u: %g + %g i\n", j, f[j].r, f[j].i);
454 static int dofip(unsigned m, mp **p, unsigned n, unsigned x)
458 for (i = n + 1; i < x; i++) {
459 *p = mp_setbit(*p, *p, i);
460 if (n ? dofip(m, p, n - 1, i) : gf_irreduciblep(*p))
462 *p = mp_clearbit(*p, *p, i);
467 static mp *fip(unsigned m)
472 p = mp_setbit(p, p, m);
474 while (!dofip(m, &p, n, m))
479 static mp *fnb(mp *p)
483 grand *r = fibrand_create(0);
485 unsigned m = mp_bits(p) - 1;
487 if ((t = onbtype(m)) == 0 || t > 2)
489 f = field_binpoly(p);
490 q = fieldpoly(m, t, r);
491 x = poly_solve(f, MP_NEW, q, r);
497 int main(int argc, char *argv[])
499 int m = atoi(argv[1]);
508 p = mp_setbit(p, p, m);
509 for (i = 2; i < argc; i++)
510 p = mp_setbit(p, p, atoi(argv[i]));
512 if (!gf_irreduciblep(p)) {
513 printf("not irreducible\n");
518 for (i = m + 1; i--; ) {
519 if (mp_testbit(p, i))
528 MP_PRINTX("beta", beta);
532 assert(mparena_count(MPARENA_GLOBAL) == 0);