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.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