55071fb754
Also the cuba-compiler-check for gcc-versions containing a certain bug (4.2, 4.4.3) has been adopted. This still needs to be tested on systems having such a gcc. If this new version of the built-in library should be installed, first make sure that the old version is completely deinstalled (including headers, pkg-config-files, etc.).
366 lines
10 KiB
C
366 lines
10 KiB
C
/*
|
|
Random.c
|
|
quasi- and pseudo-random-number generation
|
|
last modified 8 Jun 10 th
|
|
*/
|
|
|
|
/***************************************************************************
|
|
* Copyright (C) 2004-2010 by Thomas Hahn *
|
|
* hahn@feynarts.de *
|
|
* *
|
|
* This library is free software; you can redistribute it and/or *
|
|
* modify it under the terms of the GNU Lesser General Public *
|
|
* License as published by the Free Software Foundation; either *
|
|
* version 2.1 of the License, or (at your option) any later version. *
|
|
* *
|
|
* This library 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 *
|
|
* Lesser General Public License for more details. *
|
|
* *
|
|
* You should have received a copy of the GNU Lesser General Public *
|
|
* License along with this library; if not, write to the *
|
|
* Free Software Foundation, Inc., *
|
|
* 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA *
|
|
***************************************************************************/
|
|
|
|
|
|
/*
|
|
PART 1: Sobol quasi-random-number generator
|
|
adapted from ACM TOMS algorithm 659
|
|
*/
|
|
|
|
static void SobolGet(This *t, real *x)
|
|
{
|
|
number seq = t->rng.sobol.seq++;
|
|
count zerobit = 0, dim;
|
|
|
|
while( seq & 1 ) {
|
|
++zerobit;
|
|
seq >>= 1;
|
|
}
|
|
|
|
for( dim = 0; dim < t->ndim; ++dim ) {
|
|
t->rng.sobol.prev[dim] ^= t->rng.sobol.v[dim][zerobit];
|
|
x[dim] = t->rng.sobol.prev[dim]*t->rng.sobol.norm;
|
|
}
|
|
}
|
|
|
|
|
|
static void SobolSkip(This *t, number n)
|
|
{
|
|
while( n-- ) {
|
|
number seq = t->rng.sobol.seq++;
|
|
count zerobit = 0, dim;
|
|
|
|
while( seq & 1 ) {
|
|
++zerobit;
|
|
seq >>= 1;
|
|
}
|
|
|
|
for( dim = 0; dim < t->ndim; ++dim )
|
|
t->rng.sobol.prev[dim] ^= t->rng.sobol.v[dim][zerobit];
|
|
}
|
|
}
|
|
|
|
|
|
static inline void SobolIni(This *t)
|
|
{
|
|
static number ini[9*40] = {
|
|
3, 1, 0, 0, 0, 0, 0, 0, 0,
|
|
7, 1, 1, 0, 0, 0, 0, 0, 0,
|
|
11, 1, 3, 7, 0, 0, 0, 0, 0,
|
|
13, 1, 1, 5, 0, 0, 0, 0, 0,
|
|
19, 1, 3, 1, 1, 0, 0, 0, 0,
|
|
25, 1, 1, 3, 7, 0, 0, 0, 0,
|
|
37, 1, 3, 3, 9, 9, 0, 0, 0,
|
|
59, 1, 3, 7, 13, 3, 0, 0, 0,
|
|
47, 1, 1, 5, 11, 27, 0, 0, 0,
|
|
61, 1, 3, 5, 1, 15, 0, 0, 0,
|
|
55, 1, 1, 7, 3, 29, 0, 0, 0,
|
|
41, 1, 3, 7, 7, 21, 0, 0, 0,
|
|
67, 1, 1, 1, 9, 23, 37, 0, 0,
|
|
97, 1, 3, 3, 5, 19, 33, 0, 0,
|
|
91, 1, 1, 3, 13, 11, 7, 0, 0,
|
|
109, 1, 1, 7, 13, 25, 5, 0, 0,
|
|
103, 1, 3, 5, 11, 7, 11, 0, 0,
|
|
115, 1, 1, 1, 3, 13, 39, 0, 0,
|
|
131, 1, 3, 1, 15, 17, 63, 13, 0,
|
|
193, 1, 1, 5, 5, 1, 27, 33, 0,
|
|
137, 1, 3, 3, 3, 25, 17, 115, 0,
|
|
145, 1, 1, 3, 15, 29, 15, 41, 0,
|
|
143, 1, 3, 1, 7, 3, 23, 79, 0,
|
|
241, 1, 3, 7, 9, 31, 29, 17, 0,
|
|
157, 1, 1, 5, 13, 11, 3, 29, 0,
|
|
185, 1, 3, 1, 9, 5, 21, 119, 0,
|
|
167, 1, 1, 3, 1, 23, 13, 75, 0,
|
|
229, 1, 3, 3, 11, 27, 31, 73, 0,
|
|
171, 1, 1, 7, 7, 19, 25, 105, 0,
|
|
213, 1, 3, 5, 5, 21, 9, 7, 0,
|
|
191, 1, 1, 1, 15, 5, 49, 59, 0,
|
|
253, 1, 1, 1, 1, 1, 33, 65, 0,
|
|
203, 1, 3, 5, 15, 17, 19, 21, 0,
|
|
211, 1, 1, 7, 11, 13, 29, 3, 0,
|
|
239, 1, 3, 7, 5, 7, 11, 113, 0,
|
|
247, 1, 1, 5, 3, 15, 19, 61, 0,
|
|
285, 1, 3, 1, 1, 9, 27, 89, 7,
|
|
369, 1, 1, 3, 7, 31, 15, 45, 23,
|
|
299, 1, 3, 3, 9, 9, 25, 107, 39 };
|
|
|
|
count dim, bit, nbits;
|
|
number max, *pini = ini;
|
|
cnumber nmax = 2*t->maxeval;
|
|
|
|
for( nbits = 0, max = 1; max <= nmax; max <<= 1 ) ++nbits;
|
|
t->rng.sobol.norm = 1./max;
|
|
|
|
for( bit = 0; bit < nbits; ++bit )
|
|
t->rng.sobol.v[0][bit] = (max >>= 1);
|
|
|
|
for( dim = 1; dim < t->ndim; ++dim ) {
|
|
number *pv = t->rng.sobol.v[dim], *pvv = pv;
|
|
number powers = *pini++, j;
|
|
int inibits = -1, bit;
|
|
for( j = powers; j; j >>= 1 ) ++inibits;
|
|
|
|
memcpy(pv, pini, inibits*sizeof(*pini));
|
|
pini += 8;
|
|
|
|
for( bit = inibits; bit < nbits; ++bit ) {
|
|
number newv = *pvv, j = powers;
|
|
int b;
|
|
for( b = 0; b < inibits; ++b ) {
|
|
if( j & 1 ) newv ^= pvv[b] << (inibits - b);
|
|
j >>= 1;
|
|
}
|
|
pvv[inibits] = newv;
|
|
++pvv;
|
|
}
|
|
|
|
for( bit = 0; bit < nbits - 1; ++bit )
|
|
pv[bit] <<= nbits - bit - 1;
|
|
}
|
|
|
|
t->rng.sobol.seq = 0;
|
|
VecClear(t->rng.sobol.prev);
|
|
|
|
t->rng.getrandom = SobolGet;
|
|
t->rng.skiprandom = SobolSkip;
|
|
}
|
|
|
|
|
|
/*
|
|
PART 2: Mersenne Twister pseudo-random-number generator
|
|
adapted from T. Nishimura's and M. Matsumoto's C code at
|
|
http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
|
|
*/
|
|
|
|
/* 32 or 53 random bits */
|
|
#define RANDOM_BITS 32
|
|
|
|
|
|
static inline state_t Twist(state_t a, state_t b)
|
|
{
|
|
state_t mixbits = (a & 0x80000000) | (b & 0x7fffffff);
|
|
state_t matrixA = (-(b & 1)) & 0x9908b0df;
|
|
return (mixbits >> 1) ^ matrixA;
|
|
}
|
|
|
|
|
|
static inline void MersenneReload(state_t *state)
|
|
{
|
|
state_t *s = state;
|
|
int j;
|
|
|
|
for( j = MERSENNE_N - MERSENNE_M + 1; --j; ++s )
|
|
*s = s[MERSENNE_M] ^ Twist(s[0], s[1]);
|
|
for( j = MERSENNE_M; --j; ++s )
|
|
*s = s[MERSENNE_M - MERSENNE_N] ^ Twist(s[0], s[1]);
|
|
*s = s[MERSENNE_M - MERSENNE_N] ^ Twist(s[0], state[0]);
|
|
}
|
|
|
|
|
|
static inline state_t MersenneInt(state_t s)
|
|
{
|
|
s ^= s >> 11;
|
|
s ^= (s << 7) & 0x9d2c5680;
|
|
s ^= (s << 15) & 0xefc60000;
|
|
return s ^ (s >> 18);
|
|
}
|
|
|
|
|
|
static void MersenneGet(This *t, real *x)
|
|
{
|
|
count next = t->rng.mersenne.next, dim;
|
|
|
|
for( dim = 0; dim < t->ndim; ++dim ) {
|
|
#if RANDOM_BITS == 53
|
|
state_t a, b;
|
|
#endif
|
|
|
|
if( next >= MERSENNE_N ) {
|
|
MersenneReload(t->rng.mersenne.state);
|
|
next = 0;
|
|
}
|
|
|
|
#if RANDOM_BITS == 53
|
|
a = MersenneInt(t->rng.mersenne.state[next++]) >> 5;
|
|
b = MersenneInt(t->rng.mersenne.state[next++]) >> 6;
|
|
x[dim] = (67108864.*a + b)/9007199254740992.;
|
|
#else
|
|
x[dim] = MersenneInt(t->rng.mersenne.state[next++])/4294967296.;
|
|
#endif
|
|
}
|
|
|
|
t->rng.mersenne.next = next;
|
|
}
|
|
|
|
|
|
static void MersenneSkip(This *t, number n)
|
|
{
|
|
#if RANDOM_BITS == 53
|
|
n = 2*n*t->ndim + t->rng.mersenne.next;
|
|
#else
|
|
n = n*t->ndim + t->rng.mersenne.next;
|
|
#endif
|
|
t->rng.mersenne.next = n % MERSENNE_N;
|
|
n /= MERSENNE_N;
|
|
while( n-- ) MersenneReload(t->rng.mersenne.state);
|
|
}
|
|
|
|
|
|
static inline void MersenneIni(This *t)
|
|
{
|
|
state_t seed = t->seed;
|
|
state_t *next = t->rng.mersenne.state;
|
|
count j;
|
|
|
|
for( j = 1; j <= MERSENNE_N; ++j ) {
|
|
*next++ = seed;
|
|
seed = 0x6c078965*(seed ^ (seed >> 30)) + j;
|
|
/* see Knuth TAOCP Vol 2, 3rd Ed, p. 106 for multiplier */
|
|
}
|
|
|
|
MersenneReload(t->rng.mersenne.state);
|
|
t->rng.mersenne.next = 0;
|
|
|
|
t->rng.getrandom = MersenneGet;
|
|
t->rng.skiprandom = MersenneSkip;
|
|
}
|
|
|
|
|
|
/*
|
|
PART 3: Ranlux subtract-and-borrow random-number generator
|
|
proposed by Marsaglia and Zaman, implemented by F. James with
|
|
the name RCARRY in 1991, and later improved by Martin Luescher
|
|
in 1993 to produce "Luxury Pseudorandom Numbers".
|
|
Adapted from the CERNlib Fortran 77 code by F. James, 1993.
|
|
|
|
The available luxury levels are:
|
|
|
|
level 0 (p = 24): equivalent to the original RCARRY of Marsaglia
|
|
and Zaman, very long period, but fails many tests.
|
|
level 1 (p = 48): considerable improvement in quality over level 0,
|
|
now passes the gap test, but still fails spectral test.
|
|
level 2 (p = 97): passes all known tests, but theoretically still
|
|
defective.
|
|
level 3 (p = 223): DEFAULT VALUE. Any theoretically possible
|
|
correlations have very small chance of being observed.
|
|
level 4 (p = 389): highest possible luxury, all 24 bits chaotic.
|
|
*/
|
|
|
|
|
|
static inline int RanluxInt(This *t, count n)
|
|
{
|
|
int s = 0;
|
|
|
|
while( n-- ) {
|
|
s = t->rng.ranlux.state[t->rng.ranlux.j24] -
|
|
t->rng.ranlux.state[t->rng.ranlux.i24] + t->rng.ranlux.carry;
|
|
s += (t->rng.ranlux.carry = NegQ(s)) & (1 << 24);
|
|
t->rng.ranlux.state[t->rng.ranlux.i24] = s;
|
|
--t->rng.ranlux.i24;
|
|
t->rng.ranlux.i24 += NegQ(t->rng.ranlux.i24) & 24;
|
|
--t->rng.ranlux.j24;
|
|
t->rng.ranlux.j24 += NegQ(t->rng.ranlux.j24) & 24;
|
|
}
|
|
|
|
return s;
|
|
}
|
|
|
|
|
|
static void RanluxGet(This *t, real *x)
|
|
{
|
|
/* The Generator proper: "Subtract-with-borrow",
|
|
as proposed by Marsaglia and Zaman, FSU, March 1989 */
|
|
|
|
count dim;
|
|
|
|
for( dim = 0; dim < t->ndim; ++dim ) {
|
|
cint nskip = (--t->rng.ranlux.n24 >= 0) ? 0 :
|
|
(t->rng.ranlux.n24 = 24, t->rng.ranlux.nskip);
|
|
cint s = RanluxInt(t, 1 + nskip);
|
|
x[dim] = s*0x1p-24;
|
|
/* small numbers (with less than 12 significant bits) are "padded" */
|
|
if( s < (1 << 12) )
|
|
x[dim] += t->rng.ranlux.state[t->rng.ranlux.j24]*0x1p-48;
|
|
}
|
|
}
|
|
|
|
|
|
static void RanluxSkip(This *t, cnumber n)
|
|
{
|
|
RanluxInt(t, n + t->rng.ranlux.nskip*(n/24));
|
|
t->rng.ranlux.n24 = 24 - n % 24;
|
|
}
|
|
|
|
|
|
static inline void RanluxIni(This *t)
|
|
{
|
|
cint skip[] = {24, 48, 97, 223, 389,
|
|
223, 223, 223, 223, 223, 223, 223, 223, 223, 223,
|
|
223, 223, 223, 223, 223, 223, 223, 223, 223, 223};
|
|
state_t seed = t->seed;
|
|
state_t level = RNG;
|
|
count i;
|
|
|
|
if( level < sizeof skip ) level = skip[level];
|
|
t->rng.ranlux.nskip = level - 24;
|
|
|
|
t->rng.ranlux.i24 = 23;
|
|
t->rng.ranlux.j24 = 9;
|
|
t->rng.ranlux.n24 = 24;
|
|
|
|
for( i = 0; i < 24; ++i ) {
|
|
cint k = seed/53668;
|
|
seed = 40014*(seed - k*53668) - k*12211;
|
|
seed += NegQ(seed) & 2147483563;
|
|
t->rng.ranlux.state[i] = seed & ((1 << 24) - 1);
|
|
}
|
|
|
|
t->rng.ranlux.carry = ~TrueQ(t->rng.ranlux.state[23]) & (1 << 24);
|
|
|
|
t->rng.getrandom = RanluxGet;
|
|
t->rng.skiprandom = RanluxSkip;
|
|
}
|
|
|
|
|
|
/*
|
|
PART 4: User routines:
|
|
|
|
- IniRandom sets up the random-number generator to produce a
|
|
sequence of at least n ndim-dimensional random vectors.
|
|
|
|
- GetRandom retrieves one random vector.
|
|
|
|
- SkipRandom skips over n random vectors.
|
|
*/
|
|
|
|
static inline void IniRandom(This *t)
|
|
{
|
|
if( t->seed == 0 ) SobolIni(t);
|
|
else if( RNG == 0 ) MersenneIni(t);
|
|
else RanluxIni(t);
|
|
}
|
|
|