6 #include <mLib/alloc.h>
7 #include <mLib/darray.h>
10 #include <mLib/report.h>
18 #include "mpbarrett.h"
39 DA_DECL(fact_v
, fact
);
43 /*----- Factoring ---------------------------------------------------------*/
45 static void mpv_free(mp_v
*v
)
49 for (i
= 0; i
< DA_LEN(v
); i
++)
54 static mp
*smallfactors(mp
*x
, mp_v
*v
)
58 mp
*y
= MP_NEW
, *z
= MP_NEW
;
62 mp_build(&f
, &fw
, &fw
+ 1);
64 MP_PRINT("target", x
);
66 for (i
= 0; i
< NPRIME
; i
++) {
69 mp_div(&y
, &z
, x
, &f
);
72 MP_PRINT("factor", &f
);
78 MP_PRINT("remainder", x
);
80 DA_PUSH(v
, mp_fromuint(MP_NEW
, primetab
[i
]));
90 static mp
*pollardrho(mp
*x
)
93 mp
*y
= MP_NEW
, *z
= MP_NEW
, *g
= MP_NEW
, *t
= MP_NEW
;
96 mpbarrett_create(&mb
, x
);
97 y
= mprand_range(MP_NEW
, x
, rng
, 0); z
= MP_COPY(y
);
98 for (i
= 0; i
< 1024; i
++) {
99 y
= mp_sqr(y
, y
); y
= mpbarrett_reduce(&mb
, y
, y
);
100 z
= mp_sqr(z
, z
); z
= mpbarrett_reduce(&mb
, z
, z
);
101 z
= mp_sqr(z
, z
); z
= mpbarrett_reduce(&mb
, z
, z
);
103 mp_gcd(&g
, 0, 0, x
, t
);
104 if (!MP_EQ(g
, MP_ONE
))
111 mpbarrett_destroy(&mb
);
116 MP_PRINT("rho factor", g
);
122 static void addfbe(fbe_v
*v
, mp
*B
, unsigned long p
)
126 mp
*pp
= mp_fromulong(MP_NEW
, p
);
127 mp
*qq
= MP_NEW
, *rr
= MP_NEW
;
134 rr
= mp_mul(rr
, qq
, pp
);
135 if (MP_CMP(rr
, >=, B
)) break;
136 t
= rr
; rr
= qq
; qq
= t
;
146 static void makefactorbase(mp
*x
, fbe_v
*v
)
154 nb
= sqrt(mp_bits(x
)) * 1.5;
155 B
= mp_lsl(MP_NEW
, MP_ONE
, nb
);
162 for (i
= 3; i
< max
; i
++) {
163 for (j
= 0; j
< DA_LEN(v
); j
++) {
164 if (i
% DA(v
)[j
].p
== 0)
172 for (j
= 0; j
< DA_LEN(v
); j
++) {
173 printf("%lu^%u = ", DA(v
)[j
].p
, DA(v
)[j
].e
);
174 mp_writefile(DA(v
)[j
].n
, stdout
, 10);
180 static void freefactorbase(fbe_v
*v
)
184 for (i
= 0; i
< DA_LEN(v
); i
++)
189 static mp
*ecm(mp
*x
)
192 field
*f
= field_prime(x
);
195 mp
*a
= MP_NEW
, *b
= MP_NEW
;
200 makefactorbase(x
, &fb
);
202 a
= mprand_range(a
, x
, rng
, 0);
203 b
= mprand_range(b
, x
, rng
, 0);
204 c
= ec_primeproj(f
, a
, b
);
205 for (i
= 0; i
< 1024; i
++) {
207 p
.x
= mprand_range(p
.x
, x
, rng
, 0);
208 p
.y
= mprand_range(p
.y
, x
, rng
, 0);
210 for (j
= 0; j
< DA_LEN(&fb
); j
++)
211 ec_imul(c
, &p
, &p
, DA(&fb
)[j
].n
);
214 mp_gcd(&g
, 0, 0, x
, p
.z
);
215 if (!MP_EQ(g
, MP_ONE
))
223 MP_DROP(a
); MP_DROP(b
);
229 MP_PRINT("ecm factor", g
);
234 static void dofactor(mp
*x
, mp_v
*v
)
238 x
= smallfactors(x
, v
);
239 if (MP_EQ(x
, MP_ONE
)) return;
243 if (pgen_primep(x
, rng
))
254 if (pgen_primep(x
, rng
))
267 static int cmpmp(const void *a
, const void *b
)
269 mp
*const *aa
= a
, *const *bb
= b
;
270 return (mp_cmp(*aa
, *bb
));
273 static void factor(mp
*x
, fact_v
*v
)
281 qsort(DA(&fv
), DA_LEN(&fv
), sizeof(mp
*), cmpmp
);
282 for (i
= 0; i
< DA_LEN(&fv
); ) {
283 f
.p
= MP_COPY(DA(&fv
)[i
]);
285 i
< DA_LEN(&fv
) && MP_EQ(f
.p
, DA(&fv
)[i
]);
287 e
= mp_fromuint(e
, f
.e
);
288 f
.n
= mp_exp(MP_NEW
, f
.p
, e
);
290 mp_div(&f
.r
, 0, x
, f
.p
);
297 static void freefactors(fact_v
*v
)
301 for (i
= 0; i
< DA_LEN(v
); i
++) {
310 /*----- Primitive polynomials ---------------------------------------------*/
312 static int primitivep(fact_v
*f
, mp
*p
)
319 if (!gf_irreduciblep(p
))
324 gfreduce_create(&r
, p
);
325 for (i
= 0; i
< DA_LEN(f
); i
++) {
326 x
= gfreduce_exp(&r
, x
, MP_TWO
, DA(f
)[i
].r
);
328 MP_PRINT(" r", DA(f
)[i
].r
);
331 if (MP_EQ(x
, MP_ONE
)) {
336 gfreduce_destroy(&r
);
341 static int dofip(fact_v
*f
, unsigned m
, mp
**p
, unsigned n
, unsigned x
)
345 for (i
= n
+ 1; i
< x
; i
++) {
346 *p
= mp_setbit(*p
, *p
, i
);
347 if (n ?
dofip(f
, m
, p
, n
- 1, i
) : primitivep(f
, *p
))
349 *p
= mp_clearbit(*p
, *p
, i
);
354 static mp
*fip(fact_v
*f
, unsigned m
)
359 p
= mp_setbit(p
, p
, m
);
361 while (!dofip(f
, m
, &p
, n
, m
))
366 static void findprim(unsigned m
)
369 mp
*q
= mp_lsl(MP_NEW
, MP_ONE
, m
);
373 q
= mp_sub(q
, q
, MP_ONE
);
379 for (i
= 0; i
< DA_LEN(&f
); i
++) {
380 mp_writefile(DA(&f
)[i
].p
, stdout
, 10);
381 printf("^%u = ", DA(&f
)[i
].e
);
382 mp_writefile(DA(&f
)[i
].n
, stdout
, 10);
384 mp_writefile(DA(&f
)[i
].r
, stdout
, 10);
385 fputs(")\n", stdout
);
392 for (i
= m
+ 1; i
--; ) {
393 if (mp_testbit(p
, i
))
402 int main(int argc
, char *argv
[])
405 rng
= fibrand_create(0);
407 fprintf(stderr
, "Usage: %s M\n", QUIS
);
410 findprim(atoi(argv
[1]));