+#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);
+}
projects / u / mdw / catacomb / commitdiff