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

   1 #include <math.h>
   2 #include <stdio.h>
   3 #include <stdlib.h>
   4 #include <string.h>
   5
   6 #include <mLib/alloc.h>
   7 #include <mLib/darray.h>
   8 #include <mLib/dstr.h>
   9 #include <mLib/quis.h>
  10 #include <mLib/report.h>
  11
  12 #include "fibrand.h"
  13 #include "field.h"
  14 #include "ec.h"
  15 #include "gf.h"
  16 #include "gfreduce.h"
  17 #include "mp.h"
  18 #include "mpbarrett.h"
  19 #include "mpint.h"
  20 #include "mprand.h"
  21 #include "pgen.h"
  22 #include "primetab.h"
  23
  24 DA_DECL(mp_v, mp *);
  25
  26 typedef struct fbe {
  27   unsigned long p;
  28   unsigned e;
  29   mp *n;
  30 } fbe;
  31 DA_DECL(fbe_v, fbe);
  32
  33 typedef struct fact {
  34   mp *p;
  35   unsigned e;
  36   mp *n;
  37   mp *r;
  38 } fact;
  39 DA_DECL(fact_v, fact);
  40
  41 static grand *rng;
  42
  43 /*----- Factoring ---------------------------------------------------------*/
  44
  45 static void mpv_free(mp_v *v)
  46 {
  47   size_t i;
  48
  49   for (i = 0; i < DA_LEN(v); i++)
  50     mp_drop(DA(v)[i]);
  51   DA_DESTROY(v);
  52 }
  53
  54 static mp *smallfactors(mp *x, mp_v *v)
  55 {
  56   mp f;
  57   mpw fw;
  58   mp *y = MP_NEW, *z = MP_NEW;
  59   unsigned i;
  60
  61   MP_COPY(x);
  62   mp_build(&f, &fw, &fw + 1);
  63 #ifdef DEBUG
  64   MP_PRINT("target", x);
  65 #endif
  66   for (i = 0; i < NPRIME; i++) {
  67   again:
  68     fw = primetab[i];
  69     mp_div(&y, &z, x, &f);
  70     if (MP_ZEROP(z)) {
  71 #ifdef DEBUG
  72       MP_PRINT("factor", &f);
  73 #endif
  74       MP_DROP(x);
  75       x = y;
  76       y = MP_NEW;
  77 #ifdef DEBUG
  78       MP_PRINT("remainder", x);
  79 #endif
  80       DA_PUSH(v, mp_fromuint(MP_NEW, primetab[i]));
  81       goto again;
  82     }
  83   }
  84   mp_drop(y);
  85   mp_drop(z);
  86   return (x);
  87 }
  88
  89 #ifdef POLLARDRHO
  90 static mp *pollardrho(mp *x)
  91 {
  92   unsigned i;
  93   mp *y = MP_NEW, *z = MP_NEW, *g = MP_NEW, *t = MP_NEW;
  94   mpbarrett mb;
  95
  96   mpbarrett_create(&mb, x);
  97   y = mprand_range(MP_NEW, x, rng, 0); z = MP_COPY(y);
  98   for (i = 0; i < 1024; i++) {
  99     y = mp_sqr(y, y); y = mpbarrett_reduce(&mb, y, y);
 100     z = mp_sqr(z, z); z = mpbarrett_reduce(&mb, z, z);
 101     z = mp_sqr(z, z); z = mpbarrett_reduce(&mb, z, z);
 102     t = mp_sub(t, y, z);
 103     mp_gcd(&g, 0, 0, x, t);
 104     if (!MP_EQ(g, MP_ONE))
 105       goto done;
 106   }
 107   MP_DROP(g);
 108   g = 0;
 109
 110 done:
 111   mpbarrett_destroy(&mb);
 112   MP_DROP(y);
 113   MP_DROP(z);
 114   MP_DROP(t);
 115 #ifdef DEBUG
 116   MP_PRINT("rho factor", g);
 117 #endif
 118   return (g);
 119 }
 120 #endif
 121
 122 static void addfbe(fbe_v *v, mp *B, unsigned long p)
 123 {
 124   fbe e;
 125   unsigned w;
 126   mp *pp = mp_fromulong(MP_NEW, p);
 127   mp *qq = MP_NEW, *rr = MP_NEW;
 128   mp *t;
 129
 130   w = 0;
 131   qq = MP_COPY(pp);
 132   for (;;) {
 133     w++;
 134     rr = mp_mul(rr, qq, pp);
 135     if (MP_CMP(rr, >=, B)) break;
 136     t = rr; rr = qq; qq = t;
 137   }
 138   e.p = p;
 139   e.e = w;
 140   e.n = qq;
 141   MP_DROP(rr);
 142   MP_DROP(pp);
 143   DA_PUSH(v, e);
 144 }
 145
 146 static void makefactorbase(mp *x, fbe_v *v)
 147 {
 148   mp *B;
 149   unsigned nb;
 150   unsigned long max;
 151   unsigned long i;
 152   size_t j;
 153
 154   nb = sqrt(mp_bits(x)) * 1.5;
 155   B = mp_lsl(MP_NEW, MP_ONE, nb);
 156   if (nb < 32)
 157     max = 1 << nb;
 158   else
 159     max = 0xfffffff0;
 160
 161   addfbe(v, B, 2);
 162   for (i = 3; i < max; i++) {
 163     for (j = 0; j < DA_LEN(v); j++) {
 164       if (i % DA(v)[j].p == 0)
 165         goto composite;
 166     }
 167     addfbe(v, B, i);
 168   composite:;
 169   }
 170
 171 #ifdef DEBUG
 172   for (j = 0; j < DA_LEN(v); j++) {
 173     printf("%lu^%u = ", DA(v)[j].p, DA(v)[j].e);
 174     mp_writefile(DA(v)[j].n, stdout, 10);
 175     putchar('\n');
 176   }
 177 #endif
 178 }
 179
 180 static void freefactorbase(fbe_v *v)
 181 {
 182   size_t i;
 183
 184   for (i = 0; i < DA_LEN(v); i++)
 185     mp_drop(DA(v)[i].n);
 186   DA_DESTROY(v);
 187 }
 188
 189 static mp *ecm(mp *x)
 190 {
 191   fbe_v fb = DA_INIT;
 192   field *f = field_prime(x);
 193   ec_curve *c;
 194   ec p;
 195   mp *a = MP_NEW, *b = MP_NEW;
 196   mp *g = MP_NEW;
 197   unsigned i;
 198   size_t j;
 199
 200   makefactorbase(x, &fb);
 201   for (;;) {
 202     a = mprand_range(a, x, rng, 0);
 203     b = mprand_range(b, x, rng, 0);
 204     c = ec_primeproj(f, a, b);
 205     for (i = 0; i < 1024; i++) {
 206       EC_CREATE(&p);
 207       p.x = mprand_range(p.x, x, rng, 0);
 208       p.y = mprand_range(p.y, x, rng, 0);
 209       p.z = MP_ONE;
 210       for (j = 0; j < DA_LEN(&fb); j++)
 211         ec_imul(c, &p, &p, DA(&fb)[j].n);
 212       if (EC_ATINF(&p))
 213         continue;
 214       mp_gcd(&g, 0, 0, x, p.z);
 215       if (!MP_EQ(g, MP_ONE))
 216         goto done;
 217       EC_DESTROY(&p);
 218     }
 219     ec_destroycurve(c);
 220   }
 221
 222 done:
 223   MP_DROP(a); MP_DROP(b);
 224   EC_DESTROY(&p);
 225   ec_destroycurve(c);
 226   F_DESTROY(f);
 227   freefactorbase(&fb);
 228 #ifdef DEBUG
 229   MP_PRINT("ecm factor", g);
 230 #endif
 231   return (g);
 232 }
 233
 234 static void dofactor(mp *x, mp_v *v)
 235 {
 236   mp *f;
 237
 238   x = smallfactors(x, v);
 239   if (MP_EQ(x, MP_ONE)) return;
 240
 241 #ifdef POLLARDRHO
 242   for (;;) {
 243     if (pgen_primep(x, rng))
 244       goto done;
 245     f = pollardrho(x);
 246     if (!f) break;
 247     dofactor(f, v);
 248     mp_div(&x, 0, x, f);
 249     MP_DROP(f);
 250   }
 251 #endif
 252
 253   for (;;) {
 254     if (pgen_primep(x, rng))
 255       goto done;
 256     f = ecm(x);
 257     if (!f) break;
 258     dofactor(f, v);
 259     mp_div(&x, 0, x, f);
 260     MP_DROP(f);
 261   }
 262
 263 done:
 264   DA_PUSH(v, x);
 265 }
 266
 267 static int cmpmp(const void *a, const void *b)
 268 {
 269   mp *const *aa = a, *const *bb = b;
 270   return (mp_cmp(*aa, *bb));
 271 }
 272
 273 static void factor(mp *x, fact_v *v)
 274 {
 275   fact f;
 276   mp_v fv = DA_INIT;
 277   mp *e = MP_NEW;
 278   size_t i;
 279
 280   dofactor(x, &fv);
 281   qsort(DA(&fv), DA_LEN(&fv), sizeof(mp *), cmpmp);
 282   for (i = 0; i < DA_LEN(&fv); ) {
 283     f.p = MP_COPY(DA(&fv)[i]);
 284     for (f.e = 1, i++;
 285          i < DA_LEN(&fv) && MP_EQ(f.p, DA(&fv)[i]);
 286          f.e++, i++) ;
 287     e = mp_fromuint(e, f.e);
 288     f.n = mp_exp(MP_NEW, f.p, e);
 289     f.r = MP_NEW;
 290     mp_div(&f.r, 0, x, f.p);
 291     DA_PUSH(v, f);
 292   }
 293   mp_drop(e);
 294   mpv_free(&fv);
 295 }
 296
 297 static void freefactors(fact_v *v)
 298 {
 299   size_t i;
 300
 301   for (i = 0; i < DA_LEN(v); i++) {
 302     mp_drop(DA(v)[i].p);
 303     mp_drop(DA(v)[i].n);
 304     mp_drop(DA(v)[i].r);
 305   }
 306   DA_DESTROY(v);
 307 }
 308
 309
 310 /*----- Primitive polynomials ---------------------------------------------*/
 311
 312 static int primitivep(fact_v *f, mp *p)
 313 {
 314   gfreduce r;
 315   unsigned i;
 316   mp *x = MP_NEW;
 317   int rc = 1;
 318
 319   if (!gf_irreduciblep(p))
 320     return (0);
 321 #ifdef DEBUG
 322   MP_PRINTX("p", p);
 323 #endif
 324   gfreduce_create(&r, p);
 325   for (i = 0; i < DA_LEN(f); i++) {
 326     x = gfreduce_exp(&r, x, MP_TWO, DA(f)[i].r);
 327 #ifdef DEBUG
 328     MP_PRINT("  r", DA(f)[i].r);
 329     MP_PRINTX("  x", x);
 330 #endif
 331     if (MP_EQ(x, MP_ONE)) {
 332       rc = 0;
 333       break;
 334     }
 335   }
 336   gfreduce_destroy(&r);
 337   MP_DROP(x);
 338   return (rc);
 339 }
 340
 341 static int dofip(fact_v *f, unsigned m, mp **p, unsigned n, unsigned x)
 342 {
 343   unsigned i;
 344
 345   for (i = n + 1; i < x; i++) {
 346     *p = mp_setbit(*p, *p, i);
 347     if (n ? dofip(f, m, p, n - 1, i) : primitivep(f, *p))
 348       return (1);
 349     *p = mp_clearbit(*p, *p, i);
 350   }
 351   return (0);
 352 }
 353
 354 static mp *fip(fact_v *f, unsigned m)
 355 {
 356   mp *p = MP_ONE;
 357   unsigned n;
 358
 359   p = mp_setbit(p, p, m);
 360   n = 0;
 361   while (!dofip(f, m, &p, n, m))
 362     n += 2;
 363   return (p);
 364 }
 365
 366 static void findprim(unsigned m)
 367 {
 368   fact_v f = DA_INIT;
 369   mp *q = mp_lsl(MP_NEW, MP_ONE, m);
 370   mp *p;
 371   unsigned i;
 372
 373   q = mp_sub(q, q, MP_ONE);
 374   factor(q, &f);
 375
 376 #ifdef DEBUG
 377   {
 378     size_t i;
 379     for (i = 0; i < DA_LEN(&f); i++) {
 380       mp_writefile(DA(&f)[i].p, stdout, 10);
 381       printf("^%u = ", DA(&f)[i].e);
 382       mp_writefile(DA(&f)[i].n, stdout, 10);
 383       fputs(" (", stdout);
 384       mp_writefile(DA(&f)[i].r, stdout, 10);
 385       fputs(")\n", stdout);
 386     }
 387   }
 388 #endif
 389
 390   p = fip(&f, m);
 391   MP_PRINTX("p", p);
 392   for (i = m + 1; i--; ) {
 393     if (mp_testbit(p, i))
 394       printf("%u ", i);
 395   }
 396   putchar('\n');
 397   mp_drop(p);
 398   mp_drop(q);
 399   freefactors(&f);
 400 }
 401
 402 int main(int argc, char *argv[])
 403 {
 404   ego(argv[0]);
 405   rng = fibrand_create(0);
 406   if (argc != 2) {
 407     fprintf(stderr, "Usage: %s M\n", QUIS);
 408     exit(1);
 409   }
 410   findprim(atoi(argv[1]));
 411   return (0);
 412 }