chiark / gitweb /
base/dispatch.c: Recognize `CPUFEAT_ARM_NEON' as requesting ARM64 SIMD.
[catacomb] / utils / factor.c
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);
166   DA_DESTROY(v);
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
218   x = smallfactors(x, v);
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 }