1/*========================================================================
2               Copyright (C) 1996-2001 by Jorn Lind-Nielsen
3                            All rights reserved
4
5    Permission is hereby granted, without written agreement and without
6    license or royalty fees, to use, reproduce, prepare derivative
7    works, distribute, and display this software and its documentation
8    for any purpose, provided that (1) the above copyright notice and
9    the following two paragraphs appear in all copies of the source code
10    and (2) redistributions, including without limitation binaries,
11    reproduce these notices in the supporting documentation. Substantial
12    modifications to this software may be copyrighted by their authors
13    and need not follow the licensing terms described here, provided
14    that the new terms are clearly indicated in all files where they apply.
15
16    IN NO EVENT SHALL JORN LIND-NIELSEN, OR DISTRIBUTORS OF THIS
17    SOFTWARE BE LIABLE TO ANY PARTY FOR DIRECT, INDIRECT, SPECIAL,
18    INCIDENTAL, OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE USE OF THIS
19    SOFTWARE AND ITS DOCUMENTATION, EVEN IF THE AUTHORS OR ANY OF THE
20    ABOVE PARTIES HAVE BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
21
22    JORN LIND-NIELSEN SPECIFICALLY DISCLAIM ANY WARRANTIES, INCLUDING,
23    BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND
24    FITNESS FOR A PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS
25    ON AN "AS IS" BASIS, AND THE AUTHORS AND DISTRIBUTORS HAVE NO
26    OBLIGATION TO PROVIDE MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR
27    MODIFICATIONS.
28========================================================================*/
29
30/*************************************************************************
31  $Header$
32  FILE:  prime.c
33  DESCR: Prime number calculations
34  AUTH:  Jorn Lind
35  DATE:  (C) feb 2001
36*************************************************************************/
37#include <stdlib.h>
38#include <stdio.h>
39#include <time.h>
40#include "prime.h"
41
42
43#define Random(i) ( (rand() % (i)) + 1 )
44#define isEven(src) (!((src) & 0x1))
45#define hasFactor(src,n) ( (((src)!=(n)) && ((src)%(n) == 0)) )
46#define BitIsSet(src,b) ( ((src) & (1<<(b))) != 0 )
47
48#define CHECKTIMES 20
49
50#if defined(BUDDYUINT64)
51 typedef BUDDYUINT64 UINT64;
52 #define BUILTIN64
53#elif defined(__GNUC__) || defined(__KCC)
54 typedef long long UINT64;
55 #define BUILTIN64
56#elif defined(_MSV_VER)
57 typedef unsigned _int64 UINT64;
58 #define BUILTIN64
59#else
60 typedef struct __UINT64
61 {
62   unsigned int hi;
63   unsigned int lo;
64 } UINT64;
65
66 #define MAX(a,b) ((a) > (b) ? (a) : (b))
67 #define GETCARRY(a,b) ( ((a)+(b)) < MAX((a),(b)) ? 1 : 0 )
68#endif
69
70
71#ifndef BUILTIN64
72/*************************************************************************
73   64 bit unsigned int arithmetics
74*************************************************************************/
75
76static UINT64 u64_mul(unsigned int x, unsigned int y)
77{
78  UINT64 res;
79  unsigned int yh = 0;
80  unsigned int yl = y;
81  int i;
82
83  res.lo = res.hi = 0;
84
85  for (i=0 ; i<32 ; ++i)
86  {
87    if (x & 0x1)
88    {
89      unsigned int carry = GETCARRY(res.lo,yl);
90      res.lo += yl;
91      res.hi += yh + carry;
92    }
93
94    yh = (yh << 1) | (yl & 0x80000000 ? 1 : 0);
95    yl = (yl << 1);
96
97    x >>= 1;
98  }
99
100  return res;
101}
102
103
104static void u64_shl(UINT64* a, unsigned int *carryOut)
105{
106  *carryOut = (*carryOut << 1) | (a->hi & 0x80000000 ? 0x1 : 0x0);
107  a->hi     = (a->hi << 1) | (a->lo & 0x80000000 ? 0x1 : 0x0);
108  a->lo     = (a->lo << 1);
109}
110
111
112static unsigned int u64_mod(UINT64 dividend, unsigned int divisor)
113{
114  unsigned int remainder = 0;
115  int i;
116
117  u64_shl(&dividend, &remainder);
118
119  for (i=0 ; i<64 ; ++i)
120  {
121    if (remainder >= divisor)
122      remainder -= divisor;
123
124    u64_shl(&dividend, &remainder);
125  }
126
127  return remainder >> 1;
128}
129#endif /* BUILTIN64 */
130
131#ifdef BUILTIN64
132#define u64_mulmod(a,b,c) ((unsigned int)( ((UINT64)a*(UINT64)b)%(UINT64)c ));
133#else
134#define u64_mulmod(a,b,c) u64_mod( u64_mul((a),(b)), (c) );
135#endif
136
137
138/*************************************************************************
139  Miller Rabin check
140*************************************************************************/
141
142static unsigned int numberOfBits(unsigned int src)
143{
144  unsigned int b;
145
146  if (src == 0)
147    return 0;
148
149  for (b=(sizeof(unsigned int)*8)-1 ; b>0 ; --b)
150    if (BitIsSet(src,b))
151      return b+1;
152
153  return 1;
154}
155
156
157
158static int isWitness(unsigned int witness, unsigned int src)
159{
160  unsigned int bitNum = numberOfBits(src-1)-1;
161  unsigned int d = 1;
162  int i;
163
164  for (i=bitNum ; i>=0 ; --i)
165  {
166    unsigned int x = d;
167
168    d = u64_mulmod(d,d,src);
169
170    if (d == 1  &&  x != 1  &&  x != src-1)
171      return 1;
172
173    if (BitIsSet(src-1,i))
174      d = u64_mulmod(d,witness,src);
175  }
176
177  return d != 1;
178}
179
180
181static int isMillerRabinPrime(unsigned int src)
182{
183  int n;
184
185  for (n=0 ; n<CHECKTIMES ; ++n)
186  {
187    unsigned int witness = Random(src-1);
188
189    if (isWitness(witness,src))
190      return 0;
191  }
192
193  return 1;
194}
195
196
197/*************************************************************************
198  Basic prime searching stuff
199*************************************************************************/
200
201static int hasEasyFactors(unsigned int src)
202{
203  return hasFactor(src, 3)
204      || hasFactor(src, 5)
205      || hasFactor(src, 7)
206      || hasFactor(src, 11)
207      || hasFactor(src, 13);
208}
209
210
211static int isPrime(unsigned int src)
212{
213  if (hasEasyFactors(src))
214    return 0;
215
216  return isMillerRabinPrime(src);
217}
218
219
220/*************************************************************************
221  External interface
222*************************************************************************/
223
224unsigned int bdd_prime_gte(unsigned int src)
225{
226  if (isEven(src))
227    --src;
228
229  while (!isPrime(src))
230    src -= 2;
231
232  return src;
233}
234
235
236unsigned int bdd_prime_lte(unsigned int src)
237{
238  if (isEven(src))
239     --src;
240
241  while (!isPrime(src))
242     src -= 2;
243
244  return src;
245}
246
247
248
249/*************************************************************************
250   Testing
251*************************************************************************/
252
253#if 0
254int main()
255{
256  printf("Nb0 = %u\n", numberOfBits(0));
257  printf("Nb1 = %u\n", numberOfBits(1));
258  printf("Nb2 = %u\n", numberOfBits(2));
259  printf("Nb3 = %u\n", numberOfBits(3));
260  printf("Nb5 = %u\n", numberOfBits(5));
261  printf("Nb9 = %u\n", numberOfBits(9));
262  printf("Nb15 = %u\n", numberOfBits(15));
263  printf("Nb17 = %u\n", numberOfBits(17));
264  return 0;
265}
266#endif
267
268
269#if 0
270void testMul(unsigned int a, unsigned int b)
271{
272  UINT64 x = u64_mul(a,b);
273  long long z1 = (long long)a * (long long)b;
274  long long z2 = ((long long)x.hi << 32) + (long long)x.lo;
275  if (z1 != z2)
276    printf("%d * %d = %lld,%lld\n", a, b, z1, z2);
277}
278
279
280void testMod(unsigned int a, unsigned int b, unsigned int c)
281{
282  UINT64 x = u64_mul(a,b);
283
284  long long z1 = (long long)a * (long long)b;
285  long long z2 = ((long long)x.hi << 32) + (long long)x.lo;
286  unsigned int m1 = z1 % c;
287  unsigned int m2 = u64_mod(x,c);
288
289  if (z1 != z2)
290    printf("%d * %d = %lld,%lld\n", a, b, z1, z2);
291
292  if (m1 != m2)
293    printf("%llu %% %u = %u,%u\n", z1, c, m1, m2);
294}
295#endif
296
297#if 0
298int main()
299{
300  int n;
301
302  srand(time(NULL));
303
304  for (n=0 ; n<1000 ; ++n)
305  {
306    unsigned int x = Random(10000)+2;
307    int a = bdd_prime_lte(x);
308    int b=_bdd_prime_lte(x);
309    /*printf("%d: %d, %d  ", x, );*/
310    if (a != b)
311      printf("ERROR");
312    /*printf("\n");*/
313  }
314
315  return 0;
316}
317#endif
318
319
320/* EOF */
321
322