389d8222 |
1 | #include <stdio.h> |
2 | #include <stdlib.h> |
3 | |
4 | #include <mLib/alloc.h> |
5 | |
6 | #include "field.h" |
7 | #include "gf.h" |
8 | #include "mptext.h" |
9 | #include "fibrand.h" |
10 | #include "mprand.h" |
11 | #include "gfn.h" |
12 | |
13 | /*----- Polynomials over finite fields ------------------------------------*/ |
14 | |
15 | typedef struct poly { |
16 | unsigned sz, len; |
17 | mp **v; |
18 | } poly; |
19 | |
20 | #define POLY_INIT { 0, 0, 0 } |
21 | |
22 | #define POLY_ZEROP(p) (!(p)->len) |
23 | #define POLY_CONSTANTP(p) ((p)->len == 1) |
24 | #define POLY_DEGREE(p) ((p)->len - 1) |
25 | |
26 | void poly_free(poly *p) |
27 | { |
28 | unsigned i; |
29 | |
30 | for (i = 0; i < p->sz; i++) |
31 | MP_DROP(p->v[i]); |
32 | xfree(p->v); p->v = 0; |
33 | p->sz = p->len = 0; |
34 | } |
35 | |
36 | void poly_ensure(field *f, poly *p, unsigned sz) |
37 | { |
38 | mp **x; |
39 | unsigned i; |
40 | unsigned n; |
41 | |
42 | if (p->sz >= sz) |
43 | return; |
44 | n = p->sz; if (!n) n = 2; do n <<= 1; while (n < sz); |
45 | x = xmalloc(n * sizeof(mp *)); |
46 | for (i = 0; i < p->sz; i++) |
47 | x[i] = p->v[i]; |
48 | for (; i < n; i++) |
49 | x[i] = MP_COPY(f->zero); |
50 | xfree(p->v); |
51 | p->v = x; |
52 | p->sz = n; |
53 | } |
54 | |
55 | void poly_fix(field *f, poly *p, unsigned n) |
56 | { |
57 | while (n && F_ZEROP(f, p->v[n - 1])) |
58 | n--; |
59 | p->len = n; |
60 | } |
61 | |
62 | #define MAX(a, b) ((a) > (b) ? (a) : (b)) |
63 | |
64 | static void putmphex(field *f, const char *name, mp *x) |
65 | { |
66 | mp *z; |
67 | |
68 | printf("%s = 0x", name); |
69 | z = F_OUT(f, MP_NEW, x); |
70 | mp_writefile(z, stdout, 16); |
71 | MP_DROP(z); |
72 | printf("\n"); |
73 | } |
74 | |
75 | void poly_dump(field *f, const char *name, poly *p) |
76 | { |
77 | unsigned i; |
78 | mp *z = MP_NEW; |
79 | |
80 | printf("%s = (", name); |
81 | for (i = p->len; i--; ) { |
82 | printf("0x"); |
83 | z = F_OUT(f, z, p->v[i]); |
84 | mp_writefile(z, stdout, 16); |
85 | if (i) printf(", "); |
86 | } |
87 | mp_drop(z); |
88 | printf(")\n"); |
89 | } |
90 | |
91 | void poly_add(field *f, poly *d, poly *x, poly *y) |
92 | { |
93 | unsigned i; |
94 | unsigned n = MAX(x->len, y->len); |
95 | |
96 | /* poly_dump(f, "add x", x); */ |
97 | /* poly_dump(f, "add y", y); */ |
98 | poly_ensure(f, d, n); |
99 | for (i = 0; i < n; i++) { |
100 | mp *a = (i < x->len) ? x->v[i] : f->zero; |
101 | mp *b = (i < y->len) ? y->v[i] : f->zero; |
102 | d->v[i] = F_ADD(f, d->v[i], a, b); |
103 | } |
104 | poly_fix(f, d, n); |
105 | /* poly_dump(f, "add d", d); */ |
106 | } |
107 | |
108 | void poly_sub(field *f, poly *d, poly *x, poly *y) |
109 | { |
110 | unsigned i; |
111 | unsigned n = MAX(x->len, y->len); |
112 | |
113 | /* poly_dump(f, "sub x", x); */ |
114 | /* poly_dump(f, "sub y", y); */ |
115 | poly_ensure(f, d, n); |
116 | for (i = 0; i < n; i++) { |
117 | mp *a = (i < x->len) ? x->v[i] : f->zero; |
118 | mp *b = (i < y->len) ? y->v[i] : f->zero; |
119 | d->v[i] = F_SUB(f, d->v[i], a, b); |
120 | } |
121 | poly_fix(f, d, n); |
122 | /* poly_dump(f, "sub d", d); */ |
123 | } |
124 | |
125 | static void omuladd(field *f, poly *x, poly *y, mp *z, unsigned o) |
126 | { |
127 | unsigned i; |
128 | unsigned n = MAX(x->len, y->len + o); |
129 | mp *t = MP_NEW; |
130 | |
131 | /* poly_dump(f, "omuladd x", x); */ |
132 | /* poly_dump(f, "omuladd y", y); */ |
133 | /* putmphex(f, "omuladd z", z); */ |
134 | /* printf("omuladd o = %u\n", o); */ |
135 | poly_ensure(f, x, n); |
136 | for (i = o; i < n; i++) { |
137 | if (i < y->len + o) { |
138 | mp *a = (i < x->len) ? x->v[i] : f->zero; |
139 | t = F_MUL(f, t, y->v[i - o], z); |
140 | x->v[i] = F_ADD(f, x->v[i], a, t); |
141 | } |
142 | } |
143 | mp_drop(t); |
144 | poly_fix(f, x, n); |
145 | /* poly_dump(f, "omuladd d", x); */ |
146 | } |
147 | |
148 | static void omulsub(field *f, poly *x, poly *y, mp *z, unsigned o) |
149 | { |
150 | unsigned i; |
151 | unsigned n = MAX(x->len, y->len + o); |
152 | mp *t = MP_NEW; |
153 | |
154 | /* poly_dump(f, "omulsub x", x); */ |
155 | /* poly_dump(f, "omulsub y", y); */ |
156 | /* putmphex(f, "omulsub z", z); */ |
157 | /* printf("omulsub o = %u\n", o); */ |
158 | poly_ensure(f, x, n); |
159 | for (i = o; i < n; i++) { |
160 | if (i < y->len + o) { |
161 | mp *a = (i < x->len) ? x->v[i] : f->zero; |
162 | t = F_MUL(f, t, y->v[i - o], z); |
163 | x->v[i] = F_SUB(f, x->v[i], a, t); |
164 | } |
165 | } |
166 | mp_drop(t); |
167 | poly_fix(f, x, n); |
168 | /* poly_dump(f, "omulsub d", x); */ |
169 | } |
170 | |
171 | void poly_mul(field *f, poly *d, poly *x, poly *y) |
172 | { |
173 | poly dd = POLY_INIT; |
174 | unsigned i; |
175 | |
176 | /* poly_dump(f, "mul x", x); */ |
177 | /* poly_dump(f, "mul y", y); */ |
178 | poly_ensure(f, &dd, x->len + y->len); |
179 | for (i = 0; i < y->len; i++) { |
180 | if (!F_ZEROP(f, y->v[i])) |
181 | omuladd(f, &dd, x, y->v[i], i); |
182 | } |
183 | poly_free(d); |
184 | *d = dd; |
185 | /* poly_dump(f, "mul d", d); */ |
186 | } |
187 | |
188 | void poly_copy(field *f, poly *d, poly *p) |
189 | { |
190 | unsigned i; |
191 | |
192 | poly_ensure(f, d, p->len); |
193 | for (i = 0; i < p->len; i++) { |
194 | MP_DROP(d->v[i]); |
195 | d->v[i] = MP_COPY(p->v[i]); |
196 | } |
197 | d->len = p->len; |
198 | } |
199 | |
200 | void poly_div(field *f, poly *q, poly *r, poly *x, poly *y) |
201 | { |
202 | poly qq = POLY_INIT, rr = POLY_INIT; |
203 | mp *i = MP_NEW; |
204 | mp *z = MP_NEW; |
205 | unsigned j; |
206 | unsigned o, oo = 0; |
207 | |
208 | /* poly_dump(f, "div x", x); */ |
209 | /* poly_dump(f, "div y", y); */ |
210 | assert(y->len); |
211 | i = F_INV(f, MP_NEW, y->v[y->len - 1]); |
212 | poly_copy(f, &rr, x); |
213 | if (rr.len >= y->len) { |
214 | o = oo = rr.len - y->len + 1; |
215 | j = rr.len; |
216 | if (q) poly_ensure(f, &qq, o); |
217 | while (o) { |
218 | o--; j--; |
219 | if (!F_ZEROP(f, rr.v[j])) { |
220 | z = F_MUL(f, z, rr.v[j], i); |
221 | if (q) qq.v[o] = MP_COPY(z); |
222 | omulsub(f, &rr, y, z, o); |
223 | } |
224 | } |
225 | } |
226 | if (q) { |
227 | poly_fix(f, &qq, oo); |
228 | poly_free(q); |
229 | *q = qq; |
230 | /* poly_dump(f, "div q", q); */ |
231 | } |
232 | /* poly_dump(f, "div r", &rr); */ |
233 | if (!r) |
234 | poly_free(&rr); |
235 | else { |
236 | poly_free(r); |
237 | *r = rr; |
238 | } |
239 | mp_drop(i); |
240 | mp_drop(z); |
241 | } |
242 | |
243 | void poly_monic(field *f, poly *d, poly *p) |
244 | { |
245 | mp *z; |
246 | unsigned i, n; |
247 | |
248 | assert(p->len); |
249 | n = p->len; |
250 | /* poly_dump(f, "monic p", p); */ |
251 | poly_ensure(f, d, n); |
252 | z = F_INV(f, MP_NEW, p->v[n - 1]); |
253 | for (i = 0; i < n; i++) |
254 | d->v[i] = F_MUL(f, d->v[i], p->v[i], z); |
255 | poly_fix(f, d, n); |
256 | /* poly_dump(f, "monic d", d); */ |
257 | mp_drop(z); |
258 | } |
259 | |
260 | void poly_gcd(field *f, poly *g, poly *x, poly *y) |
261 | { |
262 | poly uu = POLY_INIT, vv = POLY_INIT, rr = POLY_INIT; |
263 | poly *u = &uu, *v = &vv, *r = &rr; |
264 | poly *t; |
265 | |
266 | /* poly_dump(f, "gcd x", x); */ |
267 | /* poly_dump(f, "gcd y", y); */ |
268 | poly_copy(f, u, x); poly_copy(f, v, y); |
269 | if (u->len < v->len) { t = u; u = v; v = t; } |
270 | while (!POLY_ZEROP(v)) { |
271 | poly_div(f, 0, r, u, v); |
272 | t = u; u = v; v = r; r = t; |
273 | } |
274 | poly_monic(f, g, u); |
275 | poly_free(u); |
276 | poly_free(v); |
277 | poly_free(r); |
278 | /* poly_dump(f, "gcd g", g); */ |
279 | } |
280 | |
281 | mp *poly_solve(field *f, mp *d, mp *p, grand *r) |
282 | { |
283 | poly g = POLY_INIT; |
284 | poly ut = POLY_INIT; |
285 | poly c = POLY_INIT; |
286 | poly h = POLY_INIT; |
287 | mp *z; |
288 | unsigned m = f->nbits; |
289 | unsigned i; |
290 | |
291 | poly_ensure(f, &g, m + 1); |
292 | for (i = 0; i <= m; i++) |
293 | g.v[i] = mp_testbit(p, i) ? MP_COPY(f->one) : MP_COPY(f->zero); |
294 | poly_fix(f, &g, m + 1); |
295 | poly_ensure(f, &ut, 2); |
296 | while (POLY_DEGREE(&g) > 1) { |
297 | ut.v[1] = mprand(ut.v[1], m, r, 0); |
298 | poly_fix(f, &ut, 2); |
299 | poly_copy(f, &c, &ut); |
300 | /* poly_dump(f, "c-in", &c); */ |
301 | |
302 | for (i = 1; i < m; i++) { |
303 | poly_mul(f, &c, &c, &c); |
304 | poly_add(f, &c, &c, &ut); |
305 | poly_div(f, 0, &c, &c, &g); |
306 | /* putchar('.'); fflush(stdout); */ |
307 | } |
308 | /* poly_dump(f, "c-out", &c); */ |
309 | poly_gcd(f, &h, &c, &g); |
310 | /* poly_dump(f, "h", &h); */ |
311 | if (POLY_CONSTANTP(&h) || POLY_DEGREE(&h) == POLY_DEGREE(&g)) |
312 | continue; |
313 | if (2 * POLY_DEGREE(&h) > POLY_DEGREE(&g)) |
314 | poly_div(f, &g, 0, &g, &h); |
315 | else |
316 | poly_copy(f, &g, &h); |
317 | } |
318 | z = MP_COPY(g.v[0]); |
319 | poly_free(&g); |
320 | poly_free(&ut); |
321 | poly_free(&c); |
322 | poly_free(&h); |
323 | mp_drop(d); |
324 | return (z); |
325 | } |
326 | |
327 | /*----- Other stuff -------------------------------------------------------*/ |
328 | |
329 | static int primep(unsigned x) |
330 | { |
331 | unsigned i; |
332 | |
333 | for (i = 2; i*i < x; i++) { |
334 | if (x%i == 0) |
335 | return (0); |
336 | } |
337 | return (1); |
338 | } |
339 | |
340 | static unsigned gcd(unsigned u, unsigned v) |
341 | { |
342 | unsigned r; |
343 | |
344 | if (u < v) { r = u; u = v; v = r; } |
345 | for (;;) { |
346 | r = u%v; |
347 | if (!r) return (v); |
348 | u = v; v = r; |
349 | } |
350 | } |
351 | |
352 | static int onbtype(unsigned m) |
353 | { |
354 | unsigned t; |
355 | unsigned p, h; |
356 | unsigned k, x, d; |
357 | |
358 | if (m%8 == 0) |
359 | return (0); |
360 | for (t = 1; t <= 2; t++) { |
361 | p = t*m + 1; |
362 | if (!primep(p)) |
363 | continue; |
364 | for (x = 2, k = 1; x != 1; x = (2*x)%p, k++); |
365 | h = t*m/k; |
366 | d = gcd(h, m); |
367 | if (d == 1) |
368 | return (t); |
369 | } |
370 | return (0); |
371 | } |
372 | |
373 | static mp *fieldpoly(unsigned m, int t) |
374 | { |
375 | mp *p, *q, *r, *z; |
376 | unsigned i; |
377 | |
378 | switch (t) { |
379 | case 1: |
380 | p = MP_ZERO; |
381 | for (i = 0; i <= m; i++) |
382 | p = mp_setbit(p, p, i); |
383 | break; |
384 | case 2: |
385 | q = MP_ONE; |
386 | p = MP_THREE; |
387 | r = MP_NEW; |
388 | for (i = 1; i < m; i++) { |
389 | r = mp_lsl(r, p, 1); |
390 | r = mp_xor(r, r, q); |
391 | z = r; r = q; q = p; p = z; |
392 | } |
393 | mp_drop(q); |
394 | mp_drop(r); |
395 | break; |
396 | default: |
397 | abort(); |
398 | } |
399 | return (p); |
400 | } |
401 | |
402 | static int dofip(unsigned m, mp **p, unsigned n, unsigned x) |
403 | { |
404 | unsigned i; |
405 | |
406 | for (i = n + 1; i < x; i++) { |
407 | *p = mp_setbit(*p, *p, i); |
408 | if (n ? dofip(m, p, n - 1, i) : gf_irreduciblep(*p)) |
409 | return (1); |
410 | *p = mp_clearbit(*p, *p, i); |
411 | } |
412 | return (0); |
413 | } |
414 | |
415 | static mp *fip(unsigned m) |
416 | { |
417 | mp *p = MP_ONE; |
418 | unsigned n; |
419 | |
420 | p = mp_setbit(p, p, m); |
421 | n = 0; |
422 | while (!dofip(m, &p, n, m)) |
423 | n += 2; |
424 | return (p); |
425 | } |
426 | |
427 | static mp *fnb(mp *p) |
428 | { |
429 | unsigned t; |
430 | field *f; |
431 | grand *r = fibrand_create(0); |
432 | mp *q, *x; |
433 | unsigned m = mp_bits(p) - 1; |
434 | |
435 | if ((t = onbtype(m)) == 0) |
436 | return (0); |
437 | f = field_binpoly(p); |
438 | q = fieldpoly(m, t); |
439 | x = poly_solve(f, MP_NEW, q, r); |
440 | MP_DROP(q); |
441 | F_DESTROY(f); |
442 | return (x); |
443 | } |
444 | |
445 | int main(int argc, char *argv[]) |
446 | { |
447 | int m = atoi(argv[1]); |
448 | int i; |
449 | mp *p; |
450 | mp *beta; |
451 | |
452 | if (!argv[2]) |
453 | p = fip(m); |
454 | else { |
455 | p = MP_ONE; |
456 | p = mp_setbit(p, p, m); |
457 | for (i = 2; i < argc; i++) |
458 | p = mp_setbit(p, p, atoi(argv[i])); |
459 | } |
460 | if (!gf_irreduciblep(p)) { |
461 | printf("not irreducible\n"); |
462 | return (1); |
463 | } |
464 | |
465 | MP_PRINTX("p", p); |
466 | for (i = m + 1; i--; ) { |
467 | if (mp_testbit(p, i)) |
468 | printf("%u ", i); |
469 | } |
470 | putchar('\n'); |
471 | |
472 | beta = fnb(p); |
473 | if (!beta) |
474 | printf("no onb\n"); |
475 | else |
476 | MP_PRINTX("beta", beta); |
477 | |
478 | mp_drop(p); |
479 | mp_drop(beta); |
480 | assert(mparena_count(MPARENA_GLOBAL) == 0); |
481 | return (0); |
482 | } |