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