| 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 | } |