Line data Source code
1 : /* Copyright (C) 2014 The PARI group.
2 :
3 : This file is part of the PARI/GP package.
4 :
5 : PARI/GP is free software; you can redistribute it and/or modify it under the
6 : terms of the GNU General Public License as published by the Free Software
7 : Foundation; either version 2 of the License, or (at your option) any later
8 : version. It is distributed in the hope that it will be useful, but WITHOUT
9 : ANY WARRANTY WHATSOEVER.
10 :
11 : Check the License for details. You should have received a copy of it, along
12 : with the package; see the file 'COPYING'. If not, write to the Free Software
13 : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
14 :
15 : #include "pari.h"
16 : #include "paripriv.h"
17 :
18 : #define DEBUGLEVEL DEBUGLEVEL_isprime
19 :
20 : #define dbg_mode() if (DEBUGLEVEL >= 2)
21 : #define dbg_mode1() if (DEBUGLEVEL >= 3)
22 :
23 : #define ANSI_RED "\x1b[31m"
24 : #define ANSI_GREEN "\x1b[32m"
25 : #define ANSI_YELLOW "\x1b[33m"
26 : #define ANSI_BLUE "\x1b[34m"
27 : #define ANSI_MAGENTA "\x1b[35m"
28 : #define ANSI_CYAN "\x1b[36m"
29 : #define ANSI_WHITE "\x1b[37m"
30 :
31 : #define ANSI_BRIGHT_RED "\x1b[31;1m"
32 : #define ANSI_BRIGHT_GREEN "\x1b[32;1m"
33 : #define ANSI_BRIGHT_YELLOW "\x1b[33;1m"
34 : #define ANSI_BRIGHT_BLUE "\x1b[34;1m"
35 : #define ANSI_BRIGHT_MAGENTA "\x1b[35;1m"
36 : #define ANSI_BRIGHT_CYAN "\x1b[36;1m"
37 : #define ANSI_BRIGHT_WHITE "\x1b[37;1m"
38 :
39 : #define ANSI_RESET "\x1b[0m"
40 :
41 : /* THINGS THAT DON'T BELONG */
42 :
43 : /* Assume that x is a vector such that
44 : f(x) = 1 for x <= k
45 : f(x) = 0 otherwise.
46 : Return k.
47 : */
48 : static long
49 7714 : zv_binsearch0(void *E, long (*f)(void* E, GEN x), GEN x)
50 : {
51 7714 : long lo = 1;
52 7714 : long hi = lg(x)-1;
53 25403 : while (lo < hi) {
54 17689 : long mi = lo + (hi - lo)/2 + 1;
55 17689 : if (f(E,gel(x,mi))) lo = mi;
56 10220 : else hi = mi - 1;
57 : }
58 7714 : if (f(E,gel(x,lo))) return lo;
59 3829 : return 0;
60 : }
61 :
62 : INLINE long
63 0 : time_record(GEN* X0, const char* Xx, long t)
64 : {
65 0 : long i = Xx[0]-'A'+1, j = Xx[1]-'1'+1;
66 0 : umael3(*X0, 1, i, j) += t;
67 0 : umael3(*X0, 2, i, j) ++; return t;
68 : }
69 :
70 : INLINE long
71 0 : timer_record(GEN* X0, const char* Xx, pari_timer* ti)
72 0 : { return time_record(X0, Xx, timer_delay(ti)); }
73 :
74 : INLINE long
75 5361 : FpJ_is_inf(GEN P) { return signe(gel(P, 3)) == 0; }
76 :
77 : /*****************************************************************/
78 :
79 : /* D < 0 fundamental: return the number of units in Q(sqrt(D)) */
80 : INLINE long
81 57654513 : D_get_wD(long D)
82 : {
83 57654513 : if (D == -4) return 4;
84 57652952 : if (D == -3) return 6;
85 57651139 : return 2;
86 : }
87 :
88 : /* Dinfo contains information related to the discirminant
89 : * D: the discriminant (D < 0)
90 : * h: the class number associated to D
91 : * bi: the "best invariant" for computing polclass(D, bi)
92 : * pd: the degree of polclass; equal to h except when bi is a double
93 : * eta-quotient w_p,q with p|D and q|D, where pd = h/2.
94 : * Dfac: the prime factorization of D; we have D = q0 q1* q2* ... qn*
95 : * where q0 = 1, -4, 8, -8, qi* = 1 mod 4 and |qi| is a prime.
96 : * The factorization is a vecsmall listing the indices of the qi* as
97 : * they appear in the primelist (q0 = 1 is omitted) */
98 : INLINE long
99 58191833 : Dinfo_get_D(GEN Dinfo) { return gel(Dinfo, 1)[1]; }
100 : INLINE long
101 30064335 : Dinfo_get_h(GEN Dinfo) { return gel(Dinfo, 1)[2]; }
102 : INLINE long
103 22753661 : Dinfo_get_bi(GEN Dinfo) { return gel(Dinfo, 1)[3]; }
104 : INLINE ulong
105 57799749 : Dinfo_get_pd(GEN Dinfo) { return umael(Dinfo, 1, 4); }
106 : INLINE GEN
107 32506138 : Dinfo_get_Dfac(GEN Dinfo) { return gel(Dinfo, 2); }
108 :
109 : /* primelist and indexlist
110 : *
111 : * primelist begins with 8, -4, -8; subsequent elements are the p^* for
112 : * successive odd primes <= maxsqrt, by increasing absolute value
113 : *
114 : * i | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 |
115 : * ---------+----+----+----+----+----+----+----+----+----+-----+-----+
116 : * i | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 |
117 : * Plist[i] | 8 | -4 | -8 | -3 | 5 | -7 |-11 | 13 | 17 | -19 | -23 |
118 : * p | 3 | 5 | 7 | 9 | 11 | 13 | 15 | 17 | 19 | 21 | 23 |
119 : * ---------+----+----+----+----+----+----+----+----+----+-----+-----+*/
120 :
121 : /* primelist containing primes whose absolute value is at most maxsqrt */
122 : static GEN
123 77 : ecpp_primelist_init(long maxsqrt)
124 : {
125 77 : GEN P = cgetg(maxsqrt, t_VECSMALL);
126 : forprime_t T;
127 77 : long p, j = 4;
128 77 : u_forprime_init(&T, 3, maxsqrt);
129 77 : P[1] = 8; P[2] = -4; P[3] = -8;
130 5215 : while((p = u_forprime_next(&T))) P[j++] = ((p & 3UL) == 1)? p: -p;
131 77 : setlg(P,j); return P;
132 : }
133 :
134 : static GEN
135 2555 : Dfac_to_p(GEN x, GEN P) { pari_APPLY_long(uel(P,x[i])); }
136 : static GEN
137 2555 : Dfac_to_roots(GEN x, GEN P) { pari_APPLY_type(t_VEC, gel(P,x[i])); }
138 :
139 : #if 0
140 : INLINE ulong
141 : ecpp_param_get_maxsqrt(GEN param) { return umael3(param, 1, 1, 1); }
142 : INLINE ulong
143 : ecpp_param_get_maxdisc(GEN param) { return umael3(param, 1, 1, 2); }
144 : #endif
145 : INLINE ulong
146 1120 : ecpp_param_get_maxpcdg(GEN param) { return umael3(param, 1, 1, 3); }
147 : INLINE GEN
148 1120 : ecpp_param_get_uprimelist(GEN param) { return gmael(param, 1, 2); }
149 : INLINE GEN
150 15521 : ecpp_param_get_primelist(GEN param) { return gmael(param, 1, 3); }
151 : INLINE GEN
152 1120 : ecpp_param_get_disclist(GEN param) { return gmael(param, 1, 4); }
153 : INLINE GEN
154 14345 : ecpp_param_get_primeord(GEN param) { return gmael(param, 1, 5); }
155 : INLINE GEN
156 3857 : ecpp_param_get_primorial_vec(GEN param) { return gel(param, 2); }
157 : INLINE GEN
158 4977 : ecpp_param_get_tune(GEN param) { return gel(param, 3); }
159 :
160 : /* Input: x, 20 <= x <= 35
161 : * Output: a vector whose ith entry is the product of all primes below 2^x */
162 : static GEN
163 77 : primorial_vec(ulong x)
164 : {
165 77 : pari_sp av = avma;
166 77 : long i, y = x-19;
167 77 : GEN v = primes_upto_zv(1UL << x), w = cgetg(y+1, t_VEC);
168 : /* ind[i]th prime number is the largest prime <= 2^(20+i) */
169 77 : long ind[] = { 0, 82025L, 155611L, 295947L, 564163L, 1077871L, 2063689L,
170 : 3957809L, 7603553L, 14630843L, 28192750L, 54400028L,
171 : 105097565L, 203280221L, 393615806L, 762939111L, 1480206279L};
172 77 : gel(w,1) = zv_prod_Z(vecslice(v,1,ind[1]));
173 140 : for (i = 2; i <= y; i++)
174 : {
175 63 : pari_sp av2 = avma;
176 63 : GEN q = mulii(gel(w,i-1), zv_prod_Z(vecslice(v,ind[i-1]+1,ind[i])));
177 63 : gel(w,i) = gerepileuptoint(av2, q);
178 : }
179 77 : return gerepileupto(av, w);
180 : }
181 :
182 : /* NDmqg contains
183 : N, as in the theorem in ??ecpp
184 : Dinfo, as discussed in the comment about Dinfo
185 : m, as in the theorem in ??ecpp
186 : q, as in the theorem in ??ecpp
187 : g, a quadratic nonresidue modulo N
188 : sqrt, a list of squareroots
189 : */
190 : INLINE GEN
191 1057 : NDmqg_get_N(GEN x) { return gel(x,1); }
192 : INLINE GEN
193 9786 : NDmqg_get_Dinfo(GEN x) { return gel(x,2); }
194 : INLINE GEN
195 1057 : NDmqg_get_m(GEN x) { return gel(x,3); }
196 : INLINE GEN
197 1057 : NDmqg_get_q(GEN x) { return gel(x,4); }
198 : INLINE GEN
199 1057 : NDmqg_get_g(GEN x) { return gel(x,5); }
200 : INLINE GEN
201 1057 : NDmqg_get_sqrt(GEN x) { return gel(x,6); }
202 :
203 : /* COMPARATOR FUNCTIONS */
204 :
205 : static int
206 28821478 : sort_disclist(void *data, GEN x, GEN y)
207 : {
208 : long d1, h1, g1, o1, pd1, hf1, wD1, d2, h2, g2, o2, pd2, hf2, wD2;
209 : (void)data;
210 28821478 : d1 = Dinfo_get_D(x); wD1 = D_get_wD(d1);
211 28821478 : d2 = Dinfo_get_D(y); wD2 = D_get_wD(d2);
212 : /* higher number of units means more elliptic curves to try */
213 28821478 : if (wD1 != wD2) return wD2 > wD1 ? 1 : -1;
214 : /* lower polclass degree is better: faster computation of roots modulo N */
215 28819756 : pd1 = Dinfo_get_pd(x); /* degree of polclass */
216 28819756 : pd2 = Dinfo_get_pd(y);
217 28819756 : if (pd1 != pd2) return pd1 > pd2 ? 1 : -1;
218 15029322 : g1 = lg(Dinfo_get_Dfac(x))-1; /* genus number */
219 15029322 : h1 = Dinfo_get_h(x); /* class number */
220 15029322 : o1 = h1 >> (g1-1); /* odd class number */
221 15029322 : g2 = lg(Dinfo_get_Dfac(y))-1;
222 15029322 : h2 = Dinfo_get_h(y);
223 15029322 : o2 = h2 >> (g2-1);
224 15029322 : if (o1 != o2) return g1 > g2 ? 1 : -1;
225 : /* lower class number is better: higher probability of succesful cornacchia */
226 11938752 : if (h1 != h2) return h1 > h2 ? 1 : -1;
227 : /* higher height factor is better: polclass would have lower coefficients */
228 11375567 : hf1 = modinv_height_factor(Dinfo_get_bi(x)); /* height factor for best inv. */
229 11375567 : hf2 = modinv_height_factor(Dinfo_get_bi(y));
230 11375567 : if (hf1 != hf2) return hf2 > hf1 ? 1 : -1;
231 : /* "higher" discriminant is better since its absolute value is lower */
232 5149193 : if (d1 != d2) return d2 > d1 ? 1 : -1;
233 0 : return 0;
234 : }
235 :
236 : INLINE long
237 7672 : NDmqg_get_D(GEN x) { return Dinfo_get_D(NDmqg_get_Dinfo(x)); }
238 : static int
239 3836 : sort_NDmq_by_D(void *data, GEN x, GEN y)
240 : {
241 3836 : long d1 = NDmqg_get_D(x);
242 3836 : long d2 = NDmqg_get_D(y);
243 3836 : (void)data; return d2 > d1 ? 1 : -1;
244 : }
245 :
246 : static int
247 57379 : sort_Dmq_by_q(void *data, GEN x, GEN y)
248 57379 : { (void)data; return cmpii(gel(x,3), gel(y,3)); }
249 :
250 : INLINE long
251 0 : Dmq_get_D(GEN Dmq) { return Dinfo_get_D(gel(Dmq,1)); }
252 : INLINE long
253 4634 : Dmq_get_h(GEN Dmq) { return Dinfo_get_h(gel(Dmq,1)); }
254 : static int
255 2317 : sort_Dmq_by_cnum(void *data, GEN x, GEN y)
256 : {
257 2317 : ulong h1 = Dmq_get_h(x);
258 2317 : ulong h2 = Dmq_get_h(y);
259 2317 : if (h1 != h2) return h1 > h2 ? 1 : -1;
260 924 : return sort_Dmq_by_q(data, x, y);
261 : }
262 :
263 : static void
264 392 : allh_r(GEN H, ulong minD, ulong maxD)
265 : {
266 392 : ulong a, A = usqrt(maxD/3), minD2 = (minD-1) >> 1, maxD2 = maxD >> 1;
267 158564 : for (a = 1; a <= A; a++)
268 : {
269 158172 : ulong a2 = a << 1, aa = a*a, aa4 = aa << 2, b, B;
270 : { /* b = 0 */
271 158172 : ulong D = aa << 1;
272 158172 : if (D <= minD2)
273 122990 : D += a2 * ((minD2 - D + a2) / a2);
274 35771512 : for (; D <= maxD2; D += a2) H[D]++;
275 : }
276 158172 : B = aa4 - 1; /* 4a^2 - b^2 */
277 42033320 : for (b = 1; b < a; b++)
278 : {
279 41875148 : ulong D = B >> 1; /* (4a^2 - b^2) / 2 */
280 41875148 : B -= (b << 1) + 1; if (D > maxD2) continue;
281 34463233 : if (D > minD2)
282 : {
283 3095225 : H[D]++; D += a2; /* c = a */
284 : } else
285 31368008 : D += a2 * ((minD2 - D + a2) / a2);
286 2172280887 : for (; D <= maxD2; D += a2) H[D] += 2;
287 : }
288 : { /* b = a */
289 158172 : ulong D = (aa4 - aa) >> 1;
290 158172 : if (D <= minD2)
291 142044 : D += a2 * ((minD2 - D + a2) / a2);
292 36586746 : for (; D <= maxD2; D += a2) H[D] ++;
293 : }
294 : }
295 392 : }
296 :
297 : /* return H s.t if -maxD <= D < 0 is fundamental then H[(-D)>>1] is the
298 : * ordinary class number of Q(sqrt(D)); junk at other entries. */
299 :
300 : static GEN
301 77 : allh(ulong maxD)
302 : {
303 77 : ulong maxD2 = maxD>>1, s = 1UL<<16;
304 77 : GEN H = zero_zv(maxD2);
305 : ulong a;
306 469 : for (a = 0; a < maxD; a += s)
307 392 : allh_r(H, a+1, minss(a+s,maxD));
308 77 : return H;
309 : }
310 :
311 : static GEN
312 1907423 : mkDinfo(GEN c, long D, long h)
313 : {
314 : long bi, pd, p1, p2;
315 1907423 : bi = disc_best_modinv(D);
316 1907423 : pd = (modinv_degree(&p1,&p2,bi) && (-D)%p1==0 && (-D)%p2==0)? h/2: h;
317 1907423 : return mkvec2(mkvecsmall4(D, h, bi, pd), c);
318 : }
319 :
320 : static GEN
321 77 : ecpp_disclist_init(ulong maxdisc, GEN primelist)
322 : {
323 77 : pari_sp av = avma;
324 : long t, ip, u, lp, lmerge;
325 77 : GEN merge, ev, od, Harr = allh(maxdisc); /* table of class numbers*/
326 77 : long lenv = maxdisc/4; /* max length of od/ev */
327 77 : long N = maxomegau(maxdisc) + 1;
328 :
329 : /* od[t] attached to discriminant 1-4*t, ev[t] attached to -4*t */
330 77 : od = cgetg(lenv + 1, t_VEC); /* contains 'long' at first: save memory */
331 77 : ev = cgetg(lenv + 1, t_VEC);
332 : /* first pass: squarefree sieve and restrict to maxsqrt-smooth disc */
333 5661327 : for (t = 1; t <= lenv; t++)
334 : {
335 5661250 : od[t] = 1;
336 5661250 : switch(t&7)
337 : {
338 1415309 : case 0: case 4: ev[t] = 0; break;
339 707679 : case 2: ev[t] =-8; break;
340 707630 : case 6: ev[t] = 8; break;
341 2830632 : default:ev[t] =-4; break;
342 : }
343 : }
344 77 : lp = lg(primelist);
345 5215 : for (ip = 4; ip < lp; ip++) /* skip 8,-4,-8 */
346 : { /* sieve by squares of primes */
347 5138 : long s, q = primelist[ip], p = labs(q);
348 5138 : s = (q == p)? 3: 1;
349 9597910 : for (t = (s*p+1)>>2; t <= lenv; t += p, s += 4)
350 : {
351 9592772 : long c = od[t];
352 9592772 : if (c) { if (s%p == 0) od[t] = 0; else od[t] = c*q; }
353 : }
354 5138 : s = 1;
355 9595320 : for (t = p; t <= lenv; t += p, s++)
356 : {
357 9590182 : long c = ev[t];
358 9590182 : if (c) { if (s%p == 0) ev[t] = 0; else ev[t] = c*q; }
359 : }
360 : }
361 5661327 : for (u = 0, t = 1; t <= lenv; t++)
362 : { /* restrict to maxsqrt-smooth disc */
363 5661250 : if (od[t] != -4*t+1) od[t] = 0; else u++;
364 5661250 : if (ev[t] != -4*t) ev[t] = 0; else u++;
365 : }
366 :
367 : /* second pass: sieve by primes and record factorization */
368 5661327 : for (t = 1; t <= lenv; t++)
369 : {
370 5661250 : if (od[t]) gel(od,t) = vecsmalltrunc_init(N);
371 5661250 : if (ev[t])
372 : {
373 826077 : GEN F = vecsmalltrunc_init(N);
374 : long id;
375 826077 : switch(t&7)
376 : {
377 221921 : case 2: id = 3; break;
378 221893 : case 6: id = 1; break;
379 382263 : default:id = 2; break;
380 : }
381 826077 : vecsmalltrunc_append(F, id);
382 826077 : gel(ev,t) = F;
383 : }
384 : }
385 77 : lp = lg(primelist);
386 5215 : for (ip = 4; ip < lp; ip++) /* skip 8,-4,-8 */
387 : {
388 5138 : long s, q = primelist[ip], p = labs(q);
389 5138 : s = (q == p)? 3: 1;
390 9597910 : for (t = (s*p+1)>>2; t <= lenv; t += p, s += 4)
391 : {
392 9592772 : GEN c = gel(od,t);
393 9592772 : if (c) vecsmalltrunc_append(c, ip);
394 : }
395 5138 : s = 1;
396 9595320 : for (t = p; t <= lenv; t += p, s++)
397 : {
398 9590182 : GEN c = gel(ev,t);
399 9590182 : if (c) vecsmalltrunc_append(c, ip);
400 : }
401 : }
402 : /* merging the two arrays */
403 77 : merge = cgetg(u+1, t_VEC); lmerge = 0;
404 5661327 : for (t = 1; t <= lenv; t++)
405 : {
406 : GEN c;
407 5661250 : c = gel(od,t); if (c) gel(merge, ++lmerge) = mkDinfo(c,1-4*t, Harr[2*t-1]);
408 5661250 : c = gel(ev,t); if (c) gel(merge, ++lmerge) = mkDinfo(c, -4*t, Harr[2*t]);
409 : }
410 77 : setlg(merge, lmerge);
411 77 : gen_sort_inplace(merge, NULL, &sort_disclist, NULL);
412 77 : return gerepilecopy(av, merge);
413 : }
414 :
415 : static GEN
416 77 : ecpp_primeord_init(GEN primelist, GEN disclist)
417 : {
418 77 : long i, k=1, lgdisclist = lg(disclist), lprimelist = lg(primelist);
419 77 : GEN primeord = zero_Flv(lprimelist-1);
420 1907423 : for (i=1;i < lgdisclist; i++)
421 : {
422 1907346 : GEN Dinfo = gel(disclist, i), Dfac = Dinfo_get_Dfac(Dinfo);
423 1907346 : long j, l = lg(Dfac);
424 8550640 : for (j = 1; j < l; j++)
425 : {
426 6643294 : long ip = Dfac[j];
427 6643294 : if (primeord[ip]==0)
428 5369 : primeord[ip]=k++;
429 : }
430 : }
431 77 : return perm_inv(primeord);
432 : }
433 :
434 : /* Input: a vector tune whose components are [maxsqrt,maxpcdg,tdivexp,expiN]
435 : * Output: vector param of precomputations
436 : * let x = be a component of tune then
437 : * param[1][1] = [maxsqrt, maxdisc, maxpcdg]
438 : * param[1][2] = primelist = Vecsmall of primes
439 : * param[1][3] = disclist = vector of Dinfos, sorted by quality
440 : * param[2] = primorial_vec
441 : * param[3] = tune */
442 : static GEN
443 77 : ecpp_param_set(GEN tune, GEN x)
444 : {
445 : pari_timer ti;
446 77 : pari_sp av = avma;
447 77 : ulong maxsqrt = uel(x,1), maxpcdg = uel(x,2), tdivexp = uel(x,3);
448 77 : ulong maxdisc = maxsqrt * maxsqrt;
449 : GEN Plist, U, Utree, Dlist, Olist, primorial, T;
450 77 : T = mkvecsmall3(maxsqrt, maxdisc, maxpcdg);
451 77 : dbg_mode() {
452 0 : err_printf(ANSI_BRIGHT_WHITE "\nECPP: parameters %Ps:\n"
453 : "init time: " ANSI_RESET, x);
454 0 : timer_start(&ti);
455 : }
456 77 : Plist = ecpp_primelist_init(maxsqrt);
457 77 : U = zv_abs(Plist);
458 77 : Utree = mkvec2(U, ZV_producttree(U));
459 77 : dbg_mode() err_printf("%ld",timer_delay(&ti));
460 77 : Dlist = ecpp_disclist_init(maxdisc, Plist);
461 77 : dbg_mode() err_printf(", %ld",timer_delay(&ti));
462 77 : Olist = ecpp_primeord_init(Plist, Dlist);
463 77 : dbg_mode() err_printf(", %ld",timer_delay(&ti));
464 77 : primorial = primorial_vec(tdivexp);
465 77 : dbg_mode() err_printf(", %ld\n",timer_delay(&ti));
466 77 : return gerepilecopy(av, mkvec3(mkvec5(T,Utree,Plist,Dlist,Olist),
467 : primorial, tune));
468 : }
469 :
470 : /* cert contains [N, t, s, a4, [x, y]] as documented in ??ecpp; the following
471 : * information can be obtained:
472 : * D = squarefreepart(t^2-4N)
473 : * m = (N+1-t), q = m/s, a6 = y^3 - x^2 - a4*x, J */
474 : INLINE GEN
475 3023 : cert_get_N(GEN x) { return gel(x,1); }
476 : INLINE GEN
477 1721 : cert_get_t(GEN x) { return gel(x,2); }
478 : INLINE GEN
479 168 : cert_get_D(GEN x)
480 : {
481 168 : GEN N = cert_get_N(x), t = cert_get_t(x);
482 168 : return coredisc(subii(sqri(t), shifti(N, 2)));
483 : }
484 : INLINE GEN
485 882 : cert_get_m(GEN x)
486 : {
487 882 : GEN N = cert_get_N(x), t = cert_get_t(x);
488 881 : return subii(addiu(N, 1), t);
489 : }
490 : INLINE GEN
491 882 : cert_get_s(GEN x) { return gel(x,3); }
492 : INLINE GEN
493 210 : cert_get_q(GEN x) { return diviiexact(cert_get_m(x), cert_get_s(x)); }
494 : INLINE GEN
495 7 : cert_get_qlast(GEN x) { return cert_get_q(veclast(x)); }
496 : INLINE GEN
497 1456 : cert_get_a4(GEN x) { return gel(x,4); }
498 : INLINE GEN
499 1295 : cert_get_P(GEN x) { return gel(x,5); }
500 : INLINE GEN
501 154 : cert_get_x(GEN x) { return gel(cert_get_P(x), 1); }
502 : INLINE GEN
503 469 : cert_get_a6(GEN z)
504 : {
505 469 : GEN N = cert_get_N(z), a = cert_get_a4(z), P = cert_get_P(z);
506 469 : GEN x = gel(P,1), xx = Fp_sqr(x, N);
507 469 : GEN y = gel(P,2), yy = Fp_sqr(y, N);
508 469 : return Fp_sub(yy, Fp_mul(x, Fp_add(xx,a,N), N), N);
509 : }
510 : INLINE GEN
511 168 : cert_get_E(GEN x) { return mkvec2(cert_get_a4(x), cert_get_a6(x)); }
512 : INLINE GEN
513 98 : cert_get_J(GEN x)
514 : {
515 98 : GEN a = cert_get_a4(x), b = cert_get_a6(x), N = cert_get_N(x);
516 98 : return Fp_ellj(a, b, N);
517 : }
518 :
519 : /* "Twist factor". Does not cover J = 0, 1728 */
520 : static GEN
521 49 : cert_get_lambda(GEN z)
522 : {
523 : GEN N, J, a, b, A, B;
524 49 : J = cert_get_J(z);
525 49 : N = cert_get_N(z);
526 49 : a = cert_get_a4(z);
527 49 : b = cert_get_a6(z);
528 49 : Fp_ellj_to_a4a6(J, N, &A, &B);
529 49 : return Fp_div(Fp_mul(a,B,N), Fp_mul(b,A,N), N);
530 : }
531 :
532 : /* Solves for T such that if
533 : [A, B] = [3*J*(1728-J), 2*J*(1728-J)^2]
534 : and you let
535 : L = T^3 + A*T + B, A = A*L^2, B = B*L^3
536 : then
537 : x == TL and y == L^2
538 : */
539 : static GEN
540 49 : cert_get_T(GEN z)
541 : {
542 49 : GEN N = cert_get_N(z), x = cert_get_x(z);
543 49 : GEN l = cert_get_lambda(z); /* l = 1/L */
544 49 : return Fp_mul(x, l, N);
545 : }
546 :
547 : /* database for all invariants, stolen from Hamish's code */
548 : static GEN
549 56 : polmodular_db_init_allinv(void)
550 : {
551 56 : const long DEFAULT_MODPOL_DB_LEN = 32;
552 56 : GEN a, b = cgetg(39+1, t_VEC);
553 : long i;
554 2240 : for (i = 1; i < 40; i++) gel(b,i) = zerovec_block(DEFAULT_MODPOL_DB_LEN);
555 56 : a = zerovec_block(DEFAULT_MODPOL_DB_LEN);
556 56 : return mkvec2(a, b);
557 : }
558 :
559 : /* Given D and a database of modular polynomials, return polclass(D, inv) */
560 : static GEN
561 1057 : D_polclass(long D, long inv, GEN *db)
562 : {
563 1057 : GEN HD, t = mkvec2(gel(*db, 1), inv == 0? gen_0: gmael(*db, 2, inv));
564 1057 : HD = polclass0(D, inv, 0, &t);
565 1057 : gel(*db, 1) = gel(t,1);
566 1057 : if (inv != 0) gmael(*db, 2, inv) = gel(t,2);
567 1057 : return HD;
568 : }
569 :
570 : static GEN
571 2429 : NqE_check(GEN N, GEN q, GEN a, GEN b, GEN s)
572 : {
573 2429 : GEN mP, sP, P = random_FpE(a, b, N);
574 2429 : GEN PJ = mkvec3(gel(P,1), gel(P,2), gen_1);
575 2429 : sP = FpJ_mul(PJ, s, a, N); if (FpJ_is_inf(sP)) return NULL;
576 2429 : mP = FpJ_mul(sP, q, a, N); return FpJ_is_inf(mP)? mkvec2(a, P): NULL;
577 : }
578 :
579 : /* Find an elliptic curve E modulo N and a point P on E which corresponds to
580 : * the ith element of the certificate; uses N, D, m, q, J from previous steps.
581 : * All computations are modulo N unless stated otherwise */
582 :
583 : /* g is a quadratic and cubic nonresidue modulo N */
584 : static GEN
585 203 : j0_find_g(GEN N)
586 : {
587 203 : GEN n = diviuexact(subiu(N, 1), 3);
588 : for(;;)
589 490 : {
590 693 : GEN g = randomi(N);
591 693 : if (kronecker(g, N) == -1 && !equali1(Fp_pow(g, n, N))) return g;
592 : }
593 : }
594 :
595 : static GEN
596 1057 : find_EP(GEN N, long D, GEN q, GEN g, GEN J, GEN s)
597 : {
598 : GEN A0, B0;
599 1057 : Fp_ellj_to_a4a6(J, N, &A0, &B0);
600 : for(;;)
601 0 : { /* expect one iteration: not worth saving the A's and B's */
602 1057 : GEN gg, v, A = A0, B = B0;
603 : long i;
604 1057 : if ((v = NqE_check(N, q, A, B, s))) return v;
605 693 : switch (D_get_wD(D))
606 : {
607 308 : case 2:
608 308 : gg = Fp_sqr(g, N);
609 308 : A = Fp_mul(A, gg, N);
610 308 : B = Fp_mul(Fp_mul(B, gg, N), g, N);
611 308 : if ((v = NqE_check(N, q, A, B, s))) return v;
612 0 : break;
613 140 : case 4:
614 294 : for (i = 1; i < 4; i++)
615 : {
616 294 : A = Fp_mul(A, g, N);
617 294 : if ((v = NqE_check(N, q, A, B, s))) return v;
618 : }
619 0 : break;
620 245 : default: /* 6 */
621 245 : B = Fp_mul(B, g, N);
622 245 : if ((v = NqE_check(N, q, A, B, s))) return v;
623 203 : g = j0_find_g(N);
624 630 : for (i = 1; i < 6; i++)
625 : {
626 630 : B = Fp_mul(B, g, N);
627 630 : if (i == 3) continue;
628 525 : if ((v = NqE_check(N, q, A, B, s))) return v;
629 : }
630 0 : break;
631 : }
632 : }
633 : }
634 :
635 : /* Return the genus field. If real is set only the real part */
636 : static GEN
637 413 : genusfield(GEN Dfac, GEN sq, GEN p, long real)
638 : {
639 413 : long i, j, l = lg(Dfac), dn, n = 0;
640 413 : GEN sn, s = gen_0, R = cgetg(l-1+!real, t_VECSMALL);
641 476 : for (i = 1; i < l; i++)
642 476 : if (Dfac[i] < 0) { n = i; break; }
643 413 : if (n == 0) pari_err_BUG("realgenusfield");
644 413 : dn = Dfac[n]; sn = gel(sq, n);
645 1267 : for (j = i = 1; i < l; i++)
646 854 : if (i != n)
647 : {
648 441 : long d = Dfac[i];
649 441 : GEN si = gel(sq,i);
650 441 : R[j++] = d > 0 ? d : d * dn;
651 441 : s = Fp_add(s, d > 0? si: Fp_mul(si, sn, p), p);
652 : }
653 413 : if (!real)
654 : {
655 161 : R[j++] = dn;
656 161 : s = Fp_add(s,sn, p);
657 : }
658 413 : return mkvec2(R, s);
659 : }
660 :
661 : static GEN
662 35 : realpart(GEN P, GEN R)
663 : {
664 35 : GEN G = galoisinit(rnfequation(P,R), NULL);
665 35 : GEN T = galoisfixedfield(G, gel(gal_get_gen(G),2), 1, -1);
666 35 : setvarn(T,1); return polredbest(T,0);
667 : }
668 :
669 : static GEN
670 35 : dihedralpol(long D, long l)
671 : {
672 35 : GEN P = quadpoly0(utoineg(-D), 1);
673 35 : GEN bnf = Buchall(P, 1, BIGDEFAULTPREC);
674 35 : GEN x = bnrclassfield(bnf, utoipos(l), 0, BIGDEFAULTPREC);
675 70 : pari_APPLY_same(realpart(P, gel(x,i)))
676 : }
677 :
678 : static GEN
679 1057 : FpX_classtower_oneroot(GEN P, GEN Dinfo, GEN Dfac, GEN sq, GEN p)
680 : {
681 1057 : long D = Dinfo_get_D(Dinfo), h = Dinfo_get_h(Dinfo);
682 1057 : if (degpol(P) > 1)
683 : {
684 413 : long real = !modinv_is_double_eta(Dinfo_get_bi(Dinfo));
685 413 : GEN V = genusfield(Dfac, sq, p, real), v = gel(V,1), R = gel(V,2);
686 413 : GEN N = NULL;
687 413 : long i, l = lg(v);
688 1015 : for (i = 1; i < l; i++)
689 : {
690 602 : GEN Q = deg2pol_shallow(gen_1, gen_0, stoi(-v[i]), 1);
691 602 : if (!N) N = Q;
692 : else
693 : {
694 231 : GEN cm = polcompositum0(N,Q,3);
695 231 : N = gel(cm,1); P = QXY_QXQ_evalx(P,gmael(cm,2,2),N);
696 : }
697 602 : P = liftpol_shallow(gmael(nffactor(N,P),1,1));
698 : }
699 413 : if (h%3 == 0 && (!N || degpol(N)*3<degpol(P)))
700 : {
701 35 : GEN v = dihedralpol(D,3);
702 35 : long i, l = lg(v);
703 35 : for (i = 1; i < l; i++)
704 : {
705 35 : GEN Q = gel(v,i);
706 35 : if (!N) N = Q;
707 : else
708 : {
709 0 : GEN cm = polcompositum0(N, Q, 3);
710 0 : N = gel(cm,1); P = QXY_QXQ_evalx(P,gmael(cm,2,2),N);
711 0 : R = Fp_mul(R, gel(cm,4),p);
712 : }
713 35 : R = Fp_add(R, FpX_oneroot_split(Q, p), p);
714 35 : P = liftpol_shallow(gmael(nffactor(N,P),1,1));
715 35 : if (degpol(N)*3>=degpol(P)) break;
716 : }
717 : }
718 413 : if (N)
719 406 : P = FpXY_evalx(Q_primpart(P), R, p);
720 : }
721 1057 : return P;
722 : }
723 :
724 : GEN
725 1057 : ecpp_step2_worker(GEN S, GEN HD, GEN primelist, long dbg)
726 : {
727 1057 : pari_sp av = avma;
728 : pari_timer ti;
729 : GEN J, t, s, EP, rt, res;
730 1057 : GEN N = NDmqg_get_N(S), Dinfo = NDmqg_get_Dinfo(S);
731 1057 : GEN m = NDmqg_get_m(S), q = NDmqg_get_q(S);
732 1057 : GEN g = NDmqg_get_g(S), sq = NDmqg_get_sqrt(S);
733 1057 : long D = Dinfo_get_D(Dinfo), inv = Dinfo_get_bi(Dinfo);
734 1057 : GEN Dfacp = Dfac_to_p(Dinfo_get_Dfac(Dinfo), primelist);
735 1057 : long C2 = 0, C3 = 0, C4 = 0, D1 = 0;
736 1057 : GEN P, c = getrand();
737 1057 : setrand(gen_1); /* for reproducibility */
738 : /* C2: Find a root modulo N of polclass(D,inv) */
739 1057 : if (dbg >= 2) timer_start(&ti);
740 1057 : P = FpX_classtower_oneroot(HD, Dinfo, Dfacp, sq, N);
741 1057 : if (dbg >= 2) C2 = timer_delay(&ti);
742 1057 : rt = FpX_oneroot_split(P, N);
743 1057 : if (dbg >= 2) C3 = timer_delay(&ti);
744 : /* C3: Convert root from previous step into the appropriate j-invariant */
745 1057 : J = Fp_modinv_to_j(rt, inv, N); /* root of polclass(D) */
746 1057 : if (dbg >= 2) C4 = timer_delay(&ti);
747 : /* D1: Find an elliptic curve E with a point P satisfying the theorem */
748 1057 : s = diviiexact(m, q);
749 1057 : EP = find_EP(N, D, q, g, J, s);
750 1056 : if (dbg >= 2) D1 = timer_delay(&ti);
751 :
752 : /* D2: Compute for t and s */
753 1056 : t = subii(addiu(N, 1), m); /* t = N+1-m */
754 1056 : res = mkvec2(mkvec5(N, t, s, gel(EP,1), gel(EP,2)),mkvecsmall4(C2,C3,C4,D1));
755 1056 : setrand(c);
756 1057 : return gerepilecopy(av, res);
757 : }
758 :
759 : /* This uses [N, D, m, q] from step 1 to find the appropriate j-invariants
760 : * to be used in step 3. Step 2 is divided into substeps 2a, 2b, 2c */
761 : static GEN
762 56 : ecpp_step2(GEN step1, GEN *X0, GEN primelist)
763 : {
764 56 : long j, Dprev = 0;
765 : pari_timer ti;
766 56 : GEN perm = gen_indexsort(step1, NULL, &sort_NDmq_by_D);
767 56 : GEN step2 = cgetg(lg(step1), t_VEC);
768 56 : GEN HD = NULL, db = polmodular_db_init_allinv();
769 56 : long ls = lg(step2), pending = 0;
770 : GEN worker;
771 : struct pari_mt pt;
772 56 : GEN Hi = cgetg(ls, t_VEC);
773 1113 : for (j = 1; j < ls; j++)
774 : {
775 1057 : long i = uel(perm, j);
776 1057 : GEN S = gel(step1, i);
777 1057 : GEN Dinfo = NDmqg_get_Dinfo(S);
778 1057 : long D = Dinfo_get_D(Dinfo), inv = Dinfo_get_bi(Dinfo);
779 : /* C1: Find the appropriate class polynomial modulo N */
780 1057 : dbg_mode() timer_start(&ti);
781 1057 : if (D != Dprev) HD = D_polclass(D, inv, &db);
782 1057 : dbg_mode() {
783 0 : GEN N = NDmqg_get_N(S);
784 0 : long tt = timer_record(X0, "C1", &ti);
785 0 : err_printf(ANSI_BRIGHT_GREEN "\n[ %3d | %4ld bits]" ANSI_RESET, i, expi(N));
786 0 : err_printf(ANSI_GREEN " D = %8ld poldeg = %4ld" ANSI_RESET, D, degpol(HD));
787 0 : if (D == Dprev) err_printf(" %6ld", tt);
788 0 : else err_printf(ANSI_BRIGHT_WHITE " %6ld" ANSI_RESET, tt);
789 : }
790 1057 : gel(Hi, i) = HD;
791 : }
792 56 : gunclone_deep(db);
793 56 : worker = snm_closure(is_entry("_ecpp_step2_worker"),
794 : mkvec2(primelist,stoi(DEBUGLEVEL)));
795 56 : mt_queue_start_lim(&pt, worker, ls-1);
796 1182 : for (j = 1; j < ls || pending; j++)
797 : {
798 : long jj;
799 1126 : GEN S = gel(step1, j), HD = gel(Hi, j), done;
800 1126 : mt_queue_submit(&pt, j, j < ls ? mkvec2(S, HD) : NULL);
801 1126 : done = mt_queue_get(&pt, &jj, &pending);
802 1126 : if (!done) continue;
803 1057 : dbg_mode() {
804 0 : GEN T = gel(done,2);
805 0 : GEN N = NDmqg_get_N(gel(step1, jj));
806 0 : err_printf(ANSI_BRIGHT_GREEN "\n[ %3d | %4ld bits]" ANSI_RESET, jj, expi(N));
807 0 : err_printf(" %6ld", time_record(X0, "C2", T[1]));
808 0 : err_printf(" %6ld", time_record(X0, "C3", T[2]));
809 0 : err_printf(" %6ld", time_record(X0, "C4", T[3]));
810 0 : err_printf(" %6ld", time_record(X0, "D1", T[4]));
811 : }
812 1057 : gel(step2, jj) = gel(done,1);
813 : }
814 56 : mt_queue_end(&pt);
815 56 : return step2;
816 : }
817 : /* end of functions for step 2 */
818 :
819 : GEN
820 17997 : ecpp_sqrt_worker(GEN p, GEN g, GEN N)
821 : {
822 17997 : GEN r = Fp_sqrt_i(p, g, N);
823 18003 : return r ? r: gen_0;
824 : }
825 :
826 : static void
827 14345 : ecpp_parsqrt(GEN N, GEN param, GEN kroP, GEN sqrtlist, GEN g, long nb, long *nbsqrt)
828 : {
829 14345 : GEN primelist = ecpp_param_get_primelist(param);
830 14345 : GEN ordinv = ecpp_param_get_primeord(param);
831 14345 : long i, k, l = lg(ordinv), n = *nbsqrt+1;
832 14345 : GEN worker = snm_closure(is_entry("_ecpp_sqrt_worker"), mkvec2(g, N));
833 14345 : GEN W = cgetg(nb+1,t_VEC), V;
834 50301 : for (i=n, k=1; k<=nb && i<l; i++)
835 : {
836 35956 : long pi = ordinv[i];
837 35956 : if (kroP[pi] > 0)
838 18003 : gel(W,k++) = stoi(primelist[pi]);
839 : }
840 14345 : *nbsqrt=i-1;
841 14345 : setlg(W,k);
842 14345 : V = gen_parapply(worker, W);
843 50301 : for (i=n, k=1; k<=nb && i<l; i++)
844 : {
845 35956 : long pi = ordinv[i];
846 35956 : if (kroP[pi] > 0)
847 : {
848 18003 : dbg_mode() err_printf(ANSI_MAGENTA "S" ANSI_RESET);
849 18003 : if (!signe(gel(V,k)))
850 0 : pari_err_BUG("D_find_discsqrt"); /* possible if N composite */
851 18003 : gel(sqrtlist,pi) = gel(V,k++);
852 : }
853 : }
854 14345 : }
855 :
856 : /* start of functions for step 1 */
857 :
858 : /* This finds the square root of D modulo N [given by Dfac]: find the square
859 : * root modulo N of each prime p dividing D and multiply them out */
860 : static GEN
861 40677 : D_find_discsqrt(GEN N, GEN param, GEN Dfac, GEN kroP, GEN sqrtlist, GEN g, long *nbsqrt)
862 : {
863 40677 : GEN s = NULL, sj;
864 40677 : long i, l = lg(Dfac);
865 123914 : for (i = 1; i < l; i++)
866 : {
867 83237 : long j = Dfac[i];
868 97582 : while(!signe(gel(sqrtlist,j)))
869 14345 : ecpp_parsqrt(N, param, kroP, sqrtlist, g, mt_nbthreads(), nbsqrt);
870 83237 : sj = gel(sqrtlist,j);
871 83237 : s = s? Fp_mul(s, sj, N): sj;
872 : }
873 40677 : return s;/* != NULL */
874 : }
875 :
876 : /* Given a solution U, V to U^2 + |D|V^2 = 4N, this find all the possible
877 : cardinalities of the elliptic curves over the finite field F_N
878 : whose endomorphism ring is isomorphic to the maximal order
879 : in the imaginary quadratic number field K = Q(sqrt(D)). ???
880 : */
881 : static GEN
882 10864 : NUV_find_m(GEN Dinfo, GEN N, GEN U, GEN V, long wD)
883 : {
884 10864 : GEN m, Nplus1 = addiu(N, 1), u = U, mvec = cgetg(wD+1, t_VEC);
885 : long i;
886 36358 : for (i = 0; i < wD; i++)
887 : {
888 25494 : if (odd(i)) m = subii(Nplus1, u);
889 : else
890 : {
891 12747 : if (wD == 4 && i==2) u = shifti(V, 1);
892 12250 : else if (wD == 6 && i==2) u = shifti(submuliu(U, V, 3), -1);
893 11557 : else if (wD == 6 && i==4) u = shifti(addmuliu(U, V, 3), -1);
894 12747 : m = addii(Nplus1, u);
895 : }
896 25494 : gel(mvec, i+1) = mkvec2(Dinfo, m);
897 : }
898 10864 : return mvec;
899 : }
900 :
901 : /* Populates Dmbatch with Dmvec's -- whose components are of the form [D,m],
902 : where m is a cardinalities of an elliptic curve with endomorphism ring
903 : isomorphic to the maximal order of imaginary quadratic K = Q(sqrt(D)).
904 : It returns 0 if:
905 : * Any of the p* dividing D is not a quadratic residue mod N
906 : * Cornacchia cannot find a solution U^2 + |D|V^2 = 4N.
907 : Otherwise, it returns the number of cardinalities added to Dmbatch.
908 : Finally, sqrtlist and g help compute the square root modulo N of D.
909 : */
910 : static long
911 538034 : D_collectcards(GEN N, GEN param, GEN* X0, GEN Dinfo, GEN sqrtlist, GEN g, GEN Dmbatch, GEN kroP, long *nbsqrt)
912 : {
913 538034 : long i, l, corn_succ, wD, D = Dinfo_get_D(Dinfo);
914 538034 : GEN U, V, sqrtofDmodN, Dfac = Dinfo_get_Dfac(Dinfo);
915 : pari_timer ti;
916 :
917 : /* A1: Check (p*|N) = 1 for all primes dividing D */
918 538034 : l = lg(Dfac);
919 768992 : for (i = 1; i < l; i++)
920 : {
921 728315 : long j = Dfac[i];
922 728315 : if (kroP[j] < 0) return 0;
923 : }
924 : /* A3: Get square root of D mod N */
925 40677 : dbg_mode() timer_start(&ti);
926 40677 : sqrtofDmodN = D_find_discsqrt(N, param, Dfac, kroP, sqrtlist, g, nbsqrt);
927 40677 : dbg_mode() timer_record(X0, "A3", &ti);
928 : /* A5: Use square root with Cornacchia to solve U^2 + |D|V^2 = 4N */
929 40677 : dbg_mode() timer_start(&ti);
930 40677 : corn_succ = cornacchia2_sqrt( utoi(labs(D)), N, sqrtofDmodN, &U, &V);
931 40677 : dbg_mode() timer_record(X0, "A5", &ti);
932 40677 : if (!corn_succ) {
933 29813 : dbg_mode1() err_printf(ANSI_YELLOW "c" ANSI_RESET);
934 29813 : return 0;
935 : }
936 : /* A6: Collect the w(D) cardinalities of E/(F_N) with CM by D */
937 10864 : dbg_mode() err_printf(ANSI_BRIGHT_YELLOW "D" ANSI_RESET);
938 10864 : wD = D_get_wD(D);
939 10864 : vectrunc_append_batch(Dmbatch, NUV_find_m(Dinfo,N,U,V,wD));
940 10864 : return wD;
941 : }
942 :
943 : /* Compute (S(16N, 2) + S(4096N, 4) + 4)\4 + 1, where S is the nth root
944 : * rounded down. This is at most one more than (N^1/4 + 1)^2. */
945 : static GEN
946 3857 : ecpp_qlo(GEN N)
947 : {
948 3857 : GEN a = sqrtnint(shifti(N, 4), 2);
949 3857 : GEN b = sqrtnint(shifti(N, 12), 4);
950 3857 : return addiu(shifti(addii(a, b), -2), 2);
951 : }
952 :
953 : static long
954 13356 : lessthan_qlo(void* E, GEN Dmq) { return (cmpii(gel(Dmq,3), (GEN)E) < 0); }
955 : static long
956 12047 : gained_bits(void* E, GEN Dmq) { return (expi(gel(Dmq,3)) <= (long)E); }
957 :
958 : /* Input: Dmqvec
959 : * Output: Dmqvec such that q satisfies (N^1/4 + 1)^2 < q < N/2 */
960 : static GEN
961 3857 : Dmqvec_slice(GEN N, GEN Dmqvec)
962 : {
963 : long lo, hi;
964 :
965 3857 : gen_sort_inplace(Dmqvec, NULL, &sort_Dmq_by_q, NULL); /* sort wrt q */
966 3857 : lo = zv_binsearch0((void*)ecpp_qlo(N), &lessthan_qlo, Dmqvec); lo++;
967 3857 : if (lo >= lg(Dmqvec)) return NULL;
968 :
969 3857 : hi = zv_binsearch0((void*)(expi(N)-1), &gained_bits, Dmqvec);
970 3857 : if (hi == 0) return NULL;
971 :
972 3857 : return vecslice(Dmqvec, lo, hi);
973 : }
974 :
975 : /* Given a Dmvec of [D,m]'s, simultaneously remove all prime factors of each
976 : * m less then BOUND_PRIMORIAL. This implements Franke 2004: Proving the
977 : * Primality of Very Large Numbers with fastECPP */
978 : static GEN
979 3857 : Dmvec_batchfactor(GEN Dmvec, GEN primorial)
980 : {
981 3857 : long i, l = lg(Dmvec);
982 3857 : GEN leaf, v = cgetg(l, t_VEC);
983 29351 : for (i = 1; i < l; i++)
984 : { /* cheaply remove powers of 2 */
985 25494 : GEN m = gmael(Dmvec, i, 2);
986 25494 : long v2 = vali(m);
987 25494 : gel(v,i) = v2? shifti(m,-v2): m;
988 : }
989 3857 : leaf = Z_ZV_mod_tree(primorial, v, ZV_producttree(v));
990 : /* Go through each leaf and remove small factors. */
991 29351 : for (i = 1; i < l; i++)
992 : {
993 25494 : GEN q = gel(v,i), Dm = gel(Dmvec,i), L = gel(leaf,i);
994 84315 : while (!equali1(L)) { L = gcdii(q, L); q = diviiexact(q, L); }
995 25494 : gel(v,i) = mkvec3(gel(Dm,1), gel(Dm,2), q);
996 : }
997 3857 : return v;
998 : }
999 :
1000 : /* return good parameters [maxsqrt, maxpcdg, tdivexp] for ecpp(N) */
1001 : static GEN
1002 4977 : tunevec(long expiN, GEN param)
1003 : {
1004 4977 : GEN T = ecpp_param_get_tune(param);
1005 4977 : long i, n = lg(T)-1;
1006 8106 : for (i = 1; i < n; i++) { GEN x = gel(T,i); if (expiN <= x[4]) return x; }
1007 2219 : return gel(T,n);
1008 : }
1009 : static long
1010 4977 : tunevec_tdivbd(long expiN, GEN param) { return uel(tunevec(expiN, param), 3); }
1011 : static long
1012 1120 : tunevec_batchsize(long expiN, GEN param)
1013 : {
1014 1120 : long t, b = 28 - tunevec_tdivbd(expiN, param);
1015 1120 : if (b < 0) return expiN;
1016 1120 : t = expiN >> b; return t < 1? 1: t;
1017 : }
1018 :
1019 : static GEN
1020 3857 : Dmbatch_factor_Dmqvec(GEN N, long expiN, GEN* X0, GEN Dmbatch, GEN param)
1021 : {
1022 3857 : pari_sp av = avma;
1023 : pari_timer ti;
1024 3857 : GEN Dmqvec, primorial_vec = ecpp_param_get_primorial_vec(param);
1025 3857 : GEN primorial = gel(primorial_vec, tunevec_tdivbd(expiN, param)-19);
1026 :
1027 : /* B1: Factor by batch */
1028 3857 : dbg_mode() timer_start(&ti);
1029 3857 : Dmqvec = Dmvec_batchfactor(Dmbatch, primorial);
1030 3857 : dbg_mode() timer_record(X0, "B1", &ti);
1031 :
1032 : /* B2: For each batch, remove cardinalities lower than (N^(1/4)+1)^2
1033 : * and cardinalities in which we didn't win enough bits */
1034 3857 : Dmqvec = Dmqvec_slice(N, Dmqvec);
1035 3857 : if (!Dmqvec) return gc_NULL(av); /* nothing is left */
1036 3857 : return gerepilecopy(av, Dmqvec);
1037 : }
1038 :
1039 : GEN
1040 20928 : ecpp_ispsp_worker(GEN N)
1041 20928 : { return mkvecsmall(BPSW_psp_nosmalldiv(N)); }
1042 :
1043 : static GEN
1044 1057 : mkNDmqg(GEN z, GEN N, GEN Dmq, GEN g, GEN sqrtlist)
1045 : {
1046 1057 : GEN Dinfo = gel(Dmq,1), sq = Dfac_to_roots(Dinfo_get_Dfac(Dinfo),sqrtlist);
1047 1057 : GEN NDmqg = mkcol6(N, Dinfo, gel(Dmq,2), gel(Dmq,3), g, sq);
1048 1057 : return mkvec2(NDmqg, z);
1049 : }
1050 : /* expi(N) > 64. Return [ NDmqg, N_downrun(q) ]; NDmqg is a vector of the form
1051 : * [N,D,m,q,g,sqrt]. For successive D, find a square root of D mod N (g is a
1052 : * quadratic nonresidue), solve U^2 + |D|V^2 = 4N, then find candidates for m.
1053 : * When enough m's, batch-factor them to produce the q's. If one of the q's is
1054 : * pseudoprime, recursive call with N = q. May return gen_0 at toplevel
1055 : * => N has a small prime divisor */
1056 : static GEN
1057 1120 : N_downrun(GEN N, GEN param, GEN worker, GEN *X0, long *depth, long persevere, long stopat)
1058 : {
1059 : pari_timer T, ti;
1060 1120 : long lgdisclist, lprimelist, nbsqrt = 0, t, i, j, expiN = expi(N);
1061 1120 : long persevere_next = 0, FAIL = 0;
1062 : ulong maxpcdg;
1063 : GEN primelist, uprimelist, disclist, sqrtlist, g, Dmbatch, kroP, Np;
1064 :
1065 1120 : dbg_mode() timer_start(&T);
1066 1120 : (*depth)++; /* we're going down the tree. */
1067 1120 : maxpcdg = ecpp_param_get_maxpcdg(param);
1068 1120 : primelist = ecpp_param_get_primelist(param);
1069 1120 : uprimelist = ecpp_param_get_uprimelist(param);
1070 1120 : disclist = ecpp_param_get_disclist(param);
1071 1120 : lgdisclist = lg(disclist);
1072 1120 : t = tunevec_batchsize(expiN, param);
1073 :
1074 : /* Precomputation for taking square roots, g needed for Fp_sqrt_i */
1075 1120 : g = Fp_2gener(N); if (!g) return gen_0; /* Composite if this happens. */
1076 :
1077 : /* Initialize sqrtlist for this N. */
1078 1120 : lprimelist = lg(primelist);
1079 1120 : sqrtlist = zerovec(lprimelist-1);
1080 :
1081 : /* A2: Check (p*|N) = 1 for all p */
1082 1120 : dbg_mode() timer_start(&ti);
1083 : /* kroP[i] will contain (primelist[i] | N) */
1084 1120 : kroP = cgetg(lprimelist,t_VECSMALL);
1085 1120 : switch(mod8(N))
1086 : { /* primelist[1,2,3] = [8,-4,-8]; N is odd */
1087 350 : case 1: kroP[1] = kroP[2] = kroP[3] = 1; break;
1088 259 : case 3: kroP[1] = -1; kroP[2] =-1; kroP[3] = 1; break;
1089 336 : case 5: kroP[1] = -1; kroP[2] = 1; kroP[3] =-1; break;
1090 175 : case 7: kroP[1] = 1; kroP[2] =-1; kroP[3] =-1; break;
1091 : }
1092 1120 : Np = Z_ZV_mod_tree(N, gel(uprimelist,1), gel(uprimelist,2));
1093 133161 : for(i=4; i<lprimelist; i++)
1094 : {
1095 132041 : if (Np[i]==0) return gen_0;
1096 132041 : kroP[i] = kross(Np[i],primelist[i]);
1097 : }
1098 1120 : dbg_mode() timer_record(X0, "A2", &ti);
1099 :
1100 : /* Print the start of this iteration. */
1101 1120 : dbg_mode() {
1102 0 : char o = persevere? '<': '[';
1103 0 : char c = persevere? '>': ']';
1104 0 : err_printf(ANSI_BRIGHT_CYAN "\n%c %3d | %4ld bits%c "
1105 : ANSI_RESET, o, *depth, expiN, c);
1106 : }
1107 : /* Initialize Dmbatch, populated with candidate cardinalities in Phase I
1108 : * (until #Dmbatch >= t); its elements will be factored on Phase II */
1109 1120 : Dmbatch = vectrunc_init(t+7);
1110 :
1111 : /* Number of cardinalities so far; should always be equal to lg(Dmbatch)-1. */
1112 : /* i determines which discriminant we are considering. */
1113 1120 : i = 1;
1114 3913 : while (!FAIL)
1115 : {
1116 : pari_timer F;
1117 3906 : long lgDmqlist, last_i = i, numcard = 0; /* #Dmbatch */
1118 : GEN Dmqlist;
1119 :
1120 : /* Dmbatch begins "empty", but keep the allocated memory. */
1121 3906 : setlg(Dmbatch, 1);
1122 :
1123 : /* PHASE I: Go through the D's and search for candidate m's.
1124 : * Fill up Dmbatch until there are at least t elements. */
1125 538104 : while (i < lgdisclist )
1126 : {
1127 538076 : GEN Dinfo = gel(disclist, i);
1128 : long n;
1129 538076 : if (!persevere && Dinfo_get_pd(Dinfo) > maxpcdg) { FAIL = 1; break; }
1130 538034 : n = D_collectcards(N,param, X0, Dinfo, sqrtlist, g, Dmbatch, kroP, &nbsqrt);
1131 539091 : if (n < 0) return gen_0;
1132 538034 : last_i = i++;
1133 538034 : numcard += n; if (numcard >= t) break;
1134 : }
1135 :
1136 : /* We have exhausted disclist and there are no card. to be factored */
1137 3906 : if (numcard <= 0 && (FAIL || i >= lgdisclist)) break;
1138 :
1139 : /* PHASE II: Find the corresponding q's for the m's found. Use Dmbatch */
1140 : /* Find coresponding q of the candidate m's. */
1141 3857 : dbg_mode() timer_start(&F);
1142 3857 : Dmqlist = Dmbatch_factor_Dmqvec(N, expiN, X0, Dmbatch, param);
1143 3857 : if (Dmqlist == NULL) continue; /* none left => next discriminant. */
1144 :
1145 : /* If we are cheating by adding class numbers, sort by class number */
1146 3857 : if (Dinfo_get_pd(gel(disclist, last_i)) > maxpcdg)
1147 385 : gen_sort_inplace(Dmqlist, NULL, &sort_Dmq_by_cnum, NULL);
1148 :
1149 : /* Check pseudoprimality before going down. */
1150 3857 : lgDmqlist = lg(Dmqlist);
1151 3899 : for (j = 1; j < lgDmqlist; j++)
1152 : {
1153 : GEN z, Dmq, q;
1154 : struct pari_mt pt;
1155 : pari_timer ti2;
1156 3899 : long running, stop = 0, pending = 0;
1157 3899 : mt_queue_start_lim(&pt, worker, lgDmqlist-j);
1158 3899 : dbg_mode() timer_start(&ti2);
1159 27503 : while ((running = (!stop && j < lgDmqlist)) || pending)
1160 : {
1161 : long jj;
1162 : GEN done;
1163 23604 : mt_queue_submit(&pt, j, running ? mkvec(gmael(Dmqlist,j,3)) : NULL);
1164 23604 : done = mt_queue_get(&pt, &jj, &pending);
1165 23604 : if (done)
1166 : {
1167 20930 : if (done[1] == 0)
1168 19773 : { dbg_mode() err_printf(ANSI_WHITE "." ANSI_RESET); }
1169 1157 : else if (!stop || jj < stop)
1170 1104 : stop = jj;
1171 : }
1172 23604 : j++;
1173 : }
1174 3899 : mt_queue_end(&pt);
1175 3899 : dbg_mode() timer_record(X0, "B3", &ti2);
1176 3899 : if (!stop && j >= lgDmqlist) break;
1177 1099 : j = stop; Dmq = gel(Dmqlist,j); q = gel(Dmq,3);
1178 :
1179 1099 : dbg_mode() {
1180 0 : err_printf(ANSI_BRIGHT_RED "\n %5ld bits " ANSI_RESET,
1181 0 : expi(q)-expiN);
1182 0 : err_printf(" D = %ld", Dmq_get_D(Dmq));
1183 0 : err_printf(ANSI_BRIGHT_BLUE " h = %ld" ANSI_RESET, Dmq_get_h(Dmq));
1184 : }
1185 : /* q is pseudoprime */
1186 1099 : if (expi(q) < stopat) z = gen_1; /* q is prime; sentinel */
1187 : else
1188 : {
1189 1043 : z = N_downrun(q, param, worker, X0, depth, persevere_next, stopat);
1190 1043 : if (!z) {
1191 42 : dbg_mode() { char o = persevere? '<': '[', c = persevere? '>': ']';
1192 0 : err_printf(ANSI_CYAN "\n%c %3d | %4ld bits%c "
1193 : ANSI_RESET, o, *depth, expiN, c); }
1194 42 : continue; /* downrun failed */
1195 : }
1196 2058 : if (isintzero(z)) return z;
1197 : }
1198 1057 : return mkNDmqg(z, N, Dmq, g, sqrtlist); /* SUCCESS */
1199 : }
1200 2800 : if (i >= lgdisclist) break; /* discriminants exhausted: FAIL */
1201 2793 : if (Dinfo_get_pd(gel(disclist, last_i)) > maxpcdg)
1202 : {
1203 364 : dbg_mode() err_printf(ANSI_BRIGHT_RED " !" ANSI_RESET);
1204 364 : persevere_next = 1;
1205 : }
1206 : }
1207 : /* FAILED: Out of discriminants. */
1208 63 : umael(*X0, 3, 1)++; /* FAILS++ */
1209 63 : dbg_mode() err_printf(ANSI_BRIGHT_RED " X" ANSI_RESET);
1210 63 : (*depth)--; return NULL;
1211 : }
1212 :
1213 : /* x: the output of the (recursive) downrun function; return a vector
1214 : * whose components are [N, D, m, q, g] */
1215 : static GEN
1216 56 : ecpp_flattencert(GEN x, long depth)
1217 : {
1218 56 : long i, l = depth+1;
1219 56 : GEN ret = cgetg(l, t_VEC);
1220 1113 : for (i = 1; i < l; i++) { gel(ret, i) = gel(x,1); x = gel(x,2); }
1221 56 : return ret;
1222 : }
1223 :
1224 : /* Call the first downrun then unravel the skeleton of the certificate.
1225 : * Assume expi(N) < 64. This returns one of the following:
1226 : * - a vector of the form [N, D, m, q, y]
1227 : * - gen_0 (if N is composite)
1228 : * - NULL (if FAIL) */
1229 : static GEN
1230 77 : ecpp_step1(GEN N, GEN param, GEN* X0, long stopat)
1231 : {
1232 77 : pari_sp av = avma;
1233 77 : long depth = 0;
1234 77 : GEN worker = strtofunction("_ecpp_ispsp_worker");
1235 77 : GEN downrun = N_downrun(N, param, worker, X0, &depth, 1, stopat);
1236 77 : if (downrun == NULL) return gc_NULL(av);
1237 56 : if (typ(downrun) != t_VEC) return gen_0;
1238 56 : return gerepilecopy(av, ecpp_flattencert(downrun, depth));
1239 : }
1240 :
1241 : /* The input is an integer N.
1242 : It uses the precomputation step0 done in ecpp_step0. */
1243 : static GEN
1244 77 : ecpp_param(GEN N, GEN param, long stopat)
1245 : {
1246 :
1247 : GEN step1, answer, Tv, Cv, X0;
1248 77 : pari_sp av = avma;
1249 : long i, j;
1250 :
1251 : /* Check if N is pseudoprime to begin with. */
1252 77 : if (!BPSW_psp(N)) return gen_0;
1253 :
1254 : /* Check if we should even prove it. */
1255 77 : if (expi(N) < stopat) return N;
1256 :
1257 : /* Timers and Counters */
1258 77 : Tv = mkvec4(zero_zv(5), zero_zv(3), zero_zv(4), zero_zv(1));
1259 77 : Cv = mkvec4(zero_zv(5), zero_zv(3), zero_zv(4), zero_zv(1));
1260 77 : X0 = mkvec3(Tv, Cv, zero_zv(1));
1261 :
1262 77 : step1 = ecpp_step1(N, param, &X0, stopat);
1263 77 : if (step1 == NULL) return gc_NULL(av);
1264 56 : if (typ(step1) != t_VEC) { set_avma(av); return gen_0; }
1265 :
1266 56 : answer = ecpp_step2(step1, &X0, ecpp_param_get_primelist(param));
1267 :
1268 56 : dbg_mode() {
1269 0 : for (i = 1; i < lg(Tv); i++)
1270 : {
1271 0 : GEN v = gel(Tv, i);
1272 0 : long l = lg(v);
1273 0 : for (j = 1; j < l; j++)
1274 : {
1275 0 : long t = umael3(X0,1, i,j), n = umael3(X0,2, i,j);
1276 0 : if (!t) continue;
1277 0 : err_printf("\n %c%ld: %16ld %16ld %16.3f", 'A'+i-1,j, t,n, (double)t/n);
1278 : }
1279 : }
1280 0 : err_printf("\n" ANSI_BRIGHT_RED "\nFAILS: %16ld" ANSI_RESET "\n", umael(X0, 3, 1));
1281 : }
1282 56 : return gerepilecopy(av, answer);
1283 : }
1284 :
1285 : static const long ecpp_tune[][4]=
1286 : {
1287 : { 100, 2, 20, 500 },
1288 : { 350, 14, 21, 1000 },
1289 : { 450, 18, 21, 1500 },
1290 : { 750, 22, 22, 2000 },
1291 : { 900, 26, 23, 2500 },
1292 : { 1000, 32, 23, 3000 },
1293 : { 1200, 38, 24, 3500 },
1294 : { 1400, 44, 24, 4000 },
1295 : { 1600, 50, 24, 4500 },
1296 : { 1800, 56, 25, 5000 },
1297 : { 2000, 62, 25, 5500 },
1298 : { 2200, 68, 26, 6000 },
1299 : { 2400, 74, 26, 6500 },
1300 : { 2600, 80, 26, 7000 },
1301 : { 2800, 86, 26, 7500 },
1302 : { 3000, 92, 27, 8000 },
1303 : { 3200, 98, 27, 8500 },
1304 : { 3400, 104, 28, 9000 },
1305 : { 3600, 110, 28, 9500 },
1306 : { 3800, 116, 29, 10000 },
1307 : { 4000, 122, 29, 0 }
1308 : };
1309 :
1310 : /* assume N BPSW-pseudoprime */
1311 : GEN
1312 98 : ecpp0(GEN C, long stopat)
1313 : {
1314 98 : pari_sp av = avma;
1315 : long expiN, i, tunelen;
1316 98 : GEN N = C, tune;
1317 :
1318 : /* Is it a partial certificate? */
1319 98 : if (typ(C) == t_VEC && check_ecppcert(C))
1320 7 : N = cert_get_qlast(C);
1321 98 : if (typ(N) != t_INT)
1322 0 : pari_err_TYPE("ecpp",C);
1323 :
1324 : /* Check if we should even prove it. */
1325 98 : expiN = expi(N);
1326 98 : if (stopat < 64) stopat = 64;
1327 98 : if (expiN < stopat) return C;
1328 :
1329 56 : tunelen = (expiN+499)/500;
1330 56 : tune = cgetg(tunelen+1, t_VEC);
1331 140 : for (i = 1; i <= tunelen && ecpp_tune[i-1][3]; i++)
1332 84 : gel(tune,i) = mkvecsmall4(ecpp_tune[i-1][0], ecpp_tune[i-1][1],
1333 84 : ecpp_tune[i-1][2], ecpp_tune[i-1][3]);
1334 56 : for (; i <= tunelen; i++) gel(tune,i) = mkvecsmall4(100*(i+20),4*i+42,30,500*i);
1335 : for(;;)
1336 21 : {
1337 77 : GEN C2, param, x = gel(tune, tunelen);
1338 77 : param = ecpp_param_set(tune, x);
1339 77 : C2 = ecpp_param(N, param, stopat);
1340 77 : if (C2)
1341 56 : return typ(C)==t_INT? C2: gerepilecopy(av, shallowconcat(C,C2));
1342 21 : x[1] *= 2;
1343 21 : x[2] *= 2;
1344 21 : x[3] = minss(x[3]+1, 30);
1345 : }
1346 : }
1347 :
1348 : GEN
1349 28 : ecpp(GEN N)
1350 28 : { return ecpp0(N, 0); }
1351 :
1352 : long
1353 42 : isprimeECPP(GEN N)
1354 : {
1355 42 : pari_sp av = avma;
1356 42 : if (!BPSW_psp(N)) return 0;
1357 28 : return gc_long(av, !isintzero(ecpp(N)));
1358 : }
1359 :
1360 : /* PARI ECPP Certificate -> Human-readable format */
1361 : static GEN
1362 7 : cert_out(GEN x)
1363 : {
1364 7 : long i, l = lg(x);
1365 : pari_str str;
1366 7 : str_init(&str, 1);
1367 7 : if (typ(x) == t_INT)
1368 : {
1369 0 : str_printf(&str, "%Ps is prime.\nIndeed, ispseudoprime(%Ps) = 1 and %Ps < 2^64.\n", x, x, x);
1370 0 : return strtoGENstr(str.string);
1371 : }
1372 21 : for (i = 1; i < l; i++)
1373 : {
1374 14 : GEN xi = gel(x, i);
1375 14 : str_printf(&str, "[%ld]\n", i);
1376 14 : str_printf(&str, " N = %Ps\n", cert_get_N(xi));
1377 14 : str_printf(&str, " t = %Ps\n", cert_get_t(xi));
1378 14 : str_printf(&str, " s = %Ps\n", cert_get_s(xi));
1379 14 : str_printf(&str, "a = %Ps\n", cert_get_a4(xi));
1380 14 : str_printf(&str, "D = %Ps\n", cert_get_D(xi));
1381 14 : str_printf(&str, "m = %Ps\n", cert_get_m(xi));
1382 14 : str_printf(&str, "q = %Ps\n", cert_get_q(xi));
1383 14 : str_printf(&str, "E = %Ps\n", cert_get_E(xi));
1384 14 : str_printf(&str, "P = %Ps\n", cert_get_P(xi));
1385 : }
1386 7 : return strtoGENstr(str.string);
1387 : }
1388 :
1389 : /* PARI ECPP Certificate -> Magma Certificate
1390 : * [* [*
1391 : * N, |D|, -1, m,
1392 : * [* a, b *],
1393 : * [* x, y, 1 *],
1394 : * [* [* q, 1 *] *]
1395 : * *], ... *] */
1396 : static GEN
1397 14 : magma_out(GEN x)
1398 : {
1399 14 : long i, n = lg(x)-1;
1400 : pari_str ret;
1401 14 : str_init(&ret, 1);
1402 14 : if (typ(x) == t_INT)
1403 : {
1404 0 : str_printf(&ret, "Operation not supported.");
1405 0 : return strtoGENstr(ret.string);
1406 : }
1407 14 : str_printf(&ret, "[* ");
1408 168 : for (i = 1; i<=n; i++)
1409 : {
1410 154 : GEN xi = gel(x,i), E = cert_get_E(xi), P = cert_get_P(xi);
1411 154 : str_printf(&ret, "[* %Ps, %Ps, -1, %Ps, ",
1412 : cert_get_N(xi), negi(cert_get_D(xi)), cert_get_m(xi));
1413 154 : str_printf(&ret, "[* %Ps, %Ps *], ", gel(E,1), gel(E,2));
1414 154 : str_printf(&ret, "[* %Ps, %Ps, 1 *], ", gel(P,1), gel(P,2));
1415 154 : str_printf(&ret, "[* [* %Ps, 1 *] *] *]", cert_get_q(xi));
1416 154 : if (i != n) str_printf(&ret, ", ");
1417 : }
1418 14 : str_printf(&ret, " *]");
1419 14 : return strtoGENstr(ret.string);
1420 : }
1421 :
1422 : /* Prints: label=(sign)hexvalue\n
1423 : * where sign is + or -
1424 : * hexvalue is value in hex, of the form: 0xABC123 */
1425 : static void
1426 721 : primo_printval(pari_str *ret, const char* label, GEN value)
1427 : {
1428 721 : str_printf(ret, label);
1429 721 : if (signe(value) >= 0) str_printf(ret, "=0x%P0X\n", value);
1430 175 : else str_printf(ret, "=-0x%P0X\n", negi(value));
1431 721 : }
1432 :
1433 : /* Input: PARI ECPP Certificate / Output: PRIMO Certificate
1434 : *
1435 : * Let Y^2 = X^3 + aX + b be the elliptic curve over N with the point (x,y)
1436 : * as described in the PARI certificate.
1437 : *
1438 : * If J != 0, 1728, PRIMO asks for [J,T], where T = a/A * B/b * x,
1439 : * A = 3J(1728-J) and B = 2J(1728-J)^2.
1440 : *
1441 : * If J=0 or J=1728, PRIMO asks for [A,B,T]; we let A=a and B=b => T = x */
1442 : static GEN
1443 14 : primo_out(GEN x)
1444 : {
1445 14 : long i, l = (typ(x) == t_INT) ? 1 : lg(x);
1446 : pari_str ret;
1447 14 : str_init(&ret, 1);
1448 14 : str_printf(&ret, "[PRIMO - Primality Certificate]\nFormat=4\n");
1449 14 : str_printf(&ret, "TestCount=%ld\n\n[Comments]\n", l-1);
1450 14 : str_printf(&ret, "Generated by %s\n",paricfg_version);
1451 14 : str_printf(&ret, "https://pari.math.u-bordeaux.fr/\n\n[Candidate]\n");
1452 14 : if (typ(x) == t_INT)
1453 : {
1454 0 : str_printf(&ret, "N=0x%P0X\n", x);
1455 0 : return strtoGENstr(ret.string);
1456 : }
1457 14 : str_printf(&ret, "N=0x%P0X\n\n", cert_get_N(gel(x,1)));
1458 168 : for (i = 1; i < l; i++)
1459 : {
1460 154 : GEN a4, a6, N, Nover2, xi = gel(x,i);
1461 : long Ais0, Bis0;
1462 154 : str_printf(&ret, "[%ld]\n", i);
1463 154 : N = cert_get_N(xi);
1464 154 : Nover2 = shifti(N, -1);
1465 154 : primo_printval(&ret, "S", cert_get_s(xi));
1466 154 : primo_printval(&ret, "W", cert_get_t(xi));
1467 154 : a4 = cert_get_a4(xi); Ais0 = isintzero(a4);
1468 154 : a6 = cert_get_a6(xi); Bis0 = isintzero(a6);
1469 154 : if (!Ais0 && !Bis0) { /* J != 0, 1728 */
1470 49 : primo_printval(&ret, "J", Fp_center_i(cert_get_J(xi), N, Nover2));
1471 49 : primo_printval(&ret, "T", cert_get_T(xi));
1472 : } else {
1473 105 : if (!Ais0) a4 = Fp_center_i(a4, N, Nover2);
1474 105 : if (!Bis0) a6 = Fp_center_i(a6, N, Nover2);
1475 105 : primo_printval(&ret, "A", a4);
1476 105 : primo_printval(&ret, "B", a6);
1477 105 : primo_printval(&ret, "T", cert_get_x(xi));
1478 : }
1479 154 : str_printf(&ret, "\n");
1480 : }
1481 14 : return strtoGENstr(ret.string);
1482 : }
1483 :
1484 : /* return 1 if q > (N^1/4 + 1)^2, 0 otherwise */
1485 : static long
1486 503 : Nq_isvalid(GEN N, GEN q)
1487 : {
1488 503 : GEN c = subii(sqri(subiu(q,1)), N);
1489 503 : if (signe(c) <= 0) return 0;
1490 : /* (q-1)^2 > N; check that (N - (q-1)^2)^2 > 16Nq */
1491 503 : return (cmpii(sqri(c), shifti(mulii(N,q), 4)) > 0);
1492 : }
1493 :
1494 : GEN
1495 511 : primecertisvalid_ecpp_worker(GEN certi)
1496 : {
1497 : GEN N, W, mP, sP, r, m, s, a, P, q;
1498 511 : if (lg(certi)-1 != 5) return gen_0;
1499 :
1500 511 : N = cert_get_N(certi);
1501 510 : if (typ(N) != t_INT || signe(N) <= 0) return gen_0;
1502 510 : switch(umodiu(N,6))
1503 : {
1504 503 : case 1: case 5: break; /* Check if N is not divisible by 2 or 3 */
1505 7 : default: return gen_0;
1506 : }
1507 :
1508 503 : W = cert_get_t(certi);
1509 503 : if (typ(W) != t_INT) return gen_0;
1510 : /* Check if W^2 < 4N (Hasse interval) */
1511 503 : if (!(cmpii(sqri(W), shifti(N,2)) < 0)) return gen_0;
1512 :
1513 504 : s = cert_get_s(certi);
1514 504 : if (typ(s) != t_INT || signe(s) < 0) return gen_0;
1515 :
1516 504 : m = cert_get_m(certi);
1517 504 : q = dvmdii(m, s, &r);
1518 : /* Check m%s == 0 */
1519 503 : if (!isintzero(r)) return gen_0;
1520 : /* Check q > (N^(1/4) + 1)^2 */
1521 503 : if (!Nq_isvalid(N, q)) return gen_0;
1522 :
1523 504 : a = cert_get_a4(certi);
1524 504 : if (typ(a) != t_INT) return gen_0;
1525 :
1526 504 : P = cert_get_P(certi);
1527 504 : if (typ(P) != t_VEC || lg(P) != 3) return gen_0;
1528 504 : P = FpE_to_FpJ(P);
1529 :
1530 : /* Check sP != 0 and third component is coprime to N */
1531 504 : sP = FpJ_mul(P, s, a, N);
1532 504 : if (!isint1(gcdii(gel(sP, 3), N))) return gen_0;
1533 :
1534 : /* Check mP == 0 */
1535 504 : mP = FpJ_mul(sP, q, a, N);
1536 504 : if (!FpJ_is_inf(mP)) return gen_0;
1537 :
1538 504 : return q;
1539 : }
1540 :
1541 : /* return 1 if 'cert' is a valid PARI ECPP certificate */
1542 : static long
1543 35 : ecppisvalid_i(GEN cert)
1544 : {
1545 35 : const long trustbits = 64;/* a pseudoprime < 2^trustbits is prime */
1546 35 : long i, lgcert = lg(cert);
1547 35 : GEN ql, q = gen_0, worker, check;
1548 :
1549 35 : if (typ(cert) == t_INT)
1550 : {
1551 0 : if (expi(cert) >= trustbits) return 0; /* >= 2^trustbits */
1552 0 : return BPSW_psp(cert);
1553 : }
1554 35 : if (typ(cert) != t_VEC || lgcert <= 1) return 0;
1555 35 : ql = veclast(cert); if (lg(ql)-1 != 5) return 0;
1556 35 : ql = cert_get_q(ql);
1557 35 : if (expi(ql) >= trustbits || !BPSW_psp(ql)) return 0;
1558 35 : worker = strtofunction("_primecertisvalid_ecpp_worker");
1559 35 : check = gen_parapply(worker, cert);
1560 490 : for (i = 1; i < lgcert; i++)
1561 : {
1562 462 : GEN certi = gel(cert, i);
1563 462 : GEN qq = gel(check,i), N = cert_get_N(certi);
1564 462 : if (isintzero(qq)) return 0;
1565 455 : if (i > 1 && !equalii(N, q)) return 0; /* N of this entry doesn't match q of previous */
1566 455 : q = qq;
1567 : }
1568 28 : return 1;
1569 : }
1570 : long
1571 35 : ecppisvalid(GEN cert)
1572 35 : { pari_sp av = avma; return gc_long(av, ecppisvalid_i(cert)); }
1573 :
1574 : GEN
1575 35 : ecppexport(GEN cert, long flag)
1576 : {
1577 35 : switch(flag)
1578 : {
1579 7 : case 0: return cert_out(cert);
1580 14 : case 1: return primo_out(cert);
1581 14 : case 2: return magma_out(cert);
1582 : }
1583 0 : pari_err_FLAG("primecertexport");
1584 : return NULL;/*LCOV_EXCL_LINE*/
1585 : }
|
