LCOV - code coverage report
Current view: top level - ecm/aprtcle - mpz_aprcl.c (source / functions) Hit Total Coverage
Test: unnamed Lines: 423 443 95.5 %
Date: 2019-09-21 10:55:40 Functions: 10 10 100.0 %

          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        2557 : void allocate_vars()
     218             : {
     219        2557 :   int i = 0;
     220        2557 :   aiJS = malloc(PWmax * sizeof(mpz_t));
     221        2557 :   aiJW = malloc(PWmax * sizeof(mpz_t));
     222        2557 :   aiJX = malloc(PWmax * sizeof(mpz_t));
     223        2557 :   aiJ0 = malloc(PWmax * sizeof(mpz_t));
     224        2557 :   aiJ1 = malloc(PWmax * sizeof(mpz_t));
     225        2557 :   aiJ2 = malloc(PWmax * sizeof(mpz_t));
     226        2557 :   aiJ00 = malloc(PWmax * sizeof(mpz_t));
     227        2557 :   aiJ01 = malloc(PWmax * sizeof(mpz_t));
     228       84381 :   for (i = 0 ; i < PWmax; i++)
     229             :   {
     230       81824 :     mpz_init(aiJS[i]);
     231       81824 :     mpz_init(aiJW[i]);
     232       81824 :     mpz_init(aiJX[i]);
     233       81824 :     mpz_init(aiJ0[i]);
     234       81824 :     mpz_init(aiJ1[i]);
     235       81824 :     mpz_init(aiJ2[i]);
     236       81824 :     mpz_init(aiJ00[i]);
     237       81824 :     mpz_init(aiJ01[i]);
     238             :   }
     239             : 
     240        2557 :   mpz_init(TestNbr);
     241        2557 :   mpz_init(biN);
     242        2557 :   mpz_init(biR);
     243        2557 :   mpz_init(biS);
     244        2557 :   mpz_init(biT);
     245        2557 :   mpz_init(biExp);
     246        2557 :   mpz_init(biTmp);
     247        2557 : }
     248             : 
     249             : /* ============================================================================================== */
     250             : 
     251        2557 : void free_vars()
     252             : {
     253        2557 :   int i = 0;
     254       84381 :   for (i = 0 ; i < PWmax; i++)
     255             :   {
     256       81824 :     mpz_clear(aiJS[i]);
     257       81824 :     mpz_clear(aiJW[i]);
     258       81824 :     mpz_clear(aiJX[i]);
     259       81824 :     mpz_clear(aiJ0[i]);
     260       81824 :     mpz_clear(aiJ1[i]);
     261       81824 :     mpz_clear(aiJ2[i]);
     262       81824 :     mpz_clear(aiJ00[i]);
     263       81824 :     mpz_clear(aiJ01[i]);
     264             :   }
     265        2557 :   free(aiJS);
     266        2557 :   free(aiJW);
     267        2557 :   free(aiJX);
     268        2557 :   free(aiJ0);
     269        2557 :   free(aiJ1);
     270        2557 :   free(aiJ2);
     271        2557 :   free(aiJ00);
     272        2557 :   free(aiJ01);
     273             : 
     274        2557 :   mpz_clear(TestNbr);
     275        2557 :   mpz_clear(biN);
     276        2557 :   mpz_clear(biR);
     277        2557 :   mpz_clear(biS);
     278        2557 :   mpz_clear(biT);
     279        2557 :   mpz_clear(biExp);
     280        2557 :   mpz_clear(biTmp);
     281        2557 : }
     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       67666 : int CompareSquare(mpz_t Nbr1, mpz_t Nbr2)
     290             : {
     291             :   mpz_t tmp;
     292       67666 :   int cmp = 0;
     293             : 
     294       67666 :   mpz_init(tmp);
     295       67666 :   mpz_mul(tmp, Nbr1, Nbr1);
     296             : 
     297       67666 :   cmp = mpz_cmp(tmp, Nbr2);
     298       67666 :   mpz_clear(tmp);
     299             : 
     300       67666 :   return cmp;
     301             : }
     302             : 
     303             : /* ============================================================================================== */
     304             : 
     305             : // Normalize coefficient of JS
     306    17590240 : void NormalizeJS(int PK, int PL, int PM, int P)
     307             : {
     308             :   int I, J;
     309    48495860 :   for (I = PL; I < PK; I++)
     310             :   {
     311    30905620 :     if (mpz_cmp_ui(aiJS[I], 0) != 0) /* (!BigNbrIsZero(aiJS[I])) */
     312             :     {
     313             :       /* biT = aiJS[I]; */
     314    24633060 :       mpz_set(biT, aiJS[I]);
     315   102775420 :       for (J = 1; J < P; J++)
     316             :       {
     317             :         /* SubtractBigNbrModN(aiJS[I - J * PM], biT, aiJS[I - J * PM], TestNbr, NumberLength); */
     318    78142360 :         mpz_sub(aiJS[I - J * PM], aiJS[I - J * PM], biT);
     319             :       }
     320             :       /* aiJS[I] = 0; */
     321    24633060 :       mpz_set_ui(aiJS[I], 0);
     322             :     }
     323             :   }
     324   133052960 :   for (I = 0; I < PK; I++)
     325   115462720 :     mpz_mod(aiJS[I], aiJS[I], TestNbr);
     326    17590240 : }
     327             : 
     328             : /* ============================================================================================== */
     329             : 
     330             : // Normalize coefficient of JW
     331      161607 : void NormalizeJW(int PK, int PL, int PM, int P)
     332             : {
     333             :   int I, J;
     334      431364 :   for (I = PL; I < PK; I++)
     335             :   {
     336      269757 :     if (mpz_cmp_ui(aiJW[I], 0) != 0) /* (!BigNbrIsZero(aiJW[I])) */
     337             :     {
     338             :       /* biT = aiJW[I]; */
     339      163194 :       mpz_set(biT, aiJW[I]);
     340             : 
     341      894090 :       for (J = 1; J < P; J++)
     342             :       {
     343             :         /* SubtractBigNbrModN(aiJW[I - J * PM], biT, aiJW[I - J * PM], TestNbr, NumberLength); */
     344      730896 :         mpz_sub(aiJW[I - J * PM], aiJW[I - J * PM], biT);
     345             :       }
     346             :       /* aiJW[I] = 0; */
     347      163194 :       mpz_set_ui(aiJW[I], 0);
     348             :     }
     349             :   }
     350     1351831 :   for (I = 0; I < PK; I++)
     351     1190224 :     mpz_mod(aiJW[I], aiJW[I], TestNbr);
     352      161607 : }
     353             : 
     354             : /* ============================================================================================== */
     355             : 
     356             : // Perform JS <- JS * JW
     357             : 
     358     5615658 : void JS_JW(int PK, int PL, int PM, int P)
     359             : {
     360             :   int I, J, K;
     361    33061400 :   for (I = 0; I < PL; I++)
     362             :   {
     363   226423400 :     for (J = 0; J < PL; J++)
     364             :     {
     365   198978000 :       K = (I + J) % PK;
     366             :       /* MontgomeryMult(aiJS[I], aiJW[J], biTmp); */
     367             :       /* AddBigNbrModN(aiJX[K], biTmp, aiJX[K], TestNbr, NumberLength); */
     368   198978000 :       mpz_mul(biTmp, aiJS[I], aiJW[J]);
     369   198978000 :       mpz_add(aiJX[K], aiJX[K], biTmp);
     370             :     }
     371             :   }
     372    42772180 :   for (I = 0; I < PK; I++)
     373             :   {
     374             :     /* aiJS[I] = aiJX[I]; */
     375             :     /* aiJX[I] = 0; */
     376    37156580 :     mpz_swap(aiJS[I], aiJX[I]);
     377    37156580 :     mpz_set_ui(aiJX[I], 0);
     378             :   }
     379     5615658 :   NormalizeJS(PK, PL, PM, P);
     380     5615658 : }
     381             : 
     382             : /* ============================================================================================== */
     383             : 
     384             : // Perform JS <- JS ^ 2
     385             : 
     386    11965094 : void JS_2(int PK, int PL, int PM, int P)
     387             : {
     388             :   int I, J, K;
     389    69010360 :   for (I = 0; I < PL; I++)
     390             :   {
     391    57045280 :     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    57045280 :     mpz_mul(biTmp, aiJS[I], aiJS[I]);
     396    57045280 :     mpz_add(aiJX[K], aiJX[K], biTmp);
     397    57045280 :     mpz_add(biT, aiJS[I], aiJS[I]);
     398   231277800 :     for (J = I + 1; J < PL; J++)
     399             :     {
     400   174232200 :       K = (I + J) % PK;
     401             :       /* MontgomeryMult(biT, aiJS[J], biTmp); */
     402             :       /* AddBigNbrModN(aiJX[K], biTmp, aiJX[K], TestNbr, NumberLength); */
     403   174232200 :       mpz_mul(biTmp, biT, aiJS[J]);
     404   174232200 :       mpz_add(aiJX[K], aiJX[K], biTmp);
     405             :     }
     406             :   }
     407    90139280 :   for (I = 0; I < PK; I++)
     408             :   {
     409             :     /* aiJS[I] = aiJX[I]; */
     410             :     /* aiJX[I] = 0; */
     411    78174200 :     mpz_swap(aiJS[I], aiJX[I]);
     412    78174200 :     mpz_set_ui(aiJX[I], 0);
     413             :   }
     414    11965094 :   NormalizeJS(PK, PL, PM, P);
     415    11965094 : }
     416             : 
     417             : /* ============================================================================================== */
     418             : 
     419             : // Perform JS <- JS ^ E
     420             : 
     421      161607 : void JS_E(int PK, int PL, int PM, int P)
     422             : {
     423             :   int K;
     424             :   long Mask;
     425             : 
     426      161607 :   if (mpz_cmp_ui(biExp, 1) == 0)
     427             :   {
     428       37648 :     return;
     429             :   } // Return if E == 1
     430             : 
     431             : 
     432      849256 :   for (K = 0; K < PL; K++)
     433             :   {
     434      725297 :     mpz_set(aiJW[K], aiJS[K]);
     435             :   }
     436             : 
     437             : 
     438      123959 :   Mask = mpz_sizeinbase(biExp, 2)-1;
     439             : 
     440             :   do
     441             :   {
     442    11957516 :     JS_2(PK, PL, PM, P);
     443    11957516 :     Mask--;
     444    11957516 :     if (mpz_tstbit(biExp, Mask))
     445             :     {
     446     5484850 :       JS_JW(PK, PL, PM, P);
     447             :     }
     448             :   }
     449    11957516 :   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       29044 : void JacobiSum(int mode, int P, int PL, int Q)
     460             : {
     461             :   int I, a, myP;
     462             : 
     463      153210 :   for (I = 0; I < PL; I++)
     464      124166 :     mpz_set_ui(aiJ0[I], 0);
     465             : 
     466       29044 :   myP = P; /* if (mode == 0) */
     467       29044 :   if (mode == 1) myP = 1;
     468       29044 :   if (mode == 2) myP = 4;
     469             : 
     470     1662714 :   for(a = 0; a < JPQSMAX; a++)
     471     1662714 :     if(jpqs[a].p == myP && jpqs[a].q == Q)
     472       29044 :       break;
     473             : 
     474      153210 :   for (I = 0; I < PL; I++)
     475      124166 :     mpz_set_si(aiJ0[I], sls[jpqs[a].index+I]);
     476       29044 : }
     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        3201 : 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        3201 :   int break_this = 0;
     491             : 
     492             : 
     493             :   /* make sure the input is >= 2 and odd */
     494        3201 :   if (mpz_cmp_ui(N, 2) < 0)
     495         614 :     return APRTCLE_COMPOSITE;
     496             : 
     497        2587 :   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        2567 :   if (mpz_cmp_ui(N, 3) == 0)
     508           0 :     return APRTCLE_PRIME;
     509        2567 :   if (mpz_cmp_ui(N, 7) == 0)
     510           0 :     return APRTCLE_PRIME;
     511        2567 :   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        2557 :   NumberLength = mpz_sizeinbase(N, 10);
     517        2557 :   if (NumberLength > 7000)
     518             :   {
     519           0 :     return mpz_probab_prime_p(N, 1); // mpz_probab_prime_p mpz_bpsw_prp
     520             :   }
     521             : 
     522        2557 :   allocate_vars();
     523             : 
     524        2557 :   mpz_set(TestNbr, N);
     525        2557 :   mpz_set_si(biS, 0);
     526             : 
     527        2557 :   j = PK = PL = PM = 0;
     528       84381 :   for (J = 0; J < PWmax; J++)
     529             :   {
     530             :     /* aiJX[J] = 0; */
     531       81824 :     mpz_set_ui(aiJX[J], 0);
     532             :   }
     533        2557 :   break_this = 0;
     534             : /* GetPrimes2Test : */
     535        4225 :   for (i = 0; i < LEVELmax; i++)
     536             :   {
     537             :     /* biS[0] = 2; */
     538        4225 :     mpz_set_ui(biS, 2);
     539             :     
     540       69010 :     for (j = 0; j < aiNQ[i]; j++)
     541             :     {
     542       67342 :       Q = aiQ[j];
     543       67342 :       if (aiT[i]%(Q-1) != 0) continue;
     544       67342 :       U = aiT[i] * Q;
     545             :       do
     546             :       {
     547       91188 :         U /= Q;
     548             :         /* MultBigNbrByLong(biS, Q, biS, NumberLength); */
     549       91188 :         mpz_mul_ui(biS, biS, Q);
     550             :       }
     551       91188 :       while (U % Q == 0);
     552             : 
     553             :       // Exit loop if S^2 > N.
     554       67342 :       if (CompareSquare(biS, TestNbr) > 0)
     555             :       {
     556             :         /* break GetPrimes2Test; */
     557        2557 :         break_this = 1;
     558        2557 :         break;
     559             :       }
     560             :     } /* End for j */
     561        4225 :     if (break_this) break;
     562             :   } /* End for i */
     563        2557 :   if (i == LEVELmax)
     564             :   { /* too big */
     565           0 :     free_vars();
     566           0 :     return mpz_probab_prime_p(N, 1);
     567             :   }
     568        2557 :   LEVELnow = i;
     569        2557 :   TestingQs = j;
     570        2557 :   T = aiT[LEVELnow];
     571        2557 :   NP = aiNP[LEVELnow];
     572             : 
     573        2643 : MainStart:
     574             :   for (;;)
     575             :   {
     576        9771 :     for (i = 0; i < NP; i++)
     577             :     {
     578        7652 :       P = aiP[i];
     579        7652 :       if (T%P != 0) continue;
     580             : 
     581        7652 :       SW = TestedQs = 0;
     582             :       /* Q = W = (int) BigNbrModLong(TestNbr, P * P); */
     583        7652 :       Q = W = mpz_fdiv_ui(TestNbr, P * P);
     584       20403 :       for (J = P - 2; J > 0; J--)
     585             :       {
     586       12751 :         W = (W * Q) % (P * P);
     587             :       }
     588        7652 :       if (P > 2 && W != 1)
     589             :       {
     590        4007 :         SW = 1;
     591             :       }
     592             :       for (;;)
     593             :       {
     594       80297 :         for (j = TestedQs; j <= TestingQs; j++)
     595             :         {
     596       70705 :           Q = aiQ[j] - 1;
     597             :           /* G = aiG[j]; */
     598       70705 :           K = 0;
     599      117189 :           while (Q % P == 0)
     600             :           {
     601       46484 :             K++;
     602       46484 :             Q /= P;
     603             :           }
     604       70705 :           Q = aiQ[j];
     605       70705 :           if (K == 0)
     606             :           {
     607       36578 :             continue;
     608             :           }
     609             : 
     610       34127 :           if (verbose >= APRTCLE_VERBOSE1)
     611             :           {
     612        5000 :             printf("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       34127 :           PM = 1;
     617       46484 :           for (I = 1; I < K; I++)
     618             :           {
     619       12357 :             PM = PM * P;
     620             :           }
     621       34127 :           PL = (P - 1) * PM;
     622       34127 :           PK = P * PM;
     623      195951 :           for (I = 0; I < PK; I++)
     624             :           {
     625             :             /* aiJ0[I] = aiJ1[I] = 0; */
     626      161824 :             mpz_set_ui(aiJ0[I], 0);
     627      161824 :             mpz_set_ui(aiJ1[I], 0);
     628             :           }
     629       34127 :           if (P > 2)
     630             :           {
     631       17454 :             JacobiSum(0, P, PL, Q);
     632             :           }
     633             :           else
     634             :           {
     635       16673 :             if (K != 1)
     636             :             {
     637        7576 :               JacobiSum(0, P, PL, Q);
     638       54772 :               for (I = 0; I < PK; I++)
     639             :               {
     640             :                 /* aiJW[I] = 0; */
     641       47196 :                 mpz_set_ui(aiJW[I], 0);
     642             :               }
     643        7576 :               if (K != 2)
     644             :               {
     645       14467 :                 for (I = 0; I < PM; I++)
     646             :                 {
     647             :                   /* aiJW[I] = aiJ0[I]; */
     648       12460 :                   mpz_set(aiJW[I], aiJ0[I]);
     649             :                 }
     650        2007 :                 JacobiSum(1, P, PL, Q);
     651       14467 :                 for (I = 0; I < PM; I++)
     652             :                 {
     653             :                   /* aiJS[I] = aiJ0[I]; */
     654       12460 :                   mpz_set(aiJS[I], aiJ0[I]);
     655             :                 }
     656        2007 :                 JS_JW(PK, PL, PM, P);
     657       14467 :                 for (I = 0; I < PM; I++)
     658             :                 {
     659             :                   /* aiJ1[I] = aiJS[I]; */
     660       12460 :                   mpz_set(aiJ1[I], aiJS[I]);
     661             :                 }
     662        2007 :                 JacobiSum(2, P, PL, Q);
     663       26927 :                 for (I = 0; I < PK; I++)
     664             :                 {
     665             :                   /* aiJW[I] = 0; */
     666       24920 :                   mpz_set_ui(aiJW[I], 0);
     667             :                 }
     668       14467 :                 for (I = 0; I < PM; I++)
     669             :                 {
     670             :                   /* aiJS[I] = aiJ0[I]; */
     671       12460 :                   mpz_set(aiJS[I], aiJ0[I]);
     672             :                 }
     673        2007 :                 JS_2(PK, PL, PM, P);
     674       14467 :                 for (I = 0; I < PM; I++)
     675             :                 {
     676             :                   /* aiJ2[I] = aiJS[I]; */
     677       12460 :                   mpz_set(aiJ2[I], aiJS[I]);
     678             :                 }
     679             :               }
     680             :             }
     681             :           }
     682             :           /* aiJ00[0] = aiJ01[0] = 1; */
     683       34127 :           mpz_set_ui(aiJ00[0], 1);
     684       34127 :           mpz_set_ui(aiJ01[0], 1);
     685      161824 :           for (I = 1; I < PK; I++)
     686             :           {
     687             :             /* aiJ00[I] = aiJ01[I] = 0; */
     688      127697 :             mpz_set_ui(aiJ00[I], 0);
     689      127697 :             mpz_set_ui(aiJ01[I], 0);
     690             :           }
     691             :           /* VK = (int) BigNbrModLong(TestNbr, PK); */
     692       34127 :           VK = mpz_fdiv_ui(TestNbr, PK);
     693      161824 :           for (I = 1; I < PK; I++)
     694             :           {
     695      127697 :             if (I % P != 0)
     696             :             {
     697      108343 :               U1 = 1;
     698      108343 :               U3 = I;
     699      108343 :               V1 = 0;
     700      108343 :               V3 = PK;
     701      434552 :               while (V3 != 0)
     702             :               {
     703      326209 :                 QQ = U3 / V3;
     704      326209 :                 T1 = U1 - V1 * QQ;
     705      326209 :                 T3 = U3 - V3 * QQ;
     706      326209 :                 U1 = V1;
     707      326209 :                 U3 = V3;
     708      326209 :                 V1 = T1;
     709      326209 :                 V3 = T3;
     710             :               }
     711      108343 :               aiInv[I] = (U1 + PK) % PK;
     712             :             }
     713             :             else
     714             :             {
     715       19354 :               aiInv[I] = 0;
     716             :             }
     717             :           }
     718       34127 :           if (P != 2)
     719             :           {
     720       52362 :             for (IV = 0; IV <= 1; IV++)
     721             :             {
     722      192868 :               for (X = 1; X < PK; X++)
     723             :               {
     724     1302544 :                 for (I = 0; I < PK; I++)
     725             :                 {
     726             :                   /* aiJS[I] = aiJ0[I]; */
     727     1144584 :                   mpz_set(aiJS[I], aiJ0[I]);
     728             :                 }
     729      157960 :                 if (X % P == 0)
     730             :                 {
     731        6664 :                   continue;
     732             :                 }
     733      151296 :                 if (IV == 0)
     734             :                 {
     735             :                   /* LongToBigNbr(X, biExp, NumberLength); */
     736       75648 :                   mpz_set_ui(biExp, X);
     737             :                 }
     738             :                 else
     739             :                 {
     740             :                   /* LongToBigNbr(VK * X / PK, biExp, NumberLength); */
     741       75648 :                   mpz_set_ui(biExp, (VK * X) / PK);
     742       75648 :                   if ((VK * X) / PK == 0)
     743             :                   {
     744       33320 :                     continue;
     745             :                   }
     746             :                 }
     747      117976 :                 JS_E(PK, PL, PM, P);
     748     1014376 :                 for (I = 0; I < PK; I++)
     749             :                 {
     750             :                   /* aiJW[I] = 0; */
     751      896400 :                   mpz_set_ui(aiJW[I], 0);
     752             :                 }
     753      117976 :                 InvX = aiInv[X];
     754     1014376 :                 for (I = 0; I < PK; I++)
     755             :                 {
     756      896400 :                   J = (I * InvX) % PK;
     757             :                   /* AddBigNbrModN(aiJW[J], aiJS[I], aiJW[J], TestNbr, NumberLength); */
     758      896400 :                   mpz_add(aiJW[J], aiJW[J], aiJS[I]);
     759             :                 }
     760      117976 :                 NormalizeJW(PK, PL, PM, P);
     761      117976 :                 if (IV == 0)
     762             :                 {
     763      617952 :                   for (I = 0; I < PK; I++)
     764             :                   {
     765             :                     /* aiJS[I] = aiJ00[I]; */
     766      542304 :                     mpz_set(aiJS[I], aiJ00[I]);
     767             :                   }
     768             :                 }
     769             :                 else
     770             :                 {
     771      396424 :                   for (I = 0; I < PK; I++)
     772             :                   {
     773             :                     /* aiJS[I] = aiJ01[I]; */
     774      354096 :                     mpz_set(aiJS[I], aiJ01[I]);
     775             :                   }
     776             :                 }
     777      117976 :                 JS_JW(PK, PL, PM, P);
     778      117976 :                 if (IV == 0)
     779             :                 {
     780      617952 :                   for (I = 0; I < PK; I++)
     781             :                   {
     782             :                     /* aiJ00[I] = aiJS[I]; */
     783      542304 :                     mpz_set(aiJ00[I], aiJS[I]);
     784             :                   }
     785             :                 }
     786             :                 else
     787             :                 {
     788      396424 :                   for (I = 0; I < PK; I++)
     789             :                   {
     790             :                     /* aiJ01[I] = aiJS[I]; */
     791      354096 :                     mpz_set(aiJ01[I], aiJS[I]);
     792             :                   }
     793             :                 }
     794             :               } /* end for X */
     795             :             } /* end for IV */
     796             :           }
     797             :           else
     798             :           {
     799       16673 :             if (K == 1)
     800             :             {
     801             :               /* MultBigNbrByLongModN(1, Q, aiJ00[0], TestNbr, NumberLength); */
     802        9097 :               mpz_set_ui(aiJ00[0], Q);
     803             :               /* aiJ01[0] = 1; */
     804        9097 :               mpz_set_ui(aiJ01[0], 1);
     805             :             }
     806             :             else
     807             :             {
     808        7576 :               if (K == 2)
     809             :               {
     810        5569 :                 if (VK == 1)
     811             :                 {
     812             :                   /* aiJ01[0] = 1; */
     813        2793 :                   mpz_set_ui(aiJ01[0], 1);
     814             :                 }
     815             :                 /* aiJS[0] = aiJ0[0]; */
     816             :                 /* aiJS[1] = aiJ0[1]; */
     817        5569 :                 mpz_set(aiJS[0], aiJ0[0]);
     818        5569 :                 mpz_set(aiJS[1], aiJ0[1]);
     819        5569 :                 JS_2(PK, PL, PM, P);
     820        5569 :                 if (VK == 3)
     821             :                 {
     822             :                   /* aiJ01[0] = aiJS[0]; */
     823             :                   /* aiJ01[1] = aiJS[1]; */
     824        2776 :                   mpz_set(aiJ01[0], aiJS[0]);
     825        2776 :                   mpz_set(aiJ01[1], aiJS[1]);
     826             :                 }
     827             :                 /* MultBigNbrByLongModN(aiJS[0], Q, aiJ00[0], TestNbr, NumberLength); */
     828        5569 :                 mpz_mul_ui(aiJ00[0], aiJS[0], Q);
     829             :                 /* MultBigNbrByLongModN(aiJS[1], Q, aiJ00[1], TestNbr, NumberLength); */
     830        5569 :                 mpz_mul_ui(aiJ00[1], aiJS[1], Q);
     831             :               }
     832             :               else
     833             :               {
     834        6021 :                 for (IV = 0; IV <= 1; IV++)
     835             :                 {
     836       28934 :                   for (X = 1; X < PK; X += 2)
     837             :                   {
     838      220432 :                     for (I = 0; I <= PM; I++)
     839             :                     {
     840             :                       /* aiJS[I] = aiJ1[I]; */
     841      195512 :                       mpz_set(aiJS[I], aiJ1[I]);
     842             :                     }
     843       24920 :                     if (X % 8 == 5 || X % 8 == 7)
     844             :                     {
     845       12460 :                       continue;
     846             :                     }
     847       12460 :                     if (IV == 0)
     848             :                     {
     849             :                       /* LongToBigNbr(X, biExp, NumberLength); */
     850        6230 :                       mpz_set_ui(biExp, X);
     851             :                     }
     852             :                     else
     853             :                     {
     854             :                       /* LongToBigNbr(VK * X / PK, biExp, NumberLength); */
     855        6230 :                       mpz_set_ui(biExp, VK * X / PK);
     856        6230 :                       if (VK * X / PK == 0)
     857             :                       {
     858        2956 :                         continue;
     859             :                       }
     860             :                     }
     861        9504 :                     JS_E(PK, PL, PM, P);
     862      141504 :                     for (I = 0; I < PK; I++)
     863             :                     {
     864             :                       /* aiJW[I] = 0; */
     865      132000 :                       mpz_set_ui(aiJW[I], 0);
     866             :                     }
     867        9504 :                     InvX = aiInv[X];
     868      141504 :                     for (I = 0; I < PK; I++)
     869             :                     {
     870      132000 :                       J = I * InvX % PK;
     871             :                       /* AddBigNbrModN(aiJW[J], aiJS[I], aiJW[J], TestNbr, NumberLength); */
     872      132000 :                       mpz_add(aiJW[J], aiJW[J], aiJS[I]);
     873             :                     }
     874        9504 :                     NormalizeJW(PK, PL, PM, P);
     875        9504 :                     if (IV == 0)
     876             :                     {
     877       91526 :                       for (I = 0; I < PK; I++)
     878             :                       {
     879             :                         /* aiJS[I] = aiJ00[I]; */
     880       85296 :                         mpz_set(aiJS[I], aiJ00[I]);
     881             :                       }
     882             :                     }
     883             :                     else
     884             :                     {
     885       49978 :                       for (I = 0; I < PK; I++)
     886             :                       {
     887             :                         /* aiJS[I] = aiJ01[I]; */
     888       46704 :                         mpz_set(aiJS[I], aiJ01[I]);
     889             :                       }
     890             :                     }
     891        9504 :                     NormalizeJS(PK, PL, PM, P);
     892        9504 :                     JS_JW(PK, PL, PM, P);
     893        9504 :                     if (IV == 0)
     894             :                     {
     895       91526 :                       for (I = 0; I < PK; I++)
     896             :                       {
     897             :                         /* aiJ00[I] = aiJS[I]; */
     898       85296 :                         mpz_set(aiJ00[I], aiJS[I]);
     899             :                       }
     900             :                     }
     901             :                     else
     902             :                     {
     903       49978 :                       for (I = 0; I < PK; I++)
     904             :                       {
     905             :                         /* aiJ01[I] = aiJS[I]; */
     906       46704 :                         mpz_set(aiJ01[I], aiJS[I]);
     907             :                       }
     908             :                     }
     909             :                   } /* end for X */
     910        4014 :                   if (IV == 0 || VK % 8 == 1 || VK % 8 == 3)
     911             :                   {
     912        2693 :                     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      142470 :           for (I = 0; I < PL; I++)
     938             :           {
     939             :             /* aiJS[I] = aiJ00[I]; */
     940      108343 :             mpz_set(aiJS[I], aiJ00[I]);
     941             :           }
     942       87608 :           for (; I < PK; I++)
     943             :           {
     944             :             /* aiJS[I] = 0; */
     945       53481 :             mpz_set_ui(aiJS[I], 0);
     946             :           }
     947             :           /* DivBigNbrByLong(TestNbr, PK, biExp, NumberLength); */
     948       34127 :           mpz_fdiv_q_ui(biExp, TestNbr, PK);
     949       34127 :           JS_E(PK, PL, PM, P);
     950      195951 :           for (I = 0; I < PK; I++)
     951             :           {
     952             :             /* aiJW[I] = 0; */
     953      161824 :             mpz_set_ui(aiJW[I], 0);
     954             :           }
     955      142470 :           for (I = 0; I < PL; I++)
     956             :           {
     957      671676 :             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      563333 :               mpz_mul(biTmp, aiJS[I], aiJ01[J]);
     962      563333 :               mpz_add(aiJW[(I + J) % PK], biTmp, aiJW[(I + J) % PK]);
     963             :             }
     964             :           }
     965       34127 :           NormalizeJW(PK, PL, PM, P);
     966             : /* MatchingRoot : */
     967             :           do
     968             :           {
     969       34127 :             H = -1;
     970       34127 :             W = 0;
     971      142470 :             for (I = 0; I < PL; I++)
     972             :             {
     973      108343 :               if (mpz_cmp_ui(aiJW[I], 0) != 0)/* (!BigNbrIsZero(aiJW[I])) */
     974             :               {
     975             :                 /* if (H == -1 && BigNbrAreEqual(aiJW[I], 1)) */
     976       43856 :                 if (H == -1 && (mpz_cmp_ui(aiJW[I], 1) == 0))
     977             :                 {
     978       20645 :                   H = I;
     979             :                 }
     980             :                 else
     981             :                 {
     982       23211 :                   H = -2;
     983             :                   /* AddBigNbrModN(aiJW[I], MontgomeryMultR1, biTmp, TestNbr, NumberLength); */
     984       23211 :                   mpz_add_ui(biTmp, aiJW[I], 1);
     985       23211 :                   mpz_mod(biTmp, biTmp, TestNbr);
     986       23211 :                   if (mpz_cmp_ui(biTmp, 0) == 0) /* (BigNbrIsZero(biTmp)) */
     987             :                   {
     988       22783 :                     W++;
     989             :                   }
     990             :                 }
     991             :               }
     992             :             }
     993       34127 :             if (H >= 0)
     994             :             {
     995             :               /* break MatchingRoot; */
     996       20645 :               break;
     997             :             }
     998       13482 :             if (W != P - 1)
     999             :             {
    1000             :               /* Not prime */
    1001         438 :               free_vars();
    1002         438 :               return APRTCLE_COMPOSITE;
    1003             :             }
    1004       18081 :             for (I = 0; I < PM; I++)
    1005             :             {
    1006             :               /* AddBigNbrModN(aiJW[I], 1, biTmp, TestNbr, NumberLength); */
    1007       18081 :               mpz_add_ui(biTmp, aiJW[I], 1);
    1008       18081 :               mpz_mod(biTmp, biTmp, TestNbr);
    1009       18081 :               if (mpz_cmp_ui(biTmp, 0) == 0) /* (BigNbrIsZero(biTmp)) */
    1010             :               {
    1011       13044 :                 break;
    1012             :               }
    1013             :             }
    1014       13044 :             if (I == PM)
    1015             :             {
    1016             :               /* Not prime */
    1017           0 :               free_vars();
    1018           0 :               return APRTCLE_COMPOSITE;
    1019             :             }
    1020       22783 :             for (J = 1; J <= P - 2; J++)
    1021             :             {
    1022             :               /* AddBigNbrModN(aiJW[I + J * PM], 1, biTmp, TestNbr, NumberLength); */
    1023        9739 :               mpz_add_ui(biTmp, aiJW[I + J * PM], 1);
    1024        9739 :               mpz_mod(biTmp, biTmp, TestNbr);
    1025        9739 :               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       13044 :             H = I + PL;
    1033             :           }
    1034             :           while (0);
    1035             : 
    1036       33689 :           if (SW == 1 || H % P == 0)
    1037             :           {
    1038       29120 :             continue;
    1039             :           }
    1040        4569 :           if (P != 2)
    1041             :           {
    1042         992 :             SW = 1;
    1043         992 :             continue;
    1044             :           }
    1045        3577 :           if (K == 1)
    1046             :           {
    1047        2137 :             if ((mpz_get_ui(TestNbr) & 3) == 1)
    1048             :             {
    1049         689 :               SW = 1;
    1050             :             }
    1051        2137 :             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        1440 :           mpz_set_ui(biTmp, Q);
    1058        1440 :           mpz_mod(biTmp, biTmp, TestNbr);
    1059             : 
    1060        1440 :           mpz_sub_ui(biT, TestNbr, 1); /* biT = n-1 */
    1061        1440 :           mpz_divexact_ui(biT, biT, 2); /* biT = (n-1)/2 */
    1062        1440 :           mpz_powm(biR, biTmp, biT, TestNbr); /* biR = Q^((n-1)/2) mod n */
    1063        1440 :           mpz_add_ui(biTmp, biR, 1);
    1064        1440 :           mpz_mod(biTmp, biTmp, TestNbr);
    1065             : 
    1066        1440 :           if (mpz_cmp_ui(biTmp, 0) != 0)/* (!BigNbrIsZero(biTmp)) */
    1067             :           {
    1068             :             /* Not prime */
    1069           0 :             free_vars();
    1070           0 :             return APRTCLE_COMPOSITE;
    1071             :           }
    1072        1440 :           SW = 1;
    1073             :         } /* end for j */
    1074        9592 :         if (SW == 0)
    1075             :         {
    1076        2464 :           TestedQs = TestingQs + 1;
    1077        2464 :           if (TestingQs < aiNQ[LEVELnow] - 1)
    1078             :           {
    1079        2378 :             TestingQs++;
    1080        2378 :             Q = aiQ[TestingQs];
    1081        2378 :             U = T * Q;
    1082             :             do
    1083             :             {
    1084             :               /* MultBigNbrByLong(biS, Q, biS, NumberLength); */
    1085        3068 :               mpz_mul_ui(biS, biS, Q);
    1086        3068 :               U /= Q;
    1087             :             }
    1088        3068 :             while (U % Q == 0);
    1089             : 
    1090        2378 :             continue; /* Retry */
    1091             :           }
    1092          86 :           LEVELnow++;
    1093          86 :           if (LEVELnow == LEVELmax)
    1094             :           {
    1095           0 :             free_vars();
    1096           0 :             return mpz_probab_prime_p(N, 1); /* Cannot tell */
    1097             :           }
    1098          86 :           T = aiT[LEVELnow];
    1099          86 :           NP = aiNP[LEVELnow];
    1100             :           /* biS = 2; */
    1101          86 :           mpz_set_ui(biS, 2);
    1102         324 :           for (J = 0; J <= aiNQ[LEVELnow]; J++)
    1103             :           {
    1104         324 :             Q = aiQ[J];
    1105         324 :             if (T%(Q-1) != 0) continue;
    1106         324 :             U = T * Q;
    1107             :             do
    1108             :             {
    1109             :               /* MultBigNbrByLong(biS, Q, biS, NumberLength); */
    1110         932 :               mpz_mul_ui(biS, biS, Q);
    1111         932 :               U /= Q;
    1112             :             }
    1113         932 :             while (U % Q == 0);
    1114         324 :             if (CompareSquare(biS, TestNbr) > 0)
    1115             :             {
    1116          86 :               TestingQs = J;
    1117             :               /* continue MainStart; */ /* Retry from the beginning */
    1118          86 :               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        7128 :         break;
    1125             :       } /* end for (;;) */
    1126             :     } /* end for i */
    1127             : 
    1128             :     // Final Test
    1129             :    
    1130             :     /* biR = 1 */
    1131        2119 :     mpz_set_ui(biR, 1);
    1132             :     /* biN <- TestNbr mod biS */ /* Compute N mod S */
    1133        2119 :     mpz_fdiv_r(biN, TestNbr, biS);
    1134             :    
    1135    23277900 :     for (U = 1; U <= T; U++)
    1136             :     {
    1137             :       /* biR <- (biN * biR) mod biS */
    1138    23277900 :       mpz_mul(biR, biN, biR);
    1139    23277900 :       mpz_mod(biR, biR, biS);
    1140    23277900 :       if (mpz_cmp_ui(biR, 1) == 0) /* biR == 1 */
    1141             :       {
    1142             :         /* Number is prime */
    1143        2119 :         free_vars();
    1144        2119 :         return APRTCLE_PRIME;
    1145             :       }
    1146    23275760 :       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             : /* ============================================================================================== */

Generated by: LCOV version 1.14