X-Git-Url: https://git.distorted.org.uk/u/mdw/catacomb/blobdiff_plain/9b2d01865a7dd140c3bf12432b944a9eda1ddbd5..71080a1b7e1bb95a65007aa2263562e43a477d32:/utils/prim.c diff --git a/utils/prim.c b/utils/prim.c index b279304..39b7e6f 100644 --- a/utils/prim.c +++ b/utils/prim.c @@ -1,4 +1,3 @@ -#include #include #include #include @@ -9,305 +8,11 @@ #include #include -#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);