chiark - git - mdw - catacomb/blob - utils/fnb.c

   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 }