Annotation of OpenXM_contrib/gmp/mpz/nextprime.c, Revision 1.1
1.1 ! maekawa 1: /* mpz_nextprime(p,t) - compute the next prime > t and store that in p.
! 2:
! 3: Copyright (C) 1999, 2000 Free Software Foundation, Inc.
! 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: #if __STDC__
! 27: mpz_nextprime (mpz_ptr p, mpz_srcptr t)
! 28: #else
! 29: mpz_nextprime (p, t)
! 30: mpz_ptr p;
! 31: mpz_srcptr t;
! 32: #endif
! 33: {
! 34: mpz_add_ui (p, t, 1L);
! 35: while (! mpz_probab_prime_p (p, 5))
! 36: mpz_add_ui (p, p, 1L);
! 37: }
! 38:
! 39: #if 0
! 40: /* This code is not yet tested. Will be enabled in 3.1. */
! 41:
! 42: status unsigned short primes[] =
! 43: {
! 44: 3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,
! 45: 101,103,107,109,113,127,131,137,139,149,151,157,163,167,173,179,181,
! 46: 191,193,197,199,211,223,227,229,233,239,241,251,257,263,269,271,277,
! 47: 281,283,293,307,311,313,317,331,337,347,349,353,359,367,373,379,383,
! 48: 389,397,401,409,419,421,431,433,439,443,449,457,461,463,467,479,487,
! 49: 491,499,503,509,521,523,541,547,557,563,569,571,577,587,593,599,601,
! 50: 607,613,617,619,631,641,643,647,653,659,661,673,677,683,691,701,709,
! 51: 719,727,733,739,743,751,757,761,769,773,787,797,809,811,821,823,827,
! 52: 829,839,853,857,859,863,877,881,883,887,907,911,919,929,937,941,947,
! 53: 953,967,971,977,983,991,997
! 54: };
! 55:
! 56: #define NUMBER_OF_PRIMES 167
! 57:
! 58: void
! 59: #if __STDC__
! 60: mpz_nextprime (mpz_ptr p, mpz_srcptr n)
! 61: #else
! 62: mpz_nextprime (p, n)
! 63: mpz_ptr p;
! 64: mpz_srcptr n;
! 65: #endif
! 66: {
! 67: mpz_t tmp;
! 68: unsigned short *moduli;
! 69: unsigned long difference;
! 70: int i;
! 71: int composite;
! 72:
! 73: /* First handle tiny numbers */
! 74: if (mpz_cmp_ui (n, 2) < 0)
! 75: {
! 76: mpz_set_ui (p, 2);
! 77: return;
! 78: }
! 79: mpz_add_ui (p, n, 1);
! 80: mpz_setbit (p, 0);
! 81:
! 82: if (mpz_cmp_ui (p, 7) <= 0)
! 83: return;
! 84:
! 85: prime_limit = NUMBER_OF_PRIMES - 1;
! 86: if (mpz_cmp_ui (p, primes[prime_limit]) <= 0)
! 87: /* Just use first three entries (3,5,7) of table for small numbers */
! 88: prime_limit = 3;
! 89: if (prime_limit)
! 90: {
! 91: /* Compute residues modulo small odd primes */
! 92: moduli = (unsigned short *) TMP_ALLOC (prime_limit * sizeof moduli[0]);
! 93: for (i = 0; i < prime_limit; i++)
! 94: moduli[i] = mpz_fdiv_ui (p, primes[i]);
! 95: }
! 96: for (difference = 0; ; difference += 2)
! 97: {
! 98: composite = 0;
! 99:
! 100: /* First check residues */
! 101: for (i = 0; i < prime_limit; i++)
! 102: {
! 103: int acc, pr;
! 104: composite |= (moduli[i] == 0);
! 105: acc = moduli[i] + 2;
! 106: pr = primes[i];
! 107: moduli[i] = acc >= pr ? acc - pr : acc;
! 108: }
! 109: if (composite)
! 110: continue;
! 111:
! 112: mpz_add_ui (p, p, difference);
! 113: difference = 0;
! 114:
! 115: /* Miller-Rabin test */
! 116: if (mpz_millerrabin (p, 2))
! 117: break;
! 118: }
! 119: }
! 120: #endif
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>