source/libs/gmp/gmp-src/mpz/nextprime.c

   1 /* mpz_nextprime(p,t) - compute the next prime > t and store that in p.
   2
   3 Copyright 1999-2001, 2008, 2009, 2012 Free Software Foundation, Inc.
   4
   5 Contributed to the GNU project by Niels Möller and Torbjorn Granlund.
   6
   7 This file is part of the GNU MP Library.
   8
   9 The GNU MP Library is free software; you can redistribute it and/or modify
  10 it under the terms of either:
  11
  12   * the GNU Lesser General Public License as published by the Free
  13     Software Foundation; either version 3 of the License, or (at your
  14     option) any later version.
  15
  16 or
  17
  18   * the GNU General Public License as published by the Free Software
  19     Foundation; either version 2 of the License, or (at your option) any
  20     later version.
  21
  22 or both in parallel, as here.
  23
  24 The GNU MP Library is distributed in the hope that it will be useful, but
  25 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  26 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  27 for more details.
  28
  29 You should have received copies of the GNU General Public License and the
  30 GNU Lesser General Public License along with the GNU MP Library.  If not,
  31 see https://www.gnu.org/licenses/.  */
  32
  33 #include "gmp.h"
  34 #include "gmp-impl.h"
  35 #include "longlong.h"
  36
  37 static const unsigned char primegap[] =
  38 {
  39   2,2,4,2,4,2,4,6,2,6,4,2,4,6,6,2,6,4,2,6,4,6,8,4,2,4,2,4,14,4,6,
  40   2,10,2,6,6,4,6,6,2,10,2,4,2,12,12,4,2,4,6,2,10,6,6,6,2,6,4,2,10,14,4,2,
  41   4,14,6,10,2,4,6,8,6,6,4,6,8,4,8,10,2,10,2,6,4,6,8,4,2,4,12,8,4,8,4,6,
  42   12,2,18,6,10,6,6,2,6,10,6,6,2,6,6,4,2,12,10,2,4,6,6,2,12,4,6,8,10,8,10,8,
  43   6,6,4,8,6,4,8,4,14,10,12,2,10,2,4,2,10,14,4,2,4,14,4,2,4,20,4,8,10,8,4,6,
  44   6,14,4,6,6,8,6,12
  45 };
  46
  47 #define NUMBER_OF_PRIMES 167
  48
  49 void
  50 mpz_nextprime (mpz_ptr p, mpz_srcptr n)
  51 {
  52   unsigned short *moduli;
  53   unsigned long difference;
  54   int i;
  55   unsigned prime_limit;
  56   unsigned long prime;
  57   mp_size_t pn;
  58   mp_bitcnt_t nbits;
  59   unsigned incr;
  60   TMP_SDECL;
  61
  62   /* First handle tiny numbers */
  63   if (mpz_cmp_ui (n, 2) < 0)
  64     {
  65       mpz_set_ui (p, 2);
  66       return;
  67     }
  68   mpz_add_ui (p, n, 1);
  69   mpz_setbit (p, 0);
  70
  71   if (mpz_cmp_ui (p, 7) <= 0)
  72     return;
  73
  74   pn = SIZ(p);
  75   MPN_SIZEINBASE_2EXP(nbits, PTR(p), pn, 1);
  76   if (nbits / 2 >= NUMBER_OF_PRIMES)
  77     prime_limit = NUMBER_OF_PRIMES - 1;
  78   else
  79     prime_limit = nbits / 2;
  80
  81   TMP_SMARK;
  82
  83   /* Compute residues modulo small odd primes */
  84   moduli = TMP_SALLOC_TYPE (prime_limit * sizeof moduli[0], unsigned short);
  85
  86   for (;;)
  87     {
  88       /* FIXME: Compute lazily? */
  89       prime = 3;
  90       for (i = 0; i < prime_limit; i++)
  91         {
  92           moduli[i] = mpz_fdiv_ui (p, prime);
  93           prime += primegap[i];
  94         }
  95
  96 #define INCR_LIMIT 0x10000      /* deep science */
  97
  98       for (difference = incr = 0; incr < INCR_LIMIT; difference += 2)
  99         {
 100           /* First check residues */
 101           prime = 3;
 102           for (i = 0; i < prime_limit; i++)
 103             {
 104               unsigned r;
 105               /* FIXME: Reduce moduli + incr and store back, to allow for
 106                  division-free reductions.  Alternatively, table primes[]'s
 107                  inverses (mod 2^16).  */
 108               r = (moduli[i] + incr) % prime;
 109               prime += primegap[i];
 110
 111               if (r == 0)
 112                 goto next;
 113             }
 114
 115           mpz_add_ui (p, p, difference);
 116           difference = 0;
 117
 118           /* Miller-Rabin test */
 119           if (mpz_millerrabin (p, 25))
 120             goto done;
 121         next:;
 122           incr += 2;
 123         }
 124       mpz_add_ui (p, p, difference);
 125       difference = 0;
 126     }
 127  done:
 128   TMP_SFREE;
 129 }