← Back to branch summary

~vcs-imports/openbabel/trunk

« back to all changes in this revision

Viewing changes to rand.cpp

  • Committer: ghutchis
  • Date: 2001-11-27 18:50:36 UTC
  • Revision ID: svn-v4:86f38270-e075-4da6-bf68-7dcd545c500b:openbabel/trunk:46
Initial revision

Show diffs side-by-side

added added

removed removed

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