bignum-next-prime.c

   1 /* bignum-next-prime.c
   2  *
   3  */
   4
   5 /* nettle, low-level cryptographics library
   6  *
   7  * Copyright (C) 2002 Niels Möller
   8  *
   9  * The nettle library is free software; you can redistribute it and/or modify
  10  * it under the terms of the GNU Lesser General Public License as published by
  11  * the Free Software Foundation; either version 2.1 of the License, or (at your
  12  * option) any later version.
  13  *
  14  * The nettle library is distributed in the hope that it will be useful, but
  15  * WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  16  * or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  17  * License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public License
  20  * along with the nettle library; see the file COPYING.LIB.  If not, write to
  21  * the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
  22  * MA 02111-1301, USA.
  23  */
  24
  25 #if HAVE_CONFIG_H
  26 # include "config.h"
  27 #endif
  28
  29 #include <limits.h>
  30 /* Needed for alloca on freebsd */
  31 #include <stdlib.h>
  32
  33 #include "bignum.h"
  34
  35 #include "nettle-internal.h"
  36
  37 /* From gmp.h */
  38 /* Test for gcc >= maj.min, as per __GNUC_PREREQ in glibc */
  39 #if defined (__GNUC__) && defined (__GNUC_MINOR__)
  40 #define GNUC_PREREQ(maj, min) \
  41   ((__GNUC__ << 16) + __GNUC_MINOR__ >= ((maj) << 16) + (min))
  42 #else
  43 #define GNUC_PREREQ(maj, min)  0
  44 #endif
  45
  46 #if GNUC_PREREQ (3,0)
  47 # define UNLIKELY(cond) __builtin_expect ((cond) != 0, 0)
  48 #else
  49 # define UNLIKELY(cond) cond
  50 #endif
  51
  52 /* From some benchmarking using the examples nextprime(200!) and
  53    nextprime(240!), it seems that it pays off to use a prime list up
  54    to around 5000--10000 primes. There are 6541 odd primes less than
  55    2^16. */
  56 static const uint16_t primes[] = {
  57   /* Generated by
  58
  59      ./examples/eratosthenes 65535 \
  60      | awk '{ if (NR % 10 == 2) printf ("\n"); if (NR > 1) printf("%d, ", $1); }
  61             END { printf("\n"); }' > prime-list.h
  62   */
  63   #include "prime-list.h"
  64 };
  65
  66 #define NUMBER_OF_PRIMES (sizeof(primes) / sizeof(primes[0]))
  67
  68 #ifdef mpz_millerrabin
  69 # define PRIME_P mpz_millerrabin
  70 #else
  71 # define PRIME_P mpz_probab_prime_p
  72 #endif
  73
  74 /* NOTE: The mpz_nextprime in current GMP is unoptimized. */
  75 void
  76 nettle_next_prime(mpz_t p, mpz_t n, unsigned count, unsigned prime_limit,
  77                   void *progress_ctx, nettle_progress_func *progress)
  78 {
  79   mpz_t tmp;
  80   TMP_DECL(moduli, unsigned, NUMBER_OF_PRIMES);
  81
  82   unsigned difference;
  83
  84   if (prime_limit > NUMBER_OF_PRIMES)
  85     prime_limit = NUMBER_OF_PRIMES;
  86
  87   /* First handle tiny numbers */
  88   if (mpz_cmp_ui(n, 2) <= 0)
  89     {
  90       mpz_set_ui(p, 2);
  91       return;
  92     }
  93   mpz_set(p, n);
  94   mpz_setbit(p, 0);
  95
  96   if (mpz_cmp_ui(p, 8) < 0)
  97     return;
  98
  99   mpz_init(tmp);
 100
 101   if (mpz_cmp_ui(p, primes[prime_limit-1]) <= 0)
 102     /* Use only 3, 5 and 7 */
 103     /* FIXME: Could do binary search in the table. */
 104     prime_limit = 3;
 105
 106   /* Compute residues modulo small odd primes */
 107   /* FIXME: Could be sped up by collecting limb-sized products of the
 108      primes, to reduce the calls to mpz_fdiv_ui */
 109
 110   /* FIXME: Could also handle the first few primes separately; compute
 111      the residue mod 15015 = 3 * 7 * 11 * 13, and tabulate the steps
 112      between the 5760 odd numbers in this interval that have no factor
 113      in common with 15015.
 114    */
 115   TMP_ALLOC(moduli, prime_limit);
 116   {
 117     unsigned i;
 118     for (i = 0; i < prime_limit; i++)
 119       moduli[i] = mpz_fdiv_ui(p, primes[i]);
 120   }
 121
 122   for (difference = 0; ; difference += 2)
 123     {
 124       int composite = 0;
 125       unsigned i;
 126
 127       if (difference >= UINT_MAX - 10)
 128         { /* Should not happen, at least not very often... */
 129           mpz_add_ui(p, p, difference);
 130           difference = 0;
 131         }
 132
 133       /* First check residues */
 134       for (i = 0; i < prime_limit; i++)
 135         {
 136           if (moduli[i] == 0)
 137             composite = 1;
 138
 139           moduli[i] += 2;
 140           if (UNLIKELY(moduli[i] >= primes[i]))
 141             moduli[i] -= primes[i];
 142         }
 143       if (composite)
 144         continue;
 145
 146       mpz_add_ui(p, p, difference);
 147       difference = 0;
 148
 149       if (progress)
 150         progress(progress_ctx, '.');
 151
 152       /* Miller-Rabin test */
 153       if (PRIME_P(p, count))
 154         break;
 155
 156 #if 0
 157       if (progress)
 158         progress(progress_ctx, '*');
 159 #endif
 160     }
 161   mpz_clear(tmp);
 162 }