Line data Source code
1 : /* Copyright (C) 2000 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 : /*******************************************************************/
16 : /* THE APRCL PRIMALITY TEST */
17 : /*******************************************************************/
18 : #include "pari.h"
19 : #include "paripriv.h"
20 :
21 : #define DEBUGLEVEL DEBUGLEVEL_isprime
22 :
23 : typedef struct Red {
24 : /* global data */
25 : GEN N; /* prime we are certifying */
26 : GEN N2; /* floor(N/2) */
27 : /* global data for flexible window */
28 : long k, lv;
29 : ulong mask;
30 : /* reduction data */
31 : long n;
32 : GEN C; /* polcyclo(n) */
33 : GEN (*red)(GEN x, struct Red*);
34 : } Red;
35 :
36 : #define cache_aall(C) (gel((C),1))
37 : #define cache_tall(C) (gel((C),2))
38 : #define cache_cyc(C) (gel((C),3))
39 : #define cache_E(C) (gel((C),4))
40 : #define cache_eta(C) (gel((C),5))
41 : /* the last 4 only when N = 1 mod p^k */
42 : #define cache_mat(C) (gel((C),6))
43 : #define cache_matinv(C) (gel((C),7))
44 : #define cache_a(C) (gel((C),8)) /* a has order p^k in (Z/NZ)^* */
45 : #define cache_pk(C) (gel((C),9)) /* p^k */
46 :
47 : static GEN
48 1136272 : makepoldeg1(GEN c, GEN d)
49 : {
50 : GEN z;
51 1136272 : if (signe(c)) {
52 1136272 : z = cgetg(4,t_POL); z[1] = evalsigne(1);
53 1136272 : gel(z,2) = d; gel(z,3) = c;
54 0 : } else if (signe(d)) {
55 0 : z = cgetg(3,t_POL); z[1] = evalsigne(1);
56 0 : gel(z,2) = d;
57 : } else {
58 0 : z = cgetg(2,t_POL); z[1] = evalsigne(0);
59 : }
60 1136272 : return z;
61 : }
62 :
63 : /* T mod polcyclo(p), assume deg(T) < 2p */
64 : static GEN
65 1737101 : red_cyclop(GEN T, long p)
66 : {
67 : long i, d;
68 : GEN y, z;
69 :
70 1737101 : d = degpol(T) - p; /* < p */
71 1737101 : if (d <= -2) return T;
72 :
73 : /* reduce mod (x^p - 1) */
74 1737101 : y = ZX_mod_Xnm1(T, p);
75 1737101 : z = y+2;
76 :
77 : /* reduce mod x^(p-1) + ... + 1 */
78 1737101 : d = p-1;
79 1737101 : if (degpol(y) == d)
80 10943077 : for (i=0; i<d; i++) gel(z,i) = subii(gel(z,i), gel(z,d));
81 1737101 : return normalizepol_lg(y, d+2);
82 : }
83 :
84 : /* x t_VECSMALL, as red_cyclo2n_ip */
85 : static GEN
86 7332 : u_red_cyclo2n_ip(GEN x, long n)
87 : {
88 7332 : long i, pow2 = 1L<<(n-1);
89 : GEN z;
90 :
91 37240 : for (i = lg(x)-1; i>pow2; i--) x[i-pow2] -= x[i];
92 7629 : for (; i>0; i--)
93 7629 : if (x[i]) break;
94 7332 : i += 2;
95 7332 : z = cgetg(i, t_POL); z[1] = evalsigne(1);
96 36943 : for (i--; i>=2; i--) gel(z,i) = stoi(x[i-1]);
97 7332 : return z;
98 : }
99 : /* x t_POL, n > 0. Return x mod polcyclo(2^n) = (x^(2^(n-1)) + 1). IN PLACE */
100 : static GEN
101 561014 : red_cyclo2n_ip(GEN x, long n)
102 : {
103 561014 : long i, pow2 = 1L<<(n-1);
104 5636627 : for (i = lg(x)-1; i>pow2+1; i--)
105 5075613 : if (signe(gel(x,i))) gel(x,i-pow2) = subii(gel(x,i-pow2), gel(x,i));
106 561014 : return normalizepol_lg(x, i+1);
107 : }
108 : static GEN
109 559737 : red_cyclo2n(GEN x, long n) { return red_cyclo2n_ip(leafcopy(x), n); }
110 :
111 : /* special case R->C = polcyclo(2^n) */
112 : static GEN
113 559737 : _red_cyclo2n(GEN x, Red *R)
114 559737 : { return centermod_i(red_cyclo2n(x, R->n), R->N, R->N2); }
115 : /* special case R->C = polcyclo(p) */
116 : static GEN
117 1737101 : _red_cyclop(GEN x, Red *R)
118 1737101 : { return centermod_i(red_cyclop(x, R->n), R->N, R->N2); }
119 : static GEN
120 802346 : _red(GEN x, Red *R)
121 802346 : { return centermod_i(ZX_rem(x, R->C), R->N, R->N2); }
122 : static GEN
123 15659566 : modZ(GEN x, Red *R) { return centermodii(x, R->N, R->N2); }
124 : static GEN
125 7246460 : sqrmodZ(GEN x, Red *R) { return modZ(sqri(x), R); }
126 : static GEN
127 2240142 : sqrmod(GEN x, Red *R) { return R->red(ZX_sqr(x), R); }
128 : static GEN
129 2046473 : _mul(GEN x, GEN y, Red *R)
130 1211436 : { return typ(x) == t_INT? modZ(mulii(x,y), R)
131 3257909 : : R->red(ZX_mul(x,y), R); }
132 :
133 : static GEN
134 0 : sqrconst(GEN pol, Red *R)
135 : {
136 0 : GEN z = cgetg(3,t_POL);
137 0 : gel(z,2) = modZ(sqri(gel(pol,2)), R);
138 0 : z[1] = pol[1]; return z;
139 : }
140 :
141 : /* pol^2 mod (x^2+x+1, N) */
142 : static GEN
143 662718 : sqrmod3(GEN pol, Red *R)
144 : {
145 : GEN a,b,bma,A,B;
146 662718 : long lv = lg(pol);
147 :
148 662718 : if (lv==2) return pol;
149 662718 : if (lv==3) return sqrconst(pol, R);
150 662718 : a = gel(pol,3);
151 662718 : b = gel(pol,2); bma = subii(b,a);
152 662718 : A = modZ(mulii(a,addii(b,bma)), R);
153 662718 : B = modZ(mulii(bma,addii(a,b)), R); return makepoldeg1(A,B);
154 : }
155 :
156 : /* pol^2 mod (x^2+1,N) */
157 : static GEN
158 473554 : sqrmod4(GEN pol, Red *R)
159 : {
160 : GEN a,b,A,B;
161 473554 : long lv = lg(pol);
162 :
163 473554 : if (lv==2) return pol;
164 473554 : if (lv==3) return sqrconst(pol, R);
165 473554 : a = gel(pol,3);
166 473554 : b = gel(pol,2);
167 473554 : A = modZ(mulii(a, shifti(b,1)), R);
168 473554 : B = modZ(mulii(subii(b,a),addii(b,a)), R); return makepoldeg1(A,B);
169 : }
170 :
171 : /* pol^2 mod (polcyclo(5),N) */
172 : static GEN
173 1223199 : sqrmod5(GEN pol, Red *R)
174 : {
175 : GEN c2,b,c,d,A,B,C,D;
176 1223199 : long lv = lg(pol);
177 :
178 1223199 : if (lv==2) return pol;
179 1223199 : if (lv==3) return sqrconst(pol, R);
180 1223199 : c = gel(pol,3); c2 = shifti(c,1);
181 1223199 : d = gel(pol,2);
182 1223199 : if (lv==4)
183 : {
184 0 : A = modZ(sqri(d), R);
185 0 : B = modZ(mulii(c2, d), R);
186 0 : C = modZ(sqri(c), R); return mkpoln(3,A,B,C);
187 : }
188 1223199 : b = gel(pol,4);
189 1223199 : if (lv==5)
190 : {
191 70 : A = mulii(b, subii(c2,b));
192 70 : B = addii(sqri(c), mulii(b, subii(shifti(d,1),b)));
193 70 : C = subii(mulii(c2,d), sqri(b));
194 70 : D = mulii(subii(d,b), addii(d,b));
195 : }
196 : else
197 : { /* lv == 6 */
198 1223129 : GEN a = gel(pol,5), a2 = shifti(a,1);
199 : /* 2a(d - c) + b(2c - b) */
200 1223129 : A = addii(mulii(a2, subii(d,c)), mulii(b, subii(c2,b)));
201 : /* c(c - 2a) + b(2d - b) */
202 1223129 : B = addii(mulii(c, subii(c,a2)), mulii(b, subii(shifti(d,1),b)));
203 : /* (a-b)(a+b) + 2c(d - a) */
204 1223129 : C = addii(mulii(subii(a,b),addii(a,b)), mulii(c2,subii(d,a)));
205 : /* 2a(b - c) + (d-b)(d+b) */
206 1223129 : D = addii(mulii(a2, subii(b,c)), mulii(subii(d,b), addii(d,b)));
207 : }
208 1223199 : A = modZ(A, R);
209 1223199 : B = modZ(B, R);
210 1223199 : C = modZ(C, R);
211 1223199 : D = modZ(D, R); return mkpoln(4,A,B,C,D);
212 : }
213 :
214 : /* jac^floor(N/pk) mod (N, polcyclo(pk)), flexible window */
215 : static GEN
216 43429 : _powpolmod(GEN C, GEN jac, Red *R, GEN (*_sqr)(GEN, Red *))
217 : {
218 43429 : const GEN taba = cache_aall(C);
219 43429 : const GEN tabt = cache_tall(C);
220 43429 : const long efin = lg(taba)-1, lv = R->lv;
221 43429 : GEN L, res = jac, pol2 = _sqr(res, R);
222 : long f;
223 43429 : pari_sp av0 = avma, av;
224 :
225 43429 : L = cgetg(lv+1, t_VEC); gel(L,1) = res;
226 488608 : for (f=2; f<=lv; f++) gel(L,f) = _mul(gel(L,f-1), pol2, R);
227 43429 : av = avma;
228 1688152 : for (f = efin; f >= 1; f--)
229 : {
230 1644723 : GEN t = gel(L, taba[f]);
231 1644723 : long tf = tabt[f];
232 1644723 : res = (f==efin)? t: _mul(t, res, R);
233 13442333 : while (tf--) {
234 11797610 : res = _sqr(res, R);
235 11797610 : if (gc_needed(av,1)) {
236 1150 : res = gc_GEN(av, res);
237 1150 : if(DEBUGMEM>1) pari_warn(warnmem,"powpolmod: f = %ld",f);
238 : }
239 : }
240 : }
241 43429 : return gc_GEN(av0, res);
242 : }
243 :
244 : static GEN
245 8782 : _powpolmodsimple(GEN C, Red *R, GEN jac)
246 : {
247 8782 : pari_sp av = avma;
248 8782 : GEN w = ZM_ZX_mul(cache_mat(C), jac);
249 8782 : long j, ph = lg(w);
250 :
251 8782 : R->red = &modZ;
252 34710 : for (j=1; j<ph; j++) gel(w,j) = _powpolmod(C, modZ(gel(w,j), R), R, &sqrmodZ);
253 8782 : w = centermod_i( ZM_ZC_mul(cache_matinv(C), w), R->N, R->N2);
254 8782 : w = gc_upto(av, w);
255 8782 : return RgV_to_RgX(w, 0);
256 : }
257 :
258 : static GEN
259 26283 : powpolmod(GEN C, Red *R, long p, long k, GEN jac)
260 : {
261 : GEN (*_sqr)(GEN, Red *);
262 :
263 26283 : if (!isintzero(cache_mat(C))) return _powpolmodsimple(C, R, jac);
264 17501 : if (p == 2) /* p = 2 */
265 : {
266 3283 : if (k == 2) _sqr = &sqrmod4;
267 783 : else _sqr = &sqrmod;
268 3283 : R->n = k;
269 3283 : R->red = &_red_cyclo2n;
270 14218 : } else if (k == 1)
271 : {
272 13203 : if (p == 3) _sqr = &sqrmod3;
273 10005 : else if (p == 5) _sqr = &sqrmod5;
274 4826 : else _sqr = &sqrmod;
275 13203 : R->n = p;
276 13203 : R->red = &_red_cyclop;
277 : } else {
278 1015 : R->red = &_red; _sqr = &sqrmod;
279 : }
280 17501 : return _powpolmod(C, jac, R, _sqr);
281 : }
282 :
283 : /* Return e(t) = \prod_{p-1 | t} p^{1+v_p(t)}}
284 : * faet contains the odd prime divisors of e(t) */
285 : static GEN
286 951 : compute_e(ulong t, GEN *faet)
287 : {
288 951 : GEN L, P, D = divisorsu(t);
289 951 : long l = lg(D);
290 : ulong k;
291 :
292 951 : P = vecsmalltrunc_init(l);
293 951 : L = vecsmalltrunc_init(l);
294 26976 : for (k=l-1; k>1; k--) /* k != 1: avoid d = 1 */
295 : {
296 26025 : ulong d = D[k];
297 26025 : if (uisprime(++d))
298 : { /* we want q = 1 (mod p) prime, not too large */
299 : #ifdef LONG_IS_64BIT
300 11416 : if (d > 5000000000UL) return gen_0;
301 : #endif
302 13008 : vecsmalltrunc_append(P, d);
303 13008 : vecsmalltrunc_append(L, upowuu(d, 1 + u_lval(t,d)));
304 : }
305 : }
306 951 : if (faet) *faet = P;
307 951 : return shifti(zv_prod_Z(L), 2 + u_lval(t,2));
308 : }
309 :
310 : /* Table obtained by the following script:
311 :
312 : install(compute_e, "LD&"); \\ remove 'static' first
313 : table(first = 6, step = 6, MAXT = 6983776800)=
314 : {
315 : emax = 0;
316 : forstep(t = first, MAXT, step,
317 : e = compute_e(t);
318 : if (e > 1.9*emax, emax = e;
319 : printf(" if (C < %7.2f) return %8d;\n", 2*log(e)/log(2) - 1e-2, t)
320 : );
321 : );
322 : }
323 : table(,, 147026880);
324 : table(147026880,5040, 6983776800);
325 : */
326 :
327 : /* assume C < 20003.8 */
328 : static ulong
329 951 : compute_t_small(double C)
330 : {
331 951 : if (C < 17.94) return 6;
332 951 : if (C < 31.99) return 12;
333 951 : if (C < 33.99) return 24;
334 951 : if (C < 54.07) return 36;
335 951 : if (C < 65.32) return 60;
336 951 : if (C < 68.45) return 72;
337 951 : if (C < 70.78) return 108;
338 951 : if (C < 78.04) return 120;
339 951 : if (C < 102.41) return 180;
340 951 : if (C < 127.50) return 360;
341 949 : if (C < 136.68) return 420;
342 53 : if (C < 153.44) return 540;
343 53 : if (C < 165.67) return 840;
344 53 : if (C < 169.18) return 1008;
345 53 : if (C < 178.53) return 1080;
346 53 : if (C < 192.69) return 1200;
347 53 : if (C < 206.35) return 1260;
348 53 : if (C < 211.96) return 1620;
349 53 : if (C < 222.10) return 1680;
350 53 : if (C < 225.12) return 2016;
351 53 : if (C < 244.20) return 2160;
352 53 : if (C < 270.31) return 2520;
353 53 : if (C < 279.52) return 3360;
354 53 : if (C < 293.64) return 3780;
355 53 : if (C < 346.70) return 5040;
356 53 : if (C < 348.73) return 6480;
357 53 : if (C < 383.37) return 7560;
358 53 : if (C < 396.71) return 8400;
359 53 : if (C < 426.08) return 10080;
360 53 : if (C < 458.38) return 12600;
361 53 : if (C < 527.20) return 15120;
362 53 : if (C < 595.43) return 25200;
363 53 : if (C < 636.34) return 30240;
364 7 : if (C < 672.58) return 42840;
365 7 : if (C < 684.96) return 45360;
366 7 : if (C < 708.84) return 55440;
367 7 : if (C < 771.37) return 60480;
368 7 : if (C < 775.93) return 75600;
369 7 : if (C < 859.69) return 85680;
370 7 : if (C < 893.24) return 100800;
371 7 : if (C < 912.35) return 110880;
372 7 : if (C < 966.22) return 128520;
373 7 : if (C < 1009.18) return 131040;
374 0 : if (C < 1042.04) return 166320;
375 0 : if (C < 1124.98) return 196560;
376 0 : if (C < 1251.09) return 257040;
377 0 : if (C < 1375.13) return 332640;
378 0 : if (C < 1431.11) return 393120;
379 0 : if (C < 1483.46) return 514080;
380 0 : if (C < 1546.46) return 655200;
381 0 : if (C < 1585.94) return 665280;
382 0 : if (C < 1661.44) return 786240;
383 0 : if (C < 1667.67) return 831600;
384 0 : if (C < 1677.07) return 917280;
385 0 : if (C < 1728.17) return 982800;
386 0 : if (C < 1747.57) return 1081080;
387 0 : if (C < 1773.76) return 1179360;
388 0 : if (C < 1810.81) return 1285200;
389 0 : if (C < 1924.66) return 1310400;
390 0 : if (C < 2001.27) return 1441440;
391 0 : if (C < 2096.51) return 1663200;
392 0 : if (C < 2166.02) return 1965600;
393 0 : if (C < 2321.86) return 2162160;
394 0 : if (C < 2368.45) return 2751840;
395 0 : if (C < 2377.39) return 2827440;
396 0 : if (C < 2514.97) return 3326400;
397 0 : if (C < 2588.72) return 3341520;
398 0 : if (C < 2636.84) return 3603600;
399 0 : if (C < 2667.46) return 3931200;
400 0 : if (C < 3028.92) return 4324320;
401 0 : if (C < 3045.76) return 5654880;
402 0 : if (C < 3080.78) return 6652800;
403 0 : if (C < 3121.88) return 6683040;
404 0 : if (C < 3283.38) return 7207200;
405 0 : if (C < 3514.94) return 8648640;
406 0 : if (C < 3725.71) return 10810800;
407 0 : if (C < 3817.49) return 12972960;
408 0 : if (C < 3976.57) return 14414400;
409 0 : if (C < 3980.72) return 18378360;
410 0 : if (C < 4761.70) return 21621600;
411 0 : if (C < 5067.62) return 36756720;
412 0 : if (C < 5657.30) return 43243200;
413 0 : if (C < 5959.24) return 64864800;
414 0 : if (C < 6423.60) return 73513440;
415 0 : if (C < 6497.01) return 86486400;
416 0 : if (C < 6529.89) return 113097600;
417 0 : if (C < 6899.19) return 122522400;
418 0 : if (C < 7094.26) return 129729600;
419 0 : if (C < 7494.60) return 147026880;
420 0 : if (C < 7606.21) return 172972800;
421 0 : if (C < 7785.10) return 183783600;
422 0 : if (C < 7803.68) return 216216000;
423 0 : if (C < 8024.18) return 220540320;
424 0 : if (C < 8278.12) return 245044800;
425 0 : if (C < 8316.48) return 273873600;
426 0 : if (C < 8544.02) return 294053760;
427 0 : if (C < 8634.14) return 302702400;
428 0 : if (C < 9977.69) return 367567200;
429 0 : if (C < 10053.06) return 514594080;
430 0 : if (C < 10184.29) return 551350800;
431 0 : if (C < 11798.33) return 735134400;
432 0 : if (C < 11812.60) return 821620800;
433 0 : if (C < 11935.31) return 1029188160;
434 0 : if (C < 12017.99) return 1074427200;
435 0 : if (C < 12723.99) return 1102701600;
436 0 : if (C < 13702.71) return 1470268800;
437 0 : if (C < 13748.76) return 1643241600;
438 0 : if (C < 13977.37) return 2058376320;
439 0 : if (C < 14096.03) return 2148854400UL;
440 0 : if (C < 15082.25) return 2205403200UL;
441 0 : if (C < 15344.18) return 2572970400UL;
442 0 : if (C < 15718.37) return 2940537600UL;
443 0 : if (C < 15868.65) return 3491888400UL;
444 0 : if (C < 15919.88) return 3675672000UL;
445 0 : if (C < 16217.23) return 4108104000UL;
446 : #ifdef LONG_IS_64BIT
447 0 : if (C < 17510.32) return 4410806400UL;
448 0 : if (C < 18312.87) return 5145940800UL;
449 0 : return 6983776800UL;
450 : #else
451 0 : pari_err_IMPL("APRCL for large numbers on 32bit arch");
452 0 : return 0;
453 : #endif
454 : }
455 :
456 : /* return t such that e(t) > sqrt(N), set *faet = odd prime divisors of e(t) */
457 : static ulong
458 951 : compute_t(GEN N, GEN *e, GEN *faet)
459 : { /* 2^e b <= N < 2^e (b+1), where b >= 2^52. Approximate log_2 N ~ e+log2(b),
460 : * the error is less than log(1+1/b) < 1e-15 */
461 951 : double C = dbllog2(N) + 1e-10; /* > log_2 N at least for N < 2^(2^21) */
462 : ulong t;
463 : /* Return "smallest" t such that f(t) >= C, which implies e(t) > sqrt(N) */
464 : /* For N < 2^20003.8 ~ 5.5 10^6021 */
465 951 : if (C < 20003.8)
466 : {
467 951 : t = compute_t_small(C);
468 951 : *e = compute_e(t, faet);
469 : }
470 : else
471 : {
472 : #ifdef LONG_IS_64BIT
473 0 : GEN B = sqrti(N);
474 0 : for (t = 6983776800UL+5040UL;; t+=5040)
475 0 : {
476 0 : pari_sp av = avma;
477 0 : *e = compute_e(t, faet);
478 0 : if (cmpii(*e, B) > 0) break;
479 0 : set_avma(av);
480 : }
481 : #else
482 : *e = NULL; /* LCOV_EXCL_LINE */
483 : t = 0; /* LCOV_EXCL_LINE */
484 : #endif
485 : }
486 951 : return t;
487 : }
488 :
489 : /* T[i] = discrete log of i in (Z/q)^*, q odd prime
490 : * To save on memory, compute half the table: T[q-x] = T[x] + (q-1)/2 */
491 : static GEN
492 19063 : computetabdl(ulong q)
493 : {
494 19063 : ulong g, a, i, qs2 = q>>1; /* (q-1)/2 */
495 19063 : GEN T = cgetg(qs2+2,t_VECSMALL);
496 :
497 19063 : g = pgener_Fl(q); a = 1;
498 6089505 : for (i=1; i < qs2; i++) /* g^((q-1)/2) = -1 */
499 : {
500 6070442 : a = Fl_mul(g,a,q);
501 6070442 : if (a > qs2) T[q-a] = i+qs2; else T[a] = i;
502 : }
503 19063 : T[qs2+1] = T[qs2] + qs2;
504 19063 : T[1] = 0; return T;
505 : }
506 :
507 : /* Return T: T[x] = dl of x(1-x) */
508 : static GEN
509 13008 : compute_g(ulong q)
510 : {
511 13008 : const ulong qs2 = q>>1; /* (q-1)/2 */
512 : ulong x, a;
513 13008 : GEN T = computetabdl(q); /* updated in place to save on memory */
514 13008 : a = 0; /* dl[1] */
515 3187353 : for (x=2; x<=qs2+1; x++)
516 : { /* a = dl(x) */
517 3174345 : ulong b = T[x]; /* = dl(x) */
518 3174345 : T[x] = b + a + qs2; /* dl(x) + dl(x-1) + dl(-1) */
519 3174345 : a = b;
520 : }
521 13008 : return T;
522 : }
523 :
524 : /* x a nonzero VECSMALL with >= 0 entries */
525 : static GEN
526 20228 : zv_to_ZX(GEN x)
527 : {
528 20228 : long i,j, lx = lg(x);
529 : GEN y;
530 :
531 20283 : while (lx-- && x[lx]==0) /* empty */;
532 20228 : i = lx+2; y = cgetg(i,t_POL); y[1] = evalsigne(1);
533 132329 : for (j=2; j<i; j++) gel(y,j) = utoi(x[j-1]);
534 20228 : return y;
535 : }
536 : /* p odd prime */
537 : static GEN
538 20228 : get_jac(GEN C, ulong q, long pk, GEN tabg)
539 : {
540 20228 : ulong x, qs2 = q>>1; /* (q-1)/2 */
541 20228 : GEN vpk = zero_zv(pk);
542 :
543 9929686 : for (x=2; x<=qs2; x++) vpk[ tabg[x]%pk + 1 ] += 2;
544 20228 : vpk[ tabg[x]%pk + 1 ]++; /* x = (q+1)/2 */
545 20228 : return ZX_rem(zv_to_ZX(vpk), cache_cyc(C));
546 : }
547 :
548 : /* p = 2 */
549 : static GEN
550 6055 : get_jac2(GEN N, ulong q, long k, GEN *j2q, GEN *j3q)
551 : {
552 6055 : GEN jpq, vpk, T = computetabdl(q);
553 : ulong x, pk, i, qs2;
554 :
555 : /* could store T[x+1] + T[x] + qs2 (cf compute_g).
556 : * Recompute instead, saving half the memory. */
557 6055 : pk = 1UL << k;;
558 6055 : vpk = zero_zv(pk);
559 :
560 6055 : qs2 = q>>1; /* (q-1)/2 */
561 :
562 2915160 : for (x=2; x<=qs2; x++) vpk[ (T[x]+T[x-1]+qs2)%pk + 1 ] += 2;
563 6055 : vpk[ (T[x]+T[x-1]+qs2)%pk + 1 ]++;
564 6055 : jpq = u_red_cyclo2n_ip(vpk, k);
565 :
566 6055 : if (k == 2) return jpq;
567 :
568 1021 : if (mod8(N) >= 5)
569 : {
570 256 : GEN v8 = cgetg(9,t_VECSMALL);
571 2304 : for (x=1; x<=8; x++) v8[x] = 0;
572 464072 : for (x=2; x<=qs2; x++) v8[ ((3*T[x]+T[x-1]+qs2)&7) + 1 ]++;
573 464328 : for ( ; x<=q-1; x++) v8[ ((3*T[q-x]+T[q-x+1]-3*qs2)&7) + 1 ]++;
574 256 : *j2q = RgX_inflate(ZX_sqr(u_red_cyclo2n_ip(v8,3)), pk>>3);
575 256 : *j2q = red_cyclo2n_ip(*j2q, k);
576 : }
577 19837 : for (i=1; i<=pk; i++) vpk[i] = 0;
578 2558144 : for (x=2; x<=qs2; x++) vpk[ (2*T[x]+T[x-1]+qs2)%pk + 1 ]++;
579 2559165 : for ( ; x<=q-1; x++) vpk[ (2*T[q-x]+T[q-x+1]-2*qs2)%pk + 1 ]++;
580 1021 : *j3q = ZX_mul(jpq, u_red_cyclo2n_ip(vpk,k));
581 1021 : *j3q = red_cyclo2n_ip(*j3q, k);
582 1021 : return jpq;
583 : }
584 :
585 : /* N = 1 mod p^k, return an elt of order p^k in (Z/N)^* */
586 : static GEN
587 1376 : finda(GEN Cp, GEN N, long pk, long p)
588 : {
589 : GEN a, pv;
590 1376 : if (Cp && !isintzero(cache_a(Cp))) {
591 35 : a = cache_a(Cp);
592 35 : pv = cache_pk(Cp);
593 : }
594 : else
595 : {
596 : GEN ph, b, q;
597 1341 : ulong u = 2;
598 1341 : long v = Z_lvalrem(subiu(N,1), p, &q);
599 1341 : ph = powuu(p, v-1); pv = muliu(ph, p); /* N - 1 = p^v q */
600 1341 : if (p > 2)
601 : {
602 448 : for (;;u++)
603 : {
604 1306 : a = Fp_pow(utoipos(u), q, N);
605 1306 : b = Fp_pow(a, ph, N);
606 1306 : if (!gequal1(b)) break;
607 : }
608 : }
609 : else
610 : {
611 1806 : while (krosi(u,N) >= 0) u++;
612 483 : a = Fp_pow(utoipos(u), q, N);
613 483 : b = Fp_pow(a, ph, N);
614 : }
615 : /* checking b^p = 1 mod N done economically in caller */
616 1341 : b = gcdii(subiu(b,1), N);
617 1341 : if (!gequal1(b)) return NULL;
618 :
619 1341 : if (Cp) {
620 1341 : cache_a(Cp) = a; /* a has order p^v */
621 1341 : cache_pk(Cp) = pv;
622 : }
623 : }
624 1376 : return Fp_pow(a, divis(pv, pk), N);
625 : }
626 :
627 : /* return 0: N not a prime, 1: no problem so far */
628 : static long
629 5022 : filltabs(GEN C, GEN Cp, Red *R, long p, long pk, long ltab)
630 : {
631 : pari_sp av;
632 : long i, j;
633 : long e;
634 : GEN tabt, taba, m;
635 :
636 5022 : cache_cyc(C) = polcyclo(pk,0);
637 :
638 5022 : if (p > 2)
639 : {
640 2959 : long LE = pk - pk/p + 1;
641 2959 : GEN E = cgetg(LE, t_VECSMALL), eta = cgetg(pk+1,t_VEC);
642 16079 : for (i=1,j=0; i<pk; i++)
643 13120 : if (i%p) E[++j] = i;
644 2959 : cache_E(C) = E;
645 :
646 19038 : for (i=1; i<=pk; i++)
647 : {
648 16079 : GEN z = FpX_rem(pol_xn(i-1, 0), cache_cyc(C), R->N);
649 16079 : gel(eta,i) = FpX_center_i(z, R->N, R->N2);
650 : }
651 2959 : cache_eta(C) = eta;
652 : }
653 2063 : else if (pk >= 8)
654 : {
655 161 : long LE = (pk>>2) + 1;
656 161 : GEN E = cgetg(LE, t_VECSMALL);
657 2984 : for (i=1,j=0; i<pk; i++)
658 2823 : if ((i%8)==1 || (i%8)==3) E[++j] = i;
659 161 : cache_E(C) = E;
660 : }
661 :
662 5022 : if (pk > 2 && umodiu(R->N,pk) == 1)
663 : {
664 1376 : GEN vpa, p1, p2, p3, a2 = NULL, a = finda(Cp, R->N, pk, p);
665 : long jj, ph;
666 :
667 1376 : if (!a) return 0;
668 1376 : ph = pk - pk/p;
669 1376 : vpa = cgetg(ph+1,t_COL); gel(vpa,1) = a;
670 1376 : if (pk > p) a2 = modZ(sqri(a), R);
671 1376 : jj = 1;
672 4464 : for (i=2; i<pk; i++) /* vpa = { a^i, (i,p) = 1 } */
673 3088 : if (i%p) {
674 2444 : GEN z = mulii((i%p==1) ? a2 : a, gel(vpa,jj));
675 2444 : gel(vpa,++jj) = modZ(z, R);
676 : }
677 1376 : if (!gequal1( modZ( mulii(a, gel(vpa,ph)), R) )) return 0;
678 1376 : p1 = cgetg(ph+1,t_MAT);
679 1376 : p2 = cgetg(ph+1,t_COL); gel(p1,1) = p2;
680 5196 : for (i=1; i<=ph; i++) gel(p2,i) = gen_1;
681 1376 : j = 2; gel(p1,j) = vpa; p3 = vpa;
682 2444 : for (j++; j <= ph; j++)
683 : {
684 1068 : p2 = cgetg(ph+1,t_COL); gel(p1,j) = p2;
685 7132 : for (i=1; i<=ph; i++)
686 6064 : gel(p2,i) = modZ(mulii(gel(vpa,i),gel(p3,i)), R);
687 1068 : p3 = p2;
688 : }
689 1376 : cache_mat(C) = p1;
690 1376 : cache_matinv(C) = FpM_inv(p1, R->N);
691 : }
692 :
693 5022 : tabt = cgetg(ltab+1, t_VECSMALL);
694 5022 : taba = cgetg(ltab+1, t_VECSMALL);
695 5022 : av = avma; m = divis(R->N, pk);
696 116320 : for (e=1; e<=ltab && signe(m); e++)
697 : {
698 111298 : long s = vali(m); m = shifti(m,-s);
699 111298 : tabt[e] = e==1? s: s + R->k;
700 111298 : taba[e] = signe(m)? ((mod2BIL(m) & R->mask)+1)>>1: 0;
701 111298 : m = shifti(m, -R->k);
702 : }
703 5022 : setlg(taba, e); cache_aall(C) = taba;
704 5022 : setlg(tabt, e); cache_tall(C) = tabt;
705 5022 : return gc_long(av,1);
706 : }
707 :
708 : static GEN
709 951 : calcglobs(Red *R, ulong t, long *plpC, long *pltab, GEN *pP)
710 : {
711 : GEN fat, P, E, PE;
712 : long lv, i, k, b;
713 : GEN pC;
714 :
715 951 : b = expi(R->N)+1;
716 1962 : k = 3; while (((k+1)*(k+2) << (k-1)) < b) k++;
717 951 : *pltab = (b/k)+2;
718 951 : R->k = k;
719 951 : R->lv = 1L << (k-1);
720 951 : R->mask = (1UL << k) - 1;
721 :
722 951 : fat = factoru_pow(t);
723 951 : P = gel(fat,1);
724 951 : E = gel(fat,2);
725 951 : PE= gel(fat,3);
726 951 : *plpC = lv = vecsmall_max(PE); /* max(p^e, p^e | t) */
727 951 : pC = zerovec(lv);
728 951 : gel(pC,1) = zerovec(9); /* to be used as temp in step5() */
729 7986 : for (i = 2; i <= lv; i++) gel(pC,i) = gen_0;
730 4760 : for (i=1; i<lg(P); i++)
731 : {
732 3809 : long pk, p = P[i], e = E[i];
733 3809 : pk = p;
734 8831 : for (k=1; k<=e; k++, pk*=p)
735 : {
736 5022 : gel(pC,pk) = zerovec(9);
737 5022 : if (!filltabs(gel(pC,pk), gel(pC,p), R, p,pk, *pltab)) return NULL;
738 : }
739 : }
740 951 : *pP = P; return pC;
741 : }
742 :
743 : /* sig_a^{-1}(z) for z in Q(zeta_pk) and sig_a: zeta -> zeta^a. Assume
744 : * a reduced mod pk := p^k*/
745 : static GEN
746 142015 : aut(long pk, GEN z, long a)
747 : {
748 : GEN v;
749 142015 : long b, i, dz = degpol(z);
750 142015 : if (a == 1 || dz < 0) return z;
751 120766 : v = cgetg(pk+2,t_POL); v[1] = evalvarn(0); b = 0;
752 120766 : gel(v,2) = gel(z,2); /* i = 0 */
753 1162732 : for (i = 1; i < pk; i++)
754 : {
755 1041966 : b += a; if (b > pk) b -= pk; /* b = (a*i) % pk */
756 1041966 : gel(v,i+2) = b > dz? gen_0: gel(z,b+2);
757 : }
758 120766 : return normalizepol_lg(v, pk+2);
759 : }
760 :
761 : /* z^v for v in Z[G], represented by couples [sig_x^{-1},x] */
762 : static GEN
763 21249 : autvec_TH(long pk, GEN z, GEN v, GEN C)
764 : {
765 21249 : pari_sp av = avma;
766 21249 : long i, lv = lg(v);
767 21249 : GEN s = pol_1(varn(C));
768 113091 : for (i=1; i<lv; i++)
769 : {
770 91842 : long y = v[i];
771 91842 : if (y) s = ZXQ_mul(s, ZXQ_powu(aut(pk,z, y), y, C), C);
772 : }
773 21249 : return gc_upto(av,s);
774 : }
775 :
776 : static GEN
777 21249 : autvec_AL(long pk, GEN z, GEN v, Red *R)
778 : {
779 21249 : pari_sp av = avma;
780 21249 : const long r = umodiu(R->N, pk);
781 21249 : GEN s = pol_1(varn(R->C));
782 21249 : long i, lv = lg(v);
783 113091 : for (i=1; i<lv; i++)
784 : {
785 91842 : long y = (r*v[i]) / pk;
786 91842 : if (y) s = RgXQ_mul(s, ZXQ_powu(aut(pk,z, v[i]), y, R->C), R->C);
787 : }
788 21249 : return gc_upto(av,s);
789 : }
790 :
791 : /* 0 <= i < pk, such that x^i = z mod polcyclo(pk), -1 if no such i exist */
792 : static long
793 20228 : look_eta(GEN eta, long pk, GEN z)
794 : {
795 : long i;
796 65570 : for (i=1; i<=pk; i++)
797 65570 : if (ZX_equal(z, gel(eta,i))) return i-1;
798 0 : return -1;
799 : }
800 : /* same pk = 2^k */
801 : static long
802 6055 : look_eta2(long k, GEN z)
803 : {
804 : long d, s;
805 6055 : if (typ(z) != t_POL) d = 0; /* t_INT */
806 : else
807 : {
808 6055 : if (!RgX_is_monomial(z)) return -1;
809 6055 : d = degpol(z);
810 6055 : z = gel(z,d+2); /* leading term */
811 : }
812 6055 : s = signe(z);
813 6055 : if (!s || !is_pm1(z)) return -1;
814 6055 : return (s > 0)? d: d + (1L<<(k-1));
815 : }
816 :
817 : static long
818 20228 : step4a(GEN C, Red *R, ulong q, long p, long k, GEN tabg)
819 : {
820 20228 : const long pk = upowuu(p,k);
821 : long ind;
822 : GEN jpq, s1, s2, s3;
823 :
824 20228 : if (!tabg) tabg = compute_g(q);
825 20228 : jpq = get_jac(C, q, pk, tabg);
826 20228 : s1 = autvec_TH(pk, jpq, cache_E(C), cache_cyc(C));
827 20228 : s2 = powpolmod(C,R, p,k, s1);
828 20228 : s3 = autvec_AL(pk, jpq, cache_E(C), R);
829 20228 : s3 = _red(gmul(s3,s2), R);
830 :
831 20228 : ind = look_eta(cache_eta(C), pk, s3);
832 20228 : if (ind < 0) return -1;
833 20228 : return (ind%p) != 0;
834 : }
835 :
836 : /* x == -1 mod N ? */
837 : static long
838 6257 : is_m1(GEN x, GEN N)
839 : {
840 6257 : return equalii(addiu(x,1), N);
841 : }
842 :
843 : /* p=2, k>=3 */
844 : static long
845 1021 : step4b(GEN C, Red *R, ulong q, long k)
846 : {
847 1021 : const long pk = 1L << k;
848 : long ind;
849 1021 : GEN s1, s2, s3, j2q = NULL, j3q = NULL;
850 :
851 1021 : (void)get_jac2(R->N,q,k, &j2q,&j3q);
852 :
853 1021 : s1 = autvec_TH(pk, j3q, cache_E(C), cache_cyc(C));
854 1021 : s2 = powpolmod(C,R, 2,k, s1);
855 1021 : s3 = autvec_AL(pk, j3q, cache_E(C), R);
856 1021 : s3 = _red(gmul(s3,s2), R);
857 1021 : if (j2q) s3 = _red(gmul(j2q, s3), R);
858 :
859 1021 : ind = look_eta2(k, s3);
860 1021 : if (ind < 0) return -1;
861 1021 : if ((ind&1)==0) return 0;
862 448 : s3 = Fp_pow(utoipos(q), R->N2, R->N);
863 448 : return is_m1(s3, R->N);
864 : }
865 :
866 : /* p=2, k=2 */
867 : static long
868 5034 : step4c(GEN C, Red *R, ulong q)
869 : {
870 : long ind;
871 5034 : GEN s0,s1,s3, jpq = get_jac2(R->N,q,2, NULL,NULL);
872 :
873 5034 : s0 = sqrmod4(jpq, R);
874 5034 : s1 = gmulsg(q,s0);
875 5034 : s3 = powpolmod(C,R, 2,2, s1);
876 5034 : if (mod4(R->N) == 3) s3 = _red(gmul(s0,s3), R);
877 :
878 5034 : ind = look_eta2(2, s3);
879 5034 : if (ind < 0) return -1;
880 5034 : if ((ind&1)==0) return 0;
881 2415 : s3 = Fp_pow(utoipos(q), R->N2, R->N);
882 2415 : return is_m1(s3, R->N);
883 : }
884 :
885 : /* p=2, k=1 */
886 : static long
887 6953 : step4d(Red *R, ulong q)
888 : {
889 6953 : GEN s1 = Fp_pow(utoipos(q), R->N2, R->N);
890 6953 : if (is_pm1(s1)) return 0;
891 3394 : if (is_m1(s1, R->N)) return (mod4(R->N) == 1);
892 0 : return -1;
893 : }
894 :
895 : /* return 1 [OK so far] or 0 [not a prime] */
896 : static long
897 0 : step5(GEN pC, Red *R, long p, GEN et, ulong ltab, long lpC)
898 : {
899 : pari_sp av;
900 : ulong q;
901 0 : long pk, k, fl = -1;
902 : GEN C, Cp;
903 : forprime_t T;
904 :
905 0 : (void)u_forprime_arith_init(&T, 3, ULONG_MAX, 1,p);
906 0 : while( (q = u_forprime_next(&T)) )
907 : { /* q = 1 (mod p) */
908 0 : if (umodiu(et,q) == 0) continue;
909 0 : if (umodiu(R->N,q) == 0) return 0;
910 0 : k = u_lval(q-1, p);
911 0 : pk = upowuu(p,k);
912 0 : if (pk <= lpC && !isintzero(gel(pC,pk))) {
913 0 : C = gel(pC,pk);
914 0 : Cp = gel(pC,p);
915 : } else {
916 0 : C = gel(pC,1);
917 0 : Cp = NULL;
918 0 : cache_mat(C) = gen_0; /* re-init */
919 : }
920 0 : av = avma;
921 0 : if (!filltabs(C, Cp, R, p, pk, ltab)) return 0;
922 0 : R->C = cache_cyc(C);
923 0 : if (p >= 3) fl = step4a(C,R, q,p,k, NULL);
924 0 : else if (k >= 3) fl = step4b(C,R, q,k);
925 0 : else if (k == 2) fl = step4c(C,R, q);
926 0 : else fl = step4d(R, q);
927 0 : if (fl == -1) return 0;
928 0 : if (fl == 1) return 1; /*OK*/
929 0 : set_avma(av);
930 : }
931 0 : pari_err_BUG("aprcl test fails! This is highly improbable");
932 : return 0;/*LCOV_EXCL_LINE*/
933 : }
934 :
935 : GEN
936 1180 : aprcl_step6_worker(GEN r, long t, GEN N, GEN N1, GEN et)
937 : {
938 : long i;
939 1180 : pari_sp av = avma;
940 2679349 : for (i=1; i<=t; i++)
941 : {
942 2678199 : r = remii(mulii(r,N1), et);
943 2678199 : if (equali1(r)) break;
944 2678169 : if (dvdii(N,r) && !equalii(r,N)) return gen_0; /* not prime */
945 2678169 : if ((i & 0x1f) == 0) r = gc_INT(av, r);
946 : }
947 1180 : return cgetg(1,t_VECSMALL);
948 : }
949 :
950 : /* return 1 if N prime, else 0 */
951 : static long
952 951 : step6(GEN N, ulong t, GEN et)
953 : {
954 951 : GEN worker, r, rk, N1 = remii(N, et);
955 951 : ulong k = 10000;
956 : ulong i;
957 951 : long pending = 0, res = 1;
958 : struct pari_mt pt;
959 : pari_sp btop;
960 951 : worker = snm_closure(is_entry("_aprcl_step6_worker"), mkvec3(N, N1, et));
961 951 : r = gen_1;
962 951 : rk = Fp_powu(N1, k, et);
963 951 : mt_queue_start_lim(&pt, worker, (t-1+k-1)/k);
964 951 : btop = avma;
965 2144 : for (i=1; (i<t && res) || pending; i+=k)
966 : {
967 : GEN done;
968 1193 : mt_queue_submit(&pt, i, (i<t && res)? mkvec2(r, utoi(minss(k,t-i))): NULL);
969 1193 : done = mt_queue_get(&pt, NULL, &pending);
970 1193 : if (res && done && typ(done) != t_VECSMALL) res = 0; /* not a prime */
971 1193 : r = Fp_mul(r, rk, et);
972 1193 : if (gc_needed(btop, 2))
973 : {
974 0 : if(DEBUGMEM>1) pari_warn(warnmem,"APRCL: i = %ld",i);
975 0 : r = gc_upto(btop, r);
976 : }
977 : }
978 951 : mt_queue_end(&pt);
979 951 : return res;
980 : }
981 :
982 : GEN
983 13008 : aprcl_step4_worker(ulong q, GEN pC, GEN N, GEN v)
984 : {
985 13008 : pari_sp av1 = avma, av2 = avma;
986 : long j, k;
987 : Red R;
988 13008 : GEN faq = factoru_pow(q-1), tabg = compute_g(q);
989 13008 : GEN P = gel(faq,1), E = gel(faq,2), PE = gel(faq,3);
990 13008 : long lfaq = lg(P);
991 13008 : GEN flags = cgetg(lfaq, t_VECSMALL);
992 13008 : R.N = N;
993 13008 : R.N2= shifti(N, -1);
994 13008 : R.k = v[1]; R.lv = v[2]; R.mask = uel(v,3); R.n = v[4];
995 13008 : av2 = avma;
996 46244 : for (j=1, k=1; j<lfaq; j++, set_avma(av2))
997 : {
998 33236 : long p = P[j], e = E[j], pe = PE[j], fl;
999 33236 : GEN C = gel(pC,pe);
1000 33236 : R.C = cache_cyc(C);
1001 33236 : if (p >= 3) fl = step4a(C,&R, q,p,e, tabg);
1002 13008 : else if (e >= 3) fl = step4b(C,&R, q,e);
1003 11987 : else if (e == 2) fl = step4c(C,&R, q);
1004 6953 : else fl = step4d(&R, q);
1005 33236 : if (fl == -1) return gen_0; /* not prime */
1006 33236 : if (fl == 1) flags[k++] = p;
1007 : }
1008 13008 : setlg(flags, k);
1009 13008 : return gc_leaf(av1, flags); /* OK so far */
1010 : }
1011 :
1012 : /* return 1 if prime, else 0 */
1013 : static long
1014 972 : aprcl(GEN N)
1015 : {
1016 972 : GEN pC, worker, et, fat, flaglp, faet = NULL; /*-Wall*/
1017 972 : long i, j, l, ltab, lfat, lpC, workid, pending = 0, res = 1;
1018 : ulong t;
1019 : Red R;
1020 : struct pari_mt pt;
1021 :
1022 972 : if (typ(N) != t_INT) pari_err_TYPE("aprcl",N);
1023 972 : if (cmpis(N,12) <= 0)
1024 21 : switch(itos(N))
1025 : {
1026 14 : case 2: case 3: case 5: case 7: case 11: return 1;
1027 7 : default: return 0;
1028 : }
1029 951 : if (Z_issquare(N)) return 0;
1030 951 : t = compute_t(N, &et, &faet);
1031 951 : if (DEBUGLEVEL) err_printf("Starting APRCL with t = %ld\n",t);
1032 951 : if (cmpii(sqri(et),N) < 0) pari_err_BUG("aprcl: e(t) too small");
1033 951 : if (!equali1(gcdii(N,mului(t,et)))) return 0;
1034 :
1035 951 : R.N = N;
1036 951 : R.N2= shifti(N, -1);
1037 951 : pC = calcglobs(&R, t, &lpC, <ab, &fat);
1038 951 : if (!pC) return 0;
1039 951 : lfat = lg(fat);
1040 951 : flaglp = cgetg(lfat, t_VECSMALL);
1041 951 : flaglp[1] = 0;
1042 3809 : for (i=2; i<lfat; i++)
1043 : {
1044 2858 : ulong p = fat[i];
1045 2858 : flaglp[i] = absequaliu(Fp_powu(N, p-1, sqru(p)), 1);
1046 : }
1047 951 : vecsmall_sort(faet);
1048 951 : l = lg(faet);
1049 951 : worker = snm_closure(is_entry("_aprcl_step4_worker"),
1050 951 : mkvec3(pC, N, mkvecsmall4(R.k, R.lv, R.mask, R.n)));
1051 951 : if (DEBUGLEVEL>2) err_printf("Step4: %ld q-values\n", l-1);
1052 951 : mt_queue_start_lim(&pt, worker, l-1);
1053 15382 : for (i=l-1; (i>0 && res) || pending; i--)
1054 : {
1055 : GEN done;
1056 14431 : ulong q = i>0 ? faet[i]: 0;
1057 14431 : mt_queue_submit(&pt, q, q>0? mkvec(utoi(q)): NULL);
1058 14431 : done = mt_queue_get(&pt, &workid, &pending);
1059 14431 : if (res && done)
1060 : {
1061 13008 : long lf = lg(done);
1062 13008 : if (typ(done) != t_VECSMALL) { res = 0; continue; } /* not prime */
1063 32875 : for (j=1; j<lf; j++) flaglp[ zv_search(fat, done[j]) ] = 0;
1064 13008 : if (DEBUGLEVEL>2) err_printf("testing Jacobi sums for q = %ld...OK\n", workid);
1065 : }
1066 : }
1067 951 : mt_queue_end(&pt);
1068 951 : if (!res) return 0;
1069 951 : if (DEBUGLEVEL>2) err_printf("Step5: testing conditions lp\n");
1070 3809 : for (i=2; i<lfat; i++) /*skip fat[1] = 2*/
1071 : {
1072 2858 : pari_sp av = avma;
1073 2858 : long p = fat[i];
1074 2858 : if (flaglp[i] && !step5(pC, &R, p, et, ltab, lpC)) return 0;
1075 2858 : set_avma(av);
1076 : }
1077 951 : if (DEBUGLEVEL>2)
1078 0 : err_printf("Step6: testing potential divisors\n");
1079 951 : return step6(N, t, et);
1080 : }
1081 :
1082 : long
1083 972 : isprimeAPRCL(GEN N)
1084 972 : { pari_sp av = avma; return gc_long(av, aprcl(N)); }
1085 :
1086 : /*******************************************************************/
1087 : /* DIVISORS IN RESIDUE CLASSES (LENSTRA) */
1088 : /*******************************************************************/
1089 : /* This would allow to replace e(t) > N^(1/2) by e(t) > N^(1/3), but step6
1090 : * becomes so expensive that, at least up to 6000 bits, this is useless
1091 : * in this application. */
1092 : static void
1093 385 : set_add(hashtable *H, void *d)
1094 : {
1095 385 : ulong h = H->hash(d);
1096 385 : if (!hash_search2(H, d, h)) hash_insert2(H, d, NULL, h);
1097 385 : }
1098 : static void
1099 336 : add(hashtable *H, GEN t1, GEN t2, GEN a, GEN b, GEN r, GEN s)
1100 : {
1101 336 : GEN ra, qa = dvmdii(t1, a, &ra);
1102 336 : if (!signe(ra) && dvdii(t2, b) && equalii(modii(qa, s), r))
1103 266 : set_add(H, (void*)qa);
1104 336 : }
1105 : /* T^2 - B*T + C has integer roots ? */
1106 : static void
1107 301 : check_t(hashtable *H, GEN B, GEN C4, GEN a, GEN b, GEN r, GEN s)
1108 : {
1109 301 : GEN d, t1, t2, D = subii(sqri(B), C4);
1110 301 : if (!Z_issquareall(D, &d)) return;
1111 245 : t1 = shifti(addii(B, d), -1); /* >= 0 */
1112 245 : t2 = subii(B, t1);
1113 245 : add(H, t1,t2, a,b,r,s);
1114 245 : if (signe(t2) >= 0) add(H, t2,t1, a,b,r,s);
1115 : }
1116 : /* N > s > r >= 0, (r,s) = 1 */
1117 : GEN
1118 63 : divisorslenstra(GEN N, GEN r, GEN s)
1119 : {
1120 63 : pari_sp av = avma;
1121 : GEN u, Ns2, rp, a0, a1, b0, b1, c0, c1, s2;
1122 63 : hashtable *H = hash_create_INT(11, 1);
1123 : long j;
1124 63 : if (typ(N) != t_INT) pari_err_TYPE("divisorslenstra", N);
1125 63 : if (typ(r) != t_INT) pari_err_TYPE("divisorslenstra", r);
1126 63 : if (typ(s) != t_INT) pari_err_TYPE("divisorslenstra", s);
1127 63 : u = Fp_inv(r, s);
1128 63 : rp = Fp_mul(u, N, s); /* r' */
1129 63 : s2 = sqri(s);
1130 63 : a0 = s;
1131 63 : b0 = gen_0;
1132 63 : c0 = gen_0;
1133 63 : if (dvdii(N, r)) set_add(H, (void*)r); /* case i = 0 */
1134 63 : a1 = Fp_mul(u, rp, s); if (!signe(a1)) a1 = s; /* 0 < a1 <= s */
1135 63 : b1 = gen_1;
1136 63 : c1 = Fp_mul(u, diviiexact(subii(N,mulii(r,rp)), s), s);
1137 63 : Ns2 = divii(N, s2);
1138 63 : j = 1;
1139 : for (;;)
1140 273 : {
1141 336 : GEN Cs, q, c, ab = mulii(a1,b1);
1142 : long i, lC;
1143 336 : if (j == 0) /* i even, a1 >= 0 */
1144 : {
1145 168 : if (!signe(c1)) Cs = mkvec(gen_0);
1146 : else
1147 : {
1148 112 : GEN cs = mulii(c1, s);
1149 112 : Cs = mkvec2(subii(cs,s2), cs);
1150 : }
1151 : }
1152 : else
1153 : { /* i odd, a1 > 0 */
1154 168 : GEN X = shifti(ab,1);
1155 168 : c = c1;
1156 : /* smallest c >= 2ab, c = c1 (mod s) */
1157 168 : if (cmpii(c, X) < 0)
1158 : {
1159 147 : GEN rX, qX = dvmdii(subii(X,c), s, &rX);
1160 147 : if (signe(rX)) qX = addiu(qX,1); /* ceil((X-c)/s) */
1161 147 : c = addii(c, mulii(s, qX));
1162 : }
1163 168 : Cs = (cmpii(c, addii(Ns2,ab)) <= 0)? mkvec(mulii(c,s)): cgetg(1,t_VEC);
1164 : }
1165 336 : lC = lg(Cs);
1166 336 : if (signe(a1))
1167 : {
1168 273 : GEN abN4 = shifti(mulii(ab, N), 2);
1169 273 : GEN B = addii(mulii(a1,r), mulii(b1,rp));
1170 574 : for (i = 1; i < lC; i++)
1171 301 : check_t(H, addii(B, gel(Cs,i)), abN4, a1, b1, r, s);
1172 : }
1173 : else
1174 : { /* a1 = 0, last batch */
1175 133 : for (i = 1; i < lC; i++)
1176 : {
1177 70 : GEN ry, ys = dvmdii(gel(Cs,i), b1, &ry);
1178 70 : if (!signe(ry))
1179 : {
1180 70 : GEN d = addii(ys, rp);
1181 70 : if (signe(d) > 0)
1182 : {
1183 56 : d = dvmdii(N, d, &ry);
1184 56 : if (!signe(ry)) set_add(H, (void*)d);
1185 : }
1186 : }
1187 : }
1188 63 : break; /* DONE */
1189 : }
1190 273 : j = 1-j;
1191 273 : q = dvmdii(a0, a1, &c);
1192 273 : if (j == 1 && !signe(c)) { q = subiu(q,1); c = a1; }
1193 273 : a0 = a1; a1 = c;
1194 273 : if (equali1(q)) /* frequent */
1195 : {
1196 119 : c = subii(b0, b1); b0 = b1; b1 = c;
1197 119 : c = Fp_sub(c0, c1, s); c0 = c1; c1 = c;
1198 : }
1199 : else
1200 : {
1201 154 : c = subii(b0, mulii(q,b1)); b0 = b1; b1 = c;
1202 154 : c = modii(subii(c0, mulii(q,c1)), s); c0 = c1; c1 = c;
1203 : }
1204 : }
1205 63 : return gc_upto(av, ZV_sort(hash_keys_GEN(H)));
1206 : }
|
