rts/gmp/mpz/pprime_p.c

   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 (C) 1991, 1993, 1994, 1996, 1997, 1998, 1999, 2000 Free Software
  10 Foundation, Inc.  Miller-Rabin code contributed by John Amanatides.
  11
  12 This file is part of the GNU MP Library.
  13
  14 The GNU MP Library is free software; you can redistribute it and/or modify
  15 it under the terms of the GNU Lesser General Public License as published by
  16 the Free Software Foundation; either version 2.1 of the License, or (at your
  17 option) any later version.
  18
  19 The GNU MP Library is distributed in the hope that it will be useful, but
  20 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  21 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  22 License for more details.
  23
  24 You should have received a copy of the GNU Lesser General Public License
  25 along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
  26 the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
  27 MA 02111-1307, USA. */
  28
  29 #include "gmp.h"
  30 #include "gmp-impl.h"
  31 #include "longlong.h"
  32
  33 static int isprime _PROTO ((unsigned long int t));
  34 static int mpz_millerrabin _PROTO ((mpz_srcptr n, int reps));
  35
  36 int
  37 #if __STDC__
  38 mpz_probab_prime_p (mpz_srcptr n, int reps)
  39 #else
  40 mpz_probab_prime_p (n, reps)
  41      mpz_srcptr n;
  42      int reps;
  43 #endif
  44 {
  45   mp_limb_t r;
  46
  47   /* Handle small and negative n.  */
  48   if (mpz_cmp_ui (n, 1000000L) <= 0)
  49     {
  50       int is_prime;
  51       if (mpz_sgn (n) < 0)
  52         {
  53           /* Negative number.  Negate and call ourselves.  */
  54           mpz_t n2;
  55           mpz_init (n2);
  56           mpz_neg (n2, n);
  57           is_prime = mpz_probab_prime_p (n2, reps);
  58           mpz_clear (n2);
  59           return is_prime;
  60         }
  61       is_prime = isprime (mpz_get_ui (n));
  62       return is_prime ? 2 : 0;
  63     }
  64
  65   /* If n is now even, it is not a prime.  */
  66   if ((mpz_get_ui (n) & 1) == 0)
  67     return 0;
  68
  69   /* Check if n has small factors.  */
  70   if (UDIV_TIME > (2 * UMUL_TIME + 6))
  71     r = mpn_preinv_mod_1 (PTR(n), SIZ(n), (mp_limb_t) PP, (mp_limb_t) PP_INVERTED);
  72   else
  73     r = mpn_mod_1 (PTR(n), SIZ(n), (mp_limb_t) PP);
  74   if (r % 3 == 0 || r % 5 == 0 || r % 7 == 0 || r % 11 == 0 || r % 13 == 0
  75       || r % 17 == 0 || r % 19 == 0 || r % 23 == 0 || r % 29 == 0
  76 #if BITS_PER_MP_LIMB == 64
  77       || r % 31 == 0 || r % 37 == 0 || r % 41 == 0 || r % 43 == 0
  78       || r % 47 == 0 || r % 53 == 0
  79 #endif
  80       )
  81     {
  82       return 0;
  83     }
  84
  85   /* Do more dividing.  We collect small primes, using umul_ppmm, until we
  86      overflow a single limb.  We divide our number by the small primes product,
  87      and look for factors in the remainder.  */
  88   {
  89     unsigned long int ln2;
  90     unsigned long int q;
  91     mp_limb_t p1, p0, p;
  92     unsigned int primes[15];
  93     int nprimes;
  94
  95     nprimes = 0;
  96     p = 1;
  97     ln2 = mpz_sizeinbase (n, 2) / 30; ln2 = ln2 * ln2;
  98     for (q = BITS_PER_MP_LIMB == 64 ? 59 : 31; q < ln2; q += 2)
  99       {
 100         if (isprime (q))
 101           {
 102             umul_ppmm (p1, p0, p, q);
 103             if (p1 != 0)
 104               {
 105                 r = mpn_mod_1 (PTR(n), SIZ(n), p);
 106                 while (--nprimes >= 0)
 107                   if (r % primes[nprimes] == 0)
 108                     {
 109                       if (mpn_mod_1 (PTR(n), SIZ(n), (mp_limb_t) primes[nprimes]) != 0)
 110                         abort ();
 111                       return 0;
 112                     }
 113                 p = q;
 114                 nprimes = 0;
 115               }
 116             else
 117               {
 118                 p = p0;
 119               }
 120             primes[nprimes++] = q;
 121           }
 122       }
 123   }
 124
 125   /* Perform a number of Miller-Rabin tests.  */
 126   return mpz_millerrabin (n, reps);
 127 }
 128
 129 static int
 130 #if __STDC__
 131 isprime (unsigned long int t)
 132 #else
 133 isprime (t)
 134      unsigned long int t;
 135 #endif
 136 {
 137   unsigned long int q, r, d;
 138
 139   if (t < 3 || (t & 1) == 0)
 140     return t == 2;
 141
 142   for (d = 3, r = 1; r != 0; d += 2)
 143     {
 144       q = t / d;
 145       r = t - q * d;
 146       if (q < d)
 147         return 1;
 148     }
 149   return 0;
 150 }
 151
 152 static int millerrabin _PROTO ((mpz_srcptr n, mpz_srcptr nm1,
 153                                 mpz_ptr x, mpz_ptr y,
 154                                 mpz_srcptr q, unsigned long int k));
 155
 156 static int
 157 #if __STDC__
 158 mpz_millerrabin (mpz_srcptr n, int reps)
 159 #else
 160 mpz_millerrabin (n, reps)
 161      mpz_srcptr n;
 162      int reps;
 163 #endif
 164 {
 165   int r;
 166   mpz_t nm1, x, y, q;
 167   unsigned long int k;
 168   gmp_randstate_t rstate;
 169   int is_prime;
 170   TMP_DECL (marker);
 171   TMP_MARK (marker);
 172
 173   MPZ_TMP_INIT (nm1, SIZ (n) + 1);
 174   mpz_sub_ui (nm1, n, 1L);
 175
 176   MPZ_TMP_INIT (x, SIZ (n));
 177   MPZ_TMP_INIT (y, 2 * SIZ (n)); /* mpz_powm_ui needs excessive memory!!! */
 178
 179   /* Perform a Fermat test.  */
 180   mpz_set_ui (x, 210L);
 181   mpz_powm (y, x, nm1, n);
 182   if (mpz_cmp_ui (y, 1L) != 0)
 183     {
 184       TMP_FREE (marker);
 185       return 0;
 186     }
 187
 188   MPZ_TMP_INIT (q, SIZ (n));
 189
 190   /* Find q and k, where q is odd and n = 1 + 2**k * q.  */
 191   k = mpz_scan1 (nm1, 0L);
 192   mpz_tdiv_q_2exp (q, nm1, k);
 193
 194   gmp_randinit (rstate, GMP_RAND_ALG_DEFAULT, 32L);
 195
 196   is_prime = 1;
 197   for (r = 0; r < reps && is_prime; r++)
 198     {
 199       do
 200         mpz_urandomb (x, rstate, mpz_sizeinbase (n, 2) - 1);
 201       while (mpz_cmp_ui (x, 1L) <= 0);
 202
 203       is_prime = millerrabin (n, nm1, x, y, q, k);
 204     }
 205
 206   gmp_randclear (rstate);
 207
 208   TMP_FREE (marker);
 209   return is_prime;
 210 }
 211
 212 static int
 213 #if __STDC__
 214 millerrabin (mpz_srcptr n, mpz_srcptr nm1, mpz_ptr x, mpz_ptr y,
 215              mpz_srcptr q, unsigned long int k)
 216 #else
 217 millerrabin (n, nm1, x, y, q, k)
 218      mpz_srcptr n;
 219      mpz_srcptr nm1;
 220      mpz_ptr x;
 221      mpz_ptr y;
 222      mpz_srcptr q;
 223      unsigned long int k;
 224 #endif
 225 {
 226   unsigned long int i;
 227
 228   mpz_powm (y, x, q, n);
 229
 230   if (mpz_cmp_ui (y, 1L) == 0 || mpz_cmp (y, nm1) == 0)
 231     return 1;
 232
 233   for (i = 1; i < k; i++)
 234     {
 235       mpz_powm_ui (y, y, 2L, n);
 236       if (mpz_cmp (y, nm1) == 0)
 237         return 1;
 238       if (mpz_cmp_ui (y, 1L) == 0)
 239         return 0;
 240     }
 241   return 0;
 242 }