bignum-next-prime.c 4.08 KB
Newer Older
1
2
3
4
5
6
/* bignum-next-prime.c
 *
 */

/* nettle, low-level cryptographics library
 *
Niels Möller's avatar
Niels Möller committed
7
 * Copyright (C) 2002 Niels Möller
8
9
10
11
12
13
14
15
16
17
18
19
20
 *  
 * The nettle library is free software; you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License as published by
 * the Free Software Foundation; either version 2.1 of the License, or (at your
 * option) any later version.
 * 
 * The nettle library is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
 * or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
 * License for more details.
 * 
 * You should have received a copy of the GNU Lesser General Public License
 * along with the nettle library; see the file COPYING.LIB.  If not, write to
21
22
 * the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
 * MA 02111-1301, USA.
23
24
25
26
27
28
29
 */

#if HAVE_CONFIG_H
# include "config.h"
#endif

#include <limits.h>
30
31
/* Needed for alloca on freebsd */
#include <stdlib.h>
32
33

#include "bignum.h"
34
#include "gmp-glue.h"
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50

/* From gmp.h */
/* Test for gcc >= maj.min, as per __GNUC_PREREQ in glibc */
#if defined (__GNUC__) && defined (__GNUC_MINOR__)
#define GNUC_PREREQ(maj, min) \
  ((__GNUC__ << 16) + __GNUC_MINOR__ >= ((maj) << 16) + (min))
#else
#define GNUC_PREREQ(maj, min)  0
#endif

#if GNUC_PREREQ (3,0)
# define UNLIKELY(cond) __builtin_expect ((cond) != 0, 0)
#else
# define UNLIKELY(cond) cond
#endif

Niels Möller's avatar
Niels Möller committed
51
52
/* From some benchmarking using the examples nextprime(200!) and
   nextprime(240!), it seems that it pays off to use a prime list up
53
54
55
56
57
58
59
60
61
   to around 5000--10000 primes. There are 6541 odd primes less than
   2^16. */
static const uint16_t primes[] = {
  /* Generated by

     ./examples/eratosthenes 65535 \
     | awk '{ if (NR % 10 == 2) printf ("\n"); if (NR > 1) printf("%d, ", $1); }
            END { printf("\n"); }' > prime-list.h
  */
62
63
64
65
66
  #include "prime-list.h"
};

#define NUMBER_OF_PRIMES (sizeof(primes) / sizeof(primes[0]))

67
68
69
70
71
#ifdef mpz_millerrabin
# define PRIME_P mpz_millerrabin
#else
# define PRIME_P mpz_probab_prime_p
#endif
72
73
74

/* NOTE: The mpz_nextprime in current GMP is unoptimized. */
void
75
nettle_next_prime(mpz_t p, mpz_t n, unsigned count, unsigned prime_limit,
76
		  void *progress_ctx, nettle_progress_func *progress)
77
78
79
{
  mpz_t tmp;
  unsigned difference;
80
  TMP_GMP_DECL(moduli, unsigned);
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106

  if (prime_limit > NUMBER_OF_PRIMES)
    prime_limit = NUMBER_OF_PRIMES;
  
  /* First handle tiny numbers */
  if (mpz_cmp_ui(n, 2) <= 0)
    {
      mpz_set_ui(p, 2);
      return;
    }
  mpz_set(p, n);
  mpz_setbit(p, 0);

  if (mpz_cmp_ui(p, 8) < 0)
    return;

  mpz_init(tmp);

  if (mpz_cmp_ui(p, primes[prime_limit-1]) <= 0)
    /* Use only 3, 5 and 7 */
    /* FIXME: Could do binary search in the table. */
    prime_limit = 3;
  
  /* Compute residues modulo small odd primes */
  /* FIXME: Could be sped up by collecting limb-sized products of the
     primes, to reduce the calls to mpz_fdiv_ui */
Niels Möller's avatar
Niels Möller committed
107
108
109
110
111
112

  /* FIXME: Could also handle the first few primes separately; compute
     the residue mod 15015 = 3 * 7 * 11 * 13, and tabulate the steps
     between the 5760 odd numbers in this interval that have no factor
     in common with 15015.
   */
113
114
  TMP_GMP_ALLOC(moduli, prime_limit);

115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
  {
    unsigned i;
    for (i = 0; i < prime_limit; i++)
      moduli[i] = mpz_fdiv_ui(p, primes[i]);
  }
  
  for (difference = 0; ; difference += 2)
    {
      int composite = 0;
      unsigned i;
      
      if (difference >= UINT_MAX - 10)
	{ /* Should not happen, at least not very often... */
	  mpz_add_ui(p, p, difference);
	  difference = 0;
	}

      /* First check residues */
      for (i = 0; i < prime_limit; i++)
	{
	  if (moduli[i] == 0)
	    composite = 1;
137

138
	  moduli[i] += 2;
139
	  if (UNLIKELY(moduli[i] >= primes[i]))
140
141
142
143
144
145
146
147
148
149
150
151
	    moduli[i] -= primes[i];
	}
      if (composite)
	continue;
      
      mpz_add_ui(p, p, difference);
      difference = 0;

      if (progress)
	progress(progress_ctx, '.');

      /* Miller-Rabin test */
152
      if (PRIME_P(p, count))
153
154
	break;

155
#if 0
156
157
      if (progress)
	progress(progress_ctx, '*');
158
#endif
159
160
    }
  mpz_clear(tmp);
161
  TMP_GMP_FREE(moduli);
162
}