From: Mark Wooding Date: Fri, 17 Feb 2006 12:01:17 +0000 (+0000) Subject: utils: Make very bad ECM factoring program. X-Git-Url: https://git.distorted.org.uk/u/mdw/catacomb/commitdiff_plain/71080a1b7e1bb95a65007aa2263562e43a477d32 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 index 0000000..908ce26 --- /dev/null +++ b/utils/fact.c @@ -0,0 +1,28 @@ +#include +#include +#include + +#include + +#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 index 0000000..5f29c1a --- /dev/null +++ b/utils/factor.c @@ -0,0 +1,289 @@ +#include + +#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 index 0000000..7181505 --- /dev/null +++ b/utils/factor.h @@ -0,0 +1,19 @@ +#ifndef FACTOR_H +#define FACTOR_H + +#include + +#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 --- 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);