utils: Make very bad ECM factoring program.

  * Extract factoring code from existing `prim' program.

  * Write driver front-end.

diff --git a/utils/fact.c b/utils/fact.c
new file mode 100644 (file)
index 0000000..908ce26
--- /dev/null
+++ b/utils/fact.c
@@ -0,0 +1,28 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include <mLib/quis.h>
+
+#include "factor.h"
+
+int main(int argc, char *argv[])
+{
+  mp *n;
+  fact_v fv = DA_INIT;
+  unsigned i;
+
+  ego(argv[0]);
+  if (argc != 2) {
+    fprintf(stderr, "Usage: %s N\n", QUIS);
+    exit(1);
+  }
+  n = mp_readstring(MP_NEW, argv[1], 0, 10);
+  factor(n, &fv);
+  for (i = 0; i < DA_LEN(&fv); i++) {
+    mp_writefile(DA(&fv)[i].p, stdout, 10);
+    printf("^%u\n", DA(&fv)[i].e);
+  }
+  freefactors(&fv);
+  return (0);
+}

diff --git a/utils/factor.c b/utils/factor.c
new file mode 100644 (file)
index 0000000..5f29c1a
--- /dev/null
+++ b/utils/factor.c
@@ -0,0 +1,289 @@
+#include <math.h>
+
+#include "ec.h"
+#include "fibrand.h"
+#include "field.h"
+#include "mpbarrett.h"
+#include "mpint.h"
+#include "mprand.h"
+#include "pgen.h"
+#include "primetab.h"
+
+#include "factor.h"
+
+static grand *rng = 0;
+
+DA_DECL(mp_v, mp *);
+
+typedef struct fbe {
+  unsigned long p;
+  unsigned e;
+  mp *n;
+} fbe;
+DA_DECL(fbe_v, fbe);
+
+static void mpv_free(mp_v *v)
+{
+  size_t i;
+
+  for (i = 0; i < DA_LEN(v); i++)
+    mp_drop(DA(v)[i]);
+  DA_DESTROY(v);
+}
+
+static mp *smallfactors(mp *x, mp_v *v)
+{
+  mp f;
+  mpw fw;
+  mp *y = MP_NEW, *z = MP_NEW;
+  unsigned i;
+
+  MP_COPY(x);
+  mp_build(&f, &fw, &fw + 1);
+#ifdef DEBUG
+  MP_PRINT("target", x);
+#endif
+  for (i = 0; i < NPRIME; i++) {
+  again:
+    fw = primetab[i];
+    mp_div(&y, &z, x, &f);
+    if (MP_ZEROP(z)) {
+#ifdef DEBUG
+      MP_PRINT("factor", &f);
+#endif
+      MP_DROP(x);
+      x = y;
+      y = MP_NEW;
+#ifdef DEBUG
+      MP_PRINT("remainder", x);
+#endif
+      DA_PUSH(v, mp_fromuint(MP_NEW, primetab[i]));
+      goto again;
+    }
+  }
+  mp_drop(y);
+  mp_drop(z);
+  return (x);
+}
+
+#ifdef POLLARDRHO
+static mp *pollardrho(mp *x)
+{
+  unsigned i;
+  mp *y = MP_NEW, *z = MP_NEW, *g = MP_NEW, *t = MP_NEW;
+  mpbarrett mb;
+
+  mpbarrett_create(&mb, x);
+  y = mprand_range(MP_NEW, x, rng, 0); z = MP_COPY(y);
+  for (i = 0; i < 1024; i++) {
+    y = mp_sqr(y, y); y = mpbarrett_reduce(&mb, y, y);
+    z = mp_sqr(z, z); z = mpbarrett_reduce(&mb, z, z);
+    z = mp_sqr(z, z); z = mpbarrett_reduce(&mb, z, z);
+    t = mp_sub(t, y, z);
+    mp_gcd(&g, 0, 0, x, t);
+    if (!MP_EQ(g, MP_ONE))
+      goto done;
+  }
+  MP_DROP(g);
+  g = 0;
+
+done:
+  mpbarrett_destroy(&mb);
+  MP_DROP(y);
+  MP_DROP(z);
+  MP_DROP(t);
+#ifdef DEBUG
+  MP_PRINT("rho factor", g);
+#endif
+  return (g);
+}
+#endif
+
+static void addfbe(fbe_v *v, mp *B, unsigned long p)
+{
+  fbe e;
+  unsigned w;
+  mp *pp = mp_fromulong(MP_NEW, p);
+  mp *qq = MP_NEW, *rr = MP_NEW;
+  mp *t;
+
+  w = 0;
+  qq = MP_COPY(pp);
+  for (;;) {
+    w++;
+    rr = mp_mul(rr, qq, pp);
+    if (MP_CMP(rr, >=, B)) break;
+    t = rr; rr = qq; qq = t;
+  }
+  e.p = p;
+  e.e = w;
+  e.n = qq;
+  MP_DROP(rr);
+  MP_DROP(pp);
+  DA_PUSH(v, e);
+}
+
+static void makefactorbase(mp *x, fbe_v *v)
+{
+  mp *B;
+  unsigned nb;
+  unsigned long max;
+  unsigned long i;
+  size_t j;
+
+  nb = sqrt(mp_bits(x)) * 1.5;
+  B = mp_lsl(MP_NEW, MP_ONE, nb);
+  if (nb < 32)
+    max = 1 << nb;
+  else
+    max = 0xfffffff0;
+
+  addfbe(v, B, 2);
+  for (i = 3; i < max; i++) {
+    for (j = 0; j < DA_LEN(v); j++) {
+      if (i % DA(v)[j].p == 0)
+       goto composite;
+    }
+    addfbe(v, B, i);
+  composite:;
+  }
+
+#ifdef DEBUG
+  for (j = 0; j < DA_LEN(v); j++) {
+    printf("%lu^%u = ", DA(v)[j].p, DA(v)[j].e);
+    mp_writefile(DA(v)[j].n, stdout, 10);
+    putchar('\n');
+  }
+#endif
+}
+
+static void freefactorbase(fbe_v *v)
+{
+  size_t i;
+
+  for (i = 0; i < DA_LEN(v); i++)
+    mp_drop(DA(v)[i].n);
+  DA_DESTROY(v);    
+}
+
+static mp *ecm(mp *x)
+{
+  fbe_v fb = DA_INIT;
+  field *f = field_prime(x);
+  ec_curve *c;
+  ec p;
+  mp *a = MP_NEW, *b = MP_NEW;
+  mp *g = MP_NEW;
+  unsigned i;
+  size_t j;
+
+  makefactorbase(x, &fb);
+  for (;;) {
+    a = mprand_range(a, x, rng, 0);
+    b = mprand_range(b, x, rng, 0);
+    c = ec_primeproj(f, a, b);
+    for (i = 0; i < 1024; i++) {
+      EC_CREATE(&p);
+      p.x = mprand_range(p.x, x, rng, 0);
+      p.y = mprand_range(p.y, x, rng, 0);
+      p.z = MP_ONE;
+      for (j = 0; j < DA_LEN(&fb); j++)
+       ec_imul(c, &p, &p, DA(&fb)[j].n);
+      if (EC_ATINF(&p))
+       continue;
+      mp_gcd(&g, 0, 0, x, p.z);
+      if (!MP_EQ(g, MP_ONE))
+       goto done;
+      EC_DESTROY(&p);
+    }
+    ec_destroycurve(c);
+  }
+
+done:
+  MP_DROP(a); MP_DROP(b);
+  EC_DESTROY(&p);
+  ec_destroycurve(c);
+  F_DESTROY(f);
+  freefactorbase(&fb);
+#ifdef DEBUG
+  MP_PRINT("ecm factor", g);
+#endif
+  return (g);
+}
+
+static void dofactor(mp *x, mp_v *v)
+{
+  mp *f;
+
+  x = smallfactors(x, v); 
+  if (MP_EQ(x, MP_ONE)) return;
+
+#ifdef POLLARDRHO
+  for (;;) {
+    if (pgen_primep(x, rng))
+      goto done;
+    f = pollardrho(x);
+    if (!f) break;
+    dofactor(f, v);
+    mp_div(&x, 0, x, f);
+    MP_DROP(f);
+  }
+#endif
+
+  for (;;) {
+    if (pgen_primep(x, rng))
+      goto done;
+    f = ecm(x);
+    if (!f) break;
+    dofactor(f, v);
+    mp_div(&x, 0, x, f);
+    MP_DROP(f);
+  }
+
+done:
+  DA_PUSH(v, x);
+}
+
+static int cmpmp(const void *a, const void *b)
+{
+  mp *const *aa = a, *const *bb = b;
+  return (mp_cmp(*aa, *bb));
+}
+
+void factor(mp *x, fact_v *v)
+{
+  fact f;
+  mp_v fv = DA_INIT;
+  mp *e = MP_NEW;
+  size_t i;
+
+  if (!rng) rng = fibrand_create(0);
+  dofactor(x, &fv);
+  qsort(DA(&fv), DA_LEN(&fv), sizeof(mp *), cmpmp);
+  for (i = 0; i < DA_LEN(&fv); ) {
+    f.p = MP_COPY(DA(&fv)[i]);
+    for (f.e = 1, i++;
+        i < DA_LEN(&fv) && MP_EQ(f.p, DA(&fv)[i]);
+        f.e++, i++) ;
+    e = mp_fromuint(e, f.e);
+    f.n = mp_exp(MP_NEW, f.p, e);
+    f.r = MP_NEW;
+    mp_div(&f.r, 0, x, f.p);
+    DA_PUSH(v, f);
+  }
+  mp_drop(e);
+  mpv_free(&fv);
+}
+
+void freefactors(fact_v *v)
+{
+  size_t i;
+
+  for (i = 0; i < DA_LEN(v); i++) {
+    mp_drop(DA(v)[i].p);
+    mp_drop(DA(v)[i].n);
+    mp_drop(DA(v)[i].r);
+  }
+  DA_DESTROY(v);
+}
+

diff --git a/utils/factor.h b/utils/factor.h
new file mode 100644 (file)
index 0000000..7181505
--- /dev/null
+++ b/utils/factor.h
@@ -0,0 +1,19 @@
+#ifndef FACTOR_H
+#define FACTOR_H
+
+#include <mLib/darray.h>
+
+#include "mp.h"
+
+typedef struct fact {
+  mp *p;
+  unsigned e;
+  mp *n;
+  mp *r;
+} fact;
+DA_DECL(fact_v, fact);
+
+extern void factor(mp *, fact_v *);
+extern void freefactors(fact_v *);
+
+#endif

diff --git a/utils/prim.c b/utils/prim.c
index b279304..39b7e6f 100644 (file)
--- a/utils/prim.c
+++ b/utils/prim.c
@@ -1,4 +1,3 @@
-#include <math.h>
 #include <stdio.h>
 #include <stdlib.h>
 #include <string.h>
@@ -9,305 +8,11 @@
 #include <mLib/quis.h>
 #include <mLib/report.h>
 
-#include "fibrand.h"
-#include "field.h"
-#include "ec.h"
 #include "gf.h"
 #include "gfreduce.h"
 #include "mp.h"
-#include "mpbarrett.h"
-#include "mpint.h"
-#include "mprand.h"
-#include "pgen.h"
-#include "primetab.h"
 
-DA_DECL(mp_v, mp *);
-
-typedef struct fbe {
-  unsigned long p;
-  unsigned e;
-  mp *n;
-} fbe;
-DA_DECL(fbe_v, fbe);
-
-typedef struct fact {
-  mp *p;
-  unsigned e;
-  mp *n;
-  mp *r;
-} fact;
-DA_DECL(fact_v, fact);
-
-static grand *rng;
-
-/*----- Factoring ---------------------------------------------------------*/
-
-static void mpv_free(mp_v *v)
-{
-  size_t i;
-
-  for (i = 0; i < DA_LEN(v); i++)
-    mp_drop(DA(v)[i]);
-  DA_DESTROY(v);
-}
-
-static mp *smallfactors(mp *x, mp_v *v)
-{
-  mp f;
-  mpw fw;
-  mp *y = MP_NEW, *z = MP_NEW;
-  unsigned i;
-
-  MP_COPY(x);
-  mp_build(&f, &fw, &fw + 1);
-#ifdef DEBUG
-  MP_PRINT("target", x);
-#endif
-  for (i = 0; i < NPRIME; i++) {
-  again:
-    fw = primetab[i];
-    mp_div(&y, &z, x, &f);
-    if (MP_ZEROP(z)) {
-#ifdef DEBUG
-      MP_PRINT("factor", &f);
-#endif
-      MP_DROP(x);
-      x = y;
-      y = MP_NEW;
-#ifdef DEBUG
-      MP_PRINT("remainder", x);
-#endif
-      DA_PUSH(v, mp_fromuint(MP_NEW, primetab[i]));
-      goto again;
-    }
-  }
-  mp_drop(y);
-  mp_drop(z);
-  return (x);
-}
-
-#ifdef POLLARDRHO
-static mp *pollardrho(mp *x)
-{
-  unsigned i;
-  mp *y = MP_NEW, *z = MP_NEW, *g = MP_NEW, *t = MP_NEW;
-  mpbarrett mb;
-
-  mpbarrett_create(&mb, x);
-  y = mprand_range(MP_NEW, x, rng, 0); z = MP_COPY(y);
-  for (i = 0; i < 1024; i++) {
-    y = mp_sqr(y, y); y = mpbarrett_reduce(&mb, y, y);
-    z = mp_sqr(z, z); z = mpbarrett_reduce(&mb, z, z);
-    z = mp_sqr(z, z); z = mpbarrett_reduce(&mb, z, z);
-    t = mp_sub(t, y, z);
-    mp_gcd(&g, 0, 0, x, t);
-    if (!MP_EQ(g, MP_ONE))
-      goto done;
-  }
-  MP_DROP(g);
-  g = 0;
-
-done:
-  mpbarrett_destroy(&mb);
-  MP_DROP(y);
-  MP_DROP(z);
-  MP_DROP(t);
-#ifdef DEBUG
-  MP_PRINT("rho factor", g);
-#endif
-  return (g);
-}
-#endif
-
-static void addfbe(fbe_v *v, mp *B, unsigned long p)
-{
-  fbe e;
-  unsigned w;
-  mp *pp = mp_fromulong(MP_NEW, p);
-  mp *qq = MP_NEW, *rr = MP_NEW;
-  mp *t;
-
-  w = 0;
-  qq = MP_COPY(pp);
-  for (;;) {
-    w++;
-    rr = mp_mul(rr, qq, pp);
-    if (MP_CMP(rr, >=, B)) break;
-    t = rr; rr = qq; qq = t;
-  }
-  e.p = p;
-  e.e = w;
-  e.n = qq;
-  MP_DROP(rr);
-  MP_DROP(pp);
-  DA_PUSH(v, e);
-}
-
-static void makefactorbase(mp *x, fbe_v *v)
-{
-  mp *B;
-  unsigned nb;
-  unsigned long max;
-  unsigned long i;
-  size_t j;
-
-  nb = sqrt(mp_bits(x)) * 1.5;
-  B = mp_lsl(MP_NEW, MP_ONE, nb);
-  if (nb < 32)
-    max = 1 << nb;
-  else
-    max = 0xfffffff0;
-
-  addfbe(v, B, 2);
-  for (i = 3; i < max; i++) {
-    for (j = 0; j < DA_LEN(v); j++) {
-      if (i % DA(v)[j].p == 0)
-       goto composite;
-    }
-    addfbe(v, B, i);
-  composite:;
-  }
-
-#ifdef DEBUG
-  for (j = 0; j < DA_LEN(v); j++) {
-    printf("%lu^%u = ", DA(v)[j].p, DA(v)[j].e);
-    mp_writefile(DA(v)[j].n, stdout, 10);
-    putchar('\n');
-  }
-#endif
-}
-
-static void freefactorbase(fbe_v *v)
-{
-  size_t i;
-
-  for (i = 0; i < DA_LEN(v); i++)
-    mp_drop(DA(v)[i].n);
-  DA_DESTROY(v);    
-}
-
-static mp *ecm(mp *x)
-{
-  fbe_v fb = DA_INIT;
-  field *f = field_prime(x);
-  ec_curve *c;
-  ec p;
-  mp *a = MP_NEW, *b = MP_NEW;
-  mp *g = MP_NEW;
-  unsigned i;
-  size_t j;
-
-  makefactorbase(x, &fb);
-  for (;;) {
-    a = mprand_range(a, x, rng, 0);
-    b = mprand_range(b, x, rng, 0);
-    c = ec_primeproj(f, a, b);
-    for (i = 0; i < 1024; i++) {
-      EC_CREATE(&p);
-      p.x = mprand_range(p.x, x, rng, 0);
-      p.y = mprand_range(p.y, x, rng, 0);
-      p.z = MP_ONE;
-      for (j = 0; j < DA_LEN(&fb); j++)
-       ec_imul(c, &p, &p, DA(&fb)[j].n);
-      if (EC_ATINF(&p))
-       continue;
-      mp_gcd(&g, 0, 0, x, p.z);
-      if (!MP_EQ(g, MP_ONE))
-       goto done;
-      EC_DESTROY(&p);
-    }
-    ec_destroycurve(c);
-  }
-
-done:
-  MP_DROP(a); MP_DROP(b);
-  EC_DESTROY(&p);
-  ec_destroycurve(c);
-  F_DESTROY(f);
-  freefactorbase(&fb);
-#ifdef DEBUG
-  MP_PRINT("ecm factor", g);
-#endif
-  return (g);
-}
-
-static void dofactor(mp *x, mp_v *v)
-{
-  mp *f;
-
-  x = smallfactors(x, v); 
-  if (MP_EQ(x, MP_ONE)) return;
-
-#ifdef POLLARDRHO
-  for (;;) {
-    if (pgen_primep(x, rng))
-      goto done;
-    f = pollardrho(x);
-    if (!f) break;
-    dofactor(f, v);
-    mp_div(&x, 0, x, f);
-    MP_DROP(f);
-  }
-#endif
-
-  for (;;) {
-    if (pgen_primep(x, rng))
-      goto done;
-    f = ecm(x);
-    if (!f) break;
-    dofactor(f, v);
-    mp_div(&x, 0, x, f);
-    MP_DROP(f);
-  }
-
-done:
-  DA_PUSH(v, x);
-}
-
-static int cmpmp(const void *a, const void *b)
-{
-  mp *const *aa = a, *const *bb = b;
-  return (mp_cmp(*aa, *bb));
-}
-
-static void factor(mp *x, fact_v *v)
-{
-  fact f;
-  mp_v fv = DA_INIT;
-  mp *e = MP_NEW;
-  size_t i;
-
-  dofactor(x, &fv);
-  qsort(DA(&fv), DA_LEN(&fv), sizeof(mp *), cmpmp);
-  for (i = 0; i < DA_LEN(&fv); ) {
-    f.p = MP_COPY(DA(&fv)[i]);
-    for (f.e = 1, i++;
-        i < DA_LEN(&fv) && MP_EQ(f.p, DA(&fv)[i]);
-        f.e++, i++) ;
-    e = mp_fromuint(e, f.e);
-    f.n = mp_exp(MP_NEW, f.p, e);
-    f.r = MP_NEW;
-    mp_div(&f.r, 0, x, f.p);
-    DA_PUSH(v, f);
-  }
-  mp_drop(e);
-  mpv_free(&fv);
-}
-
-static void freefactors(fact_v *v)
-{
-  size_t i;
-
-  for (i = 0; i < DA_LEN(v); i++) {
-    mp_drop(DA(v)[i].p);
-    mp_drop(DA(v)[i].n);
-    mp_drop(DA(v)[i].r);
-  }
-  DA_DESTROY(v);
-}
-
-
-/*----- Primitive polynomials ---------------------------------------------*/
+#include "factor.h"
 
 static int primitivep(fact_v *f, mp *p)
 {
@@ -402,7 +107,6 @@ static void findprim(unsigned m)
 int main(int argc, char *argv[])
 {
   ego(argv[0]);
-  rng = fibrand_create(0);
   if (argc != 2) {
     fprintf(stderr, "Usage: %s M\n", QUIS);
     exit(1);