f4535c64 |
1 | #include <math.h> |
2 | #include <stdio.h> |
3 | #include <stdlib.h> |
4 | #include <string.h> |
5 | |
6 | #include <mLib/alloc.h> |
7 | #include <mLib/darray.h> |
8 | #include <mLib/dstr.h> |
9 | #include <mLib/quis.h> |
10 | #include <mLib/report.h> |
11 | |
12 | #include "fibrand.h" |
13 | #include "field.h" |
14 | #include "ec.h" |
15 | #include "gf.h" |
16 | #include "gfreduce.h" |
17 | #include "mp.h" |
18 | #include "mpbarrett.h" |
19 | #include "mpint.h" |
20 | #include "mprand.h" |
21 | #include "pgen.h" |
22 | #include "primetab.h" |
23 | |
24 | DA_DECL(mp_v, mp *); |
25 | |
26 | typedef struct fbe { |
27 | unsigned long p; |
28 | unsigned e; |
29 | mp *n; |
30 | } fbe; |
31 | DA_DECL(fbe_v, fbe); |
32 | |
33 | typedef struct fact { |
34 | mp *p; |
35 | unsigned e; |
36 | mp *n; |
37 | mp *r; |
38 | } fact; |
39 | DA_DECL(fact_v, fact); |
40 | |
41 | static grand *rng; |
42 | |
43 | /*----- Factoring ---------------------------------------------------------*/ |
44 | |
45 | static void mpv_free(mp_v *v) |
46 | { |
47 | size_t i; |
48 | |
49 | for (i = 0; i < DA_LEN(v); i++) |
50 | mp_drop(DA(v)[i]); |
51 | DA_DESTROY(v); |
52 | } |
53 | |
54 | static mp *smallfactors(mp *x, mp_v *v) |
55 | { |
56 | mp f; |
57 | mpw fw; |
58 | mp *y = MP_NEW, *z = MP_NEW; |
59 | unsigned i; |
60 | |
61 | MP_COPY(x); |
62 | mp_build(&f, &fw, &fw + 1); |
63 | #ifdef DEBUG |
64 | MP_PRINT("target", x); |
65 | #endif |
66 | for (i = 0; i < NPRIME; i++) { |
67 | again: |
68 | fw = primetab[i]; |
69 | mp_div(&y, &z, x, &f); |
70 | if (MP_ZEROP(z)) { |
71 | #ifdef DEBUG |
72 | MP_PRINT("factor", &f); |
73 | #endif |
74 | MP_DROP(x); |
75 | x = y; |
76 | y = MP_NEW; |
77 | #ifdef DEBUG |
78 | MP_PRINT("remainder", x); |
79 | #endif |
80 | DA_PUSH(v, mp_fromuint(MP_NEW, primetab[i])); |
81 | goto again; |
82 | } |
83 | } |
84 | mp_drop(y); |
85 | mp_drop(z); |
86 | return (x); |
87 | } |
88 | |
89 | #ifdef POLLARDRHO |
90 | static mp *pollardrho(mp *x) |
91 | { |
92 | unsigned i; |
93 | mp *y = MP_NEW, *z = MP_NEW, *g = MP_NEW, *t = MP_NEW; |
94 | mpbarrett mb; |
95 | |
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); |
102 | t = mp_sub(t, y, z); |
103 | mp_gcd(&g, 0, 0, x, t); |
104 | if (!MP_EQ(g, MP_ONE)) |
105 | goto done; |
106 | } |
107 | MP_DROP(g); |
108 | g = 0; |
109 | |
110 | done: |
111 | mpbarrett_destroy(&mb); |
112 | MP_DROP(y); |
113 | MP_DROP(z); |
114 | MP_DROP(t); |
115 | #ifdef DEBUG |
116 | MP_PRINT("rho factor", g); |
117 | #endif |
118 | return (g); |
119 | } |
120 | #endif |
121 | |
122 | static void addfbe(fbe_v *v, mp *B, unsigned long p) |
123 | { |
124 | fbe e; |
125 | unsigned w; |
126 | mp *pp = mp_fromulong(MP_NEW, p); |
127 | mp *qq = MP_NEW, *rr = MP_NEW; |
128 | mp *t; |
129 | |
130 | w = 0; |
131 | qq = MP_COPY(pp); |
132 | for (;;) { |
133 | w++; |
134 | rr = mp_mul(rr, qq, pp); |
135 | if (MP_CMP(rr, >=, B)) break; |
136 | t = rr; rr = qq; qq = t; |
137 | } |
138 | e.p = p; |
139 | e.e = w; |
140 | e.n = qq; |
141 | MP_DROP(rr); |
142 | MP_DROP(pp); |
143 | DA_PUSH(v, e); |
144 | } |
145 | |
146 | static void makefactorbase(mp *x, fbe_v *v) |
147 | { |
148 | mp *B; |
149 | unsigned nb; |
150 | unsigned long max; |
151 | unsigned long i; |
152 | size_t j; |
153 | |
154 | nb = sqrt(mp_bits(x)) * 1.5; |
155 | B = mp_lsl(MP_NEW, MP_ONE, nb); |
156 | if (nb < 32) |
157 | max = 1 << nb; |
158 | else |
159 | max = 0xfffffff0; |
160 | |
161 | addfbe(v, B, 2); |
162 | for (i = 3; i < max; i++) { |
163 | for (j = 0; j < DA_LEN(v); j++) { |
164 | if (i % DA(v)[j].p == 0) |
165 | goto composite; |
166 | } |
167 | addfbe(v, B, i); |
168 | composite:; |
169 | } |
170 | |
171 | #ifdef DEBUG |
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); |
175 | putchar('\n'); |
176 | } |
177 | #endif |
178 | } |
179 | |
180 | static void freefactorbase(fbe_v *v) |
181 | { |
182 | size_t i; |
183 | |
184 | for (i = 0; i < DA_LEN(v); i++) |
185 | mp_drop(DA(v)[i].n); |
186 | DA_DESTROY(v); |
187 | } |
188 | |
189 | static mp *ecm(mp *x) |
190 | { |
191 | fbe_v fb = DA_INIT; |
192 | field *f = field_prime(x); |
193 | ec_curve *c; |
194 | ec p; |
195 | mp *a = MP_NEW, *b = MP_NEW; |
196 | mp *g = MP_NEW; |
197 | unsigned i; |
198 | size_t j; |
199 | |
200 | makefactorbase(x, &fb); |
201 | for (;;) { |
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++) { |
206 | EC_CREATE(&p); |
207 | p.x = mprand_range(p.x, x, rng, 0); |
208 | p.y = mprand_range(p.y, x, rng, 0); |
209 | p.z = MP_ONE; |
210 | for (j = 0; j < DA_LEN(&fb); j++) |
211 | ec_imul(c, &p, &p, DA(&fb)[j].n); |
212 | if (EC_ATINF(&p)) |
213 | continue; |
214 | mp_gcd(&g, 0, 0, x, p.z); |
215 | if (!MP_EQ(g, MP_ONE)) |
216 | goto done; |
217 | EC_DESTROY(&p); |
218 | } |
219 | ec_destroycurve(c); |
220 | } |
221 | |
222 | done: |
223 | MP_DROP(a); MP_DROP(b); |
224 | EC_DESTROY(&p); |
225 | ec_destroycurve(c); |
226 | F_DESTROY(f); |
227 | freefactorbase(&fb); |
228 | #ifdef DEBUG |
229 | MP_PRINT("ecm factor", g); |
230 | #endif |
231 | return (g); |
232 | } |
233 | |
234 | static void dofactor(mp *x, mp_v *v) |
235 | { |
236 | mp *f; |
237 | |
238 | x = smallfactors(x, v); |
239 | if (MP_EQ(x, MP_ONE)) return; |
240 | |
241 | #ifdef POLLARDRHO |
242 | for (;;) { |
243 | if (pgen_primep(x, rng)) |
244 | goto done; |
245 | f = pollardrho(x); |
246 | if (!f) break; |
247 | dofactor(f, v); |
248 | mp_div(&x, 0, x, f); |
249 | MP_DROP(f); |
250 | } |
251 | #endif |
252 | |
253 | for (;;) { |
254 | if (pgen_primep(x, rng)) |
255 | goto done; |
256 | f = ecm(x); |
257 | if (!f) break; |
258 | dofactor(f, v); |
259 | mp_div(&x, 0, x, f); |
260 | MP_DROP(f); |
261 | } |
262 | |
263 | done: |
264 | DA_PUSH(v, x); |
265 | } |
266 | |
267 | static int cmpmp(const void *a, const void *b) |
268 | { |
269 | mp *const *aa = a, *const *bb = b; |
270 | return (mp_cmp(*aa, *bb)); |
271 | } |
272 | |
273 | static void factor(mp *x, fact_v *v) |
274 | { |
275 | fact f; |
276 | mp_v fv = DA_INIT; |
277 | mp *e = MP_NEW; |
278 | size_t i; |
279 | |
280 | dofactor(x, &fv); |
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]); |
284 | for (f.e = 1, i++; |
285 | i < DA_LEN(&fv) && MP_EQ(f.p, DA(&fv)[i]); |
286 | f.e++, i++) ; |
287 | e = mp_fromuint(e, f.e); |
288 | f.n = mp_exp(MP_NEW, f.p, e); |
289 | f.r = MP_NEW; |
290 | mp_div(&f.r, 0, x, f.p); |
291 | DA_PUSH(v, f); |
292 | } |
293 | mp_drop(e); |
294 | mpv_free(&fv); |
295 | } |
296 | |
297 | static void freefactors(fact_v *v) |
298 | { |
299 | size_t i; |
300 | |
301 | for (i = 0; i < DA_LEN(v); i++) { |
302 | mp_drop(DA(v)[i].p); |
303 | mp_drop(DA(v)[i].n); |
304 | mp_drop(DA(v)[i].r); |
305 | } |
306 | DA_DESTROY(v); |
307 | } |
308 | |
309 | |
310 | /*----- Primitive polynomials ---------------------------------------------*/ |
311 | |
312 | static int primitivep(fact_v *f, mp *p) |
313 | { |
314 | gfreduce r; |
315 | unsigned i; |
316 | mp *x = MP_NEW; |
317 | int rc = 1; |
318 | |
319 | if (!gf_irreduciblep(p)) |
320 | return (0); |
321 | #ifdef DEBUG |
322 | MP_PRINTX("p", p); |
323 | #endif |
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); |
327 | #ifdef DEBUG |
328 | MP_PRINT(" r", DA(f)[i].r); |
329 | MP_PRINTX(" x", x); |
330 | #endif |
331 | if (MP_EQ(x, MP_ONE)) { |
332 | rc = 0; |
333 | break; |
334 | } |
335 | } |
336 | gfreduce_destroy(&r); |
337 | MP_DROP(x); |
338 | return (rc); |
339 | } |
340 | |
341 | static int dofip(fact_v *f, unsigned m, mp **p, unsigned n, unsigned x) |
342 | { |
343 | unsigned i; |
344 | |
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)) |
348 | return (1); |
349 | *p = mp_clearbit(*p, *p, i); |
350 | } |
351 | return (0); |
352 | } |
353 | |
354 | static mp *fip(fact_v *f, unsigned m) |
355 | { |
356 | mp *p = MP_ONE; |
357 | unsigned n; |
358 | |
359 | p = mp_setbit(p, p, m); |
360 | n = 0; |
361 | while (!dofip(f, m, &p, n, m)) |
362 | n += 2; |
363 | return (p); |
364 | } |
365 | |
366 | static void findprim(unsigned m) |
367 | { |
368 | fact_v f = DA_INIT; |
369 | mp *q = mp_lsl(MP_NEW, MP_ONE, m); |
370 | mp *p; |
371 | unsigned i; |
372 | |
373 | q = mp_sub(q, q, MP_ONE); |
374 | factor(q, &f); |
375 | |
376 | #ifdef DEBUG |
377 | { |
378 | size_t i; |
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); |
383 | fputs(" (", stdout); |
384 | mp_writefile(DA(&f)[i].r, stdout, 10); |
385 | fputs(")\n", stdout); |
386 | } |
387 | } |
388 | #endif |
389 | |
390 | p = fip(&f, m); |
391 | MP_PRINTX("p", p); |
392 | for (i = m + 1; i--; ) { |
393 | if (mp_testbit(p, i)) |
394 | printf("%u ", i); |
395 | } |
396 | putchar('\n'); |
397 | mp_drop(p); |
398 | mp_drop(q); |
399 | freefactors(&f); |
400 | } |
401 | |
402 | int main(int argc, char *argv[]) |
403 | { |
404 | ego(argv[0]); |
405 | rng = fibrand_create(0); |
406 | if (argc != 2) { |
407 | fprintf(stderr, "Usage: %s M\n", QUIS); |
408 | exit(1); |
409 | } |
410 | findprim(atoi(argv[1])); |
411 | return (0); |
412 | } |