view release on metacpan or  search on metacpan

spec/bpsw/trn.c  view on Meta::CPAN

mpz_powm(mpzRem, mpzB, mpzd, mpzN);
if(mpz_cmp_ui(mpzRem, 1)==0)RETURN(1);
if(mpz_cmp(mpzRem, mpzNm1)==0)RETURN(1);

/* Now calculate B^2d, B^4d, B^8d, ..., B^((N-1)/2) by successive
   squaring. If any of these is congruent to -1 mod N, N is a
   sprp base B. */

for(r=1; r < s; r++)
  {
  mpz_mul(mpzRem, mpzRem, mpzRem);
  mpz_mod(mpzRem, mpzRem, mpzN);
  if(mpz_cmp(mpzRem, mpzNm1)==0)RETURN(1);
  }
RETURN(0);
}
/**********************************************************************/
int iLucasSelfridge(mpz_t mpzN)
{
/* Test mpzN for primality using Lucas's test with Selfridge's parameters.
   Returns 1 if mpzN is prime or a Lucas-Selfridge pseudoprime. Returns
   0 if mpzN is definitely composite. Note that a Lucas-Selfridge test
   typically requires three to seven times as many bit operations as a
   single Miller-Rabin test. The frequency of Lucas-Selfridge pseudoprimes
   appears to be roughly four times that of base-2 strong pseudoprimes;
   the Baillie-PSW test is based on the hope (verified by the author,
   May, 2005, for all N < 10^13; and by Martin Fuller, January, 2007,
   for all N < 10^15) that the two tests have no common pseudoprimes. */

int iComp2, iP, iJ, iSign;
long lDabs, lD, lQ;
unsigned long ulMaxBits, ulNbits, ul, ulGCD;
mpz_t mpzU, mpzV, mpzNplus1, mpzU2m, mpzV2m, mpzQm, mpz2Qm,
      mpzT1, mpzT2, mpzT3, mpzT4, mpzD;

#undef RETURN
#define RETURN(n)           \
  {                         \
  mpz_clear(mpzU);          \
  mpz_clear(mpzV);          \
  mpz_clear(mpzNplus1);     \
  mpz_clear(mpzU2m);        \
  mpz_clear(mpzV2m);        \
  mpz_clear(mpzQm);         \
  mpz_clear(mpz2Qm);        \
  mpz_clear(mpzT1);         \
  mpz_clear(mpzT2);         \
  mpz_clear(mpzT3);         \
  mpz_clear(mpzT4);         \
  mpz_clear(mpzD);          \
  return(n);                \
  }

/* This implementation of the algorithm assumes N is an odd integer > 2,
   so we first eliminate all N < 3 and all even N. As a practical matter,
   we also need to filter out all perfect square values of N, such as
   1093^2 (a base-2 strong pseudoprime); this is because we will later
   require an integer D for which Jacobi(D,N) = -1, and no such integer
   exists if N is a perfect square. The algorithm as written would
   still eventually return zero in this case, but would require
   nearly sqrt(N)/2 iterations. */

iComp2=mpz_cmp_si(mpzN, 2);
if(iComp2 < 0)return(0);
if(iComp2==0)return(1);
if(mpz_even_p(mpzN))return(0);
if(mpz_perfect_square_p(mpzN))return(0);

/* Allocate storage for the mpz_t variables. Most require twice
   the storage of mpzN, since multiplications of order O(mpzN)*O(mpzN)
   will be performed. */

ulMaxBits=2*mpz_sizeinbase(mpzN, 2) + mp_bits_per_limb;
mpz_init2(mpzU, ulMaxBits);
mpz_init2(mpzV, ulMaxBits);
mpz_init2(mpzNplus1, ulMaxBits);
mpz_init2(mpzU2m, ulMaxBits);
mpz_init2(mpzV2m, ulMaxBits);
mpz_init2(mpzQm, ulMaxBits);
mpz_init2(mpz2Qm, ulMaxBits);
mpz_init2(mpzT1, ulMaxBits);
mpz_init2(mpzT2, ulMaxBits);
mpz_init2(mpzT3, ulMaxBits);
mpz_init2(mpzT4, ulMaxBits);
mpz_init(mpzD);

/* Find the first element D in the sequence {5, -7, 9, -11, 13, ...}
   such that Jacobi(D,N) = -1 (Selfridge's algorithm). Although
   D will nearly always be "small" (perfect square N's having
   been eliminated), an overflow trap for D is present. */

lDabs=5;
iSign=1;
while(1)
  {
  lD=iSign*lDabs;
  iSign = -iSign;
  ulGCD=mpz_gcd_ui(NULL, mpzN, lDabs);
  /* if 1 < GCD < N then N is composite with factor lDabs, and
     Jacobi(D,N) is technically undefined (but often returned
     as zero). */
  if((ulGCD > 1) && mpz_cmp_ui(mpzN, ulGCD) > 0)RETURN(0);
  mpz_set_si(mpzD, lD);
  iJ=mpz_jacobi(mpzD, mpzN);
  if(iJ==-1)break;
  lDabs += 2;
  if(lDabs > ulDmax)ulDmax=lDabs;  /* tracks global max of |D| */
  if(lDabs > LONG_MAX-2)
    {
    fprintf(stderr,
      "\n ERROR: D overflows signed long in Lucas-Selfridge test.");
    fprintf(stderr, "\n N=");
    mpz_out_str(stderr, 10, mpzN);
    fprintf(stderr, "\n |D|=%ld\n\n", lDabs);
    exit(EXIT_FAILURE);
    }
  }

iP=1;         /* Selfridge's choice */
lQ=(1-lD)/4;  /* Required so D = P*P - 4*Q */

spec/bpsw/trn.c  view on Meta::CPAN

  }

