/**********************************************************************

rand.cpp - Pseudo random number generator.

Copyright (C) 1998-2001 by OpenEye Scientific Software, Inc.

Some portions Copyright (C) 2001-2006 by Geoffrey R. Hutchison

This file is part of the Open Babel project.

For more information, see <http://openbabel.org/>

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 <openbabel/babelconfig.h>

#include <cstdlib>

#include <cstdio>

#include <cmath>

#include <limits>

#include "rand.h"

#if TIME_WITH_SYS_TIME

#include <sys/time.h>

#include <time.h>

#else

#if HAVE_SYS_TIME_H

#include <sys/time.h>

#else

#include <time.h>

#endif

#endif

#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 )

  {

    unsigned int temp;

    unsigned int rem;

    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 )

  {

    unsigned int root;

    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 );

  }

  static void DoubleAdd( DoubleType *x, unsigned int y )

  {

    x->lo += y;

    if( x->lo < y )

      x->hi++;

  }

  static void DoubleMultiply( unsigned int x, unsigned int y, DoubleType *z )

  {

    unsigned int x0, x1, x2, x3;

    unsigned int hx, lx;

    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];

  }

  static unsigned int DoubleModulus( DoubleType *n,  unsigned int d )

  {

    unsigned int d1, d0;

    unsigned int r1, r0;

    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 )

  {

    DoubleType d;

    unsigned int k;

    unsigned int b;

    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 )

  {

    unsigned int *ptr;

    unsigned int half;

    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 )

  {

    unsigned int hi,lo;

    unsigned int half;

    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;

  }

  static int DetermineSequence( unsigned int m, unsigned int *pm,

                         unsigned int *pa,

                         unsigned int *pc )

  {

    unsigned int fact[MAXFACT];

    unsigned int a=0, c;

    unsigned int b;

    int pot,best;

    int count;

    int flag;

    int i;

    do

      {

        best = 0;

        count = DetermineFactors(m,fact);

        if( (m&3) == 0 )

          fact[0] = 4;

        if( count )

          {

            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;

                      }

                  }

              }

          }

        m++;

      }

    while( best < 3 );

    m--;

    c = DetermineIncrement(m);

    *pm = m;

    *pa = a;

    *pc = c;

    return best;

  }

  static void GenerateSequence( unsigned int p, unsigned int m,

                         unsigned int a, unsigned int c )

  {

    unsigned int i;

    unsigned int x;

    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 */

    if (useSysRand)

      this->TimeSeed();

  }

  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 *)nullptr);

    x = (time.tv_usec%p);

    srand( x );

#ifdef HAVE_SRANDDEV

    sranddev();

#endif

#endif

  }

  double OBRandom::Normal(double mu, double sigma) {

    const int n = 12;

    double z = -0.5 * n;

    for (int _ = 0; _ < n; _++) {

      z += NextFloat();

    }

    return mu + z * sigma;

  }

  bool OBRandom::Bernoulli(double p) {

    return UniformReal(0.0, 1.0) <= p;

  }

  uint_fast64_t OBRandomMT::randomSeed() const {

    auto ob_random_seed = std::getenv("OB_RANDOM_SEED");

    if (ob_random_seed && ob_random_seed[0] != '\0') {

      return std::atol(ob_random_seed);

    }

#if defined(WIN32) || defined(__MINGW32__)

    // for VC++ do it this way

    time_t ltime;

    time(&ltime);

    return ltime;

#else

    return std::random_device{}();

#endif

  }

  int OBRandomMT::NextInt() const {

    return UniformInt(0, std::numeric_limits<int>::max());

  }

  double OBRandomMT::NextFloat() const {

    return UniformReal(0.0, 1.0);

  }

  double OBRandomMT::Normal(double mu, double sigma) const {

    std::normal_distribution<double> n{mu, sigma};

    return n(*prng);

  }

  bool OBRandomMT::Bernoulli(double p) const {

    std::bernoulli_distribution b{p};

    return b(*prng);

  }

} //end namespace OpenBabel

//! \file rand.cpp

//! \brief Pseudo random number generator.