-#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.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 ---------------------------------------------*/
+#include "factor.h"
static int primitivep(fact_v *f, mp *p)
{
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);
projects / u / mdw / catacomb / blobdiff