[BACK]Return to pprime_p.c CVS log [TXT][DIR] Up to [local] / OpenXM_contrib / gmp / mpz

Diff for /OpenXM_contrib/gmp/mpz/Attic/pprime_p.c between version 1.1 and 1.1.1.4

version 1.1, 2000/01/10 15:35:27 version 1.1.1.4, 2003/08/25 16:06:33
Line 2 
Line 2 
    An implementation of the probabilistic primality test found in Knuth's     An implementation of the probabilistic primality test found in Knuth's
    Seminumerical Algorithms book.  If the function mpz_probab_prime_p()     Seminumerical Algorithms book.  If the function mpz_probab_prime_p()
    returns 0 then n is not prime.  If it returns 1, then n is 'probably'     returns 0 then n is not prime.  If it returns 1, then n is 'probably'
    prime.  The probability of a false positive is (1/4)**reps, where     prime.  If it returns 2, n is surely prime.  The probability of a false
    reps is the number of internal passes of the probabilistic algorithm.     positive is (1/4)**reps, where reps is the number of internal passes of the
    Knuth indicates that 25 passes are reasonable.     probabilistic algorithm.  Knuth indicates that 25 passes are reasonable.
   
 Copyright (C) 1991, 1993, 1994 Free Software Foundation, Inc.  Copyright 1991, 1993, 1994, 1996, 1997, 1998, 1999, 2000, 2001 Free Software
 Contributed by John Amanatides.  Foundation, Inc.  Miller-Rabin code contributed by John Amanatides.
   
 This file is part of the GNU MP Library.  This file is part of the GNU MP Library.
   
 The GNU MP Library is free software; you can redistribute it and/or modify  The GNU MP Library is free software; you can redistribute it and/or modify
 it under the terms of the GNU Library General Public License as published by  it under the terms of the GNU Lesser General Public License as published by
 the Free Software Foundation; either version 2 of the License, or (at your  the Free Software Foundation; either version 2.1 of the License, or (at your
 option) any later version.  option) any later version.
   
 The GNU MP Library is distributed in the hope that it will be useful, but  The GNU MP Library is distributed in the hope that it will be useful, but
 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public  or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
 License for more details.  License for more details.
   
 You should have received a copy of the GNU Library General Public License  You should have received a copy of the GNU Lesser General Public License
 along with the GNU MP Library; see the file COPYING.LIB.  If not, write to  along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
 the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,  the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
 MA 02111-1307, USA. */  MA 02111-1307, USA. */
   
 #include "gmp.h"  #include "gmp.h"
   #include "gmp-impl.h"
   #include "longlong.h"
   
 static int  static int isprime _PROTO ((unsigned long int t));
 possibly_prime (n, n_minus_1, x, y, q, k)  
      mpz_srcptr n;  int
      mpz_srcptr n_minus_1;  mpz_probab_prime_p (mpz_srcptr n, int reps)
      mpz_ptr x;  
      mpz_ptr y;  
      mpz_srcptr q;  
      unsigned long int k;  
 {  {
   unsigned long int i;    mp_limb_t r;
   
   /* find random x s.t. 1 < x < n */    /* Handle small and negative n.  */
   do    if (mpz_cmp_ui (n, 1000000L) <= 0)
     {      {
       mpz_random (x, mpz_size (n));        int is_prime;
       mpz_mmod (x, x, n);        if (mpz_sgn (n) < 0)
           {
             /* Negative number.  Negate and call ourselves.  */
             mpz_t n2;
             mpz_init (n2);
             mpz_neg (n2, n);
             is_prime = mpz_probab_prime_p (n2, reps);
             mpz_clear (n2);
             return is_prime;
           }
         is_prime = isprime (mpz_get_ui (n));
         return is_prime ? 2 : 0;
     }      }
   while (mpz_cmp_ui (x, 1L) <= 0);  
   
   mpz_powm (y, x, q, n);    /* If n is now even, it is not a prime.  */
     if ((mpz_get_ui (n) & 1) == 0)
       return 0;
   
   if (mpz_cmp_ui (y, 1L) == 0 || mpz_cmp (y, n_minus_1) == 0)  #if defined (PP)
     return 1;    /* Check if n has small factors.  */
   #if defined (PP_INVERTED)
   for (i = 1; i < k; i++)    r = MPN_MOD_OR_PREINV_MOD_1 (PTR(n), SIZ(n), (mp_limb_t) PP,
     {                                 (mp_limb_t) PP_INVERTED);
       mpz_powm_ui (y, y, 2L, n);  
       if (mpz_cmp (y, n_minus_1) == 0)  
         return 1;  
       if (mpz_cmp_ui (y, 1L) == 0)  
         return 0;  
     }  
   return 0;  
 }  
   
 int  
 #if __STDC__  
 mpz_probab_prime_p (mpz_srcptr m, int reps)  
 #else  #else
 mpz_probab_prime_p (m, reps)    r = mpn_mod_1 (PTR(n), SIZ(n), (mp_limb_t) PP);
      mpz_srcptr m;  
      int reps;  
 #endif  #endif
 {    if (r % 3 == 0
   mpz_t n, n_minus_1, x, y, q;  #if BITS_PER_MP_LIMB >= 4
   int i, is_prime;        || r % 5 == 0
   unsigned long int k;  #endif
   #if BITS_PER_MP_LIMB >= 8
   mpz_init (n);        || r % 7 == 0
   /* Take the absolute value of M, to handle positive and negative primes.  */  #endif
   mpz_abs (n, m);  #if BITS_PER_MP_LIMB >= 16
         || r % 11 == 0 || r % 13 == 0
   if (mpz_cmp_ui (n, 3L) <= 0)  #endif
   #if BITS_PER_MP_LIMB >= 32
         || r % 17 == 0 || r % 19 == 0 || r % 23 == 0 || r % 29 == 0
   #endif
   #if BITS_PER_MP_LIMB >= 64
         || r % 31 == 0 || r % 37 == 0 || r % 41 == 0 || r % 43 == 0
         || r % 47 == 0 || r % 53 == 0
   #endif
         )
     {      {
       mpz_clear (n);        return 0;
       return mpz_cmp_ui (n, 1L) > 0;  
     }      }
   #endif /* PP */
   
   if ((mpz_get_ui (n) & 1) == 0)    /* Do more dividing.  We collect small primes, using umul_ppmm, until we
     {       overflow a single limb.  We divide our number by the small primes product,
       mpz_clear (n);       and look for factors in the remainder.  */
       return 0;                 /* even */    {
     }      unsigned long int ln2;
       unsigned long int q;
       mp_limb_t p1, p0, p;
       unsigned int primes[15];
       int nprimes;
   
   mpz_init (n_minus_1);      nprimes = 0;
   mpz_sub_ui (n_minus_1, n, 1L);      p = 1;
   mpz_init (x);      ln2 = mpz_sizeinbase (n, 2) / 30; ln2 = ln2 * ln2;
   mpz_init (y);      for (q = PP_FIRST_OMITTED; q < ln2; q += 2)
         {
           if (isprime (q))
             {
               umul_ppmm (p1, p0, p, q);
               if (p1 != 0)
                 {
                   r = mpn_mod_1 (PTR(n), SIZ(n), p);
                   while (--nprimes >= 0)
                     if (r % primes[nprimes] == 0)
                       {
                         ASSERT_ALWAYS (mpn_mod_1 (PTR(n), SIZ(n), (mp_limb_t) primes[nprimes]) == 0);
                         return 0;
                       }
                   p = q;
                   nprimes = 0;
                 }
               else
                 {
                   p = p0;
                 }
               primes[nprimes++] = q;
             }
         }
     }
   
   /* find q and k, s.t.  n = 1 + 2**k * q */    /* Perform a number of Miller-Rabin tests.  */
   mpz_init_set (q, n_minus_1);    return mpz_millerrabin (n, reps);
   k = mpz_scan1 (q, 0);  }
   mpz_tdiv_q_2exp (q, q, k);  
   
   is_prime = 1;  static int
   for (i = 0; i < reps && is_prime; i++)  isprime (unsigned long int t)
     is_prime &= possibly_prime (n, n_minus_1, x, y, q, k);  {
     unsigned long int q, r, d;
   
   mpz_clear (n_minus_1);    if (t < 3 || (t & 1) == 0)
   mpz_clear (n);      return t == 2;
   mpz_clear (x);  
   mpz_clear (y);    for (d = 3, r = 1; r != 0; d += 2)
   mpz_clear (q);      {
   return is_prime;        q = t / d;
         r = t - q * d;
         if (q < d)
           return 1;
       }
     return 0;
 }  }
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.1.1.4
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>