Return to nextprime.c CVS log

Up to [local] / OpenXM_contrib / gmp / mpz

Annotation of OpenXM_contrib/gmp/mpz/nextprime.c, Revision 1.1.1.2

1.1       maekawa     1: /* mpz_nextprime(p,t) - compute the next prime > t and store that in p.
                      2:
1.1.1.2 ! ohara       3: Copyright 1999, 2000, 2001 Free Software Foundation, Inc.
1.1       maekawa     4:
                      5: This file is part of the GNU MP Library.
                      6:
                      7: The GNU MP Library is free software; you can redistribute it and/or modify
                      8: it under the terms of the GNU Lesser General Public License as published by
                      9: the Free Software Foundation; either version 2.1 of the License, or (at your
                     10: option) any later version.
                     11:
                     12: The GNU MP Library is distributed in the hope that it will be useful, but
                     13: WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
                     14: or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
                     15: License for more details.
                     16:
                     17: You should have received a copy of the GNU Lesser General Public License
                     18: along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
                     19: the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
                     20: MA 02111-1307, USA. */
                     21:
                     22: #include "gmp.h"
                     23: #include "gmp-impl.h"
                     24:
                     25: void
                     26: mpz_nextprime (mpz_ptr p, mpz_srcptr t)
                     27: {
                     28:   mpz_add_ui (p, t, 1L);
                     29:   while (! mpz_probab_prime_p (p, 5))
                     30:     mpz_add_ui (p, p, 1L);
                     31: }
                     32:
                     33: #if 0
1.1.1.2 ! ohara      34: /* This code is not yet tested.  Will be enabled some time. */
1.1       maekawa    35:
                     36: status unsigned short primes[] =
                     37: {
                     38: 3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,
                     39: 101,103,107,109,113,127,131,137,139,149,151,157,163,167,173,179,181,
                     40: 191,193,197,199,211,223,227,229,233,239,241,251,257,263,269,271,277,
                     41: 281,283,293,307,311,313,317,331,337,347,349,353,359,367,373,379,383,
                     42: 389,397,401,409,419,421,431,433,439,443,449,457,461,463,467,479,487,
                     43: 491,499,503,509,521,523,541,547,557,563,569,571,577,587,593,599,601,
                     44: 607,613,617,619,631,641,643,647,653,659,661,673,677,683,691,701,709,
                     45: 719,727,733,739,743,751,757,761,769,773,787,797,809,811,821,823,827,
                     46: 829,839,853,857,859,863,877,881,883,887,907,911,919,929,937,941,947,
                     47: 953,967,971,977,983,991,997
                     48: };
                     49:
                     50: #define NUMBER_OF_PRIMES 167
                     51:
                     52: void
                     53: mpz_nextprime (mpz_ptr p, mpz_srcptr n)
                     54: {
                     55:   mpz_t tmp;
                     56:   unsigned short *moduli;
                     57:   unsigned long difference;
                     58:   int i;
                     59:   int composite;
                     60:
                     61:   /* First handle tiny numbers */
                     62:   if (mpz_cmp_ui (n, 2) < 0)
                     63:     {
                     64:       mpz_set_ui (p, 2);
                     65:       return;
                     66:     }
                     67:   mpz_add_ui (p, n, 1);
                     68:   mpz_setbit (p, 0);
                     69:
                     70:   if (mpz_cmp_ui (p, 7) <= 0)
                     71:     return;
                     72:
                     73:   prime_limit = NUMBER_OF_PRIMES - 1;
                     74:   if (mpz_cmp_ui (p, primes[prime_limit]) <= 0)
                     75:     /* Just use first three entries (3,5,7) of table for small numbers */
                     76:     prime_limit = 3;
                     77:   if (prime_limit)
                     78:     {
                     79:       /* Compute residues modulo small odd primes */
                     80:       moduli = (unsigned short *) TMP_ALLOC (prime_limit * sizeof moduli[0]);
                     81:       for (i = 0; i < prime_limit; i++)
                     82:        moduli[i] = mpz_fdiv_ui (p, primes[i]);
                     83:     }
                     84:   for (difference = 0; ; difference += 2)
                     85:     {
                     86:       composite = 0;
                     87:
                     88:       /* First check residues */
                     89:       for (i = 0; i < prime_limit; i++)
                     90:        {
                     91:          int acc, pr;
                     92:          composite |= (moduli[i] == 0);
                     93:          acc = moduli[i] + 2;
                     94:          pr = primes[i];
                     95:          moduli[i] = acc >= pr ? acc - pr : acc;
                     96:        }
                     97:       if (composite)
                     98:        continue;
                     99:
                    100:       mpz_add_ui (p, p, difference);
                    101:       difference = 0;
                    102:
                    103:       /* Miller-Rabin test */
                    104:       if (mpz_millerrabin (p, 2))
                    105:        break;
                    106:     }
                    107: }
                    108: #endif