src/openbabel/rand.cpp

/**********************************************************************
rand.cpp - Pseudo random number generator.
 
Copyright (C) 1998-2001 by OpenEye Scientific Software, Inc.
Some portions Copyright (C) 2001-2005 by Geoffrey R. Hutchison
 
This file is part of the Open Babel project.
For more information, see <http://openbabel.sourceforge.net/>
 
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation version 2 of the License.
 
This program 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 General Public License for more details.
***********************************************************************/

/* derived from sample.c
 * Pseudo-Random Number Generator Generator
 * Roger Sayle, Metaphorics LLC
 * Version 1.2, September 1998
 */

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "mol.hpp"
#include "vector3.hpp"
#include "obutil.hpp"

#ifndef True
#define True   1
#define False  0
#endif

#define IsEven(x) (((x)&1)==0)
#define IsOdd(x)  (((x)&1)!=0)

#define BothEven(x,y) IsEven((x)|(y))
#define IsPrime(x)    (!IsEven((x))&&IsOddPrime((x)))

#define HiPart(x)   (x>>16)
#define LoPart(x)   ((x)&0xffff)

/* The maximum number of unique prime factors of a 32 bit number */
/* 2*3*5*7*11*13*17*19*23*29 = 6469693230 > 2^32 = 4294967296    */
#define MAXFACT    10

