fb8db84d |
1 | #include <math.h> |
389d8222 |
2 | #include <stdio.h> |
3 | #include <stdlib.h> |
4 | |
5 | #include <mLib/alloc.h> |
6 | |
7 | #include "field.h" |
8 | #include "gf.h" |
9 | #include "mptext.h" |
10 | #include "fibrand.h" |
11 | #include "mprand.h" |
12 | #include "gfn.h" |
13 | |
14 | /*----- Polynomials over finite fields ------------------------------------*/ |
15 | |
16 | typedef struct poly { |
17 | unsigned sz, len; |
18 | mp **v; |
19 | } poly; |
20 | |
21 | #define POLY_INIT { 0, 0, 0 } |
22 | |
23 | #define POLY_ZEROP(p) (!(p)->len) |
24 | #define POLY_CONSTANTP(p) ((p)->len == 1) |
25 | #define POLY_DEGREE(p) ((p)->len - 1) |
26 | |
27 | void poly_free(poly *p) |
28 | { |
29 | unsigned i; |
30 | |
31 | for (i = 0; i < p->sz; i++) |
32 | MP_DROP(p->v[i]); |
33 | xfree(p->v); p->v = 0; |
34 | p->sz = p->len = 0; |
35 | } |
36 | |
37 | void poly_ensure(field *f, poly *p, unsigned sz) |
38 | { |
39 | mp **x; |
40 | unsigned i; |
41 | unsigned n; |
42 | |
43 | if (p->sz >= sz) |
44 | return; |
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++) |
48 | x[i] = p->v[i]; |
49 | for (; i < n; i++) |
50 | x[i] = MP_COPY(f->zero); |
51 | xfree(p->v); |
52 | p->v = x; |
53 | p->sz = n; |
54 | } |
55 | |
56 | void poly_fix(field *f, poly *p, unsigned n) |
57 | { |
58 | while (n && F_ZEROP(f, p->v[n - 1])) |
59 | n--; |
60 | p->len = n; |
61 | } |
62 | |
63 | #define MAX(a, b) ((a) > (b) ? (a) : (b)) |
64 | |
65 | static void putmphex(field *f, const char *name, mp *x) |
66 | { |
67 | mp *z; |
68 | |
69 | printf("%s = 0x", name); |
70 | z = F_OUT(f, MP_NEW, x); |
71 | mp_writefile(z, stdout, 16); |
72 | MP_DROP(z); |
73 | printf("\n"); |
74 | } |
75 | |
76 | void poly_dump(field *f, const char *name, poly *p) |
77 | { |
78 | unsigned i; |
79 | mp *z = MP_NEW; |
80 | |
81 | printf("%s = (", name); |
82 | for (i = p->len; i--; ) { |
83 | printf("0x"); |
84 | z = F_OUT(f, z, p->v[i]); |
85 | mp_writefile(z, stdout, 16); |
86 | if (i) printf(", "); |
87 | } |
88 | mp_drop(z); |
89 | printf(")\n"); |
90 | } |
91 | |
92 | void poly_add(field *f, poly *d, poly *x, poly *y) |
93 | { |
94 | unsigned i; |
95 | unsigned n = MAX(x->len, y->len); |
96 | |
97 | /* poly_dump(f, "add x", x); */ |
98 | /* poly_dump(f, "add y", y); */ |
99 | poly_ensure(f, d, n); |
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); |
104 | } |
105 | poly_fix(f, d, n); |
106 | /* poly_dump(f, "add d", d); */ |
107 | } |
108 | |
109 | void poly_sub(field *f, poly *d, poly *x, poly *y) |
110 | { |
111 | unsigned i; |
112 | unsigned n = MAX(x->len, y->len); |
113 | |
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); |
121 | } |
122 | poly_fix(f, d, n); |
123 | /* poly_dump(f, "sub d", d); */ |
124 | } |
125 | |
126 | static void omuladd(field *f, poly *x, poly *y, mp *z, unsigned o) |
127 | { |
128 | unsigned i; |
129 | unsigned n = MAX(x->len, y->len + o); |
130 | mp *t = MP_NEW; |
131 | |
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); |
142 | } |
143 | } |
144 | mp_drop(t); |
145 | poly_fix(f, x, n); |
146 | /* poly_dump(f, "omuladd d", x); */ |
147 | } |
148 | |
149 | static void omulsub(field *f, poly *x, poly *y, mp *z, unsigned o) |
150 | { |
151 | unsigned i; |
152 | unsigned n = MAX(x->len, y->len + o); |
153 | mp *t = MP_NEW; |
154 | |
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); |
165 | } |
166 | } |
167 | mp_drop(t); |
168 | poly_fix(f, x, n); |
169 | /* poly_dump(f, "omulsub d", x); */ |
170 | } |
171 | |
172 | void poly_mul(field *f, poly *d, poly *x, poly *y) |
173 | { |
174 | poly dd = POLY_INIT; |
175 | unsigned i; |
176 | |
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); |
183 | } |
184 | poly_free(d); |
185 | *d = dd; |
186 | /* poly_dump(f, "mul d", d); */ |
187 | } |
188 | |
189 | void poly_copy(field *f, poly *d, poly *p) |
190 | { |
191 | unsigned i; |
192 | |
193 | poly_ensure(f, d, p->len); |
194 | for (i = 0; i < p->len; i++) { |
195 | MP_DROP(d->v[i]); |
196 | d->v[i] = MP_COPY(p->v[i]); |
197 | } |
198 | d->len = p->len; |
199 | } |
200 | |
201 | void poly_div(field *f, poly *q, poly *r, poly *x, poly *y) |
202 | { |
203 | poly qq = POLY_INIT, rr = POLY_INIT; |
204 | mp *i = MP_NEW; |
205 | mp *z = MP_NEW; |
206 | unsigned j; |
207 | unsigned o, oo = 0; |
208 | |
209 | /* poly_dump(f, "div x", x); */ |
210 | /* poly_dump(f, "div y", y); */ |
211 | assert(y->len); |
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; |
216 | j = rr.len; |
217 | if (q) poly_ensure(f, &qq, o); |
218 | while (o) { |
219 | o--; j--; |
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); |
224 | } |
225 | } |
226 | } |
227 | if (q) { |
228 | poly_fix(f, &qq, oo); |
229 | poly_free(q); |
230 | *q = qq; |
231 | /* poly_dump(f, "div q", q); */ |
232 | } |
233 | /* poly_dump(f, "div r", &rr); */ |
234 | if (!r) |
235 | poly_free(&rr); |
236 | else { |
237 | poly_free(r); |
238 | *r = rr; |
239 | } |
240 | mp_drop(i); |
241 | mp_drop(z); |
242 | } |
243 | |
244 | void poly_monic(field *f, poly *d, poly *p) |
245 | { |
246 | mp *z; |
247 | unsigned i, n; |
248 | |
249 | assert(p->len); |
250 | n = p->len; |
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); |
256 | poly_fix(f, d, n); |
257 | /* poly_dump(f, "monic d", d); */ |
258 | mp_drop(z); |
259 | } |
260 | |
261 | void poly_gcd(field *f, poly *g, poly *x, poly *y) |
262 | { |
263 | poly uu = POLY_INIT, vv = POLY_INIT, rr = POLY_INIT; |
264 | poly *u = &uu, *v = &vv, *r = &rr; |
265 | poly *t; |
266 | |
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; |
274 | } |
275 | poly_monic(f, g, u); |
276 | poly_free(u); |
277 | poly_free(v); |
278 | poly_free(r); |
279 | /* poly_dump(f, "gcd g", g); */ |
280 | } |
281 | |
282 | mp *poly_solve(field *f, mp *d, mp *p, grand *r) |
283 | { |
284 | poly g = POLY_INIT; |
285 | poly ut = POLY_INIT; |
286 | poly c = POLY_INIT; |
287 | poly h = POLY_INIT; |
288 | mp *z; |
289 | unsigned m = f->nbits; |
290 | unsigned i; |
291 | |
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); |
299 | poly_fix(f, &ut, 2); |
300 | poly_copy(f, &c, &ut); |
301 | /* poly_dump(f, "c-in", &c); */ |
302 | |
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); |
45c0fd36 |
307 | /* putchar('.'); fflush(stdout); */ |
389d8222 |
308 | } |
45c0fd36 |
309 | /* poly_dump(f, "c-out", &c); */ |
389d8222 |
310 | poly_gcd(f, &h, &c, &g); |
45c0fd36 |
311 | /* poly_dump(f, "h", &h); */ |
389d8222 |
312 | if (POLY_CONSTANTP(&h) || POLY_DEGREE(&h) == POLY_DEGREE(&g)) |
313 | continue; |
314 | if (2 * POLY_DEGREE(&h) > POLY_DEGREE(&g)) |
315 | poly_div(f, &g, 0, &g, &h); |
316 | else |
317 | poly_copy(f, &g, &h); |
318 | } |
319 | z = MP_COPY(g.v[0]); |
320 | poly_free(&g); |
321 | poly_free(&ut); |
322 | poly_free(&c); |
323 | poly_free(&h); |
324 | mp_drop(d); |
325 | return (z); |
326 | } |
327 | |
328 | /*----- Other stuff -------------------------------------------------------*/ |
329 | |
330 | static int primep(unsigned x) |
331 | { |
332 | unsigned i; |
333 | |
334 | for (i = 2; i*i < x; i++) { |
335 | if (x%i == 0) |
336 | return (0); |
337 | } |
338 | return (1); |
339 | } |
340 | |
341 | static unsigned gcd(unsigned u, unsigned v) |
342 | { |
343 | unsigned r; |
344 | |
345 | if (u < v) { r = u; u = v; v = r; } |
346 | for (;;) { |
347 | r = u%v; |
348 | if (!r) return (v); |
349 | u = v; v = r; |
350 | } |
351 | } |
352 | |
fb8db84d |
353 | static unsigned order(unsigned x, unsigned p) |
354 | { |
355 | unsigned y, k; |
356 | |
357 | if (!x || x == 1) return (0); |
358 | for (y = x, k = 1; y != 1; y = (y*x)%p, k++); |
359 | return (k); |
360 | } |
361 | |
389d8222 |
362 | static int onbtype(unsigned m) |
363 | { |
364 | unsigned t; |
365 | unsigned p, h; |
fb8db84d |
366 | unsigned k, d; |
389d8222 |
367 | |
368 | if (m%8 == 0) |
369 | return (0); |
fb8db84d |
370 | for (t = 1;; t++) { |
389d8222 |
371 | p = t*m + 1; |
372 | if (!primep(p)) |
373 | continue; |
fb8db84d |
374 | k = order(2, p); |
389d8222 |
375 | h = t*m/k; |
376 | d = gcd(h, m); |
377 | if (d == 1) |
378 | return (t); |
379 | } |
380 | return (0); |
381 | } |
382 | |
fb8db84d |
383 | #define PI 3.1415926535897932384626433832795028841971693993751 |
384 | |
385 | static mp *fieldpoly(unsigned m, int t, grand *rr) |
389d8222 |
386 | { |
387 | mp *p, *q, *r, *z; |
388 | unsigned i; |
389 | |
390 | switch (t) { |
391 | case 1: |
392 | p = MP_ZERO; |
393 | for (i = 0; i <= m; i++) |
394 | p = mp_setbit(p, p, i); |
395 | break; |
396 | case 2: |
397 | q = MP_ONE; |
398 | p = MP_THREE; |
399 | r = MP_NEW; |
400 | for (i = 1; i < m; i++) { |
401 | r = mp_lsl(r, p, 1); |
402 | r = mp_xor(r, r, q); |
403 | z = r; r = q; q = p; p = z; |
404 | } |
405 | mp_drop(q); |
406 | mp_drop(r); |
407 | break; |
fb8db84d |
408 | default: { |
409 | #ifdef notdef |
410 | unsigned pp = t*m + 1; |
411 | unsigned uu; |
412 | unsigned j; |
413 | struct cplx { double r, i; }; |
414 | struct cplx e, n; |
415 | struct cplx *f; |
416 | |
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; |
420 | f[0].r = 1; |
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); |
424 | putchar('\n'); |
425 | for (i = 1; i <= m; i++) { |
426 | e.r = e.i = 0; |
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); |
430 | } |
431 | printf("new factor: %g + %g i\n", e.r, e.i); |
432 | for (j = i; j--; ) { |
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; |
437 | f[j] = n; |
438 | } |
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); |
442 | putchar('\n'); |
443 | } |
444 | xfree(f); |
445 | p = MP_ZERO; |
446 | #else |
447 | abort(); |
448 | #endif |
45c0fd36 |
449 | } break; |
389d8222 |
450 | } |
451 | return (p); |
452 | } |
453 | |
454 | static int dofip(unsigned m, mp **p, unsigned n, unsigned x) |
455 | { |
456 | unsigned i; |
457 | |
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)) |
461 | return (1); |
462 | *p = mp_clearbit(*p, *p, i); |
463 | } |
464 | return (0); |
465 | } |
466 | |
467 | static mp *fip(unsigned m) |
468 | { |
469 | mp *p = MP_ONE; |
470 | unsigned n; |
471 | |
472 | p = mp_setbit(p, p, m); |
473 | n = 0; |
474 | while (!dofip(m, &p, n, m)) |
475 | n += 2; |
45c0fd36 |
476 | return (p); |
389d8222 |
477 | } |
478 | |
479 | static mp *fnb(mp *p) |
480 | { |
481 | unsigned t; |
482 | field *f; |
483 | grand *r = fibrand_create(0); |
484 | mp *q, *x; |
485 | unsigned m = mp_bits(p) - 1; |
486 | |
fb8db84d |
487 | if ((t = onbtype(m)) == 0 || t > 2) |
389d8222 |
488 | return (0); |
489 | f = field_binpoly(p); |
fb8db84d |
490 | q = fieldpoly(m, t, r); |
389d8222 |
491 | x = poly_solve(f, MP_NEW, q, r); |
492 | MP_DROP(q); |
493 | F_DESTROY(f); |
494 | return (x); |
495 | } |
496 | |
497 | int main(int argc, char *argv[]) |
498 | { |
499 | int m = atoi(argv[1]); |
500 | int i; |
501 | mp *p; |
502 | mp *beta; |
503 | |
504 | if (!argv[2]) |
505 | p = fip(m); |
506 | else { |
507 | p = MP_ONE; |
508 | p = mp_setbit(p, p, m); |
509 | for (i = 2; i < argc; i++) |
510 | p = mp_setbit(p, p, atoi(argv[i])); |
511 | } |
512 | if (!gf_irreduciblep(p)) { |
513 | printf("not irreducible\n"); |
514 | return (1); |
515 | } |
516 | |
517 | MP_PRINTX("p", p); |
518 | for (i = m + 1; i--; ) { |
519 | if (mp_testbit(p, i)) |
520 | printf("%u ", i); |
521 | } |
522 | putchar('\n'); |
523 | |
524 | beta = fnb(p); |
525 | if (!beta) |
526 | printf("no onb\n"); |
527 | else |
528 | MP_PRINTX("beta", beta); |
529 | |
530 | mp_drop(p); |
531 | mp_drop(beta); |
532 | assert(mparena_count(MPARENA_GLOBAL) == 0); |
533 | return (0); |
534 | } |