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
++) {