Line data Source code
1 : /* Copyright 2011-2015 David Cleaver
2 : *
3 : * This program is free software: you can redistribute it and/or modify
4 : * it under the terms of the GNU Lesser General Public License as published by
5 : * the Free Software Foundation, either version 3 of the License, or
6 : * (at your option) any later version.
7 : *
8 : * This program is distributed in the hope that it will be useful,
9 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
10 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11 : * GNU Lesser General Public License for more details.
12 : *
13 : * You should have received a copy of the GNU Lesser General Public License
14 : * along with this program. If not, see <http://www.gnu.org/licenses/>.
15 : */
16 :
17 : /*
18 : * v1.0 Posted to SourceForge on 2013/07/04
19 : *
20 : * v1.1 Posted to SourceForge on 2013/12/27
21 : * [The following improvements/fixes were recommended by Laurent Desnogues in 2013/08]
22 : * - Speed improvement 1: Removed extraneous NormalizeJS calls in ARPCL
23 : * - Speed improvement 2: Removed/consolidated calls to mpz_mod in APRCL
24 : * (these improvements make the APRCL code about 1.5-2.2x faster)
25 : * - Bug fix: Final test in APRCL routine is now correct
26 : */
27 :
28 :
29 : #include <stdio.h>
30 : #include <stdlib.h>
31 : #include <gmp.h>
32 : #include "mpz_aprcl.h"
33 :
34 : #ifndef HAVE_U64_T
35 : #define HAVE_U64_T
36 : typedef long long s64_t;
37 : typedef unsigned long long u64_t;
38 : #endif
39 :
40 :
41 : /* **********************************************************************************
42 : * APR-CL (also known as APRT-CLE) is a prime proving algorithm developed by:
43 : * L. Adleman, C. Pomerance, R. Rumely, H. Cohen, and H. W. Lenstra
44 : * APRT-CLE = Adleman-Pomerance-Rumely Test Cohen-Lenstra Extended version
45 : * You can find all the details of this implementation in the Cohen & Lenstra paper:
46 : * H. Cohen and A. K. Lenstra, "Implementation of a new primality test",
47 : * Math. Comp., 48 (1987) 103--121
48 : *
49 : * ----------------------------------------------------------------------------------
50 : *
51 : * This C/GMP version is a conversion of Dario Alpern's Java based APRT-CLE code
52 : * His code was based on Yuji Kida's UBASIC code
53 : *
54 : * Based on APRT-CLE Written by Dario Alejandro Alpern (Buenos Aires - Argentina)
55 : * From: Last updated September 10th, 2011. See http://www.alpertron.com.ar/ECM.HTM
56 : *
57 : * On 2012/11/12 Dario Alpern has approved the conversion, from Java to C/GMP, of
58 : * his implementation of the APR-CL algorithm, and that it be licensed under the LGPL.
59 : *
60 : * ----------------------------------------------------------------------------------
61 : *
62 : * With improvements based on Jason Moxham's APRCL v1.15 code, from 2003/01/01
63 : *
64 : * On 2013/04/14 Toby Moxham has approved the APR-CL code and data tables,
65 : * originally written by his brother Jason Moxham on 2003/01/01, to be released
66 : * under the LGPL.
67 : *
68 : * *********************************************************************************/
69 :
70 : /* int PWmax = 32, Qmax = 55441, LEVELmax = 11; */
71 : /* QMax is the largest Q in the aiQ array */
72 : /* PWmax is the largest p^k that divides any (Q-1) */
73 : /* #define Qmax 55441 */
74 : /* #define PWmax 32 */
75 : /*
76 : 2-max = 2^5 = 32 from t[ 4]= 4324320
77 : 3-max = 3^3 = 27 from t[ 4]= 4324320
78 : 5-max = 5^2 = 25 from t[ 6]= 367567200
79 : 7-max = 7^1 = 7 from t[ 1]= 5040
80 : 11-max = 11^1 = 11 from t[ 2]= 55440
81 : 13-max = 13^1 = 13 from t[ 3]= 720720
82 : 17-max = 17^1 = 17 from t[ 5]= 73513440
83 : 19-max = 19^1 = 19 from t[ 7]= 1396755360
84 : */
85 :
86 : /* largest value of p^k that divides any t value */
87 : #define PWmax 32
88 :
89 : /* Qmax[]: list of largest q-prime for each t value */
90 : int Qmax[] = {61,2521,55441,180181,4324321,10501921,367567201,232792561,
91 : 1745944201};
92 :
93 : /* number of values in Qmax[], aiNP[], aiNQ[] */
94 : int LEVELmax = 9;
95 :
96 : /* list of primes that divide our t values */
97 : int aiP[] = {2,3,5,7,11,13,17,19};
98 :
99 : /* list of q-primes from all t values */
100 : int aiQ[] = {2,3,5,7,11,13,31,61,17,19,29,37,41,43,71,73,113,127,181,211,241,
101 : 281,337,421,631,1009,2521,23,67,89,199,331,397,463,617,661,881,991,1321,2311,
102 : 3697,4621,9241,18481,55441,53,79,131,157,313,521,547,859,911,937,1093,1171,
103 : 1873,2003,2341,2731,2861,3121,3433,6007,6553,8009,8191,8581,16381,20021,20593,
104 : 21841,25741,36037,48049,51481,65521,72073,120121,180181,97,109,271,353,379,433,
105 : 541,673,757,1249,2017,2081,2161,2377,2971,3169,3361,3511,4159,5281,7393,7561,
106 : 7723,8317,8737,9829,13729,14561,15121,16633,23761,24571,26209,28081,30241,
107 : 38611,39313,47521,66529,96097,108109,110881,123553,131041,196561,216217,270271,
108 : 332641,393121,432433,540541,617761,4324321,103,137,239,307,409,443,613,919,953,
109 : 1021,1123,1327,1361,1429,1531,1871,2143,2381,2857,3061,3571,3673,4421,4591,
110 : 5237,6121,6427,6733,7481,8161,9181,9283,9521,10099,10711,12241,12377,12853,
111 : 14281,15913,16831,17137,17681,19891,22441,23563,23869,24481,27847,29173,29921,
112 : 30941,34273,36721,42841,43759,46411,47737,52361,53857,59671,63649,70687,72931,
113 : 74257,78541,79561,87517,92821,97241,100981,102103,116689,117811,128521,145861,
114 : 148513,157081,161569,167077,185641,201961,209441,235621,238681,269281,291721,
115 : 314161,371281,388961,417691,445537,471241,477361,514081,565489,612613,656371,
116 : 680681,700129,816817,1633633,1670761,1837837,2625481,4084081,5250961,5654881,
117 : 8168161,9189181,10501921,101,151,401,601,701,1051,1201,1301,1801,1951,2551,
118 : 2801,3301,3851,4201,4951,5101,5851,6301,7151,9901,11551,11701,12601,14851,
119 : 15401,15601,17551,17851,18701,19801,21601,23801,28051,33151,34651,40801,42901,
120 : 44201,50051,53551,54601,56101,66301,70201,77351,79201,81901,91801,92401,93601,
121 : 103951,107101,109201,118801,122401,140401,150151,151201,160651,193051,198901,
122 : 200201,218401,224401,232051,243101,257401,300301,321301,367201,415801,428401,
123 : 448801,450451,504901,530401,600601,673201,729301,795601,800801,982801,1029601,
124 : 1093951,1178101,1201201,1458601,2088451,2187901,2402401,2570401,2702701,
125 : 3088801,3141601,3712801,5105101,5834401,6806801,7068601,8353801,17503201,
126 : 22972951,52509601,183783601,367567201,191,229,419,457,571,647,761,1483,1597,
127 : 2053,2129,2281,2927,3041,3877,4447,4523,4561,4789,6271,6689,6841,6917,7411,
128 : 7753,8209,8779,8893,10337,11287,11971,12541,13339,13567,13681,14821,16417,
129 : 17291,17443,18089,19381,20521,20749,21319,21737,22573,25841,27361,28729,29641,
130 : 30097,31123,35531,35569,35911,38039,39521,40699,43891,46817,47881,48907,51871,
131 : 53353,56431,57457,58787,59281,63841,71821,72353,75583,77521,87211,90289,97813,
132 : 105337,106591,108529,114913,117041,124489,131671,134369,135661,139537,140449,
133 : 146719,163021,177841,186733,207481,213181,217057,217361,225721,251941,279073,
134 : 287281,300961,302329,342343,351121,377911,391249,406981,451441,456457,461891,
135 : 489061,511633,526681,554269,568481,608609,651169,652081,697681,733591,782497,
136 : 790021,813961,895357,1027027,1053361,1058149,1108537,1133731,1264033,1279081,
137 : 1369369,1492261,1580041,1790713,1813969,1867321,1939939,2217073,2238391,
138 : 2282281,2351441,2489761,2645371,2771341,2934361,2984521,3233231,3627937,
139 : 3837241,3912481,3979361,4157011,4232593,4476781,5135131,5372137,5868721,
140 : 6046561,6348889,6651217,6715171,6846841,7162849,7674481,9767521,11737441,
141 : 12471031,12697777,17907121,24942061,27387361,31744441,35814241,41081041,
142 : 46558513,53721361,107442721,174594421,232792561,1901,2851,5701,39901,41801,
143 : 53201,62701,64601,74101,79801,98801,113051,119701,135851,148201,205201,219451,
144 : 290701,292601,319201,333451,339151,359101,410401,452201,478801,501601,532951,
145 : 564301,658351,666901,778051,839801,957601,1037401,1065901,1128601,1222651,
146 : 1259701,1504801,1808801,1889551,2074801,2173601,2445301,2667601,3052351,
147 : 3511201,3730651,3779101,3950101,4069801,4149601,4408951,5038801,6104701,
148 : 6224401,8558551,9781201,11191951,11411401,14922601,16279201,17117101,17635801,
149 : 19186201,19562401,22383901,22822801,23514401,25581601,25675651,31600801,
150 : 35271601,37346401,38372401,45349201,59690401,67151701,83140201,129329201,
151 : 134303401,193993801,249420601,436486051,634888801,1163962801,1745944201};
152 : /* number of q primes in the above array: 618 */
153 :
154 : /* a primitive root for each q-prime in the above array */
155 : int aiG[] = {1,2,2,3,2,2,3,2,3,2,2,2,6,3,7,5,3,3,2,2,7,3,10,2,3,11,17,5,2,3,3,
156 : 3,5,3,3,2,3,6,13,3,5,2,13,13,38,2,3,2,5,10,3,2,2,17,5,5,2,10,5,7,3,2,7,5,3,10,
157 : 3,17,6,2,3,5,11,6,2,17,17,17,5,29,6,5,6,6,3,2,5,2,5,2,7,5,3,23,5,10,7,22,7,3,7,
158 : 5,13,3,6,5,10,23,6,11,15,7,7,11,19,11,3,10,17,19,5,2,69,5,17,11,15,7,19,23,5,
159 : 14,19,17,5,3,7,5,21,2,2,7,3,10,2,3,3,6,2,14,3,3,11,6,2,5,3,11,3,7,3,2,6,7,2,2,
160 : 3,2,3,7,6,5,19,5,6,5,3,2,14,2,2,11,6,2,3,2,5,37,23,3,3,5,6,5,12,7,5,10,5,2,7,2,
161 : 6,3,7,5,7,2,22,2,5,26,13,5,41,13,3,10,29,7,14,37,19,3,12,29,19,14,33,13,2,2,6,
162 : 28,5,5,19,5,7,19,29,13,23,6,7,2,6,3,7,2,7,11,2,11,3,6,3,6,2,11,6,6,2,10,7,2,7,
163 : 6,11,2,6,23,3,2,2,13,7,3,2,3,2,13,6,6,2,3,22,7,22,14,14,13,2,13,34,22,6,6,17,
164 : 17,7,11,6,17,2,2,2,3,22,31,3,2,29,2,2,11,31,19,13,2,2,23,37,23,7,19,3,17,7,3,6,
165 : 31,23,7,2,19,26,6,31,26,29,10,14,3,19,29,37,6,37,41,23,19,6,2,13,3,5,6,2,11,2,
166 : 3,7,5,3,2,3,5,11,2,11,3,22,2,2,10,7,11,5,3,3,10,14,2,3,22,2,10,6,2,3,7,11,2,14,
167 : 6,6,3,7,22,7,10,3,6,11,12,7,3,2,3,3,29,2,7,7,3,5,2,7,17,2,3,5,7,13,23,2,10,3,
168 : 23,5,3,11,3,3,10,5,17,6,6,7,5,31,10,10,6,17,6,10,13,7,7,3,29,3,7,6,29,5,18,17,
169 : 13,29,6,3,3,22,14,14,6,10,17,13,6,7,34,2,5,2,10,31,43,6,13,13,21,29,2,5,7,17,3,
170 : 22,7,7,7,29,14,5,13,21,6,10,15,6,2,5,14,14,11,5,7,23,13,7,37,29,11,5,13,22,37,
171 : 58,26,29,5,43,23,2,71,2,2,2,2,3,3,2,3,2,23,3,6,2,2,17,14,2,2,3,23,26,3,2,11,11,
172 : 29,31,15,2,7,10,3,3,22,11,6,14,3,6,31,6,3,47,3,10,7,6,13,10,6,6,13,17,3,7,2,11,
173 : 6,42,3,23,13,37,26,11,21,7,6,37,6,7,2,13,29,59,26,22,59,10,31,3,23,53,42,19,11,
174 : 46,23};
175 :
176 :
177 : /* number of primes from aiP, not necessarily in order, that divides each t */
178 : int aiNP[] = {3,4,5,6,6,7,7,8,8};
179 :
180 : /* number of q-primes for each t */
181 : int aiNQ[] = {8,27,45,81,134,245,351,424,618};
182 :
183 : /* t | e(t) | #Qp | Qmax | divisors of t */
184 : /* --------------|--------------|-----|------------|------------------------ */
185 : s64_t aiT[] = {
186 : 60, /* | 6.8144 E 9 | 8 | 61 | p={2,3,5} */
187 : 5040, /* | 1.5321 E 52 | 27 | 2521 | p={2,3,5,7} */
188 : 55440, /* | 4.9209 E 106 | 45 | 55441 | p={2,3,5,7,11} */
189 : 720720, /* | 2.5992 E 237 | 81 | 180181 | p={2,3,5,7,11,13} */
190 : 4324320, /* | 7.9285 E 455 | 134 | 4324321 | p={2,3,5,7,11,13} */
191 : 73513440, /* | 7.0821 E 966 | 245 | 10501921 | p={2,3,5,7,11,13,17} */
192 : 367567200, /* | 6.2087 E1501 | 351 | 367567201 | p={2,3,5,7,11,13,17} */
193 : 1396755360, /* | 4.0165 E1913 | 424 | 232792561 | p={2,3,5,7,11,13,17,19} */
194 : 6983776800};/* | 7.4712 E3010 | 618 | 1745944201 | p={2,3,5,7,11,13,17,19} */
195 :
196 :
197 : int aiInv[PWmax];
198 : mpz_t biTmp;
199 : mpz_t biExp;
200 : mpz_t biN;
201 : mpz_t biR;
202 : mpz_t biS;
203 : mpz_t biT;
204 : mpz_t *aiJS; /* [PWmax] */
205 : mpz_t *aiJW; /* [PWmax] */
206 : mpz_t *aiJX; /* [PWmax] */
207 : mpz_t *aiJ0; /* [PWmax] */
208 : mpz_t *aiJ1; /* [PWmax] */
209 : mpz_t *aiJ2; /* [PWmax] */
210 : mpz_t *aiJ00; /* [PWmax] */
211 : mpz_t *aiJ01; /* [PWmax] */
212 : int NumberLength; /* Length of multiple precision nbrs */
213 : mpz_t TestNbr;
214 :
215 : /* ============================================================================================== */
216 :
217 2587 : void allocate_vars()
218 : {
219 2587 : int i = 0;
220 2587 : aiJS = malloc(PWmax * sizeof(mpz_t));
221 2587 : aiJW = malloc(PWmax * sizeof(mpz_t));
222 2587 : aiJX = malloc(PWmax * sizeof(mpz_t));
223 2587 : aiJ0 = malloc(PWmax * sizeof(mpz_t));
224 2587 : aiJ1 = malloc(PWmax * sizeof(mpz_t));
225 2587 : aiJ2 = malloc(PWmax * sizeof(mpz_t));
226 2587 : aiJ00 = malloc(PWmax * sizeof(mpz_t));
227 2587 : aiJ01 = malloc(PWmax * sizeof(mpz_t));
228 85371 : for (i = 0 ; i < PWmax; i++)
229 : {
230 82784 : mpz_init(aiJS[i]);
231 82784 : mpz_init(aiJW[i]);
232 82784 : mpz_init(aiJX[i]);
233 82784 : mpz_init(aiJ0[i]);
234 82784 : mpz_init(aiJ1[i]);
235 82784 : mpz_init(aiJ2[i]);
236 82784 : mpz_init(aiJ00[i]);
237 82784 : mpz_init(aiJ01[i]);
238 : }
239 :
240 2587 : mpz_init(TestNbr);
241 2587 : mpz_init(biN);
242 2587 : mpz_init(biR);
243 2587 : mpz_init(biS);
244 2587 : mpz_init(biT);
245 2587 : mpz_init(biExp);
246 2587 : mpz_init(biTmp);
247 2587 : }
248 :
249 : /* ============================================================================================== */
250 :
251 2587 : void free_vars()
252 : {
253 2587 : int i = 0;
254 85371 : for (i = 0 ; i < PWmax; i++)
255 : {
256 82784 : mpz_clear(aiJS[i]);
257 82784 : mpz_clear(aiJW[i]);
258 82784 : mpz_clear(aiJX[i]);
259 82784 : mpz_clear(aiJ0[i]);
260 82784 : mpz_clear(aiJ1[i]);
261 82784 : mpz_clear(aiJ2[i]);
262 82784 : mpz_clear(aiJ00[i]);
263 82784 : mpz_clear(aiJ01[i]);
264 : }
265 2587 : free(aiJS);
266 2587 : free(aiJW);
267 2587 : free(aiJX);
268 2587 : free(aiJ0);
269 2587 : free(aiJ1);
270 2587 : free(aiJ2);
271 2587 : free(aiJ00);
272 2587 : free(aiJ01);
273 :
274 2587 : mpz_clear(TestNbr);
275 2587 : mpz_clear(biN);
276 2587 : mpz_clear(biR);
277 2587 : mpz_clear(biS);
278 2587 : mpz_clear(biT);
279 2587 : mpz_clear(biExp);
280 2587 : mpz_clear(biTmp);
281 2587 : }
282 :
283 : /* ============================================================================================== */
284 :
285 : // Compare Nbr1^2 vs. Nbr2
286 : // -1 -> Nbr1^2 < Nbr2
287 : // 0 -> Nbr1^2 == Nbr2
288 : // 1 -> Nbr1^2 > Nbr2
289 68636 : int CompareSquare(mpz_t Nbr1, mpz_t Nbr2)
290 : {
291 : mpz_t tmp;
292 68636 : int cmp = 0;
293 :
294 68636 : mpz_init(tmp);
295 68636 : mpz_mul(tmp, Nbr1, Nbr1);
296 :
297 68636 : cmp = mpz_cmp(tmp, Nbr2);
298 68636 : mpz_clear(tmp);
299 :
300 68636 : return cmp;
301 : }
302 :
303 : /* ============================================================================================== */
304 :
305 : // Normalize coefficient of JS
306 18143805 : void NormalizeJS(int PK, int PL, int PM, int P)
307 : {
308 : int I, J;
309 50106484 : for (I = PL; I < PK; I++)
310 : {
311 31962679 : if (mpz_cmp_ui(aiJS[I], 0) != 0) /* (!BigNbrIsZero(aiJS[I])) */
312 : {
313 : /* biT = aiJS[I]; */
314 25467342 : mpz_set(biT, aiJS[I]);
315 105659183 : for (J = 1; J < P; J++)
316 : {
317 : /* SubtractBigNbrModN(aiJS[I - J * PM], biT, aiJS[I - J * PM], TestNbr, NumberLength); */
318 80191841 : mpz_sub(aiJS[I - J * PM], aiJS[I - J * PM], biT);
319 : }
320 : /* aiJS[I] = 0; */
321 25467342 : mpz_set_ui(aiJS[I], 0);
322 : }
323 : }
324 136940207 : for (I = 0; I < PK; I++)
325 118796402 : mpz_mod(aiJS[I], aiJS[I], TestNbr);
326 18143805 : }
327 :
328 : /* ============================================================================================== */
329 :
330 : // Normalize coefficient of JW
331 167857 : void NormalizeJW(int PK, int PL, int PM, int P)
332 : {
333 : int I, J;
334 450114 : for (I = PL; I < PK; I++)
335 : {
336 282257 : if (mpz_cmp_ui(aiJW[I], 0) != 0) /* (!BigNbrIsZero(aiJW[I])) */
337 : {
338 : /* biT = aiJW[I]; */
339 170644 : mpz_set(biT, aiJW[I]);
340 :
341 926880 : for (J = 1; J < P; J++)
342 : {
343 : /* SubtractBigNbrModN(aiJW[I - J * PM], biT, aiJW[I - J * PM], TestNbr, NumberLength); */
344 756236 : mpz_sub(aiJW[I - J * PM], aiJW[I - J * PM], biT);
345 : }
346 : /* aiJW[I] = 0; */
347 170644 : mpz_set_ui(aiJW[I], 0);
348 : }
349 : }
350 1404971 : for (I = 0; I < PK; I++)
351 1237114 : mpz_mod(aiJW[I], aiJW[I], TestNbr);
352 167857 : }
353 :
354 : /* ============================================================================================== */
355 :
356 : // Perform JS <- JS * JW
357 :
358 5801557 : void JS_JW(int PK, int PL, int PM, int P)
359 : {
360 : int I, J, K;
361 34022339 : for (I = 0; I < PL; I++)
362 : {
363 231679974 : for (J = 0; J < PL; J++)
364 : {
365 203459192 : K = (I + J) % PK;
366 : /* MontgomeryMult(aiJS[I], aiJW[J], biTmp); */
367 : /* AddBigNbrModN(aiJX[K], biTmp, aiJX[K], TestNbr, NumberLength); */
368 203459192 : mpz_mul(biTmp, aiJS[I], aiJW[J]);
369 203459192 : mpz_add(aiJX[K], aiJX[K], biTmp);
370 : }
371 : }
372 44089135 : for (I = 0; I < PK; I++)
373 : {
374 : /* aiJS[I] = aiJX[I]; */
375 : /* aiJX[I] = 0; */
376 38287578 : mpz_swap(aiJS[I], aiJX[I]);
377 38287578 : mpz_set_ui(aiJX[I], 0);
378 : }
379 5801557 : NormalizeJS(PK, PL, PM, P);
380 5801557 : }
381 :
382 : /* ============================================================================================== */
383 :
384 : // Perform JS <- JS ^ 2
385 :
386 12332224 : void JS_2(int PK, int PL, int PM, int P)
387 : {
388 : int I, J, K;
389 70875605 : for (I = 0; I < PL; I++)
390 : {
391 58543381 : K = 2 * I % PK;
392 : /* MontgomeryMult(aiJS[I], aiJS[I], biTmp); */
393 : /* AddBigNbrModN(aiJX[K], biTmp, aiJX[K], TestNbr, NumberLength); */
394 : /* AddBigNbrModN(aiJS[I], aiJS[I], biT, TestNbr, NumberLength); */
395 58543381 : mpz_mul(biTmp, aiJS[I], aiJS[I]);
396 58543381 : mpz_add(aiJX[K], aiJX[K], biTmp);
397 58543381 : mpz_add(biT, aiJS[I], aiJS[I]);
398 236302927 : for (J = I + 1; J < PL; J++)
399 : {
400 177759546 : K = (I + J) % PK;
401 : /* MontgomeryMult(biT, aiJS[J], biTmp); */
402 : /* AddBigNbrModN(aiJX[K], biTmp, aiJX[K], TestNbr, NumberLength); */
403 177759546 : mpz_mul(biTmp, biT, aiJS[J]);
404 177759546 : mpz_add(aiJX[K], aiJX[K], biTmp);
405 : }
406 : }
407 92701928 : for (I = 0; I < PK; I++)
408 : {
409 : /* aiJS[I] = aiJX[I]; */
410 : /* aiJX[I] = 0; */
411 80369704 : mpz_swap(aiJS[I], aiJX[I]);
412 80369704 : mpz_set_ui(aiJX[I], 0);
413 : }
414 12332224 : NormalizeJS(PK, PL, PM, P);
415 12332224 : }
416 :
417 : /* ============================================================================================== */
418 :
419 : // Perform JS <- JS ^ E
420 :
421 167857 : void JS_E(int PK, int PL, int PM, int P)
422 : {
423 : int K;
424 : long Mask;
425 :
426 167857 : if (mpz_cmp_ui(biExp, 1) == 0)
427 : {
428 39008 : return;
429 : } // Return if E == 1
430 :
431 :
432 881176 : for (K = 0; K < PL; K++)
433 : {
434 752327 : mpz_set(aiJW[K], aiJS[K]);
435 : }
436 :
437 :
438 128849 : Mask = mpz_sizeinbase(biExp, 2)-1;
439 :
440 : do
441 : {
442 12324358 : JS_2(PK, PL, PM, P);
443 12324358 : Mask--;
444 12324358 : if (mpz_tstbit(biExp, Mask))
445 : {
446 5665619 : JS_JW(PK, PL, PM, P);
447 : }
448 : }
449 12324358 : while (Mask > 0);
450 : }
451 :
452 : /* ============================================================================================== */
453 :
454 : // if mode==0 then store J(p,q) ie x+f(x)
455 : // if mode==1 then p=2, look for p==1, and stores Jstar(q) ie 2x+f(x)
456 : // if mode==2 then p=2, look for p==4, and stores Jhash(q) ie 3x+f(x)
457 : // This is based on ideas and code from Jason Moxham
458 :
459 30224 : void JacobiSum(int mode, int P, int PL, int Q)
460 : {
461 : int I, a, myP;
462 :
463 159870 : for (I = 0; I < PL; I++)
464 129646 : mpz_set_ui(aiJ0[I], 0);
465 :
466 30224 : myP = P; /* if (mode == 0) */
467 30224 : if (mode == 1) myP = 1;
468 30224 : if (mode == 2) myP = 4;
469 :
470 1708232 : for(a = 0; a < JPQSMAX; a++)
471 1708232 : if(jpqs[a].p == myP && jpqs[a].q == Q)
472 30224 : break;
473 :
474 159870 : for (I = 0; I < PL; I++)
475 129646 : mpz_set_si(aiJ0[I], sls[jpqs[a].index+I]);
476 30224 : }
477 :
478 : /* ============================================================================================== */
479 :
480 : /* Prime checking routine */
481 : /* Return codes: 0 = N is composite. */
482 : /* 1 = N is a bpsw probable prime */
483 : /* 2 = N is prime. */
484 3230 : int mpz_aprtcle(mpz_t N, int verbose)
485 : {
486 : s64_t T, U;
487 : int i, j, H, I, J, K, P, Q, W, X;
488 : int IV, InvX, LEVELnow, NP, PK, PL, PM, SW, VK, TestedQs, TestingQs;
489 : int QQ, T1, T3, U1, U3, V1, V3;
490 3230 : int break_this = 0;
491 :
492 :
493 : /* make sure the input is >= 2 and odd */
494 3230 : if (mpz_cmp_ui(N, 2) < 0)
495 613 : return APRTCLE_COMPOSITE;
496 :
497 2617 : if (mpz_divisible_ui_p(N, 2))
498 : {
499 20 : if (mpz_cmp_ui(N, 2) == 0)
500 20 : return APRTCLE_PRIME;
501 : else
502 0 : return APRTCLE_COMPOSITE;
503 : }
504 :
505 : /* only three small exceptions for this implementation */
506 : /* with this set of P and Q primes */
507 2597 : if (mpz_cmp_ui(N, 3) == 0)
508 0 : return APRTCLE_PRIME;
509 2597 : if (mpz_cmp_ui(N, 7) == 0)
510 0 : return APRTCLE_PRIME;
511 2597 : if (mpz_cmp_ui(N, 11) == 0)
512 10 : return APRTCLE_PRIME;
513 :
514 : /* If the input number is larger than 7000 decimal digits
515 : we will just return whether it is a BPSW (probable) prime */
516 2587 : NumberLength = mpz_sizeinbase(N, 10);
517 2587 : if (NumberLength > 7000)
518 : {
519 0 : return mpz_probab_prime_p(N, 1); // mpz_probab_prime_p mpz_bpsw_prp
520 : }
521 :
522 2587 : allocate_vars();
523 :
524 2587 : mpz_set(TestNbr, N);
525 2587 : mpz_set_si(biS, 0);
526 :
527 2587 : j = PK = PL = PM = 0;
528 85371 : for (J = 0; J < PWmax; J++)
529 : {
530 : /* aiJX[J] = 0; */
531 82784 : mpz_set_ui(aiJX[J], 0);
532 : }
533 2587 : break_this = 0;
534 : /* GetPrimes2Test : */
535 4285 : for (i = 0; i < LEVELmax; i++)
536 : {
537 : /* biS[0] = 2; */
538 4285 : mpz_set_ui(biS, 2);
539 :
540 69950 : for (j = 0; j < aiNQ[i]; j++)
541 : {
542 68252 : Q = aiQ[j];
543 68252 : if (aiT[i]%(Q-1) != 0) continue;
544 68252 : U = aiT[i] * Q;
545 : do
546 : {
547 92468 : U /= Q;
548 : /* MultBigNbrByLong(biS, Q, biS, NumberLength); */
549 92468 : mpz_mul_ui(biS, biS, Q);
550 : }
551 92468 : while (U % Q == 0);
552 :
553 : // Exit loop if S^2 > N.
554 68252 : if (CompareSquare(biS, TestNbr) > 0)
555 : {
556 : /* break GetPrimes2Test; */
557 2587 : break_this = 1;
558 2587 : break;
559 : }
560 : } /* End for j */
561 4285 : if (break_this) break;
562 : } /* End for i */
563 2587 : if (i == LEVELmax)
564 : { /* too big */
565 0 : free_vars();
566 0 : return mpz_probab_prime_p(N, 1);
567 : }
568 2587 : LEVELnow = i;
569 2587 : TestingQs = j;
570 2587 : T = aiT[LEVELnow];
571 2587 : NP = aiNP[LEVELnow];
572 :
573 2683 : MainStart:
574 : for (;;)
575 : {
576 9941 : for (i = 0; i < NP; i++)
577 : {
578 7792 : P = aiP[i];
579 7792 : if (T%P != 0) continue;
580 :
581 7792 : SW = TestedQs = 0;
582 : /* Q = W = (int) BigNbrModLong(TestNbr, P * P); */
583 7792 : Q = W = mpz_fdiv_ui(TestNbr, P * P);
584 20903 : for (J = P - 2; J > 0; J--)
585 : {
586 13111 : W = (W * Q) % (P * P);
587 : }
588 7792 : if (P > 2 && W != 1)
589 : {
590 4077 : SW = 1;
591 : }
592 : for (;;)
593 : {
594 82977 : for (j = TestedQs; j <= TestingQs; j++)
595 : {
596 73185 : Q = aiQ[j] - 1;
597 : /* G = aiG[j]; */
598 73185 : K = 0;
599 121469 : while (Q % P == 0)
600 : {
601 48284 : K++;
602 48284 : Q /= P;
603 : }
604 73185 : Q = aiQ[j];
605 73185 : if (K == 0)
606 : {
607 37828 : continue;
608 : }
609 :
610 35357 : if (verbose >= APRTCLE_VERBOSE1)
611 : {
612 5000 : printf("APR primality test: P = %2d, Q = %12d (%3.2f%%)\r", P, Q, (i * (TestingQs + 1) + j) * 100.0 / (NP * (TestingQs + 1)));
613 5000 : fflush(stdout);
614 : }
615 :
616 35357 : PM = 1;
617 48284 : for (I = 1; I < K; I++)
618 : {
619 12927 : PM = PM * P;
620 : }
621 35357 : PL = (P - 1) * PM;
622 35357 : PK = P * PM;
623 203771 : for (I = 0; I < PK; I++)
624 : {
625 : /* aiJ0[I] = aiJ1[I] = 0; */
626 168414 : mpz_set_ui(aiJ0[I], 0);
627 168414 : mpz_set_ui(aiJ1[I], 0);
628 : }
629 35357 : if (P > 2)
630 : {
631 18124 : JacobiSum(0, P, PL, Q);
632 : }
633 : else
634 : {
635 17233 : if (K != 1)
636 : {
637 7866 : JacobiSum(0, P, PL, Q);
638 57142 : for (I = 0; I < PK; I++)
639 : {
640 : /* aiJW[I] = 0; */
641 49276 : mpz_set_ui(aiJW[I], 0);
642 : }
643 7866 : if (K != 2)
644 : {
645 15257 : for (I = 0; I < PM; I++)
646 : {
647 : /* aiJW[I] = aiJ0[I]; */
648 13140 : mpz_set(aiJW[I], aiJ0[I]);
649 : }
650 2117 : JacobiSum(1, P, PL, Q);
651 15257 : for (I = 0; I < PM; I++)
652 : {
653 : /* aiJS[I] = aiJ0[I]; */
654 13140 : mpz_set(aiJS[I], aiJ0[I]);
655 : }
656 2117 : JS_JW(PK, PL, PM, P);
657 15257 : for (I = 0; I < PM; I++)
658 : {
659 : /* aiJ1[I] = aiJS[I]; */
660 13140 : mpz_set(aiJ1[I], aiJS[I]);
661 : }
662 2117 : JacobiSum(2, P, PL, Q);
663 28397 : for (I = 0; I < PK; I++)
664 : {
665 : /* aiJW[I] = 0; */
666 26280 : mpz_set_ui(aiJW[I], 0);
667 : }
668 15257 : for (I = 0; I < PM; I++)
669 : {
670 : /* aiJS[I] = aiJ0[I]; */
671 13140 : mpz_set(aiJS[I], aiJ0[I]);
672 : }
673 2117 : JS_2(PK, PL, PM, P);
674 15257 : for (I = 0; I < PM; I++)
675 : {
676 : /* aiJ2[I] = aiJS[I]; */
677 13140 : mpz_set(aiJ2[I], aiJS[I]);
678 : }
679 : }
680 : }
681 : }
682 : /* aiJ00[0] = aiJ01[0] = 1; */
683 35357 : mpz_set_ui(aiJ00[0], 1);
684 35357 : mpz_set_ui(aiJ01[0], 1);
685 168414 : for (I = 1; I < PK; I++)
686 : {
687 : /* aiJ00[I] = aiJ01[I] = 0; */
688 133057 : mpz_set_ui(aiJ00[I], 0);
689 133057 : mpz_set_ui(aiJ01[I], 0);
690 : }
691 : /* VK = (int) BigNbrModLong(TestNbr, PK); */
692 35357 : VK = mpz_fdiv_ui(TestNbr, PK);
693 168414 : for (I = 1; I < PK; I++)
694 : {
695 133057 : if (I % P != 0)
696 : {
697 112733 : U1 = 1;
698 112733 : U3 = I;
699 112733 : V1 = 0;
700 112733 : V3 = PK;
701 452402 : while (V3 != 0)
702 : {
703 339669 : QQ = U3 / V3;
704 339669 : T1 = U1 - V1 * QQ;
705 339669 : T3 = U3 - V3 * QQ;
706 339669 : U1 = V1;
707 339669 : U3 = V3;
708 339669 : V1 = T1;
709 339669 : V3 = T3;
710 : }
711 112733 : aiInv[I] = (U1 + PK) % PK;
712 : }
713 : else
714 : {
715 20324 : aiInv[I] = 0;
716 : }
717 : }
718 35357 : if (P != 2)
719 : {
720 54372 : for (IV = 0; IV <= 1; IV++)
721 : {
722 200808 : for (X = 1; X < PK; X++)
723 : {
724 1356784 : for (I = 0; I < PK; I++)
725 : {
726 : /* aiJS[I] = aiJ0[I]; */
727 1192224 : mpz_set(aiJS[I], aiJ0[I]);
728 : }
729 164560 : if (X % P == 0)
730 : {
731 7104 : continue;
732 : }
733 157456 : if (IV == 0)
734 : {
735 : /* LongToBigNbr(X, biExp, NumberLength); */
736 78728 : mpz_set_ui(biExp, X);
737 : }
738 : else
739 : {
740 : /* LongToBigNbr(VK * X / PK, biExp, NumberLength); */
741 78728 : mpz_set_ui(biExp, (VK * X) / PK);
742 78728 : if ((VK * X) / PK == 0)
743 : {
744 34980 : continue;
745 : }
746 : }
747 122476 : JS_E(PK, PL, PM, P);
748 1052056 : for (I = 0; I < PK; I++)
749 : {
750 : /* aiJW[I] = 0; */
751 929580 : mpz_set_ui(aiJW[I], 0);
752 : }
753 122476 : InvX = aiInv[X];
754 1052056 : for (I = 0; I < PK; I++)
755 : {
756 929580 : J = (I * InvX) % PK;
757 : /* AddBigNbrModN(aiJW[J], aiJS[I], aiJW[J], TestNbr, NumberLength); */
758 929580 : mpz_add(aiJW[J], aiJW[J], aiJS[I]);
759 : }
760 122476 : NormalizeJW(PK, PL, PM, P);
761 122476 : if (IV == 0)
762 : {
763 642872 : for (I = 0; I < PK; I++)
764 : {
765 : /* aiJS[I] = aiJ00[I]; */
766 564144 : mpz_set(aiJS[I], aiJ00[I]);
767 : }
768 : }
769 : else
770 : {
771 409184 : for (I = 0; I < PK; I++)
772 : {
773 : /* aiJS[I] = aiJ01[I]; */
774 365436 : mpz_set(aiJS[I], aiJ01[I]);
775 : }
776 : }
777 122476 : JS_JW(PK, PL, PM, P);
778 122476 : if (IV == 0)
779 : {
780 642872 : for (I = 0; I < PK; I++)
781 : {
782 : /* aiJ00[I] = aiJS[I]; */
783 564144 : mpz_set(aiJ00[I], aiJS[I]);
784 : }
785 : }
786 : else
787 : {
788 409184 : for (I = 0; I < PK; I++)
789 : {
790 : /* aiJ01[I] = aiJS[I]; */
791 365436 : mpz_set(aiJ01[I], aiJS[I]);
792 : }
793 : }
794 : } /* end for X */
795 : } /* end for IV */
796 : }
797 : else
798 : {
799 17233 : if (K == 1)
800 : {
801 : /* MultBigNbrByLongModN(1, Q, aiJ00[0], TestNbr, NumberLength); */
802 9367 : mpz_set_ui(aiJ00[0], Q);
803 : /* aiJ01[0] = 1; */
804 9367 : mpz_set_ui(aiJ01[0], 1);
805 : }
806 : else
807 : {
808 7866 : if (K == 2)
809 : {
810 5749 : if (VK == 1)
811 : {
812 : /* aiJ01[0] = 1; */
813 2823 : mpz_set_ui(aiJ01[0], 1);
814 : }
815 : /* aiJS[0] = aiJ0[0]; */
816 : /* aiJS[1] = aiJ0[1]; */
817 5749 : mpz_set(aiJS[0], aiJ0[0]);
818 5749 : mpz_set(aiJS[1], aiJ0[1]);
819 5749 : JS_2(PK, PL, PM, P);
820 5749 : if (VK == 3)
821 : {
822 : /* aiJ01[0] = aiJS[0]; */
823 : /* aiJ01[1] = aiJS[1]; */
824 2926 : mpz_set(aiJ01[0], aiJS[0]);
825 2926 : mpz_set(aiJ01[1], aiJS[1]);
826 : }
827 : /* MultBigNbrByLongModN(aiJS[0], Q, aiJ00[0], TestNbr, NumberLength); */
828 5749 : mpz_mul_ui(aiJ00[0], aiJS[0], Q);
829 : /* MultBigNbrByLongModN(aiJS[1], Q, aiJ00[1], TestNbr, NumberLength); */
830 5749 : mpz_mul_ui(aiJ00[1], aiJS[1], Q);
831 : }
832 : else
833 : {
834 6351 : for (IV = 0; IV <= 1; IV++)
835 : {
836 30514 : for (X = 1; X < PK; X += 2)
837 : {
838 232432 : for (I = 0; I <= PM; I++)
839 : {
840 : /* aiJS[I] = aiJ1[I]; */
841 206152 : mpz_set(aiJS[I], aiJ1[I]);
842 : }
843 26280 : if (X % 8 == 5 || X % 8 == 7)
844 : {
845 13140 : continue;
846 : }
847 13140 : if (IV == 0)
848 : {
849 : /* LongToBigNbr(X, biExp, NumberLength); */
850 6570 : mpz_set_ui(biExp, X);
851 : }
852 : else
853 : {
854 : /* LongToBigNbr(VK * X / PK, biExp, NumberLength); */
855 6570 : mpz_set_ui(biExp, VK * X / PK);
856 6570 : if (VK * X / PK == 0)
857 : {
858 3116 : continue;
859 : }
860 : }
861 10024 : JS_E(PK, PL, PM, P);
862 149144 : for (I = 0; I < PK; I++)
863 : {
864 : /* aiJW[I] = 0; */
865 139120 : mpz_set_ui(aiJW[I], 0);
866 : }
867 10024 : InvX = aiInv[X];
868 149144 : for (I = 0; I < PK; I++)
869 : {
870 139120 : J = I * InvX % PK;
871 : /* AddBigNbrModN(aiJW[J], aiJS[I], aiJW[J], TestNbr, NumberLength); */
872 139120 : mpz_add(aiJW[J], aiJW[J], aiJS[I]);
873 : }
874 10024 : NormalizeJW(PK, PL, PM, P);
875 10024 : if (IV == 0)
876 : {
877 96506 : for (I = 0; I < PK; I++)
878 : {
879 : /* aiJS[I] = aiJ00[I]; */
880 89936 : mpz_set(aiJS[I], aiJ00[I]);
881 : }
882 : }
883 : else
884 : {
885 52638 : for (I = 0; I < PK; I++)
886 : {
887 : /* aiJS[I] = aiJ01[I]; */
888 49184 : mpz_set(aiJS[I], aiJ01[I]);
889 : }
890 : }
891 10024 : NormalizeJS(PK, PL, PM, P);
892 10024 : JS_JW(PK, PL, PM, P);
893 10024 : if (IV == 0)
894 : {
895 96506 : for (I = 0; I < PK; I++)
896 : {
897 : /* aiJ00[I] = aiJS[I]; */
898 89936 : mpz_set(aiJ00[I], aiJS[I]);
899 : }
900 : }
901 : else
902 : {
903 52638 : for (I = 0; I < PK; I++)
904 : {
905 : /* aiJ01[I] = aiJS[I]; */
906 49184 : mpz_set(aiJ01[I], aiJS[I]);
907 : }
908 : }
909 : } /* end for X */
910 4234 : if (IV == 0 || VK % 8 == 1 || VK % 8 == 3)
911 : {
912 2913 : continue;
913 : }
914 9505 : for (I = 0; I < PM; I++)
915 : {
916 : /* aiJW[I] = aiJ2[I]; */
917 : /* aiJS[I] = aiJ01[I]; */
918 8184 : mpz_set(aiJW[I], aiJ2[I]);
919 8184 : mpz_set(aiJS[I], aiJ01[I]);
920 : }
921 9505 : for (; I < PK; I++)
922 : {
923 : /* aiJW[I] = aiJS[I] = 0; */
924 8184 : mpz_set_ui(aiJW[I], 0);
925 8184 : mpz_set_ui(aiJS[I], 0);
926 : }
927 1321 : JS_JW(PK, PL, PM, P);
928 9505 : for (I = 0; I < PM; I++)
929 : {
930 : /* aiJ01[I] = aiJS[I]; */
931 8184 : mpz_set(aiJ01[I], aiJS[I]);
932 : }
933 : } /* end for IV */
934 : }
935 : }
936 : }
937 148090 : for (I = 0; I < PL; I++)
938 : {
939 : /* aiJS[I] = aiJ00[I]; */
940 112733 : mpz_set(aiJS[I], aiJ00[I]);
941 : }
942 91038 : for (; I < PK; I++)
943 : {
944 : /* aiJS[I] = 0; */
945 55681 : mpz_set_ui(aiJS[I], 0);
946 : }
947 : /* DivBigNbrByLong(TestNbr, PK, biExp, NumberLength); */
948 35357 : mpz_fdiv_q_ui(biExp, TestNbr, PK);
949 35357 : JS_E(PK, PL, PM, P);
950 203771 : for (I = 0; I < PK; I++)
951 : {
952 : /* aiJW[I] = 0; */
953 168414 : mpz_set_ui(aiJW[I], 0);
954 : }
955 148090 : for (I = 0; I < PL; I++)
956 : {
957 699136 : for (J = 0; J < PL; J++)
958 : {
959 : /* MontgomeryMult(aiJS[I], aiJ01[J], biTmp); */
960 : /* AddBigNbrModN(biTmp, aiJW[(I + J) % PK], aiJW[(I + J) % PK], TestNbr, NumberLength); */
961 586403 : mpz_mul(biTmp, aiJS[I], aiJ01[J]);
962 586403 : mpz_add(aiJW[(I + J) % PK], biTmp, aiJW[(I + J) % PK]);
963 : }
964 : }
965 35357 : NormalizeJW(PK, PL, PM, P);
966 : /* MatchingRoot : */
967 : do
968 : {
969 35357 : H = -1;
970 35357 : W = 0;
971 148090 : for (I = 0; I < PL; I++)
972 : {
973 112733 : if (mpz_cmp_ui(aiJW[I], 0) != 0)/* (!BigNbrIsZero(aiJW[I])) */
974 : {
975 : /* if (H == -1 && BigNbrAreEqual(aiJW[I], 1)) */
976 45706 : if (H == -1 && (mpz_cmp_ui(aiJW[I], 1) == 0))
977 : {
978 21355 : H = I;
979 : }
980 : else
981 : {
982 24351 : H = -2;
983 : /* AddBigNbrModN(aiJW[I], MontgomeryMultR1, biTmp, TestNbr, NumberLength); */
984 24351 : mpz_add_ui(biTmp, aiJW[I], 1);
985 24351 : mpz_mod(biTmp, biTmp, TestNbr);
986 24351 : if (mpz_cmp_ui(biTmp, 0) == 0) /* (BigNbrIsZero(biTmp)) */
987 : {
988 23923 : W++;
989 : }
990 : }
991 : }
992 : }
993 35357 : if (H >= 0)
994 : {
995 : /* break MatchingRoot; */
996 21355 : break;
997 : }
998 14002 : if (W != P - 1)
999 : {
1000 : /* Not prime */
1001 438 : free_vars();
1002 438 : return APRTCLE_COMPOSITE;
1003 : }
1004 18831 : for (I = 0; I < PM; I++)
1005 : {
1006 : /* AddBigNbrModN(aiJW[I], 1, biTmp, TestNbr, NumberLength); */
1007 18831 : mpz_add_ui(biTmp, aiJW[I], 1);
1008 18831 : mpz_mod(biTmp, biTmp, TestNbr);
1009 18831 : if (mpz_cmp_ui(biTmp, 0) == 0) /* (BigNbrIsZero(biTmp)) */
1010 : {
1011 13564 : break;
1012 : }
1013 : }
1014 13564 : if (I == PM)
1015 : {
1016 : /* Not prime */
1017 0 : free_vars();
1018 0 : return APRTCLE_COMPOSITE;
1019 : }
1020 23923 : for (J = 1; J <= P - 2; J++)
1021 : {
1022 : /* AddBigNbrModN(aiJW[I + J * PM], 1, biTmp, TestNbr, NumberLength); */
1023 10359 : mpz_add_ui(biTmp, aiJW[I + J * PM], 1);
1024 10359 : mpz_mod(biTmp, biTmp, TestNbr);
1025 10359 : if (mpz_cmp_ui(biTmp, 0) != 0)/* (!BigNbrIsZero(biTmp)) */
1026 : {
1027 : /* Not prime */
1028 0 : free_vars();
1029 0 : return APRTCLE_COMPOSITE;
1030 : }
1031 : }
1032 13564 : H = I + PL;
1033 : }
1034 : while (0);
1035 :
1036 34919 : if (SW == 1 || H % P == 0)
1037 : {
1038 30230 : continue;
1039 : }
1040 4689 : if (P != 2)
1041 : {
1042 1022 : SW = 1;
1043 1022 : continue;
1044 : }
1045 3667 : if (K == 1)
1046 : {
1047 2207 : if ((mpz_get_ui(TestNbr) & 3) == 1)
1048 : {
1049 699 : SW = 1;
1050 : }
1051 2207 : continue;
1052 : }
1053 :
1054 : // if (Q^((N-1)/2) mod N != N-1), N is not prime.
1055 :
1056 : /* MultBigNbrByLongModN(1, Q, biTmp, TestNbr, NumberLength); */
1057 1460 : mpz_set_ui(biTmp, Q);
1058 1460 : mpz_mod(biTmp, biTmp, TestNbr);
1059 :
1060 1460 : mpz_sub_ui(biT, TestNbr, 1); /* biT = n-1 */
1061 1460 : mpz_divexact_ui(biT, biT, 2); /* biT = (n-1)/2 */
1062 1460 : mpz_powm(biR, biTmp, biT, TestNbr); /* biR = Q^((n-1)/2) mod n */
1063 1460 : mpz_add_ui(biTmp, biR, 1);
1064 1460 : mpz_mod(biTmp, biTmp, TestNbr);
1065 :
1066 1460 : if (mpz_cmp_ui(biTmp, 0) != 0)/* (!BigNbrIsZero(biTmp)) */
1067 : {
1068 : /* Not prime */
1069 0 : free_vars();
1070 0 : return APRTCLE_COMPOSITE;
1071 : }
1072 1460 : SW = 1;
1073 : } /* end for j */
1074 9792 : if (SW == 0)
1075 : {
1076 2534 : TestedQs = TestingQs + 1;
1077 2534 : if (TestingQs < aiNQ[LEVELnow] - 1)
1078 : {
1079 2438 : TestingQs++;
1080 2438 : Q = aiQ[TestingQs];
1081 2438 : U = T * Q;
1082 : do
1083 : {
1084 : /* MultBigNbrByLong(biS, Q, biS, NumberLength); */
1085 3128 : mpz_mul_ui(biS, biS, Q);
1086 3128 : U /= Q;
1087 : }
1088 3128 : while (U % Q == 0);
1089 :
1090 2438 : continue; /* Retry */
1091 : }
1092 96 : LEVELnow++;
1093 96 : if (LEVELnow == LEVELmax)
1094 : {
1095 0 : free_vars();
1096 0 : return mpz_probab_prime_p(N, 1); /* Cannot tell */
1097 : }
1098 96 : T = aiT[LEVELnow];
1099 96 : NP = aiNP[LEVELnow];
1100 : /* biS = 2; */
1101 96 : mpz_set_ui(biS, 2);
1102 384 : for (J = 0; J <= aiNQ[LEVELnow]; J++)
1103 : {
1104 384 : Q = aiQ[J];
1105 384 : if (T%(Q-1) != 0) continue;
1106 384 : U = T * Q;
1107 : do
1108 : {
1109 : /* MultBigNbrByLong(biS, Q, biS, NumberLength); */
1110 1072 : mpz_mul_ui(biS, biS, Q);
1111 1072 : U /= Q;
1112 : }
1113 1072 : while (U % Q == 0);
1114 384 : if (CompareSquare(biS, TestNbr) > 0)
1115 : {
1116 96 : TestingQs = J;
1117 : /* continue MainStart; */ /* Retry from the beginning */
1118 96 : goto MainStart;
1119 : }
1120 : } /* end for J */
1121 0 : free_vars();
1122 0 : return mpz_probab_prime_p(N, 1); /* Program error */
1123 : } /* end if */
1124 7258 : break;
1125 : } /* end for (;;) */
1126 : } /* end for i */
1127 :
1128 : // Final Test
1129 :
1130 : /* biR = 1 */
1131 2149 : mpz_set_ui(biR, 1);
1132 : /* biN <- TestNbr mod biS */ /* Compute N mod S */
1133 2149 : mpz_fdiv_r(biN, TestNbr, biS);
1134 :
1135 23891100 : for (U = 1; U <= T; U++)
1136 : {
1137 : /* biR <- (biN * biR) mod biS */
1138 23891100 : mpz_mul(biR, biN, biR);
1139 23891100 : mpz_mod(biR, biR, biS);
1140 23891100 : if (mpz_cmp_ui(biR, 1) == 0) /* biR == 1 */
1141 : {
1142 : /* Number is prime */
1143 2149 : free_vars();
1144 2149 : return APRTCLE_PRIME;
1145 : }
1146 23888951 : if (mpz_divisible_p(TestNbr, biR) && mpz_cmp(biR, TestNbr) < 0) /* biR < N and biR | TestNbr */
1147 : {
1148 : /* Number is composite */
1149 0 : free_vars();
1150 0 : return APRTCLE_COMPOSITE;
1151 : }
1152 : } /* End for U */
1153 : /* This should never be reached. */
1154 0 : free_vars();
1155 0 : return mpz_probab_prime_p(N, 1);
1156 : }
1157 : }
1158 :
1159 : /* ============================================================================================== */
|
