| 1 | #include <math.h> |
| 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); |
| 307 | /* putchar('.'); fflush(stdout); */ |
| 308 | } |
| 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)) |
| 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 | |
| 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 | |
| 362 | static int onbtype(unsigned m) |
| 363 | { |
| 364 | unsigned t; |
| 365 | unsigned p, h; |
| 366 | unsigned k, d; |
| 367 | |
| 368 | if (m%8 == 0) |
| 369 | return (0); |
| 370 | for (t = 1;; t++) { |
| 371 | p = t*m + 1; |
| 372 | if (!primep(p)) |
| 373 | continue; |
| 374 | k = order(2, p); |
| 375 | h = t*m/k; |
| 376 | d = gcd(h, m); |
| 377 | if (d == 1) |
| 378 | return (t); |
| 379 | } |
| 380 | return (0); |
| 381 | } |
| 382 | |
| 383 | #define PI 3.1415926535897932384626433832795028841971693993751 |
| 384 | |
| 385 | static mp *fieldpoly(unsigned m, int t, grand *rr) |
| 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; |
| 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 |
| 449 | } break; |
| 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; |
| 476 | return (p); |
| 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 | |
| 487 | if ((t = onbtype(m)) == 0 || t > 2) |
| 488 | return (0); |
| 489 | f = field_binpoly(p); |
| 490 | q = fieldpoly(m, t, r); |
| 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 | } |