/* If U_(N - Jacobi(D,N)) is congruent to 0 mod N, then N is
   a prime or a Lucas pseudoprime; otherwise it is definitely
   composite. */

if(mpz_sgn(mpzU)==0)RETURN(1);
RETURN(0);
}
/**********************************************************************/
int iStrongLucasSelfridge(mpz_t mpzN)
{
/* Test N for primality using the strong Lucas test with Selfridge's
   parameters. Returns 1 if N is prime or a strong Lucas-Selfridge
   pseudoprime (in which case N is also a pseudoprime to the standard
   Lucas-Selfridge test). Returns 0 if N is definitely composite.

   The running time of the strong Lucas-Selfridge test is, on average,
   roughly 10 % greater than the running time for the standard
   Lucas-Selfridge test (3 to 7 times that of a single Miller-Rabin
   test). However, the frequency of strong Lucas pseudoprimes appears
   to be only (roughly) 30 % that of (standard) Lucas pseudoprimes, and
   only slightly greater than the frequency of base-2 strong pseudoprimes,
   indicating that the strong Lucas-Selfridge test is more computationally
   effective than the standard version. */

int iComp2, iP, iJ, iSign;
long lDabs, lD, lQ;
unsigned long ulMaxBits, uldbits, ul, ulGCD, r, s;
mpz_t mpzU, mpzV, mpzNplus1, mpzU2m, mpzV2m, mpzQm, mpz2Qm,
      mpzT1, mpzT2, mpzT3, mpzT4, mpzD, mpzd, mpzQkd, mpz2Qkd;

#undef RETURN
#define RETURN(n)           \
  {                         \
  mpz_clear(mpzU);          \
  mpz_clear(mpzV);          \
  mpz_clear(mpzNplus1);     \
  mpz_clear(mpzU2m);        \
  mpz_clear(mpzV2m);        \
  mpz_clear(mpzQm);         \
  mpz_clear(mpz2Qm);        \
  mpz_clear(mpzT1);         \
  mpz_clear(mpzT2);         \
  mpz_clear(mpzT3);         \
  mpz_clear(mpzT4);         \
  mpz_clear(mpzD);          \
  mpz_clear(mpzd);          \
  mpz_clear(mpzQkd);        \
  mpz_clear(mpz2Qkd);       \
  return(n);                \
  }

/* This implementation of the algorithm assumes N is an odd integer > 2,
   so we first eliminate all N < 3 and all even N. As a practical matter,
   we also need to filter out all perfect square values of N, such as
   1093^2 (a base-2 strong pseudoprime); this is because we will later
   require an integer D for which Jacobi(D,N) = -1, and no such integer
   exists if N is a perfect square. The algorithm as written would
   still eventually return zero in this case, but would require
   nearly sqrt(N)/2 iterations. */

iComp2=mpz_cmp_si(mpzN, 2);
if(iComp2 < 0)return(0);
if(iComp2==0)return(1);
if(mpz_even_p(mpzN))return(0);
if(mpz_perfect_square_p(mpzN))return(0);

/* Allocate storage for the mpz_t variables. Most require twice
   the storage of mpzN, since multiplications of order O(mpzN)*O(mpzN)
   will be performed. */

ulMaxBits=2*mpz_sizeinbase(mpzN, 2) + mp_bits_per_limb;
mpz_init2(mpzU, ulMaxBits);
mpz_init2(mpzV, ulMaxBits);
mpz_init2(mpzNplus1, ulMaxBits);
mpz_init2(mpzU2m, ulMaxBits);
mpz_init2(mpzV2m, ulMaxBits);
mpz_init2(mpzQm, ulMaxBits);
mpz_init2(mpz2Qm, ulMaxBits);
mpz_init2(mpzT1, ulMaxBits);
mpz_init2(mpzT2, ulMaxBits);
mpz_init2(mpzT3, ulMaxBits);
mpz_init2(mpzT4, ulMaxBits);
mpz_init(mpzD);
mpz_init2(mpzd, ulMaxBits);
mpz_init2(mpzQkd, ulMaxBits);
mpz_init2(mpz2Qkd, ulMaxBits);

/* Find the first element D in the sequence {5, -7, 9, -11, 13, ...}
   such that Jacobi(D,N) = -1 (Selfridge's algorithm). Theory
   indicates that, if N is not a perfect square, D will "nearly
   always" be "small." Just in case, an overflow trap for D is
   included. */

lDabs=5;
iSign=1;
while(1)
  {
  lD=iSign*lDabs;
  iSign = -iSign;
  ulGCD=mpz_gcd_ui(NULL, mpzN, lDabs);
  /* if 1 < GCD < N then N is composite with factor lDabs, and
     Jacobi(D,N) is technically undefined (but often returned
     as zero). */
  if((ulGCD > 1) && mpz_cmp_ui(mpzN, ulGCD) > 0)RETURN(0);
  mpz_set_si(mpzD, lD);
  iJ=mpz_jacobi(mpzD, mpzN);
  if(iJ==-1)break;
  lDabs += 2;
  if(lDabs > ulDmax)ulDmax=lDabs;  /* tracks global max of |D| */
  if(lDabs > LONG_MAX-2)
    {
    fprintf(stderr,
      "\n ERROR: D overflows signed long in Lucas-Selfridge test.");
    fprintf(stderr, "\n N=");
    mpz_out_str(stderr, 10, mpzN);
    fprintf(stderr, "\n |D|=%ld\n\n", lDabs);
    exit(EXIT_FAILURE);
    }
  }