14 static grand *rng = 0;
25 static void mpv_free(mp_v *v)
29 for (i = 0; i < DA_LEN(v); i++)
34 static mp *smallfactors(mp *x, mp_v *v)
38 mp *y = MP_NEW, *z = MP_NEW;
42 mp_build(&f, &fw, &fw + 1);
44 MP_PRINT("target", x);
46 for (i = 0; i < NPRIME; i++) {
49 mp_div(&y, &z, x, &f);
52 MP_PRINT("factor", &f);
58 MP_PRINT("remainder", x);
60 DA_PUSH(v, mp_fromuint(MP_NEW, primetab[i]));
70 static mp *pollardrho(mp *x)
73 mp *y = MP_NEW, *z = MP_NEW, *g = MP_NEW, *t = MP_NEW;
76 mpbarrett_create(&mb, x);
77 y = mprand_range(MP_NEW, x, rng, 0); z = MP_COPY(y);
78 for (i = 0; i < 1024; i++) {
79 y = mp_sqr(y, y); y = mpbarrett_reduce(&mb, y, y);
80 z = mp_sqr(z, z); z = mpbarrett_reduce(&mb, z, z);
81 z = mp_sqr(z, z); z = mpbarrett_reduce(&mb, z, z);
83 mp_gcd(&g, 0, 0, x, t);
84 if (!MP_EQ(g, MP_ONE))
91 mpbarrett_destroy(&mb);
96 MP_PRINT("rho factor", g);
102 static void addfbe(fbe_v *v, mp *B, unsigned long p)
106 mp *pp = mp_fromulong(MP_NEW, p);
107 mp *qq = MP_NEW, *rr = MP_NEW;
114 rr = mp_mul(rr, qq, pp);
115 if (MP_CMP(rr, >=, B)) break;
116 t = rr; rr = qq; qq = t;
126 static void makefactorbase(mp *x, fbe_v *v)
134 nb = sqrt(mp_bits(x)) * 1.5;
135 B = mp_lsl(MP_NEW, MP_ONE, nb);
142 for (i = 3; i < max; i++) {
143 for (j = 0; j < DA_LEN(v); j++) {
144 if (i % DA(v)[j].p == 0)
152 for (j = 0; j < DA_LEN(v); j++) {
153 printf("%lu^%u = ", DA(v)[j].p, DA(v)[j].e);
154 mp_writefile(DA(v)[j].n, stdout, 10);
160 static void freefactorbase(fbe_v *v)
164 for (i = 0; i < DA_LEN(v); i++)
169 static mp *ecm(mp *x)
172 field *f = field_prime(x);
175 mp *a = MP_NEW, *b = MP_NEW;
180 makefactorbase(x, &fb);
182 a = mprand_range(a, x, rng, 0);
183 b = mprand_range(b, x, rng, 0);
184 c = ec_primeproj(f, a, b);
185 for (i = 0; i < 1024; i++) {
187 p.x = mprand_range(p.x, x, rng, 0);
188 p.y = mprand_range(p.y, x, rng, 0);
190 for (j = 0; j < DA_LEN(&fb); j++)
191 ec_imul(c, &p, &p, DA(&fb)[j].n);
194 mp_gcd(&g, 0, 0, x, p.z);
195 if (!MP_EQ(g, MP_ONE))
203 MP_DROP(a); MP_DROP(b);
209 MP_PRINT("ecm factor", g);
214 static void dofactor(mp *x, mp_v *v)
218 x = smallfactors(x, v);
219 if (MP_EQ(x, MP_ONE)) return;
223 if (pgen_primep(x, rng))
234 if (pgen_primep(x, rng))
247 static int cmpmp(const void *a, const void *b)
249 mp *const *aa = a, *const *bb = b;
250 return (mp_cmp(*aa, *bb));
253 void factor(mp *x, fact_v *v)
260 if (!rng) rng = fibrand_create(0);
262 qsort(DA(&fv), DA_LEN(&fv), sizeof(mp *), cmpmp);
263 for (i = 0; i < DA_LEN(&fv); ) {
264 f.p = MP_COPY(DA(&fv)[i]);
266 i < DA_LEN(&fv) && MP_EQ(f.p, DA(&fv)[i]);
268 e = mp_fromuint(e, f.e);
269 f.n = mp_exp(MP_NEW, f.p, e);
271 mp_div(&f.r, 0, x, f.p);
278 void freefactors(fact_v *v)
282 for (i = 0; i < DA_LEN(v); i++) {