git.distorted.org.uk Git - u/mdw/catacomb/blob - utils/fnb.c

   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 }