[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.3 and 1.1.1.4

version 1.1.1.3, 2000/12/01 05:45:11 version 1.1.1.4, 2003/08/25 16:06:33
Line 6 
Line 6 
    positive is (1/4)**reps, where reps is the number of internal passes of the     positive is (1/4)**reps, where reps is the number of internal passes of the
    probabilistic algorithm.  Knuth indicates that 25 passes are reasonable.     probabilistic algorithm.  Knuth indicates that 25 passes are reasonable.
   
 Copyright (C) 1991, 1993, 1994, 1996, 1997, 1998, 1999, 2000 Free Software  Copyright 1991, 1993, 1994, 1996, 1997, 1998, 1999, 2000, 2001 Free Software
 Foundation, Inc.  Miller-Rabin code 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.
Line 31  MA 02111-1307, USA. */
Line 31  MA 02111-1307, USA. */
 #include "longlong.h"  #include "longlong.h"
   
 static int isprime _PROTO ((unsigned long int t));  static int isprime _PROTO ((unsigned long int t));
 static int mpz_millerrabin _PROTO ((mpz_srcptr n, int reps));  
   
 int  int
 #if __STDC__  
 mpz_probab_prime_p (mpz_srcptr n, int reps)  mpz_probab_prime_p (mpz_srcptr n, int reps)
 #else  
 mpz_probab_prime_p (n, reps)  
      mpz_srcptr n;  
      int reps;  
 #endif  
 {  {
   mp_limb_t r;    mp_limb_t r;
   
Line 66  mpz_probab_prime_p (n, reps)
Line 59  mpz_probab_prime_p (n, reps)
   if ((mpz_get_ui (n) & 1) == 0)    if ((mpz_get_ui (n) & 1) == 0)
     return 0;      return 0;
   
   #if defined (PP)
   /* Check if n has small factors.  */    /* Check if n has small factors.  */
   if (UDIV_TIME > (2 * UMUL_TIME + 6))  #if defined (PP_INVERTED)
     r = mpn_preinv_mod_1 (PTR(n), SIZ(n), (mp_limb_t) PP, (mp_limb_t) PP_INVERTED);    r = MPN_MOD_OR_PREINV_MOD_1 (PTR(n), SIZ(n), (mp_limb_t) PP,
   else                                 (mp_limb_t) PP_INVERTED);
     r = mpn_mod_1 (PTR(n), SIZ(n), (mp_limb_t) PP);  #else
   if (r % 3 == 0 || r % 5 == 0 || r % 7 == 0 || r % 11 == 0 || r % 13 == 0    r = mpn_mod_1 (PTR(n), SIZ(n), (mp_limb_t) PP);
   #endif
     if (r % 3 == 0
   #if BITS_PER_MP_LIMB >= 4
         || r % 5 == 0
   #endif
   #if BITS_PER_MP_LIMB >= 8
         || r % 7 == 0
   #endif
   #if BITS_PER_MP_LIMB >= 16
         || r % 11 == 0 || r % 13 == 0
   #endif
   #if BITS_PER_MP_LIMB >= 32
       || r % 17 == 0 || r % 19 == 0 || r % 23 == 0 || r % 29 == 0        || r % 17 == 0 || r % 19 == 0 || r % 23 == 0 || r % 29 == 0
 #if BITS_PER_MP_LIMB == 64  #endif
   #if BITS_PER_MP_LIMB >= 64
       || r % 31 == 0 || r % 37 == 0 || r % 41 == 0 || r % 43 == 0        || r % 31 == 0 || r % 37 == 0 || r % 41 == 0 || r % 43 == 0
       || r % 47 == 0 || r % 53 == 0        || r % 47 == 0 || r % 53 == 0
 #endif  #endif
Line 81  mpz_probab_prime_p (n, reps)
Line 88  mpz_probab_prime_p (n, reps)
     {      {
       return 0;        return 0;
     }      }
   #endif /* PP */
   
   /* Do more dividing.  We collect small primes, using umul_ppmm, until we    /* 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,       overflow a single limb.  We divide our number by the small primes product,
Line 95  mpz_probab_prime_p (n, reps)
Line 103  mpz_probab_prime_p (n, reps)
     nprimes = 0;      nprimes = 0;
     p = 1;      p = 1;
     ln2 = mpz_sizeinbase (n, 2) / 30; ln2 = ln2 * ln2;      ln2 = mpz_sizeinbase (n, 2) / 30; ln2 = ln2 * ln2;
     for (q = BITS_PER_MP_LIMB == 64 ? 59 : 31; q < ln2; q += 2)      for (q = PP_FIRST_OMITTED; q < ln2; q += 2)
       {        {
         if (isprime (q))          if (isprime (q))
           {            {
Line 106  mpz_probab_prime_p (n, reps)
Line 114  mpz_probab_prime_p (n, reps)
                 while (--nprimes >= 0)                  while (--nprimes >= 0)
                   if (r % primes[nprimes] == 0)                    if (r % primes[nprimes] == 0)
                     {                      {
                       if (mpn_mod_1 (PTR(n), SIZ(n), (mp_limb_t) primes[nprimes]) != 0)                        ASSERT_ALWAYS (mpn_mod_1 (PTR(n), SIZ(n), (mp_limb_t) primes[nprimes]) == 0);
                         abort ();  
                       return 0;                        return 0;
                     }                      }
                 p = q;                  p = q;
Line 127  mpz_probab_prime_p (n, reps)
Line 134  mpz_probab_prime_p (n, reps)
 }  }
   
 static int  static int
 #if __STDC__  
 isprime (unsigned long int t)  isprime (unsigned long int t)
 #else  
 isprime (t)  
      unsigned long int t;  
 #endif  
 {  {
   unsigned long int q, r, d;    unsigned long int q, r, d;
   
Line 145  isprime (t)
Line 147  isprime (t)
       r = t - q * d;        r = t - q * d;
       if (q < d)        if (q < d)
         return 1;          return 1;
     }  
   return 0;  
 }  
   
 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__  
 mpz_millerrabin (mpz_srcptr n, int reps)  
 #else  
 mpz_millerrabin (n, reps)  
      mpz_srcptr n;  
      int reps;  
 #endif  
 {  
   int r;  
   mpz_t nm1, x, y, q;  
   unsigned long int k;  
   gmp_randstate_t rstate;  
   int is_prime;  
   TMP_DECL (marker);  
   TMP_MARK (marker);  
   
   MPZ_TMP_INIT (nm1, SIZ (n) + 1);  
   mpz_sub_ui (nm1, n, 1L);  
   
   MPZ_TMP_INIT (x, SIZ (n));  
   MPZ_TMP_INIT (y, 2 * SIZ (n)); /* mpz_powm_ui needs excessive memory!!! */  
   
   /* Perform a Fermat test.  */  
   mpz_set_ui (x, 210L);  
   mpz_powm (y, x, nm1, n);  
   if (mpz_cmp_ui (y, 1L) != 0)  
     {  
       TMP_FREE (marker);  
       return 0;  
     }  
   
   MPZ_TMP_INIT (q, SIZ (n));  
   
   /* Find q and k, where q is odd and n = 1 + 2**k * q.  */  
   k = mpz_scan1 (nm1, 0L);  
   mpz_tdiv_q_2exp (q, nm1, k);  
   
   gmp_randinit (rstate, GMP_RAND_ALG_DEFAULT, 32L);  
   
   is_prime = 1;  
   for (r = 0; r < reps && is_prime; r++)  
     {  
       do  
         mpz_urandomb (x, rstate, mpz_sizeinbase (n, 2) - 1);  
       while (mpz_cmp_ui (x, 1L) <= 0);  
   
       is_prime = millerrabin (n, nm1, x, y, q, k);  
     }  
   
   gmp_randclear (rstate);  
   
   TMP_FREE (marker);  
   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;    return 0;
 }  }
Legend:
Removed from v.1.1.1.3  
changed lines
  Added in v.1.1.1.4
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>