[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.1.1 and 1.1.1.2

version 1.1.1.1, 2000/01/10 15:35:27 version 1.1.1.2, 2000/09/09 14:12:56
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 (C) 1991, 1993, 1994, 1996, 1997, 1998, 1999, 2000 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)  static int mpz_millerrabin _PROTO ((mpz_srcptr n, int reps));
   
   int
   #if __STDC__
   mpz_probab_prime_p (mpz_srcptr n, int reps)
   #else
   mpz_probab_prime_p (n, reps)
      mpz_srcptr n;       mpz_srcptr n;
      mpz_srcptr n_minus_1;       int reps;
      mpz_ptr x;  #endif
      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)    /* Check if n has small factors.  */
     return 1;    if (UDIV_TIME > (2 * UMUL_TIME + 6))
       r = mpn_preinv_mod_1 (PTR(n), SIZ(n), (mp_limb_t) PP, (mp_limb_t) PP_INVERTED);
     else
       r = mpn_mod_1 (PTR(n), SIZ(n), (mp_limb_t) PP);
     if (r % 3 == 0 || r % 5 == 0 || r % 7 == 0 || r % 11 == 0 || r % 13 == 0
         || r % 17 == 0 || r % 19 == 0 || r % 23 == 0 || r % 29 == 0
   #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
         )
       {
         return 0;
       }
   
   for (i = 1; i < k; i++)    /* 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,
        and look for factors in the remainder.  */
     {
       unsigned long int ln2;
       unsigned long int q;
       mp_limb_t p1, p0, p;
       unsigned int primes[15];
       int nprimes;
   
       nprimes = 0;
       p = 1;
       ln2 = mpz_sizeinbase (n, 2) / 30; ln2 = ln2 * ln2;
       for (q = BITS_PER_MP_LIMB == 64 ? 59 : 31; 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)
                       {
                         if (mpn_mod_1 (PTR(n), SIZ(n), (mp_limb_t) primes[nprimes]) != 0)
                           abort ();
                         return 0;
                       }
                   p = q;
                   nprimes = 0;
                 }
               else
                 {
                   p = p0;
                 }
               primes[nprimes++] = q;
             }
         }
     }
   
     /* Perform a number of Miller-Rabin tests.  */
     return mpz_millerrabin (n, reps);
   }
   
   static int
   #if __STDC__
   isprime (unsigned long int t)
   #else
   isprime (t)
        unsigned long int t;
   #endif
   {
     unsigned long int q, r, d;
   
     if (t < 3 || (t & 1) == 0)
       return t == 2;
   
     for (d = 3, r = 1; r != 0; d += 2)
     {      {
       mpz_powm_ui (y, y, 2L, n);        q = t / d;
       if (mpz_cmp (y, n_minus_1) == 0)        r = t - q * d;
         if (q < d)
         return 1;          return 1;
       if (mpz_cmp_ui (y, 1L) == 0)  
         return 0;  
     }      }
   return 0;    return 0;
 }  }
   
 int  static int millerrabin _PROTO ((mpz_srcptr n, mpz_srcptr nm1,
                                   mpz_ptr x, mpz_ptr y,
                                   mpz_srcptr q, unsigned long int k));
   
   static int
 #if __STDC__  #if __STDC__
 mpz_probab_prime_p (mpz_srcptr m, int reps)  mpz_millerrabin (mpz_srcptr n, int reps)
 #else  #else
 mpz_probab_prime_p (m, reps)  mpz_millerrabin (n, reps)
      mpz_srcptr m;       mpz_srcptr n;
      int reps;       int reps;
 #endif  #endif
 {  {
   mpz_t n, n_minus_1, x, y, q;    int r;
   int i, is_prime;    mpz_t nm1, x, y, q;
   unsigned long int k;    unsigned long int k;
     gmp_randstate_t rstate;
     int is_prime;
     TMP_DECL (marker);
     TMP_MARK (marker);
   
   mpz_init (n);    MPZ_TMP_INIT (nm1, SIZ (n) + 1);
   /* Take the absolute value of M, to handle positive and negative primes.  */    mpz_sub_ui (nm1, n, 1L);
   mpz_abs (n, m);  
   
   if (mpz_cmp_ui (n, 3L) <= 0)    MPZ_TMP_INIT (x, SIZ (n));
     {    MPZ_TMP_INIT (y, 2 * SIZ (n)); /* mpz_powm_ui needs excessive memory!!! */
       mpz_clear (n);  
       return mpz_cmp_ui (n, 1L) > 0;  
     }  
   
   if ((mpz_get_ui (n) & 1) == 0)    /* Perform a Fermat test.  */
     mpz_set_ui (x, 210L);
     mpz_powm (y, x, nm1, n);
     if (mpz_cmp_ui (y, 1L) != 0)
     {      {
       mpz_clear (n);        return 0;
       return 0;                 /* even */        TMP_FREE (marker);
     }      }
   
   mpz_init (n_minus_1);    MPZ_TMP_INIT (q, SIZ (n));
   mpz_sub_ui (n_minus_1, n, 1L);  
   mpz_init (x);  
   mpz_init (y);  
   
   /* find q and k, s.t.  n = 1 + 2**k * q */    /* Find q and k, where q is odd and n = 1 + 2**k * q.  */
   mpz_init_set (q, n_minus_1);    k = mpz_scan1 (nm1, 0L);
   k = mpz_scan1 (q, 0);    mpz_tdiv_q_2exp (q, nm1, k);
   mpz_tdiv_q_2exp (q, q, k);  
   
     gmp_randinit (rstate, GMP_RAND_ALG_DEFAULT, 32L);
   
   is_prime = 1;    is_prime = 1;
   for (i = 0; i < reps && is_prime; i++)    for (r = 0; r < reps && is_prime; r++)
     is_prime &= possibly_prime (n, n_minus_1, x, y, q, k);      {
         do
           mpz_urandomb (x, rstate, mpz_sizeinbase (n, 2) - 1);
         while (mpz_cmp_ui (x, 1L) <= 0);
   
   mpz_clear (n_minus_1);        is_prime = millerrabin (n, nm1, x, y, q, k);
   mpz_clear (n);      }
   mpz_clear (x);  
   mpz_clear (y);    gmp_randclear (rstate);
   mpz_clear (q);  
     TMP_FREE (marker);
   return is_prime;    return is_prime;
   }
   
   static int
   #if __STDC__
   millerrabin (mpz_srcptr n, mpz_srcptr nm1, mpz_ptr x, mpz_ptr y,
                mpz_srcptr q, unsigned long int k)
   #else
   millerrabin (n, nm1, x, y, q, k)
        mpz_srcptr n;
        mpz_srcptr nm1;
        mpz_ptr x;
        mpz_ptr y;
        mpz_srcptr q;
        unsigned long int k;
   #endif
   {
     unsigned long int i;
   
     mpz_powm (y, x, q, n);
   
     if (mpz_cmp_ui (y, 1L) == 0 || mpz_cmp (y, nm1) == 0)
       return 1;
   
     for (i = 1; i < k; i++)
       {
         mpz_powm_ui (y, y, 2L, n);
         if (mpz_cmp (y, nm1) == 0)
           return 1;
         if (mpz_cmp_ui (y, 1L) == 0)
           return 0;
       }
     return 0;
 }  }
Legend:
Removed from v.1.1.1.1  
changed lines
  Added in v.1.1.1.2
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>