chiark / gitweb /
symm/t/chacha: Missing test from RFC8439.
[catacomb] / utils / factor.c
CommitLineData
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
14static grand *rng = 0;
15
16DA_DECL(mp_v, mp *);
17
18typedef struct fbe {
19 unsigned long p;
20 unsigned e;
21 mp *n;
22} fbe;
23DA_DECL(fbe_v, fbe);
24
25static 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
34static 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
70static 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
90done:
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
102static 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
126static 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
160static 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
169static 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
202done:
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
214static 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
243done:
244 DA_PUSH(v, x);
245}
246
247static 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
253void 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
278void 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}