text: /* Copyright 2011-2015 David Cleaver
   *
   * This program is free software: you can redistribute it and/or modify
   * it under the terms of the GNU Lesser General Public License as published by
   * the Free Software Foundation, either version 3 of the License, or
   * (at your option) any later version.
   *
   * This program is distributed in the hope that it will be useful,
   * but WITHOUT ANY WARRANTY; without even the implied warranty of
   * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   * GNU Lesser General Public License for more details.
   *
   * You should have received a copy of the GNU Lesser General Public License
   * along with this program.  If not, see <http://www.gnu.org/licenses/>.
   */
  
  /*
   * v1.0 Posted to SourceForge on 2013/07/04
   *
   * v1.1 Posted to SourceForge on 2013/12/27
   * [The following improvements/fixes were recommended by Laurent Desnogues in 2013/08]
   *      - Speed improvement 1: Removed extraneous NormalizeJS calls in ARPCL
   *      - Speed improvement 2: Removed/consolidated calls to mpz_mod in APRCL
   *        (these improvements make the APRCL code about 1.5-2.2x faster)
   *      - Bug fix: Final test in APRCL routine is now correct
   */
  
  
  #include <stdio.h>
  #include <stdlib.h>
  #include <gmp.h>
  #include "mpz_aprcl.h"
  
  #ifndef HAVE_U64_T
  #define HAVE_U64_T
  typedef long long s64_t;
  typedef unsigned long long u64_t;
  #endif
  
  
  /* **********************************************************************************
   * APR-CL (also known as APRT-CLE) is a prime proving algorithm developed by:
   * L. Adleman, C. Pomerance, R. Rumely, H. Cohen, and H. W. Lenstra
   * APRT-CLE = Adleman-Pomerance-Rumely Test Cohen-Lenstra Extended version
   * You can find all the details of this implementation in the Cohen & Lenstra paper:
   *    H. Cohen and A. K. Lenstra, "Implementation of a new primality test",
   *    Math. Comp., 48 (1987) 103--121
   *
   * ----------------------------------------------------------------------------------
   *
   * This C/GMP version is a conversion of Dario Alpern's Java based APRT-CLE code
   * His code was based on Yuji Kida's UBASIC code
   *
   * Based on APRT-CLE Written by Dario Alejandro Alpern (Buenos Aires - Argentina)
   * From: Last updated September 10th, 2011. See http://www.alpertron.com.ar/ECM.HTM
   *
   * On 2012/11/12 Dario Alpern has approved the conversion, from Java to C/GMP, of 
   * his implementation of the APR-CL algorithm, and that it be licensed under the LGPL.
   *
   * ----------------------------------------------------------------------------------
   *
   * With improvements based on Jason Moxham's APRCL v1.15 code, from 2003/01/01
   *
   * On 2013/04/14 Toby Moxham has approved the APR-CL code and data tables, 
   * originally written by his brother Jason Moxham on 2003/01/01, to be released 
   * under the LGPL.
   *
   * *********************************************************************************/
  
  /* int PWmax = 32, Qmax = 55441, LEVELmax = 11; */
  /* QMax is the largest Q in the aiQ array */
  /* PWmax is the largest p^k that divides any (Q-1) */
  /* #define Qmax 55441 */
  /* #define PWmax 32 */
  /*
   2-max =  2^5 =  32 from t[ 4]=    4324320
   3-max =  3^3 =  27 from t[ 4]=    4324320
   5-max =  5^2 =  25 from t[ 6]=  367567200
   7-max =  7^1 =   7 from t[ 1]=       5040
  11-max = 11^1 =  11 from t[ 2]=      55440
  13-max = 13^1 =  13 from t[ 3]=     720720
  17-max = 17^1 =  17 from t[ 5]=   73513440
  19-max = 19^1 =  19 from t[ 7]= 1396755360
  */
  
  /* largest value of p^k that divides any t value */
  #define PWmax 32
  
  /* Qmax[]: list of largest q-prime for each t value */
  int Qmax[] = {61,2521,55441,180181,4324321,10501921,367567201,232792561,
  1745944201};
  
  /* number of values in Qmax[], aiNP[], aiNQ[] */
  int LEVELmax = 9;
  
  /* list of primes that divide our t values */
  int aiP[] = {2,3,5,7,11,13,17,19};
  
  /* list of q-primes from all t values */
  int aiQ[] = {2,3,5,7,11,13,31,61,17,19,29,37,41,43,71,73,113,127,181,211,241,
  281,337,421,631,1009,2521,23,67,89,199,331,397,463,617,661,881,991,1321,2311,
  3697,4621,9241,18481,55441,53,79,131,157,313,521,547,859,911,937,1093,1171,
  1873,2003,2341,2731,2861,3121,3433,6007,6553,8009,8191,8581,16381,20021,20593,
  21841,25741,36037,48049,51481,65521,72073,120121,180181,97,109,271,353,379,433,
  541,673,757,1249,2017,2081,2161,2377,2971,3169,3361,3511,4159,5281,7393,7561,
  7723,8317,8737,9829,13729,14561,15121,16633,23761,24571,26209,28081,30241,
  38611,39313,47521,66529,96097,108109,110881,123553,131041,196561,216217,270271,
  332641,393121,432433,540541,617761,4324321,103,137,239,307,409,443,613,919,953,
  1021,1123,1327,1361,1429,1531,1871,2143,2381,2857,3061,3571,3673,4421,4591,
  5237,6121,6427,6733,7481,8161,9181,9283,9521,10099,10711,12241,12377,12853,
  14281,15913,16831,17137,17681,19891,22441,23563,23869,24481,27847,29173,29921,
  30941,34273,36721,42841,43759,46411,47737,52361,53857,59671,63649,70687,72931,
  74257,78541,79561,87517,92821,97241,100981,102103,116689,117811,128521,145861,
  148513,157081,161569,167077,185641,201961,209441,235621,238681,269281,291721,
  314161,371281,388961,417691,445537,471241,477361,514081,565489,612613,656371,
  680681,700129,816817,1633633,1670761,1837837,2625481,4084081,5250961,5654881,
  8168161,9189181,10501921,101,151,401,601,701,1051,1201,1301,1801,1951,2551,
  2801,3301,3851,4201,4951,5101,5851,6301,7151,9901,11551,11701,12601,14851,
  15401,15601,17551,17851,18701,19801,21601,23801,28051,33151,34651,40801,42901,
  44201,50051,53551,54601,56101,66301,70201,77351,79201,81901,91801,92401,93601,
  103951,107101,109201,118801,122401,140401,150151,151201,160651,193051,198901,
  200201,218401,224401,232051,243101,257401,300301,321301,367201,415801,428401,
  448801,450451,504901,530401,600601,673201,729301,795601,800801,982801,1029601,
  1093951,1178101,1201201,1458601,2088451,2187901,2402401,2570401,2702701,
  3088801,3141601,3712801,5105101,5834401,6806801,7068601,8353801,17503201,
  22972951,52509601,183783601,367567201,191,229,419,457,571,647,761,1483,1597,
  2053,2129,2281,2927,3041,3877,4447,4523,4561,4789,6271,6689,6841,6917,7411,
  7753,8209,8779,8893,10337,11287,11971,12541,13339,13567,13681,14821,16417,
  17291,17443,18089,19381,20521,20749,21319,21737,22573,25841,27361,28729,29641,
  30097,31123,35531,35569,35911,38039,39521,40699,43891,46817,47881,48907,51871,
  53353,56431,57457,58787,59281,63841,71821,72353,75583,77521,87211,90289,97813,
  105337,106591,108529,114913,117041,124489,131671,134369,135661,139537,140449,
  146719,163021,177841,186733,207481,213181,217057,217361,225721,251941,279073,
  287281,300961,302329,342343,351121,377911,391249,406981,451441,456457,461891,
  489061,511633,526681,554269,568481,608609,651169,652081,697681,733591,782497,
  790021,813961,895357,1027027,1053361,1058149,1108537,1133731,1264033,1279081,
  1369369,1492261,1580041,1790713,1813969,1867321,1939939,2217073,2238391,
  2282281,2351441,2489761,2645371,2771341,2934361,2984521,3233231,3627937,
  3837241,3912481,3979361,4157011,4232593,4476781,5135131,5372137,5868721,
  6046561,6348889,6651217,6715171,6846841,7162849,7674481,9767521,11737441,
  12471031,12697777,17907121,24942061,27387361,31744441,35814241,41081041,
  46558513,53721361,107442721,174594421,232792561,1901,2851,5701,39901,41801,
  53201,62701,64601,74101,79801,98801,113051,119701,135851,148201,205201,219451,
  290701,292601,319201,333451,339151,359101,410401,452201,478801,501601,532951,
  564301,658351,666901,778051,839801,957601,1037401,1065901,1128601,1222651,
  1259701,1504801,1808801,1889551,2074801,2173601,2445301,2667601,3052351,
  3511201,3730651,3779101,3950101,4069801,4149601,4408951,5038801,6104701,
  6224401,8558551,9781201,11191951,11411401,14922601,16279201,17117101,17635801,
  19186201,19562401,22383901,22822801,23514401,25581601,25675651,31600801,
  35271601,37346401,38372401,45349201,59690401,67151701,83140201,129329201,
  134303401,193993801,249420601,436486051,634888801,1163962801,1745944201};
  /* number of q primes in the above array: 618 */
  
  /* a primitive root for each q-prime in the above array */
  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,
  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,
  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,
  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,
  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,
  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,
  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,
  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,
  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,
  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,
  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,
  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,
  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,
  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,
  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,
  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,
  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,
  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,
  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,
  46,23};
  
  
  /* number of primes from aiP, not necessarily in order, that divides each t */
  int aiNP[] = {3,4,5,6,6,7,7,8,8};
  
  /* number of q-primes for each t */
  int aiNQ[] = {8,27,45,81,134,245,351,424,618};
  
  /*         t     |       e(t)   | #Qp |      Qmax  |   divisors of t         */
  /* --------------|--------------|-----|------------|------------------------ */
  s64_t aiT[] =  {
            60, /* | 6.8144 E   9 |   8 |         61 | p={2,3,5}               */
          5040, /* | 1.5321 E  52 |  27 |       2521 | p={2,3,5,7}             */
         55440, /* | 4.9209 E 106 |  45 |      55441 | p={2,3,5,7,11}          */
        720720, /* | 2.5992 E 237 |  81 |     180181 | p={2,3,5,7,11,13}       */
       4324320, /* | 7.9285 E 455 | 134 |    4324321 | p={2,3,5,7,11,13}       */
      73513440, /* | 7.0821 E 966 | 245 |   10501921 | p={2,3,5,7,11,13,17}    */
     367567200, /* | 6.2087 E1501 | 351 |  367567201 | p={2,3,5,7,11,13,17}    */
    1396755360, /* | 4.0165 E1913 | 424 |  232792561 | p={2,3,5,7,11,13,17,19} */
    6983776800};/* | 7.4712 E3010 | 618 | 1745944201 | p={2,3,5,7,11,13,17,19} */
  
  
  int aiInv[PWmax];
  mpz_t biTmp;
  mpz_t biExp;
  mpz_t biN;
  mpz_t biR;
  mpz_t biS;
  mpz_t biT;
  mpz_t *aiJS; /* [PWmax] */
  mpz_t *aiJW; /* [PWmax] */
  mpz_t *aiJX; /* [PWmax] */
  mpz_t *aiJ0; /* [PWmax] */
  mpz_t *aiJ1; /* [PWmax] */
  mpz_t *aiJ2; /* [PWmax] */
  mpz_t *aiJ00; /* [PWmax] */
  mpz_t *aiJ01; /* [PWmax] */
  int NumberLength; /* Length of multiple precision nbrs */
  mpz_t TestNbr;
  
  /* ============================================================================================== */
  
  void allocate_vars()
  {
    int i = 0;
    aiJS = malloc(PWmax * sizeof(mpz_t));
    aiJW = malloc(PWmax * sizeof(mpz_t));
    aiJX = malloc(PWmax * sizeof(mpz_t));
    aiJ0 = malloc(PWmax * sizeof(mpz_t));
    aiJ1 = malloc(PWmax * sizeof(mpz_t));
    aiJ2 = malloc(PWmax * sizeof(mpz_t));
    aiJ00 = malloc(PWmax * sizeof(mpz_t));
    aiJ01 = malloc(PWmax * sizeof(mpz_t));
    for (i = 0 ; i < PWmax; i++)
    {
      mpz_init(aiJS[i]);
      mpz_init(aiJW[i]);
      mpz_init(aiJX[i]);
      mpz_init(aiJ0[i]);
      mpz_init(aiJ1[i]);
      mpz_init(aiJ2[i]);
      mpz_init(aiJ00[i]);
      mpz_init(aiJ01[i]);
    }
  
    mpz_init(TestNbr);
    mpz_init(biN);
    mpz_init(biR);
    mpz_init(biS);
    mpz_init(biT);
    mpz_init(biExp);
    mpz_init(biTmp);
  }
  
  /* ============================================================================================== */
  
  void free_vars()
  {
    int i = 0;
    for (i = 0 ; i < PWmax; i++)
    {
      mpz_clear(aiJS[i]);
      mpz_clear(aiJW[i]);
      mpz_clear(aiJX[i]);
      mpz_clear(aiJ0[i]);
      mpz_clear(aiJ1[i]);
      mpz_clear(aiJ2[i]);
      mpz_clear(aiJ00[i]);
      mpz_clear(aiJ01[i]);
    }
    free(aiJS);
    free(aiJW);
    free(aiJX);
    free(aiJ0);
    free(aiJ1);
    free(aiJ2);
    free(aiJ00);
    free(aiJ01);
  
    mpz_clear(TestNbr);
    mpz_clear(biN);
    mpz_clear(biR);
    mpz_clear(biS);
    mpz_clear(biT);
    mpz_clear(biExp);
    mpz_clear(biTmp);
  }
  
  /* ============================================================================================== */
  
  // Compare Nbr1^2 vs. Nbr2
  // -1 -> Nbr1^2 < Nbr2
  //  0 -> Nbr1^2 == Nbr2
  //  1 -> Nbr1^2 > Nbr2
  int CompareSquare(mpz_t Nbr1, mpz_t Nbr2)
  {
    mpz_t tmp;
    int cmp = 0;
  
    mpz_init(tmp);
    mpz_mul(tmp, Nbr1, Nbr1);
  
    cmp = mpz_cmp(tmp, Nbr2);
    mpz_clear(tmp);
  
    return cmp;
  }
  
  /* ============================================================================================== */
  
  // Normalize coefficient of JS
  void NormalizeJS(int PK, int PL, int PM, int P)
  {
    int I, J;
    for (I = PL; I < PK; I++)
    {
      if (mpz_cmp_ui(aiJS[I], 0) != 0) /* (!BigNbrIsZero(aiJS[I])) */
      {
        /* biT = aiJS[I]; */
        mpz_set(biT, aiJS[I]);
        for (J = 1; J < P; J++)
        {
          /* SubtractBigNbrModN(aiJS[I - J * PM], biT, aiJS[I - J * PM], TestNbr, NumberLength); */
          mpz_sub(aiJS[I - J * PM], aiJS[I - J * PM], biT);
        }
        /* aiJS[I] = 0; */
        mpz_set_ui(aiJS[I], 0);
      }
    }
    for (I = 0; I < PK; I++)
      mpz_mod(aiJS[I], aiJS[I], TestNbr);
  }
  
  /* ============================================================================================== */
  
  // Normalize coefficient of JW
  void NormalizeJW(int PK, int PL, int PM, int P)
  {
    int I, J;
    for (I = PL; I < PK; I++)
    {
      if (mpz_cmp_ui(aiJW[I], 0) != 0) /* (!BigNbrIsZero(aiJW[I])) */
      {
        /* biT = aiJW[I]; */
        mpz_set(biT, aiJW[I]);
  
        for (J = 1; J < P; J++)
        {
          /* SubtractBigNbrModN(aiJW[I - J * PM], biT, aiJW[I - J * PM], TestNbr, NumberLength); */
          mpz_sub(aiJW[I - J * PM], aiJW[I - J * PM], biT);
        }
        /* aiJW[I] = 0; */
        mpz_set_ui(aiJW[I], 0);
      }
    }
    for (I = 0; I < PK; I++)
      mpz_mod(aiJW[I], aiJW[I], TestNbr);
  }
  
  /* ============================================================================================== */
  
  // Perform JS <- JS * JW
  
  void JS_JW(int PK, int PL, int PM, int P)
  {
    int I, J, K;
    for (I = 0; I < PL; I++)
    {
      for (J = 0; J < PL; J++)
      {
        K = (I + J) % PK;
        /* MontgomeryMult(aiJS[I], aiJW[J], biTmp); */
        /* AddBigNbrModN(aiJX[K], biTmp, aiJX[K], TestNbr, NumberLength); */
        mpz_mul(biTmp, aiJS[I], aiJW[J]);
        mpz_add(aiJX[K], aiJX[K], biTmp);
      }
    }
    for (I = 0; I < PK; I++)
    {
      /* aiJS[I] = aiJX[I]; */
      /* aiJX[I] = 0; */
      mpz_swap(aiJS[I], aiJX[I]);
      mpz_set_ui(aiJX[I], 0);
    }
    NormalizeJS(PK, PL, PM, P);
  }
  
  /* ============================================================================================== */
  
  // Perform JS <- JS ^ 2
  
  void JS_2(int PK, int PL, int PM, int P)
  {
    int I, J, K;
    for (I = 0; I < PL; I++)
    {
      K = 2 * I % PK;
      /* MontgomeryMult(aiJS[I], aiJS[I], biTmp); */
      /* AddBigNbrModN(aiJX[K], biTmp, aiJX[K], TestNbr, NumberLength); */
      /* AddBigNbrModN(aiJS[I], aiJS[I], biT, TestNbr, NumberLength); */
      mpz_mul(biTmp, aiJS[I], aiJS[I]);
      mpz_add(aiJX[K], aiJX[K], biTmp);
      mpz_add(biT, aiJS[I], aiJS[I]);
      for (J = I + 1; J < PL; J++)
      {
        K = (I + J) % PK;
        /* MontgomeryMult(biT, aiJS[J], biTmp); */
        /* AddBigNbrModN(aiJX[K], biTmp, aiJX[K], TestNbr, NumberLength); */
        mpz_mul(biTmp, biT, aiJS[J]);
        mpz_add(aiJX[K], aiJX[K], biTmp);
      }
    }
    for (I = 0; I < PK; I++)
    {
      /* aiJS[I] = aiJX[I]; */
      /* aiJX[I] = 0; */
      mpz_swap(aiJS[I], aiJX[I]);
      mpz_set_ui(aiJX[I], 0);
    }
    NormalizeJS(PK, PL, PM, P);
  }
  
  /* ============================================================================================== */
  
  // Perform JS <- JS ^ E
  
  void JS_E(int PK, int PL, int PM, int P)
  {
    int K;
    long Mask;
  
    if (mpz_cmp_ui(biExp, 1) == 0)
    {
      return;
    } // Return if E == 1
  
  
    for (K = 0; K < PL; K++)
    {
      mpz_set(aiJW[K], aiJS[K]);
    }
  
  
    Mask = mpz_sizeinbase(biExp, 2)-1;
  
    do
    {
      JS_2(PK, PL, PM, P);
      Mask--;
      if (mpz_tstbit(biExp, Mask))
      {
        JS_JW(PK, PL, PM, P);
      }
    }
    while (Mask > 0);
  }
  
  /* ============================================================================================== */
  
  // if mode==0 then store J(p,q) ie x+f(x)
  // if mode==1 then p=2, look for p==1, and stores Jstar(q) ie 2x+f(x)
  // if mode==2 then p=2, look for p==4, and stores Jhash(q) ie 3x+f(x)
  // This is based on ideas and code from Jason Moxham
  
  void JacobiSum(int mode, int P, int PL, int Q)
  {
    int I, a, myP;
  
    for (I = 0; I < PL; I++)
      mpz_set_ui(aiJ0[I], 0);
  
    myP = P; /* if (mode == 0) */
    if (mode == 1) myP = 1;
    if (mode == 2) myP = 4;
  
    for(a = 0; a < JPQSMAX; a++)
      if(jpqs[a].p == myP && jpqs[a].q == Q)
        break;
  
    for (I = 0; I < PL; I++)
      mpz_set_si(aiJ0[I], sls[jpqs[a].index+I]);
  }
  
  /* ============================================================================================== */
  
  /* Prime checking routine                       */
  /* Return codes: 0 = N is composite.            */
  /*               1 = N is a bpsw probable prime */
  /*               2 = N is prime.                */
  int mpz_aprtcle(mpz_t N, int verbose)
  {
    s64_t T, U;
    int i, j, H, I, J, K, P, Q, W, X;
    int IV, InvX, LEVELnow, NP, PK, PL, PM, SW, VK, TestedQs, TestingQs;
    int QQ, T1, T3, U1, U3, V1, V3;
    int break_this = 0;
  
  
    /* make sure the input is >= 2 and odd */
    if (mpz_cmp_ui(N, 2) < 0)
      return APRTCLE_COMPOSITE;
  
    if (mpz_divisible_ui_p(N, 2))
    {
      if (mpz_cmp_ui(N, 2) == 0)
        return APRTCLE_PRIME;
      else
        return APRTCLE_COMPOSITE;
    }
  
  /* only three small exceptions for this implementation */
  /* with this set of P and Q primes */
    if (mpz_cmp_ui(N, 3) == 0)
      return APRTCLE_PRIME;
    if (mpz_cmp_ui(N, 7) == 0)
      return APRTCLE_PRIME;
    if (mpz_cmp_ui(N, 11) == 0)
      return APRTCLE_PRIME;
  
    /* If the input number is larger than 7000 decimal digits
       we will just return whether it is a BPSW (probable) prime */
    NumberLength = mpz_sizeinbase(N, 10);
    if (NumberLength > 7000)
    {
      return mpz_probab_prime_p(N, 1); // mpz_probab_prime_p mpz_bpsw_prp
    }
  
    allocate_vars();
  
    mpz_set(TestNbr, N);
    mpz_set_si(biS, 0);
  
    j = PK = PL = PM = 0;
    for (J = 0; J < PWmax; J++)
    {
      /* aiJX[J] = 0; */
      mpz_set_ui(aiJX[J], 0);
    }
    break_this = 0;
  /* GetPrimes2Test : */
    for (i = 0; i < LEVELmax; i++)
    {
      /* biS[0] = 2; */
      mpz_set_ui(biS, 2);
      
      for (j = 0; j < aiNQ[i]; j++)
      {
        Q = aiQ[j];
        if (aiT[i]%(Q-1) != 0) continue;
        U = aiT[i] * Q;
        do
        {
          U /= Q;
          /* MultBigNbrByLong(biS, Q, biS, NumberLength); */
          mpz_mul_ui(biS, biS, Q);
        }
        while (U % Q == 0);
  
        // Exit loop if S^2 > N.
        if (CompareSquare(biS, TestNbr) > 0)
        {
          /* break GetPrimes2Test; */
          break_this = 1;
          break;
        }
      } /* End for j */
      if (break_this) break;
    } /* End for i */
    if (i == LEVELmax)
    { /* too big */
      free_vars();
      return mpz_probab_prime_p(N, 1);
    }
    LEVELnow = i;
    TestingQs = j;
    T = aiT[LEVELnow];
    NP = aiNP[LEVELnow];
  
  MainStart:
    for (;;)
    {
      for (i = 0; i < NP; i++)
      {
        P = aiP[i];
        if (T%P != 0) continue;
  
        SW = TestedQs = 0;
        /* Q = W = (int) BigNbrModLong(TestNbr, P * P); */
        Q = W = mpz_fdiv_ui(TestNbr, P * P);
        for (J = P - 2; J > 0; J--)
        {
          W = (W * Q) % (P * P);
        }
        if (P > 2 && W != 1)
        {
          SW = 1;
        }
        for (;;)
        {
          for (j = TestedQs; j <= TestingQs; j++)
          {
            Q = aiQ[j] - 1;
            /* G = aiG[j]; */
            K = 0;
            while (Q % P == 0)
            {
              K++;
              Q /= P;
            }
            Q = aiQ[j];
            if (K == 0)
            {
              continue;
            }
  
            if (verbose >= APRTCLE_VERBOSE1)
            {
              printf("P = %2d, Q = %12d  (%3.2f%%)\r", P, Q, (i * (TestingQs + 1) + j) * 100.0 / (NP * (TestingQs + 1)));
              fflush(stdout);
            }
  
            PM = 1;
            for (I = 1; I < K; I++)
            {
              PM = PM * P;
            }
            PL = (P - 1) * PM;
            PK = P * PM;
            for (I = 0; I < PK; I++)
            {
              /* aiJ0[I] = aiJ1[I] = 0; */
              mpz_set_ui(aiJ0[I], 0);
              mpz_set_ui(aiJ1[I], 0);
            }
            if (P > 2)
            {
              JacobiSum(0, P, PL, Q);
            }
            else
            {
              if (K != 1)
              {
                JacobiSum(0, P, PL, Q);
                for (I = 0; I < PK; I++)
                {
                  /* aiJW[I] = 0; */
                  mpz_set_ui(aiJW[I], 0);
                }
                if (K != 2)
                {
                  for (I = 0; I < PM; I++)
                  {
                    /* aiJW[I] = aiJ0[I]; */
                    mpz_set(aiJW[I], aiJ0[I]);
                  }
                  JacobiSum(1, P, PL, Q);
                  for (I = 0; I < PM; I++)
                  {
                    /* aiJS[I] = aiJ0[I]; */
                    mpz_set(aiJS[I], aiJ0[I]);
                  }
                  JS_JW(PK, PL, PM, P);
                  for (I = 0; I < PM; I++)
                  {
                    /* aiJ1[I] = aiJS[I]; */
                    mpz_set(aiJ1[I], aiJS[I]);
                  }
                  JacobiSum(2, P, PL, Q);
                  for (I = 0; I < PK; I++)
                  {
                    /* aiJW[I] = 0; */
                    mpz_set_ui(aiJW[I], 0);
                  }
                  for (I = 0; I < PM; I++)
                  {
                    /* aiJS[I] = aiJ0[I]; */
                    mpz_set(aiJS[I], aiJ0[I]);
                  }
                  JS_2(PK, PL, PM, P);
                  for (I = 0; I < PM; I++)
                  {
                    /* aiJ2[I] = aiJS[I]; */
                    mpz_set(aiJ2[I], aiJS[I]);
                  }
                }
              }
            }
            /* aiJ00[0] = aiJ01[0] = 1; */
            mpz_set_ui(aiJ00[0], 1);
            mpz_set_ui(aiJ01[0], 1);
            for (I = 1; I < PK; I++)
            {
              /* aiJ00[I] = aiJ01[I] = 0; */
              mpz_set_ui(aiJ00[I], 0);
              mpz_set_ui(aiJ01[I], 0);
            }
            /* VK = (int) BigNbrModLong(TestNbr, PK); */
            VK = mpz_fdiv_ui(TestNbr, PK);
            for (I = 1; I < PK; I++)
            {
              if (I % P != 0)
              {
                U1 = 1;
                U3 = I;
                V1 = 0;
                V3 = PK;
                while (V3 != 0)
                {
                  QQ = U3 / V3;
                  T1 = U1 - V1 * QQ;
                  T3 = U3 - V3 * QQ;
                  U1 = V1;
                  U3 = V3;
                  V1 = T1;
                  V3 = T3;
                }
                aiInv[I] = (U1 + PK) % PK;
              }
              else
              {
                aiInv[I] = 0;
              }
            }
            if (P != 2)
            {
              for (IV = 0; IV <= 1; IV++)
              {
                for (X = 1; X < PK; X++)
                {
                  for (I = 0; I < PK; I++)
                  {
                    /* aiJS[I] = aiJ0[I]; */
                    mpz_set(aiJS[I], aiJ0[I]);
                  }
                  if (X % P == 0)
                  {
                    continue;
                  }
                  if (IV == 0)
                  {
                    /* LongToBigNbr(X, biExp, NumberLength); */
                    mpz_set_ui(biExp, X);
                  }
                  else
                  {
                    /* LongToBigNbr(VK * X / PK, biExp, NumberLength); */
                    mpz_set_ui(biExp, (VK * X) / PK);
                    if ((VK * X) / PK == 0)
                    {
                      continue;
                    }
                  }
                  JS_E(PK, PL, PM, P);
                  for (I = 0; I < PK; I++)
                  {
                    /* aiJW[I] = 0; */
                    mpz_set_ui(aiJW[I], 0);
                  }
                  InvX = aiInv[X];
                  for (I = 0; I < PK; I++)
                  {
                    J = (I * InvX) % PK;
                    /* AddBigNbrModN(aiJW[J], aiJS[I], aiJW[J], TestNbr, NumberLength); */
                    mpz_add(aiJW[J], aiJW[J], aiJS[I]);
                  }
                  NormalizeJW(PK, PL, PM, P);
                  if (IV == 0)
                  {
                    for (I = 0; I < PK; I++)
                    {
                      /* aiJS[I] = aiJ00[I]; */
                      mpz_set(aiJS[I], aiJ00[I]);
                    }
                  }
                  else
                  {
                    for (I = 0; I < PK; I++)
                    {
                      /* aiJS[I] = aiJ01[I]; */
                      mpz_set(aiJS[I], aiJ01[I]);
                    }
                  }
                  JS_JW(PK, PL, PM, P);
                  if (IV == 0)
                  {
                    for (I = 0; I < PK; I++)
                    {
                      /* aiJ00[I] = aiJS[I]; */
                      mpz_set(aiJ00[I], aiJS[I]);
                    }
                  }
                  else
                  {
                    for (I = 0; I < PK; I++)
                    {
                      /* aiJ01[I] = aiJS[I]; */
                      mpz_set(aiJ01[I], aiJS[I]);
                    }
                  }
                } /* end for X */
              } /* end for IV */
            }
            else
            {
              if (K == 1)
              {
                /* MultBigNbrByLongModN(1, Q, aiJ00[0], TestNbr, NumberLength); */
                mpz_set_ui(aiJ00[0], Q);
                /* aiJ01[0] = 1; */
                mpz_set_ui(aiJ01[0], 1);
              }
              else
              {
                if (K == 2)
                {
                  if (VK == 1)
                  {
                    /* aiJ01[0] = 1; */
                    mpz_set_ui(aiJ01[0], 1);
                  }
                  /* aiJS[0] = aiJ0[0]; */
                  /* aiJS[1] = aiJ0[1]; */
                  mpz_set(aiJS[0], aiJ0[0]);
                  mpz_set(aiJS[1], aiJ0[1]);
                  JS_2(PK, PL, PM, P);
                  if (VK == 3)
                  {
                    /* aiJ01[0] = aiJS[0]; */
                    /* aiJ01[1] = aiJS[1]; */
                    mpz_set(aiJ01[0], aiJS[0]);
                    mpz_set(aiJ01[1], aiJS[1]);
                  }
                  /* MultBigNbrByLongModN(aiJS[0], Q, aiJ00[0], TestNbr, NumberLength); */
                  mpz_mul_ui(aiJ00[0], aiJS[0], Q);
                  /* MultBigNbrByLongModN(aiJS[1], Q, aiJ00[1], TestNbr, NumberLength); */
                  mpz_mul_ui(aiJ00[1], aiJS[1], Q);
                }
                else
                {
                  for (IV = 0; IV <= 1; IV++)
                  {
                    for (X = 1; X < PK; X += 2)
                    {
                      for (I = 0; I <= PM; I++)
                      {
                        /* aiJS[I] = aiJ1[I]; */
                        mpz_set(aiJS[I], aiJ1[I]);
                      }
                      if (X % 8 == 5 || X % 8 == 7)
                      {
                        continue;
                      }
                      if (IV == 0)
                      {
                        /* LongToBigNbr(X, biExp, NumberLength); */
                        mpz_set_ui(biExp, X);
                      }
                      else
                      {
                        /* LongToBigNbr(VK * X / PK, biExp, NumberLength); */
                        mpz_set_ui(biExp, VK * X / PK);
                        if (VK * X / PK == 0)
                        {
                          continue;
                        }
                      }
                      JS_E(PK, PL, PM, P);
                      for (I = 0; I < PK; I++)
                      {
                        /* aiJW[I] = 0; */
                        mpz_set_ui(aiJW[I], 0);
                      }
                      InvX = aiInv[X];
                      for (I = 0; I < PK; I++)
                      {
                        J = I * InvX % PK;
                        /* AddBigNbrModN(aiJW[J], aiJS[I], aiJW[J], TestNbr, NumberLength); */
                        mpz_add(aiJW[J], aiJW[J], aiJS[I]);
                      }
                      NormalizeJW(PK, PL, PM, P);
                      if (IV == 0)
                      {
                        for (I = 0; I < PK; I++)
                        {
                          /* aiJS[I] = aiJ00[I]; */
                          mpz_set(aiJS[I], aiJ00[I]);
                        }
                      }
                      else
                      {
                        for (I = 0; I < PK; I++)
                        {
                          /* aiJS[I] = aiJ01[I]; */
                          mpz_set(aiJS[I], aiJ01[I]);
                        }
                      }
                      NormalizeJS(PK, PL, PM, P);
                      JS_JW(PK, PL, PM, P);
                      if (IV == 0)
                      {
                        for (I = 0; I < PK; I++)
                        {
                          /* aiJ00[I] = aiJS[I]; */
                          mpz_set(aiJ00[I], aiJS[I]);
                        }
                      }
                      else
                      {
                        for (I = 0; I < PK; I++)
                        {
                          /* aiJ01[I] = aiJS[I]; */
                          mpz_set(aiJ01[I], aiJS[I]);
                        }
                      }
                    } /* end for X */
                    if (IV == 0 || VK % 8 == 1 || VK % 8 == 3)
                    {
                      continue;
                    }
                    for (I = 0; I < PM; I++)
                    {
                      /* aiJW[I] = aiJ2[I]; */
                      /* aiJS[I] = aiJ01[I]; */
                      mpz_set(aiJW[I], aiJ2[I]);
                      mpz_set(aiJS[I], aiJ01[I]);
                    }
                    for (; I < PK; I++)
                    {
                      /* aiJW[I] = aiJS[I] = 0; */
                      mpz_set_ui(aiJW[I], 0);
                      mpz_set_ui(aiJS[I], 0);
                    }
                    JS_JW(PK, PL, PM, P);
                    for (I = 0; I < PM; I++)
                    {
                      /* aiJ01[I] = aiJS[I]; */
                      mpz_set(aiJ01[I], aiJS[I]);
                    }
                  } /* end for IV */
                }
              }
            }
            for (I = 0; I < PL; I++)
            {
              /* aiJS[I] = aiJ00[I]; */
              mpz_set(aiJS[I], aiJ00[I]);
            }
            for (; I < PK; I++)
            {
              /* aiJS[I] = 0; */
              mpz_set_ui(aiJS[I], 0);
            }
            /* DivBigNbrByLong(TestNbr, PK, biExp, NumberLength); */
            mpz_fdiv_q_ui(biExp, TestNbr, PK);
            JS_E(PK, PL, PM, P);
            for (I = 0; I < PK; I++)
            {
              /* aiJW[I] = 0; */
              mpz_set_ui(aiJW[I], 0);
            }
            for (I = 0; I < PL; I++)
            {
              for (J = 0; J < PL; J++)
              {
                /* MontgomeryMult(aiJS[I], aiJ01[J], biTmp); */
                /* AddBigNbrModN(biTmp, aiJW[(I + J) % PK], aiJW[(I + J) % PK], TestNbr, NumberLength); */
                mpz_mul(biTmp, aiJS[I], aiJ01[J]);
                mpz_add(aiJW[(I + J) % PK], biTmp, aiJW[(I + J) % PK]);
              }
            }
            NormalizeJW(PK, PL, PM, P);
  /* MatchingRoot : */
            do
            {
              H = -1;
              W = 0;
              for (I = 0; I < PL; I++)
              {
                if (mpz_cmp_ui(aiJW[I], 0) != 0)/* (!BigNbrIsZero(aiJW[I])) */
                {
                  /* if (H == -1 && BigNbrAreEqual(aiJW[I], 1)) */
                  if (H == -1 && (mpz_cmp_ui(aiJW[I], 1) == 0))
                  {
                    H = I;
                  }
                  else
                  {
                    H = -2;
                    /* AddBigNbrModN(aiJW[I], MontgomeryMultR1, biTmp, TestNbr, NumberLength); */
                    mpz_add_ui(biTmp, aiJW[I], 1);
                    mpz_mod(biTmp, biTmp, TestNbr);
                    if (mpz_cmp_ui(biTmp, 0) == 0) /* (BigNbrIsZero(biTmp)) */
                    {
                      W++;
                    }
                  }
                }
              }
              if (H >= 0)
              {
                /* break MatchingRoot; */
                break;
              }
              if (W != P - 1)
              {
                /* Not prime */
                free_vars();
                return APRTCLE_COMPOSITE;
              }
              for (I = 0; I < PM; I++)
              {
                /* AddBigNbrModN(aiJW[I], 1, biTmp, TestNbr, NumberLength); */
                mpz_add_ui(biTmp, aiJW[I], 1);
                mpz_mod(biTmp, biTmp, TestNbr);
                if (mpz_cmp_ui(biTmp, 0) == 0) /* (BigNbrIsZero(biTmp)) */
                {
                  break;
                }
              }
              if (I == PM)
              {
                /* Not prime */
                free_vars();
                return APRTCLE_COMPOSITE;
              }
              for (J = 1; J <= P - 2; J++)
              {
                /* AddBigNbrModN(aiJW[I + J * PM], 1, biTmp, TestNbr, NumberLength); */
                mpz_add_ui(biTmp, aiJW[I + J * PM], 1);
                mpz_mod(biTmp, biTmp, TestNbr);
                if (mpz_cmp_ui(biTmp, 0) != 0)/* (!BigNbrIsZero(biTmp)) */
                {
                  /* Not prime */
                  free_vars();
                  return APRTCLE_COMPOSITE;
                }
              }
              H = I + PL;
            }
            while (0);
  
            if (SW == 1 || H % P == 0)
            {
              continue;
            }
            if (P != 2)
            {
              SW = 1;
              continue;
            }
            if (K == 1)
            {
              if ((mpz_get_ui(TestNbr) & 3) == 1)
              {
                SW = 1;
              }
              continue;
            }
  
            // if (Q^((N-1)/2) mod N != N-1), N is not prime.
  
            /* MultBigNbrByLongModN(1, Q, biTmp, TestNbr, NumberLength); */
            mpz_set_ui(biTmp, Q);
            mpz_mod(biTmp, biTmp, TestNbr);
  
            mpz_sub_ui(biT, TestNbr, 1); /* biT = n-1 */
            mpz_divexact_ui(biT, biT, 2); /* biT = (n-1)/2 */
            mpz_powm(biR, biTmp, biT, TestNbr); /* biR = Q^((n-1)/2) mod n */
            mpz_add_ui(biTmp, biR, 1);
            mpz_mod(biTmp, biTmp, TestNbr);
  
            if (mpz_cmp_ui(biTmp, 0) != 0)/* (!BigNbrIsZero(biTmp)) */
            {
              /* Not prime */
              free_vars();
              return APRTCLE_COMPOSITE;
            }
            SW = 1;
          } /* end for j */
          if (SW == 0)
          {
            TestedQs = TestingQs + 1;
            if (TestingQs < aiNQ[LEVELnow] - 1)
            {
              TestingQs++;
              Q = aiQ[TestingQs];
              U = T * Q;
              do
              {
                /* MultBigNbrByLong(biS, Q, biS, NumberLength); */
                mpz_mul_ui(biS, biS, Q);
                U /= Q;
              }
              while (U % Q == 0);
  
              continue; /* Retry */
            }
            LEVELnow++;
            if (LEVELnow == LEVELmax)
            {
              free_vars();
              return mpz_probab_prime_p(N, 1); /* Cannot tell */
            }
            T = aiT[LEVELnow];
            NP = aiNP[LEVELnow];
            /* biS = 2; */
            mpz_set_ui(biS, 2);
            for (J = 0; J <= aiNQ[LEVELnow]; J++)
            {
              Q = aiQ[J];
              if (T%(Q-1) != 0) continue;
              U = T * Q;
              do
              {
                /* MultBigNbrByLong(biS, Q, biS, NumberLength); */
                mpz_mul_ui(biS, biS, Q);
                U /= Q;
              }
              while (U % Q == 0);
              if (CompareSquare(biS, TestNbr) > 0)
              {
                TestingQs = J;
                /* continue MainStart; */ /* Retry from the beginning */
                goto MainStart;
              }
            } /* end for J */
            free_vars();
            return mpz_probab_prime_p(N, 1); /* Program error */
          } /* end if */
          break;
        } /* end for (;;) */
      } /* end for i */
  
      // Final Test
     
      /* biR = 1 */
      mpz_set_ui(biR, 1);
      /* biN <- TestNbr mod biS */ /* Compute N mod S */
      mpz_fdiv_r(biN, TestNbr, biS);
     
      for (U = 1; U <= T; U++)
      {
        /* biR <- (biN * biR) mod biS */
        mpz_mul(biR, biN, biR);
        mpz_mod(biR, biR, biS);
        if (mpz_cmp_ui(biR, 1) == 0) /* biR == 1 */
        {
          /* Number is prime */
          free_vars();
          return APRTCLE_PRIME;
        }
        if (mpz_divisible_p(TestNbr, biR) && mpz_cmp(biR, TestNbr) < 0) /* biR < N and biR | TestNbr */
        {
          /* Number is composite */
          free_vars();
          return APRTCLE_COMPOSITE;
        }
      } /* End for U */
      /* This should never be reached. */
      free_vars();
      return mpz_probab_prime_p(N, 1);
    }
  }
  
  /* ============================================================================================== */
  