X-Git-Url: https://git.distorted.org.uk/u/mdw/catacomb/blobdiff_plain/e40955b56456a6191ebff48b00c667a0674b3bdb..f4535c6454395e6d56ce0091a07b6d4f7d54a47f:/utils/prim.c?ds=sidebyside

diff --git a/utils/prim.c b/utils/prim.c
new file mode 100644
index 0000000..b279304
--- /dev/null
+++ b/utils/prim.c
@@ -0,0 +1,412 @@
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include <mLib/alloc.h>
+#include <mLib/darray.h>
+#include <mLib/dstr.h>
+#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 ---------------------------------------------*/
+
+static int primitivep(fact_v *f, mp *p)
+{
+  gfreduce r;
+  unsigned i;
+  mp *x = MP_NEW;
+  int rc = 1;
+
+  if (!gf_irreduciblep(p))
+    return (0);
+#ifdef DEBUG
+  MP_PRINTX("p", p);
+#endif
+  gfreduce_create(&r, p);
+  for (i = 0; i < DA_LEN(f); i++) {
+    x = gfreduce_exp(&r, x, MP_TWO, DA(f)[i].r);
+#ifdef DEBUG
+    MP_PRINT("  r", DA(f)[i].r);
+    MP_PRINTX("  x", x);
+#endif
+    if (MP_EQ(x, MP_ONE)) {
+      rc = 0;
+      break;
+    }
+  }
+  gfreduce_destroy(&r);
+  MP_DROP(x);
+  return (rc);
+}
+
+static int dofip(fact_v *f, unsigned m, mp **p, unsigned n, unsigned x)
+{
+  unsigned i;
+
+  for (i = n + 1; i < x; i++) {
+    *p = mp_setbit(*p, *p, i);
+    if (n ? dofip(f, m, p, n - 1, i) : primitivep(f, *p))
+      return (1);
+    *p = mp_clearbit(*p, *p, i);
+  }
+  return (0);
+}
+
+static mp *fip(fact_v *f, unsigned m)
+{
+  mp *p = MP_ONE;
+  unsigned n;
+
+  p = mp_setbit(p, p, m);
+  n = 0;
+  while (!dofip(f, m, &p, n, m))
+    n += 2;
+  return (p);  
+}
+
+static void findprim(unsigned m)
+{
+  fact_v f = DA_INIT;
+  mp *q = mp_lsl(MP_NEW, MP_ONE, m);
+  mp *p;
+  unsigned i;
+
+  q = mp_sub(q, q, MP_ONE);
+  factor(q, &f);
+
+#ifdef DEBUG
+  {
+    size_t i;
+    for (i = 0; i < DA_LEN(&f); i++) {
+      mp_writefile(DA(&f)[i].p, stdout, 10);
+      printf("^%u = ", DA(&f)[i].e);
+      mp_writefile(DA(&f)[i].n, stdout, 10);
+      fputs(" (", stdout);
+      mp_writefile(DA(&f)[i].r, stdout, 10);
+      fputs(")\n", stdout);
+    }
+  }
+#endif
+
+  p = fip(&f, m);
+  MP_PRINTX("p", p);
+  for (i = m + 1; i--; ) {
+    if (mp_testbit(p, i))
+      printf("%u ", i);
+  }
+  putchar('\n');
+  mp_drop(p);
+  mp_drop(q);
+  freefactors(&f);
+}
+
+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);
+  }
+  findprim(atoi(argv[1]));
+  return (0);
+}