namespace OpenBabel
{
#define MAXPRIMES  256
static int primes[MAXPRIMES] = {
                                   1,    2,    3,    5,    7,   11,   13,   17,   19,   23,
                                   29,   31,   37,   41,   43,   47,   53,   59,   61,   67,
                                   71,   73,   79,   83,   89,   97,  101,  103,  107,  109,
                                   113,  127,  131,  137,  139,  149,  151,  157,  163,  167,
                                   173,  179,  181,  191,  193,  197,  199,  211,  223,  227,
                                   229,  233,  239,  241,  251,  257,  263,  269,  271,  277,
                                   281,  283,  293,  307,  311,  313,  317,  331,  337,  347,
                                   349,  353,  359,  367,  373,  379,  383,  389,  397,  401,
                                   409,  419,  421,  431,  433,  439,  443,  449,  457,  461,
                                   463,  467,  479,  487,  491,  499,  503,  509,  521,  523,
                                   541,  547,  557,  563,  569,  571,  577,  587,  593,  599,
                                   601,  607,  613,  617,  619,  631,  641,  643,  647,  653,
                                   659,  661,  673,  677,  683,  691,  701,  709,  719,  727,
                                   733,  739,  743,  751,  757,  761,  769,  773,  787,  797,
                                   809,  811,  821,  823,  827,  829,  839,  853,  857,  859,
                                   863,  877,  881,  883,  887,  907,  911,  919,  929,  937,
                                   941,  947,  953,  967,  971,  977,  983,  991,  997, 1009,
                                   1013, 1019, 1021, 1031, 1033, 1039, 1049, 1051, 1061, 1063,
                                   1069, 1087, 1091, 1093, 1097, 1103, 1109, 1117, 1123, 1129,
                                   1151, 1153, 1163, 1171, 1181, 1187, 1193, 1201, 1213, 1217,
                                   1223, 1229, 1231, 1237, 1249, 1259, 1277, 1279, 1283, 1289,
                                   1291, 1297, 1301, 1303, 1307, 1319, 1321, 1327, 1361, 1367,
                                   1373, 1381, 1399, 1409, 1423, 1427, 1429, 1433, 1439, 1447,
                                   1451, 1453, 1459, 1471, 1481, 1483, 1487, 1489, 1493, 1499,
                                   1511, 1523, 1531, 1543, 1549, 1553, 1559, 1567, 1571, 1579,
                                   1583, 1597, 1601, 1607, 1609, 1613
                               };

static unsigned int isqrt( unsigned int val )
{
    register unsigned int temp;
    register unsigned int rem;
    register int i,result;

    i = 16;
    while( !(val&((unsigned int)3<<30)) && i )
    {
        val <<= 2;
        i--;
    }

    if( i )
    {
        rem = (val>>30)-1;
        val <<= 2;
        result = 1;
        i--;

        while( i )
        {
            rem = (rem<<2) | (val>>30);
            result <<= 1;
            val <<= 2;

            temp = result<<1;
            if( rem > temp )
            {
                rem -= temp|1;
                result |= 1;
            }
            i--;
        }
        return result;
    }
    else
        return 0;
}

static int IsOddPrime( unsigned int x )
{
    register unsigned int root;
    register unsigned int i;

    root = isqrt(x);
    for( i=2; i<MAXPRIMES-1; i++ )
    {
        if( (x%primes[i]) == 0 )
            return False;
        if( (unsigned int) primes[i] >= root )
            return True;
    }

    for( i=primes[MAXPRIMES-1]; i<=root; i+=2 )
        if( (x%i) == 0 )
            return False;
    return True;
}

static int RelativelyPrime( unsigned int x, unsigned int y )
{
    if( BothEven(x,y) )
        return False;

    if( IsEven(x) )
    {
        do
        {
            x >>= 1;
        }
        while( IsEven(x) );
    }
    else
        while( IsEven(y) )
            y >>= 1;

    while( x != y )
    {
        if( x > y )
        {
            x -= y;
            do
            {
                x >>= 1;
            }
            while( IsEven(x) );
        }
        else if( x < y )
        {
            y -= x;
            do
            {
                y >>= 1;
            }
            while( IsEven(y) );
        }
    }
    return( x == 1 );
}

void DoubleAdd( DoubleType *x, unsigned int y )
{
    x->lo += y;
    if( x->lo < y )
        x->hi++;
}

void DoubleMultiply( unsigned int x, unsigned int y, DoubleType *z )
{
    register unsigned int x0, x1, x2, x3;
    register unsigned int hx, lx;
    register unsigned int hy, ly;

    hx = HiPart(x);
    lx = LoPart(x);
    hy = HiPart(y);
    ly = LoPart(y);

    x0 = lx*ly;
    x1 = lx*hy;
    x2 = hx*ly;
    x3 = hx*hy;

    x1 += HiPart(x0);
    x1 += x2;
    if( x1 < x2 )
        x3 += (1<<16);

    z->hi = HiPart(x1) + x3;
    z->lo = (LoPart(x1)<<16) + LoPart(x0);
}

static int LeadingZeros( unsigned int x )
{
    static int table[256] = {
                                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,
                                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,
                                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,
                                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,
                                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,
                                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,
                                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,
                                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
                            };

    if( x >= (1<<16) )
    {
        if( x >= (1<<24) )
        {
            return  8-table[x>>24];
        }
        else
            return 16-table[x>>16];
    }
    else if( x >= (1<<8) )
    {
        return 24-table[x>>8];
    }
    else
        return 32-table[x];
}

unsigned int DoubleModulus( DoubleType *n,  unsigned int d )
{
    register unsigned int d1, d0;
    register unsigned int r1, r0;
    register unsigned int m,s;

    s = LeadingZeros(d);
    if( s > 0 )
    {
        d = d<<s;
        n->hi = (n->hi<<s) | (n->lo>>(32-s));
        n->lo = n->lo << s;
    }

    d1 = HiPart(d);
    d0 = LoPart(d);

    m = ((unsigned int)(n->hi/d1)) * d0;
    r1 = ((n->hi%d1)<<16) + HiPart(n->lo);
    if( r1 < m )
    {
        r1 += d;
        if( (r1>=d) && (r1<m) )
            r1 += d;
    }
    r1 -= m;

    m = ((unsigned int)(r1/d1)) * d0;
    r0 = ((r1%d1)<<16) + LoPart(n->lo);
    if( r0 < m )
    {
        r0 += d;
        if( (r0>=d) && (r0<m) )
            r0 += d;
    }
    r0 -= m;

    return r0 >> s;
}

static int DeterminePotency( unsigned int m, unsigned int a )
{
    auto DoubleType d;
    register unsigned int k;
    register unsigned int b;
    register int s;

    b = a-1;
    k = b;
    s = 1;
    while( (k!=0) && (s<100) )
    {
        DoubleMultiply(k,b,&d);
        k = DoubleModulus(&d,m);
        s++;
    }
    return s;
}

static int DetermineFactors( unsigned int x, unsigned int *factors )
{
    register unsigned int *ptr;
    register unsigned int half;
    register unsigned int i;

    half = x/2;
    ptr = factors;
    for( i=1; i<MAXPRIMES; i++ )
    {
        if( (x%primes[i]) == 0 )
            *ptr++ = primes[i];
        if( (unsigned int)(primes[i]) >= half )
            return ptr-factors;
    }

    for( i=primes[MAXPRIMES-1]+2; i<=half; i+=2 )
        if( IsOddPrime(i) && ((x%i)==0) )
            *ptr++ = i;
    return ptr-factors;
}

static unsigned int DetermineIncrement( unsigned int m )
{
    register unsigned int hi,lo;
    register unsigned int half;
    register unsigned int i;

    /* 1/2 + sqrt(3)/6 */
    hi = (int)floor(0.7886751345948*m+0.5);
    if( RelativelyPrime(m,hi) )
        return hi;

    /* 1/2 - sqrt(3)/6 */
    lo = (int)floor(0.2113248654052*m+0.5);
    if( RelativelyPrime(m,lo) )
        return lo;

    half = m/2;
    for( i=1; i<half; i++ )
    {
        if( RelativelyPrime(m,hi+i) )
            return hi+i;
        if( RelativelyPrime(m,hi-i) )
            return hi-i;
        if( RelativelyPrime(m,lo+i) )
            return lo+i;
        if( RelativelyPrime(m,lo-i) )
            return lo-i;
    }
    return 1;
}

int DetermineSequence( unsigned int m, unsigned int *pm,
                       unsigned int *pa,
                       unsigned int *pc )
{
    auto unsigned int fact[MAXFACT];
    register unsigned int a=0, c;
    register unsigned int b;
    register int pot,best;
    register int count;
    register int flag;
    register int i;

    do
    {
        best = 0;

        count = DetermineFactors(m,fact);
        if( (m&3) == 0 )
            fact[0] = 4;

#ifdef DEBUG

        printf("factors(%d):",m);
#endif

        if( count )
        {
#ifdef DEBUG
            printf("%u",fact[0]);
            for( i=1; i<count; i++ )
                printf(",%u",fact[i]);
#endif

            for( b=m-2; b>0; b-- )
            {
                flag = True;
                for( i=0; i<count; i++ )
                    if( b%fact[i] )
                    {
                        flag = False;
                        break;
                    }

                if( flag )
                {
                    pot = DeterminePotency(m,b+1);
                    if( pot > best )
                    {
                        best = pot;
                        a = b+1;
                    }
                }
            }

#ifdef DEBUG
            if( best )
                printf(" [%d]",best);
#endif

        }
#ifdef DEBUG
        printf("\n");
#endif

        m++;
    }
    while( best < 3 );
    m--;

    c = DetermineIncrement(m);

    *pm = m;
    *pa = a;
    *pc = c;
    return best;
}


void GenerateSequence( unsigned int p, unsigned int m,
                       unsigned int a, unsigned int c )
{
    register unsigned int i;
    register unsigned int x;
    auto DoubleType d;

    x = 0;  /* seed */
    for( i=0; i<p; i++ )
    {
        printf("%u\n",x);

        do
        {
            DoubleMultiply(a,x,&d);
            DoubleAdd(&d,c);
            x = DoubleModulus(&d,m);
        }
        while( x >= p );
    }
}

//**********************************************
//***** Member functions from Random class *****
//**********************************************

OBRandom::OBRandom(bool useSysRand)
{

    this->OBRandomUseSysRand= useSysRand;
    p = 70092;
    DetermineSequence(p,&m,&a,&c);
    x = 0;  /* seed */
}

int OBRandom::NextInt()
{
    if (OBRandomUseSysRand)
    {
        return(rand());
    }
    do
    {
        DoubleMultiply(a,x,&d);
        DoubleAdd(&d,c);
        x = DoubleModulus(&d,m);
    }
    while( x >= p );

    return(x);
}

double OBRandom::NextFloat()
{

    if (OBRandomUseSysRand)
    {
        return(double(rand())/double(RAND_MAX));
    }
    do
    {
        DoubleMultiply(a,x,&d);
        DoubleAdd(&d,c);
        x = DoubleModulus(&d,m);
    }
    while( x >= p );

    return((double)x/p);
}

void OBRandom::TimeSeed()
{
#ifdef WIN32
    // for VC++ do it this way
    time_t ltime;
    time(&ltime);
    const long secs= long(ltime);
    x= secs % p;
    srand( (unsigned)time( NULL ) );
#else

    timeval time;
    gettimeofday(&time,(struct timezone *)NULL);
    x = (time.tv_usec%p);
    srand( x );
#ifdef HAVE_SRANDDEV
    sranddev();
#endif

#endif
}

} //end namespace OpenBabel

//! \file rand.cpp
//! \brief Pseudo random number generator.
Revision:	741
Committed:	Wed Nov 16 19:42:11 2005 UTC (19 years, 5 months ago) by tim
File size:	13194 byte(s)
Log Message:	adding openbabel
#	User	Rev	Content
1	tim	741	/**********************************************************************
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			/* 2357111317192329 = 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(&ltime);
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.