github.com/aergoio/aergo@v1.3.1/libtool/src/gmp-6.1.2/mpz/pprime_p.c (about)

     1  /* mpz_probab_prime_p --
     2     An implementation of the probabilistic primality test found in Knuth's
     3     Seminumerical Algorithms book.  If the function mpz_probab_prime_p()
     4     returns 0 then n is not prime.  If it returns 1, then n is 'probably'
     5     prime.  If it returns 2, n is surely prime.  The probability of a false
     6     positive is (1/4)**reps, where reps is the number of internal passes of the
     7     probabilistic algorithm.  Knuth indicates that 25 passes are reasonable.
     8  
     9  Copyright 1991, 1993, 1994, 1996-2002, 2005 Free Software Foundation, Inc.
    10  
    11  This file is part of the GNU MP Library.
    12  
    13  The GNU MP Library is free software; you can redistribute it and/or modify
    14  it under the terms of either:
    15  
    16    * the GNU Lesser General Public License as published by the Free
    17      Software Foundation; either version 3 of the License, or (at your
    18      option) any later version.
    19  
    20  or
    21  
    22    * the GNU General Public License as published by the Free Software
    23      Foundation; either version 2 of the License, or (at your option) any
    24      later version.
    25  
    26  or both in parallel, as here.
    27  
    28  The GNU MP Library is distributed in the hope that it will be useful, but
    29  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
    30  or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    31  for more details.
    32  
    33  You should have received copies of the GNU General Public License and the
    34  GNU Lesser General Public License along with the GNU MP Library.  If not,
    35  see https://www.gnu.org/licenses/.  */
    36  
    37  #include "gmp.h"
    38  #include "gmp-impl.h"
    39  #include "longlong.h"
    40  
    41  static int isprime (unsigned long int);
    42  
    43  
    44  /* MPN_MOD_OR_MODEXACT_1_ODD can be used instead of mpn_mod_1 for the trial
    45     division.  It gives a result which is not the actual remainder r but a
    46     value congruent to r*2^n mod d.  Since all the primes being tested are
    47     odd, r*2^n mod p will be 0 if and only if r mod p is 0.  */
    48  
    49  int
    50  mpz_probab_prime_p (mpz_srcptr n, int reps)
    51  {
    52    mp_limb_t r;
    53    mpz_t n2;
    54  
    55    /* Handle small and negative n.  */
    56    if (mpz_cmp_ui (n, 1000000L) <= 0)
    57      {
    58        int is_prime;
    59        if (mpz_cmpabs_ui (n, 1000000L) <= 0)
    60  	{
    61  	  is_prime = isprime (mpz_get_ui (n));
    62  	  return is_prime ? 2 : 0;
    63  	}
    64        /* Negative number.  Negate and fall out.  */
    65        PTR(n2) = PTR(n);
    66        SIZ(n2) = -SIZ(n);
    67        n = n2;
    68      }
    69  
    70    /* If n is now even, it is not a prime.  */
    71    if ((mpz_get_ui (n) & 1) == 0)
    72      return 0;
    73  
    74  #if defined (PP)
    75    /* Check if n has small factors.  */
    76  #if defined (PP_INVERTED)
    77    r = MPN_MOD_OR_PREINV_MOD_1 (PTR(n), (mp_size_t) SIZ(n), (mp_limb_t) PP,
    78  			       (mp_limb_t) PP_INVERTED);
    79  #else
    80    r = mpn_mod_1 (PTR(n), (mp_size_t) SIZ(n), (mp_limb_t) PP);
    81  #endif
    82    if (r % 3 == 0
    83  #if GMP_LIMB_BITS >= 4
    84        || r % 5 == 0
    85  #endif
    86  #if GMP_LIMB_BITS >= 8
    87        || r % 7 == 0
    88  #endif
    89  #if GMP_LIMB_BITS >= 16
    90        || r % 11 == 0 || r % 13 == 0
    91  #endif
    92  #if GMP_LIMB_BITS >= 32
    93        || r % 17 == 0 || r % 19 == 0 || r % 23 == 0 || r % 29 == 0
    94  #endif
    95  #if GMP_LIMB_BITS >= 64
    96        || r % 31 == 0 || r % 37 == 0 || r % 41 == 0 || r % 43 == 0
    97        || r % 47 == 0 || r % 53 == 0
    98  #endif
    99        )
   100      {
   101        return 0;
   102      }
   103  #endif /* PP */
   104  
   105    /* Do more dividing.  We collect small primes, using umul_ppmm, until we
   106       overflow a single limb.  We divide our number by the small primes product,
   107       and look for factors in the remainder.  */
   108    {
   109      unsigned long int ln2;
   110      unsigned long int q;
   111      mp_limb_t p1, p0, p;
   112      unsigned int primes[15];
   113      int nprimes;
   114  
   115      nprimes = 0;
   116      p = 1;
   117      ln2 = mpz_sizeinbase (n, 2);	/* FIXME: tune this limit */
   118      for (q = PP_FIRST_OMITTED; q < ln2; q += 2)
   119        {
   120  	if (isprime (q))
   121  	  {
   122  	    umul_ppmm (p1, p0, p, q);
   123  	    if (p1 != 0)
   124  	      {
   125  		r = MPN_MOD_OR_MODEXACT_1_ODD (PTR(n), (mp_size_t) SIZ(n), p);
   126  		while (--nprimes >= 0)
   127  		  if (r % primes[nprimes] == 0)
   128  		    {
   129  		      ASSERT_ALWAYS (mpn_mod_1 (PTR(n), (mp_size_t) SIZ(n), (mp_limb_t) primes[nprimes]) == 0);
   130  		      return 0;
   131  		    }
   132  		p = q;
   133  		nprimes = 0;
   134  	      }
   135  	    else
   136  	      {
   137  		p = p0;
   138  	      }
   139  	    primes[nprimes++] = q;
   140  	  }
   141        }
   142    }
   143  
   144    /* Perform a number of Miller-Rabin tests.  */
   145    return mpz_millerrabin (n, reps);
   146  }
   147  
   148  static int
   149  isprime (unsigned long int t)
   150  {
   151    unsigned long int q, r, d;
   152  
   153    if (t < 3 || (t & 1) == 0)
   154      return t == 2;
   155  
   156    for (d = 3, r = 1; r != 0; d += 2)
   157      {
   158        q = t / d;
   159        r = t - q * d;
   160        if (q < d)
   161  	return 1;
   162      }
   163    return 0;
   164  }