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

   1 #include <math.h>
   2
   3 #include "ec.h"
   4 #include "fibrand.h"
   5 #include "field.h"
   6 #include "mpbarrett.h"
   7 #include "mpint.h"
   8 #include "mprand.h"
   9 #include "pgen.h"
  10 #include "primetab.h"
  11
  12 #include "factor.h"
  13
  14 static grand *rng = 0;
  15
  16 DA_DECL(mp_v, mp *);
  17
  18 typedef struct fbe {
  19   unsigned long p;
  20   unsigned e;
  21   mp *n;
  22 } fbe;
  23 DA_DECL(fbe_v, fbe);
  24
  25 static void mpv_free(mp_v *v)
  26 {
  27   size_t i;
  28
  29   for (i = 0; i < DA_LEN(v); i++)
  30     mp_drop(DA(v)[i]);
  31   DA_DESTROY(v);
  32 }
  33
  34 static mp *smallfactors(mp *x, mp_v *v)
  35 {
  36   mp f;
  37   mpw fw;
  38   mp *y = MP_NEW, *z = MP_NEW;
  39   unsigned i;
  40
  41   MP_COPY(x);
  42   mp_build(&f, &fw, &fw + 1);
  43 #ifdef DEBUG
  44   MP_PRINT("target", x);
  45 #endif
  46   for (i = 0; i < NPRIME; i++) {
  47   again:
  48     fw = primetab[i];
  49     mp_div(&y, &z, x, &f);
  50     if (MP_ZEROP(z)) {
  51 #ifdef DEBUG
  52       MP_PRINT("factor", &f);
  53 #endif
  54       MP_DROP(x);
  55       x = y;
  56       y = MP_NEW;
  57 #ifdef DEBUG
  58       MP_PRINT("remainder", x);
  59 #endif
  60       DA_PUSH(v, mp_fromuint(MP_NEW, primetab[i]));
  61       goto again;
  62     }
  63   }
  64   mp_drop(y);
  65   mp_drop(z);
  66   return (x);
  67 }
  68
  69 #ifdef POLLARDRHO
  70 static mp *pollardrho(mp *x)
  71 {
  72   unsigned i;
  73   mp *y = MP_NEW, *z = MP_NEW, *g = MP_NEW, *t = MP_NEW;
  74   mpbarrett mb;
  75
  76   mpbarrett_create(&mb, x);
  77   y = mprand_range(MP_NEW, x, rng, 0); z = MP_COPY(y);
  78   for (i = 0; i < 1024; i++) {
  79     y = mp_sqr(y, y); y = mpbarrett_reduce(&mb, y, y);
  80     z = mp_sqr(z, z); z = mpbarrett_reduce(&mb, z, z);
  81     z = mp_sqr(z, z); z = mpbarrett_reduce(&mb, z, z);
  82     t = mp_sub(t, y, z);
  83     mp_gcd(&g, 0, 0, x, t);
  84     if (!MP_EQ(g, MP_ONE))
  85       goto done;
  86   }
  87   MP_DROP(g);
  88   g = 0;
  89
  90 done:
  91   mpbarrett_destroy(&mb);
  92   MP_DROP(y);
  93   MP_DROP(z);
  94   MP_DROP(t);
  95 #ifdef DEBUG
  96   MP_PRINT("rho factor", g);
  97 #endif
  98   return (g);
  99 }
 100 #endif
 101
 102 static void addfbe(fbe_v *v, mp *B, unsigned long p)
 103 {
 104   fbe e;
 105   unsigned w;
 106   mp *pp = mp_fromulong(MP_NEW, p);
 107   mp *qq = MP_NEW, *rr = MP_NEW;
 108   mp *t;
 109
 110   w = 0;
 111   qq = MP_COPY(pp);
 112   for (;;) {
 113     w++;
 114     rr = mp_mul(rr, qq, pp);
 115     if (MP_CMP(rr, >=, B)) break;
 116     t = rr; rr = qq; qq = t;
 117   }
 118   e.p = p;
 119   e.e = w;
 120   e.n = qq;
 121   MP_DROP(rr);
 122   MP_DROP(pp);
 123   DA_PUSH(v, e);
 124 }
 125
 126 static void makefactorbase(mp *x, fbe_v *v)
 127 {
 128   mp *B;
 129   unsigned nb;
 130   unsigned long max;
 131   unsigned long i;
 132   size_t j;
 133
 134   nb = sqrt(mp_bits(x)) * 1.5;
 135   B = mp_lsl(MP_NEW, MP_ONE, nb);
 136   if (nb < 32)
 137     max = 1 << nb;
 138   else
 139     max = 0xfffffff0;
 140
 141   addfbe(v, B, 2);
 142   for (i = 3; i < max; i++) {
 143     for (j = 0; j < DA_LEN(v); j++) {
 144       if (i % DA(v)[j].p == 0)
 145         goto composite;
 146     }
 147     addfbe(v, B, i);
 148   composite:;
 149   }
 150
 151 #ifdef DEBUG
 152   for (j = 0; j < DA_LEN(v); j++) {
 153     printf("%lu^%u = ", DA(v)[j].p, DA(v)[j].e);
 154     mp_writefile(DA(v)[j].n, stdout, 10);
 155     putchar('\n');
 156   }
 157 #endif
 158 }
 159
 160 static void freefactorbase(fbe_v *v)
 161 {
 162   size_t i;
 163
 164   for (i = 0; i < DA_LEN(v); i++)
 165     mp_drop(DA(v)[i].n);
 166   DA_DESTROY(v);
 167 }
 168
 169 static mp *ecm(mp *x)
 170 {
 171   fbe_v fb = DA_INIT;
 172   field *f = field_prime(x);
 173   ec_curve *c;
 174   ec p;
 175   mp *a = MP_NEW, *b = MP_NEW;
 176   mp *g = MP_NEW;
 177   unsigned i;
 178   size_t j;
 179
 180   makefactorbase(x, &fb);
 181   for (;;) {
 182     a = mprand_range(a, x, rng, 0);
 183     b = mprand_range(b, x, rng, 0);
 184     c = ec_primeproj(f, a, b);
 185     for (i = 0; i < 1024; i++) {
 186       EC_CREATE(&p);
 187       p.x = mprand_range(p.x, x, rng, 0);
 188       p.y = mprand_range(p.y, x, rng, 0);
 189       p.z = MP_ONE;
 190       for (j = 0; j < DA_LEN(&fb); j++)
 191         ec_imul(c, &p, &p, DA(&fb)[j].n);
 192       if (EC_ATINF(&p))
 193         continue;
 194       mp_gcd(&g, 0, 0, x, p.z);
 195       if (!MP_EQ(g, MP_ONE))
 196         goto done;
 197       EC_DESTROY(&p);
 198     }
 199     ec_destroycurve(c);
 200   }
 201
 202 done:
 203   MP_DROP(a); MP_DROP(b);
 204   EC_DESTROY(&p);
 205   ec_destroycurve(c);
 206   F_DESTROY(f);
 207   freefactorbase(&fb);
 208 #ifdef DEBUG
 209   MP_PRINT("ecm factor", g);
 210 #endif
 211   return (g);
 212 }
 213
 214 static void dofactor(mp *x, mp_v *v)
 215 {
 216   mp *f;
 217
 218   x = smallfactors(x, v);
 219   if (MP_EQ(x, MP_ONE)) return;
 220
 221 #ifdef POLLARDRHO
 222   for (;;) {
 223     if (pgen_primep(x, rng))
 224       goto done;
 225     f = pollardrho(x);
 226     if (!f) break;
 227     dofactor(f, v);
 228     mp_div(&x, 0, x, f);
 229     MP_DROP(f);
 230   }
 231 #endif
 232
 233   for (;;) {
 234     if (pgen_primep(x, rng))
 235       goto done;
 236     f = ecm(x);
 237     if (!f) break;
 238     dofactor(f, v);
 239     mp_div(&x, 0, x, f);
 240     MP_DROP(f);
 241   }
 242
 243 done:
 244   DA_PUSH(v, x);
 245 }
 246
 247 static int cmpmp(const void *a, const void *b)
 248 {
 249   mp *const *aa = a, *const *bb = b;
 250   return (mp_cmp(*aa, *bb));
 251 }
 252
 253 void factor(mp *x, fact_v *v)
 254 {
 255   fact f;
 256   mp_v fv = DA_INIT;
 257   mp *e = MP_NEW;
 258   size_t i;
 259
 260   if (!rng) rng = fibrand_create(0);
 261   dofactor(x, &fv);
 262   qsort(DA(&fv), DA_LEN(&fv), sizeof(mp *), cmpmp);
 263   for (i = 0; i < DA_LEN(&fv); ) {
 264     f.p = MP_COPY(DA(&fv)[i]);
 265     for (f.e = 1, i++;
 266          i < DA_LEN(&fv) && MP_EQ(f.p, DA(&fv)[i]);
 267          f.e++, i++) ;
 268     e = mp_fromuint(e, f.e);
 269     f.n = mp_exp(MP_NEW, f.p, e);
 270     f.r = MP_NEW;
 271     mp_div(&f.r, 0, x, f.p);
 272     DA_PUSH(v, f);
 273   }
 274   mp_drop(e);
 275   mpv_free(&fv);
 276 }
 277
 278 void freefactors(fact_v *v)
 279 {
 280   size_t i;
 281
 282   for (i = 0; i < DA_LEN(v); i++) {
 283     mp_drop(DA(v)[i].p);
 284     mp_drop(DA(v)[i].n);
 285     mp_drop(DA(v)[i].r);
 286   }
 287   DA_DESTROY(v);
 288 }