| 1 |
/********************************************************************** |
| 2 |
rand.cpp - Pseudo random number generator. |
| 3 |
|
| 4 |
Copyright (C) 1998-2001 by OpenEye Scientific Software, Inc. |
| 5 |
Some portions Copyright (C) 2001-2005 by Geoffrey R. Hutchison |
| 6 |
|
| 7 |
This file is part of the Open Babel project. |
| 8 |
For more information, see <http://openbabel.sourceforge.net/> |
| 9 |
|
| 10 |
This program is free software; you can redistribute it and/or modify |
| 11 |
it under the terms of the GNU General Public License as published by |
| 12 |
the Free Software Foundation version 2 of the License. |
| 13 |
|
| 14 |
This program is distributed in the hope that it will be useful, |
| 15 |
but WITHOUT ANY WARRANTY; without even the implied warranty of |
| 16 |
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
| 17 |
GNU General Public License for more details. |
| 18 |
***********************************************************************/ |
| 19 |
|
| 20 |
/* derived from sample.c |
| 21 |
* Pseudo-Random Number Generator Generator |
| 22 |
* Roger Sayle, Metaphorics LLC |
| 23 |
* Version 1.2, September 1998 |
| 24 |
*/ |
| 25 |
|
| 26 |
#include <stdlib.h> |
| 27 |
#include <stdio.h> |
| 28 |
#include <math.h> |
| 29 |
#include "mol.hpp" |
| 30 |
#include "vector3.hpp" |
| 31 |
#include "obutil.hpp" |
| 32 |
|
| 33 |
#ifndef True |
| 34 |
#define True 1 |
| 35 |
#define False 0 |
| 36 |
#endif |
| 37 |
|
| 38 |
#define IsEven(x) (((x)&1)==0) |
| 39 |
#define IsOdd(x) (((x)&1)!=0) |
| 40 |
|
| 41 |
#define BothEven(x,y) IsEven((x)|(y)) |
| 42 |
#define IsPrime(x) (!IsEven((x))&&IsOddPrime((x))) |
| 43 |
|
| 44 |
#define HiPart(x) (x>>16) |
| 45 |
#define LoPart(x) ((x)&0xffff) |
| 46 |
|
| 47 |
/* The maximum number of unique prime factors of a 32 bit number */ |
| 48 |
/* 2*3*5*7*11*13*17*19*23*29 = 6469693230 > 2^32 = 4294967296 */ |
| 49 |
#define MAXFACT 10 |
| 50 |
|
| 51 |
namespace OpenBabel |
| 52 |
{ |
| 53 |
#define MAXPRIMES 256 |
| 54 |
static int primes[MAXPRIMES] = { |
| 55 |
1, 2, 3, 5, 7, 11, 13, 17, 19, 23, |
| 56 |
29, 31, 37, 41, 43, 47, 53, 59, 61, 67, |
| 57 |
71, 73, 79, 83, 89, 97, 101, 103, 107, 109, |
| 58 |
113, 127, 131, 137, 139, 149, 151, 157, 163, 167, |
| 59 |
173, 179, 181, 191, 193, 197, 199, 211, 223, 227, |
| 60 |
229, 233, 239, 241, 251, 257, 263, 269, 271, 277, |
| 61 |
281, 283, 293, 307, 311, 313, 317, 331, 337, 347, |
| 62 |
349, 353, 359, 367, 373, 379, 383, 389, 397, 401, |
| 63 |
409, 419, 421, 431, 433, 439, 443, 449, 457, 461, |
| 64 |
463, 467, 479, 487, 491, 499, 503, 509, 521, 523, |
| 65 |
541, 547, 557, 563, 569, 571, 577, 587, 593, 599, |
| 66 |
601, 607, 613, 617, 619, 631, 641, 643, 647, 653, |
| 67 |
659, 661, 673, 677, 683, 691, 701, 709, 719, 727, |
| 68 |
733, 739, 743, 751, 757, 761, 769, 773, 787, 797, |
| 69 |
809, 811, 821, 823, 827, 829, 839, 853, 857, 859, |
| 70 |
863, 877, 881, 883, 887, 907, 911, 919, 929, 937, |
| 71 |
941, 947, 953, 967, 971, 977, 983, 991, 997, 1009, |
| 72 |
1013, 1019, 1021, 1031, 1033, 1039, 1049, 1051, 1061, 1063, |
| 73 |
1069, 1087, 1091, 1093, 1097, 1103, 1109, 1117, 1123, 1129, |
| 74 |
1151, 1153, 1163, 1171, 1181, 1187, 1193, 1201, 1213, 1217, |
| 75 |
1223, 1229, 1231, 1237, 1249, 1259, 1277, 1279, 1283, 1289, |
| 76 |
1291, 1297, 1301, 1303, 1307, 1319, 1321, 1327, 1361, 1367, |
| 77 |
1373, 1381, 1399, 1409, 1423, 1427, 1429, 1433, 1439, 1447, |
| 78 |
1451, 1453, 1459, 1471, 1481, 1483, 1487, 1489, 1493, 1499, |
| 79 |
1511, 1523, 1531, 1543, 1549, 1553, 1559, 1567, 1571, 1579, |
| 80 |
1583, 1597, 1601, 1607, 1609, 1613 |
| 81 |
}; |
| 82 |
|
| 83 |
static unsigned int isqrt( unsigned int val ) |
| 84 |
{ |
| 85 |
register unsigned int temp; |
| 86 |
register unsigned int rem; |
| 87 |
register int i,result; |
| 88 |
|
| 89 |
i = 16; |
| 90 |
while( !(val&((unsigned int)3<<30)) && i ) |
| 91 |
{ |
| 92 |
val <<= 2; |
| 93 |
i--; |
| 94 |
} |
| 95 |
|
| 96 |
if( i ) |
| 97 |
{ |
| 98 |
rem = (val>>30)-1; |
| 99 |
val <<= 2; |
| 100 |
result = 1; |
| 101 |
i--; |
| 102 |
|
| 103 |
while( i ) |
| 104 |
{ |
| 105 |
rem = (rem<<2) | (val>>30); |
| 106 |
result <<= 1; |
| 107 |
val <<= 2; |
| 108 |
|
| 109 |
temp = result<<1; |
| 110 |
if( rem > temp ) |
| 111 |
{ |
| 112 |
rem -= temp|1; |
| 113 |
result |= 1; |
| 114 |
} |
| 115 |
i--; |
| 116 |
} |
| 117 |
return result; |
| 118 |
} |
| 119 |
else |
| 120 |
return 0; |
| 121 |
} |
| 122 |
|
| 123 |
static int IsOddPrime( unsigned int x ) |
| 124 |
{ |
| 125 |
register unsigned int root; |
| 126 |
register unsigned int i; |
| 127 |
|
| 128 |
root = isqrt(x); |
| 129 |
for( i=2; i<MAXPRIMES-1; i++ ) |
| 130 |
{ |
| 131 |
if( (x%primes[i]) == 0 ) |
| 132 |
return False; |
| 133 |
if( (unsigned int) primes[i] >= root ) |
| 134 |
return True; |
| 135 |
} |
| 136 |
|
| 137 |
for( i=primes[MAXPRIMES-1]; i<=root; i+=2 ) |
| 138 |
if( (x%i) == 0 ) |
| 139 |
return False; |
| 140 |
return True; |
| 141 |
} |
| 142 |
|
| 143 |
static int RelativelyPrime( unsigned int x, unsigned int y ) |
| 144 |
{ |
| 145 |
if( BothEven(x,y) ) |
| 146 |
return False; |
| 147 |
|
| 148 |
if( IsEven(x) ) |
| 149 |
{ |
| 150 |
do |
| 151 |
{ |
| 152 |
x >>= 1; |
| 153 |
} |
| 154 |
while( IsEven(x) ); |
| 155 |
} |
| 156 |
else |
| 157 |
while( IsEven(y) ) |
| 158 |
y >>= 1; |
| 159 |
|
| 160 |
while( x != y ) |
| 161 |
{ |
| 162 |
if( x > y ) |
| 163 |
{ |
| 164 |
x -= y; |
| 165 |
do |
| 166 |
{ |
| 167 |
x >>= 1; |
| 168 |
} |
| 169 |
while( IsEven(x) ); |
| 170 |
} |
| 171 |
else if( x < y ) |
| 172 |
{ |
| 173 |
y -= x; |
| 174 |
do |
| 175 |
{ |
| 176 |
y >>= 1; |
| 177 |
} |
| 178 |
while( IsEven(y) ); |
| 179 |
} |
| 180 |
} |
| 181 |
return( x == 1 ); |
| 182 |
} |
| 183 |
|
| 184 |
void DoubleAdd( DoubleType *x, unsigned int y ) |
| 185 |
{ |
| 186 |
x->lo += y; |
| 187 |
if( x->lo < y ) |
| 188 |
x->hi++; |
| 189 |
} |
| 190 |
|
| 191 |
void DoubleMultiply( unsigned int x, unsigned int y, DoubleType *z ) |
| 192 |
{ |
| 193 |
register unsigned int x0, x1, x2, x3; |
| 194 |
register unsigned int hx, lx; |
| 195 |
register unsigned int hy, ly; |
| 196 |
|
| 197 |
hx = HiPart(x); |
| 198 |
lx = LoPart(x); |
| 199 |
hy = HiPart(y); |
| 200 |
ly = LoPart(y); |
| 201 |
|
| 202 |
x0 = lx*ly; |
| 203 |
x1 = lx*hy; |
| 204 |
x2 = hx*ly; |
| 205 |
x3 = hx*hy; |
| 206 |
|
| 207 |
x1 += HiPart(x0); |
| 208 |
x1 += x2; |
| 209 |
if( x1 < x2 ) |
| 210 |
x3 += (1<<16); |
| 211 |
|
| 212 |
z->hi = HiPart(x1) + x3; |
| 213 |
z->lo = (LoPart(x1)<<16) + LoPart(x0); |
| 214 |
} |
| 215 |
|
| 216 |
static int LeadingZeros( unsigned int x ) |
| 217 |
{ |
| 218 |
static int table[256] = { |
| 219 |
0,1,2,2,3,3,3,3,4,4,4,4,4,4,4,4,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5, |
| 220 |
6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6, |
| 221 |
7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7, |
| 222 |
7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7, |
| 223 |
8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8, |
| 224 |
8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8, |
| 225 |
8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8, |
| 226 |
8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8 |
| 227 |
}; |
| 228 |
|
| 229 |
if( x >= (1<<16) ) |
| 230 |
{ |
| 231 |
if( x >= (1<<24) ) |
| 232 |
{ |
| 233 |
return 8-table[x>>24]; |
| 234 |
} |
| 235 |
else |
| 236 |
return 16-table[x>>16]; |
| 237 |
} |
| 238 |
else if( x >= (1<<8) ) |
| 239 |
{ |
| 240 |
return 24-table[x>>8]; |
| 241 |
} |
| 242 |
else |
| 243 |
return 32-table[x]; |
| 244 |
} |
| 245 |
|
| 246 |
unsigned int DoubleModulus( DoubleType *n, unsigned int d ) |
| 247 |
{ |
| 248 |
register unsigned int d1, d0; |
| 249 |
register unsigned int r1, r0; |
| 250 |
register unsigned int m,s; |
| 251 |
|
| 252 |
s = LeadingZeros(d); |
| 253 |
if( s > 0 ) |
| 254 |
{ |
| 255 |
d = d<<s; |
| 256 |
n->hi = (n->hi<<s) | (n->lo>>(32-s)); |
| 257 |
n->lo = n->lo << s; |
| 258 |
} |
| 259 |
|
| 260 |
d1 = HiPart(d); |
| 261 |
d0 = LoPart(d); |
| 262 |
|
| 263 |
m = ((unsigned int)(n->hi/d1)) * d0; |
| 264 |
r1 = ((n->hi%d1)<<16) + HiPart(n->lo); |
| 265 |
if( r1 < m ) |
| 266 |
{ |
| 267 |
r1 += d; |
| 268 |
if( (r1>=d) && (r1<m) ) |
| 269 |
r1 += d; |
| 270 |
} |
| 271 |
r1 -= m; |
| 272 |
|
| 273 |
m = ((unsigned int)(r1/d1)) * d0; |
| 274 |
r0 = ((r1%d1)<<16) + LoPart(n->lo); |
| 275 |
if( r0 < m ) |
| 276 |
{ |
| 277 |
r0 += d; |
| 278 |
if( (r0>=d) && (r0<m) ) |
| 279 |
r0 += d; |
| 280 |
} |
| 281 |
r0 -= m; |
| 282 |
|
| 283 |
return r0 >> s; |
| 284 |
} |
| 285 |
|
| 286 |
static int DeterminePotency( unsigned int m, unsigned int a ) |
| 287 |
{ |
| 288 |
auto DoubleType d; |
| 289 |
register unsigned int k; |
| 290 |
register unsigned int b; |
| 291 |
register int s; |
| 292 |
|
| 293 |
b = a-1; |
| 294 |
k = b; |
| 295 |
s = 1; |
| 296 |
while( (k!=0) && (s<100) ) |
| 297 |
{ |
| 298 |
DoubleMultiply(k,b,&d); |
| 299 |
k = DoubleModulus(&d,m); |
| 300 |
s++; |
| 301 |
} |
| 302 |
return s; |
| 303 |
} |
| 304 |
|
| 305 |
static int DetermineFactors( unsigned int x, unsigned int *factors ) |
| 306 |
{ |
| 307 |
register unsigned int *ptr; |
| 308 |
register unsigned int half; |
| 309 |
register unsigned int i; |
| 310 |
|
| 311 |
half = x/2; |
| 312 |
ptr = factors; |
| 313 |
for( i=1; i<MAXPRIMES; i++ ) |
| 314 |
{ |
| 315 |
if( (x%primes[i]) == 0 ) |
| 316 |
*ptr++ = primes[i]; |
| 317 |
if( (unsigned int)(primes[i]) >= half ) |
| 318 |
return ptr-factors; |
| 319 |
} |
| 320 |
|
| 321 |
for( i=primes[MAXPRIMES-1]+2; i<=half; i+=2 ) |
| 322 |
if( IsOddPrime(i) && ((x%i)==0) ) |
| 323 |
*ptr++ = i; |
| 324 |
return ptr-factors; |
| 325 |
} |
| 326 |
|
| 327 |
static unsigned int DetermineIncrement( unsigned int m ) |
| 328 |
{ |
| 329 |
register unsigned int hi,lo; |
| 330 |
register unsigned int half; |
| 331 |
register unsigned int i; |
| 332 |
|
| 333 |
/* 1/2 + sqrt(3)/6 */ |
| 334 |
hi = (int)floor(0.7886751345948*m+0.5); |
| 335 |
if( RelativelyPrime(m,hi) ) |
| 336 |
return hi; |
| 337 |
|
| 338 |
/* 1/2 - sqrt(3)/6 */ |
| 339 |
lo = (int)floor(0.2113248654052*m+0.5); |
| 340 |
if( RelativelyPrime(m,lo) ) |
| 341 |
return lo; |
| 342 |
|
| 343 |
half = m/2; |
| 344 |
for( i=1; i<half; i++ ) |
| 345 |
{ |
| 346 |
if( RelativelyPrime(m,hi+i) ) |
| 347 |
return hi+i; |
| 348 |
if( RelativelyPrime(m,hi-i) ) |
| 349 |
return hi-i; |
| 350 |
if( RelativelyPrime(m,lo+i) ) |
| 351 |
return lo+i; |
| 352 |
if( RelativelyPrime(m,lo-i) ) |
| 353 |
return lo-i; |
| 354 |
} |
| 355 |
return 1; |
| 356 |
} |
| 357 |
|
| 358 |
int DetermineSequence( unsigned int m, unsigned int *pm, |
| 359 |
unsigned int *pa, |
| 360 |
unsigned int *pc ) |
| 361 |
{ |
| 362 |
auto unsigned int fact[MAXFACT]; |
| 363 |
register unsigned int a=0, c; |
| 364 |
register unsigned int b; |
| 365 |
register int pot,best; |
| 366 |
register int count; |
| 367 |
register int flag; |
| 368 |
register int i; |
| 369 |
|
| 370 |
do |
| 371 |
{ |
| 372 |
best = 0; |
| 373 |
|
| 374 |
count = DetermineFactors(m,fact); |
| 375 |
if( (m&3) == 0 ) |
| 376 |
fact[0] = 4; |
| 377 |
|
| 378 |
#ifdef DEBUG |
| 379 |
|
| 380 |
printf("factors(%d):",m); |
| 381 |
#endif |
| 382 |
|
| 383 |
if( count ) |
| 384 |
{ |
| 385 |
#ifdef DEBUG |
| 386 |
printf("%u",fact[0]); |
| 387 |
for( i=1; i<count; i++ ) |
| 388 |
printf(",%u",fact[i]); |
| 389 |
#endif |
| 390 |
|
| 391 |
for( b=m-2; b>0; b-- ) |
| 392 |
{ |
| 393 |
flag = True; |
| 394 |
for( i=0; i<count; i++ ) |
| 395 |
if( b%fact[i] ) |
| 396 |
{ |
| 397 |
flag = False; |
| 398 |
break; |
| 399 |
} |
| 400 |
|
| 401 |
if( flag ) |
| 402 |
{ |
| 403 |
pot = DeterminePotency(m,b+1); |
| 404 |
if( pot > best ) |
| 405 |
{ |
| 406 |
best = pot; |
| 407 |
a = b+1; |
| 408 |
} |
| 409 |
} |
| 410 |
} |
| 411 |
|
| 412 |
#ifdef DEBUG |
| 413 |
if( best ) |
| 414 |
printf(" [%d]",best); |
| 415 |
#endif |
| 416 |
|
| 417 |
} |
| 418 |
#ifdef DEBUG |
| 419 |
printf("\n"); |
| 420 |
#endif |
| 421 |
|
| 422 |
m++; |
| 423 |
} |
| 424 |
while( best < 3 ); |
| 425 |
m--; |
| 426 |
|
| 427 |
c = DetermineIncrement(m); |
| 428 |
|
| 429 |
*pm = m; |
| 430 |
*pa = a; |
| 431 |
*pc = c; |
| 432 |
return best; |
| 433 |
} |
| 434 |
|
| 435 |
|
| 436 |
void GenerateSequence( unsigned int p, unsigned int m, |
| 437 |
unsigned int a, unsigned int c ) |
| 438 |
{ |
| 439 |
register unsigned int i; |
| 440 |
register unsigned int x; |
| 441 |
auto DoubleType d; |
| 442 |
|
| 443 |
x = 0; /* seed */ |
| 444 |
for( i=0; i<p; i++ ) |
| 445 |
{ |
| 446 |
printf("%u\n",x); |
| 447 |
|
| 448 |
do |
| 449 |
{ |
| 450 |
DoubleMultiply(a,x,&d); |
| 451 |
DoubleAdd(&d,c); |
| 452 |
x = DoubleModulus(&d,m); |
| 453 |
} |
| 454 |
while( x >= p ); |
| 455 |
} |
| 456 |
} |
| 457 |
|
| 458 |
//********************************************** |
| 459 |
//***** Member functions from Random class ***** |
| 460 |
//********************************************** |
| 461 |
|
| 462 |
OBRandom::OBRandom(bool useSysRand) |
| 463 |
{ |
| 464 |
|
| 465 |
this->OBRandomUseSysRand= useSysRand; |
| 466 |
p = 70092; |
| 467 |
DetermineSequence(p,&m,&a,&c); |
| 468 |
x = 0; /* seed */ |
| 469 |
} |
| 470 |
|
| 471 |
int OBRandom::NextInt() |
| 472 |
{ |
| 473 |
if (OBRandomUseSysRand) |
| 474 |
{ |
| 475 |
return(rand()); |
| 476 |
} |
| 477 |
do |
| 478 |
{ |
| 479 |
DoubleMultiply(a,x,&d); |
| 480 |
DoubleAdd(&d,c); |
| 481 |
x = DoubleModulus(&d,m); |
| 482 |
} |
| 483 |
while( x >= p ); |
| 484 |
|
| 485 |
return(x); |
| 486 |
} |
| 487 |
|
| 488 |
double OBRandom::NextFloat() |
| 489 |
{ |
| 490 |
|
| 491 |
if (OBRandomUseSysRand) |
| 492 |
{ |
| 493 |
return(double(rand())/double(RAND_MAX)); |
| 494 |
} |
| 495 |
do |
| 496 |
{ |
| 497 |
DoubleMultiply(a,x,&d); |
| 498 |
DoubleAdd(&d,c); |
| 499 |
x = DoubleModulus(&d,m); |
| 500 |
} |
| 501 |
while( x >= p ); |
| 502 |
|
| 503 |
return((double)x/p); |
| 504 |
} |
| 505 |
|
| 506 |
void OBRandom::TimeSeed() |
| 507 |
{ |
| 508 |
#ifdef WIN32 |
| 509 |
// for VC++ do it this way |
| 510 |
time_t ltime; |
| 511 |
time(<ime); |
| 512 |
const long secs= long(ltime); |
| 513 |
x= secs % p; |
| 514 |
srand( (unsigned)time( NULL ) ); |
| 515 |
#else |
| 516 |
|
| 517 |
timeval time; |
| 518 |
gettimeofday(&time,(struct timezone *)NULL); |
| 519 |
x = (time.tv_usec%p); |
| 520 |
srand( x ); |
| 521 |
#ifdef HAVE_SRANDDEV |
| 522 |
sranddev(); |
| 523 |
#endif |
| 524 |
|
| 525 |
#endif |
| 526 |
} |
| 527 |
|
| 528 |
} //end namespace OpenBabel |
| 529 |
|
| 530 |
//! \file rand.cpp |
| 531 |
//! \brief Pseudo random number generator. |
