Commit | Line | Data |
---|---|---|
71080a1b MW |
1 | #include <math.h> |
2 | ||
3 | #include "ec.h" | |
4 | #include "fibrand.h" | |
5 | #include "field.h" | |
6 | #include "mpbarrett.h" | |
7 | #include "mpint.h" | |
8 | #include "mprand.h" | |
9 | #include "pgen.h" | |
10 | #include "primetab.h" | |
11 | ||
12 | #include "factor.h" | |
13 | ||
14 | static grand *rng = 0; | |
15 | ||
16 | DA_DECL(mp_v, mp *); | |
17 | ||
18 | typedef struct fbe { | |
19 | unsigned long p; | |
20 | unsigned e; | |
21 | mp *n; | |
22 | } fbe; | |
23 | DA_DECL(fbe_v, fbe); | |
24 | ||
25 | static void mpv_free(mp_v *v) | |
26 | { | |
27 | size_t i; | |
28 | ||
29 | for (i = 0; i < DA_LEN(v); i++) | |
30 | mp_drop(DA(v)[i]); | |
31 | DA_DESTROY(v); | |
32 | } | |
33 | ||
34 | static mp *smallfactors(mp *x, mp_v *v) | |
35 | { | |
36 | mp f; | |
37 | mpw fw; | |
38 | mp *y = MP_NEW, *z = MP_NEW; | |
39 | unsigned i; | |
40 | ||
41 | MP_COPY(x); | |
42 | mp_build(&f, &fw, &fw + 1); | |
43 | #ifdef DEBUG | |
44 | MP_PRINT("target", x); | |
45 | #endif | |
46 | for (i = 0; i < NPRIME; i++) { | |
47 | again: | |
48 | fw = primetab[i]; | |
49 | mp_div(&y, &z, x, &f); | |
50 | if (MP_ZEROP(z)) { | |
51 | #ifdef DEBUG | |
52 | MP_PRINT("factor", &f); | |
53 | #endif | |
54 | MP_DROP(x); | |
55 | x = y; | |
56 | y = MP_NEW; | |
57 | #ifdef DEBUG | |
58 | MP_PRINT("remainder", x); | |
59 | #endif | |
60 | DA_PUSH(v, mp_fromuint(MP_NEW, primetab[i])); | |
61 | goto again; | |
62 | } | |
63 | } | |
64 | mp_drop(y); | |
65 | mp_drop(z); | |
66 | return (x); | |
67 | } | |
68 | ||
69 | #ifdef POLLARDRHO | |
70 | static mp *pollardrho(mp *x) | |
71 | { | |
72 | unsigned i; | |
73 | mp *y = MP_NEW, *z = MP_NEW, *g = MP_NEW, *t = MP_NEW; | |
74 | mpbarrett mb; | |
75 | ||
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); | |
82 | t = mp_sub(t, y, z); | |
83 | mp_gcd(&g, 0, 0, x, t); | |
84 | if (!MP_EQ(g, MP_ONE)) | |
85 | goto done; | |
86 | } | |
87 | MP_DROP(g); | |
88 | g = 0; | |
89 | ||
90 | done: | |
91 | mpbarrett_destroy(&mb); | |
92 | MP_DROP(y); | |
93 | MP_DROP(z); | |
94 | MP_DROP(t); | |
95 | #ifdef DEBUG | |
96 | MP_PRINT("rho factor", g); | |
97 | #endif | |
98 | return (g); | |
99 | } | |
100 | #endif | |
101 | ||
102 | static void addfbe(fbe_v *v, mp *B, unsigned long p) | |
103 | { | |
104 | fbe e; | |
105 | unsigned w; | |
106 | mp *pp = mp_fromulong(MP_NEW, p); | |
107 | mp *qq = MP_NEW, *rr = MP_NEW; | |
108 | mp *t; | |
109 | ||
110 | w = 0; | |
111 | qq = MP_COPY(pp); | |
112 | for (;;) { | |
113 | w++; | |
114 | rr = mp_mul(rr, qq, pp); | |
115 | if (MP_CMP(rr, >=, B)) break; | |
116 | t = rr; rr = qq; qq = t; | |
117 | } | |
118 | e.p = p; | |
119 | e.e = w; | |
120 | e.n = qq; | |
121 | MP_DROP(rr); | |
122 | MP_DROP(pp); | |
123 | DA_PUSH(v, e); | |
124 | } | |
125 | ||
126 | static void makefactorbase(mp *x, fbe_v *v) | |
127 | { | |
128 | mp *B; | |
129 | unsigned nb; | |
130 | unsigned long max; | |
131 | unsigned long i; | |
132 | size_t j; | |
133 | ||
134 | nb = sqrt(mp_bits(x)) * 1.5; | |
135 | B = mp_lsl(MP_NEW, MP_ONE, nb); | |
136 | if (nb < 32) | |
137 | max = 1 << nb; | |
138 | else | |
139 | max = 0xfffffff0; | |
140 | ||
141 | addfbe(v, B, 2); | |
142 | for (i = 3; i < max; i++) { | |
143 | for (j = 0; j < DA_LEN(v); j++) { | |
144 | if (i % DA(v)[j].p == 0) | |
145 | goto composite; | |
146 | } | |
147 | addfbe(v, B, i); | |
148 | composite:; | |
149 | } | |
150 | ||
151 | #ifdef DEBUG | |
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); | |
155 | putchar('\n'); | |
156 | } | |
157 | #endif | |
158 | } | |
159 | ||
160 | static void freefactorbase(fbe_v *v) | |
161 | { | |
162 | size_t i; | |
163 | ||
164 | for (i = 0; i < DA_LEN(v); i++) | |
165 | mp_drop(DA(v)[i].n); | |
45c0fd36 | 166 | DA_DESTROY(v); |
71080a1b MW |
167 | } |
168 | ||
169 | static mp *ecm(mp *x) | |
170 | { | |
171 | fbe_v fb = DA_INIT; | |
172 | field *f = field_prime(x); | |
173 | ec_curve *c; | |
174 | ec p; | |
175 | mp *a = MP_NEW, *b = MP_NEW; | |
176 | mp *g = MP_NEW; | |
177 | unsigned i; | |
178 | size_t j; | |
179 | ||
180 | makefactorbase(x, &fb); | |
181 | for (;;) { | |
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++) { | |
186 | EC_CREATE(&p); | |
187 | p.x = mprand_range(p.x, x, rng, 0); | |
188 | p.y = mprand_range(p.y, x, rng, 0); | |
189 | p.z = MP_ONE; | |
190 | for (j = 0; j < DA_LEN(&fb); j++) | |
191 | ec_imul(c, &p, &p, DA(&fb)[j].n); | |
192 | if (EC_ATINF(&p)) | |
193 | continue; | |
194 | mp_gcd(&g, 0, 0, x, p.z); | |
195 | if (!MP_EQ(g, MP_ONE)) | |
196 | goto done; | |
197 | EC_DESTROY(&p); | |
198 | } | |
199 | ec_destroycurve(c); | |
200 | } | |
201 | ||
202 | done: | |
203 | MP_DROP(a); MP_DROP(b); | |
204 | EC_DESTROY(&p); | |
205 | ec_destroycurve(c); | |
206 | F_DESTROY(f); | |
207 | freefactorbase(&fb); | |
208 | #ifdef DEBUG | |
209 | MP_PRINT("ecm factor", g); | |
210 | #endif | |
211 | return (g); | |
212 | } | |
213 | ||
214 | static void dofactor(mp *x, mp_v *v) | |
215 | { | |
216 | mp *f; | |
217 | ||
45c0fd36 | 218 | x = smallfactors(x, v); |
71080a1b MW |
219 | if (MP_EQ(x, MP_ONE)) return; |
220 | ||
221 | #ifdef POLLARDRHO | |
222 | for (;;) { | |
223 | if (pgen_primep(x, rng)) | |
224 | goto done; | |
225 | f = pollardrho(x); | |
226 | if (!f) break; | |
227 | dofactor(f, v); | |
228 | mp_div(&x, 0, x, f); | |
229 | MP_DROP(f); | |
230 | } | |
231 | #endif | |
232 | ||
233 | for (;;) { | |
234 | if (pgen_primep(x, rng)) | |
235 | goto done; | |
236 | f = ecm(x); | |
237 | if (!f) break; | |
238 | dofactor(f, v); | |
239 | mp_div(&x, 0, x, f); | |
240 | MP_DROP(f); | |
241 | } | |
242 | ||
243 | done: | |
244 | DA_PUSH(v, x); | |
245 | } | |
246 | ||
247 | static int cmpmp(const void *a, const void *b) | |
248 | { | |
249 | mp *const *aa = a, *const *bb = b; | |
250 | return (mp_cmp(*aa, *bb)); | |
251 | } | |
252 | ||
253 | void factor(mp *x, fact_v *v) | |
254 | { | |
255 | fact f; | |
256 | mp_v fv = DA_INIT; | |
257 | mp *e = MP_NEW; | |
258 | size_t i; | |
259 | ||
260 | if (!rng) rng = fibrand_create(0); | |
261 | dofactor(x, &fv); | |
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]); | |
265 | for (f.e = 1, i++; | |
266 | i < DA_LEN(&fv) && MP_EQ(f.p, DA(&fv)[i]); | |
267 | f.e++, i++) ; | |
268 | e = mp_fromuint(e, f.e); | |
269 | f.n = mp_exp(MP_NEW, f.p, e); | |
270 | f.r = MP_NEW; | |
271 | mp_div(&f.r, 0, x, f.p); | |
272 | DA_PUSH(v, f); | |
273 | } | |
274 | mp_drop(e); | |
275 | mpv_free(&fv); | |
276 | } | |
277 | ||
278 | void freefactors(fact_v *v) | |
279 | { | |
280 | size_t i; | |
281 | ||
282 | for (i = 0; i < DA_LEN(v); i++) { | |
283 | mp_drop(DA(v)[i].p); | |
284 | mp_drop(DA(v)[i].n); | |
285 | mp_drop(DA(v)[i].r); | |
286 | } | |
287 | DA_DESTROY(v); | |
288 | } |