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(÷nd, &remainder); 118 119 for (i=0 ; i<64 ; ++i) 120 { 121 if (remainder >= divisor) 122 remainder -= divisor; 123 124 u64_shl(÷nd, &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