View Single Post
Old 12-17-2003, 01:43 PM   #25 (permalink)
dimbulb
Riiiiight........
 
PHP Code:
#ifndef RNG_H
#define RNG_H

/***********************************************************************\
*
* File: RngStream.cpp for multiple streams of Random Numbers
* Language: C++
* Copyright: Pierre L’Ecuyer, University of Montreal
* Notice: This code can be used freely for personnal, academic,
* or non-commercial purposes. For commercial purposes,
* please contact P. L’Ecuyer at: [email]lecuyer@iro.umontreal.ca[/email]
* Date: 14 August 2001
*
\***********************************************************************/
#include <cstdlib>
#include <iostream>
#include "RngStream.h"
using namespace std;
namespace
{
const 
double m1 4294967087.0;
const 
double m2 4294944443.0;
const 
double norm 1.0 / (m1 1.0);
const 
double a12 1403580.0;
const 
double a13n 810728.0;
const 
double a21 527612.0;
const 
double a23n 1370589.0;
const 
double two17 131072.0;
const 
double two53 9007199254740992.0;
const 
double fact 5.9604644775390625e-8/* 1 / 2^24 */
// The following are the transition matrices of the two MRG components
// (in matrix form), raised to the powers -1, 1, 2^76, and 2^127, resp.
const double InvA1[3][3] = { // Inverse of A1p0
184888585.00.01945170933.0 },
1.00.00.0 },
0.01.00.0 }
};
const 
double InvA2[3][3] = { // Inverse of A2p0
0.0360363334.04225571728.0 },
1.00.00.0 },
0.01.00.0 }
};

const 
double A1p0[3][3] = {
0.01.00.0 },
0.00.01.0 },
{ -
810728.01403580.00.0 }
};
const 
double A2p0[3][3] = {
0.01.00.0 },
0.00.01.0 },
{ -
1370589.00.0527612.0 }
};
const 
double A1p76[3][3] = {
82758667.01871391091.04127413238.0 },
3672831523.069195019.01871391091.0 },
3672091415.03528743235.069195019.0 }
};
const 
double A2p76[3][3] = {
1511326704.03759209742.01610795712.0 },
4292754251.01511326704.03889917532.0 },
3859662829.04292754251.03708466080.0 }
};
const 
double A1p127[3][3] = {
2427906178.03580155704.0949770784.0 },
226153695.01230515664.03580155704.0 },
1988835001.0986791581.01230515664.0 }
};
const 
double A2p127[3][3] = {
1464411153.0277697599.01610723613.0 },
32183930.01464411153.01022607788.0 },
2824425944.032183930.02093834863.0 }
};
//-------------------------------------------------------------------------
// Return (a*s + c) MOD m; a, s, c and m must be < 2^35
//
double MultModM (double adouble sdouble cdouble m)
{
double v;
long a1;
=*+c;


if (
>= two53 || <= -two53) {
a1 static_cast<long> (two17); -= a1 two17;
a1 *s;
a1 static_cast<long> (m); -= a1 m;
two17 +*+c;
}
a1 static_cast<long> (m);
/* in case v < 0)*/
if ((-= a1 m) < 0.0) return += m; else return v;
}
//-------------------------------------------------------------------------
// Compute the vector v = A*s MOD m. Assume that -m < s[i] < m.
// Works also when v = s.
//
void MatVecModM (const double A[3][3], const double s[3], double v[3],
double m)
{
int i;
double x[3]; // Necessary if v =s
for (=0;<3; ++i) {
x[i] = MultModM (A[i][0], s[0], 0.0m);
x[i] = MultModM (A[i][1], s[1], x[i], m);
x[i] = MultModM (A[i][2], s[2], x[i], m);
}
for (
=0;<3; ++i)
v[i] = x[i];
}
//-------------------------------------------------------------------------
// Compute the matrix C = A*B MOD m. Assume that -m < s[i] < m.
// Note: works also if A =C or B =C or A =B =C.
//
void MatMatModM (const double A[3][3], const double B[3][3],
double C[3][3], double m)
{
int ij;
double V[3], W[3][3];
for (
=0;<3; ++i) {
for (
=0;<3; ++j)
V[j] = B[j][i];
MatVecModM (AVVm);
for (
=0;<3; ++j)

W[j][i] = V[j];
}
for (
=0;<3; ++i)
for (
=0;<3; ++j)
C[i][j] = W[i][j];
}
//-------------------------------------------------------------------------
// Compute the matrix B = (A^(2^e) Mod m); works also if A =B.
//
void MatTwoPowModM (const double A[3][3], double B[3][3], double mlong e)
{
int ij;
/* initialize: B =A */
if (!= B) {
for (
=0;<3; ++i)
for (
=0;<3; ++j)
B[i][j] = A[i][j];
}
/* Compute B = A^(2^e) mod m */
for (0ei++)
MatMatModM (BBBm);
}
//-------------------------------------------------------------------------
// Compute the matrix B = (A^n Mod m); works even if A =B.
//
void MatPowModM (const double A[3][3], double B[3][3], double mlong n)
{
int ij;
double W[3][3];
/* initialize: W = A; B = I */
for (=0;<3; ++i)
for (
=0;<3; ++j) {
W[i][j] = A[i][j];
B[i][j] = 0.0;
}
for (
=0;<3; ++j)
B[j][j] = 1.0;
/* Compute B = A^n mod m using the binary decomposition of n */
while (0) {
if (
2MatMatModM (WBBm);
MatMatModM (WWWm);

/=2;
}
}
//-------------------------------------------------------------------------
// Check that the seeds are legitimate values. Returns 0 if legal seeds,
// -1 otherwise.
//
int CheckSeed (const unsigned long seed[6])
{
int i;
for (
=0;<3; ++i) {
if (
seed[i] >= m1) {
cerr << "****************************************\n\n"
<< "ERROR: Seed[" << << "] >= 4294967087, Seed is not set."
<< "\n\n****************************************\n\n";
return (-
1);
}
}
for (
=3;<6; ++i) {
if (
seed[i] >= m2) {
cerr << "*****************************************\n\n"
<< "ERROR: Seed[" << << "] >= 4294944443, Seed is not set."
<< "\n\n*****************************************\n\n";
return (-
1);
}
}
if (
seed[0] == && seed[1] == && seed[2] == 0) {
cerr << "****************************\n\n"
<< "ERROR: First 3 seeds = 0.\n\n"
<< "****************************\n\n";
return (-
1);
}
if (
seed[3] == && seed[4] == && seed[5] == 0) {
cerr << "****************************\n\n"
<< "ERROR: Last 3 seeds = 0.\n\n"
<< "****************************\n\n";
return (-
1);
}
return 
0;
}
// end of anonymous namespace


//-------------------------------------------------------------------------
// Generate the next random number.
//
double RngStream::U01 ()
{
long k;
double p1p2u;
/* Component 1 */
p1 a12 Cg[1] - a13n Cg[0];
static_cast<long> (p1 m1);
p1 -= m1;
if (
p1 0.0p1 += m1;
Cg[0] = Cg[1]; Cg[1] = Cg[2]; Cg[2] = p1;
/* Component 2 */
p2 a21 Cg[5] - a23n Cg[3];
static_cast<long> (p2 m2);
p2 -= m2;
if (
p2 0.0p2 += m2;
Cg[3] = Cg[4]; Cg[4] = Cg[5]; Cg[5] = p2;
/* Combination */
= ((p1 p2) ? (p1 p2) * norm : (p1 p2 m1) * norm);
return (
anti == false) ?:(-u);
}
//-------------------------------------------------------------------------
// Generate the next random number with extended (53 bits) precision.
//
double RngStream::U01d ()
{
double u;
U01();
if (
anti) {
// Don’t forget that U01() returns 1 -u in the antithetic case
+= (U01() - 1.0) * fact;
return (
0.0) ?+1.0 :u;
} else {
+= U01() * fact;
return (
1.0) ?:(1.0);
}
}
//*************************************************************************

// Public members of the class start here
//-------------------------------------------------------------------------
// The default seed of the package; will be the seed of the first
// declared RngStream, unless SetPackageSeed is called.
//
double RngStream::nextSeed[6] =
{
12345.012345.012345.012345.012345.012345.0
};
//-------------------------------------------------------------------------
// constructor
//
RngStream::RngStream (const char *s) : name (s)
{
anti false;
incPrec false;
/* Information on a stream. The arrays {Cg, Bg, Ig} contain the current
state of the stream, the starting state of the current SubStream, and the
starting state of the stream. This stream generates antithetic variates
if anti = true. It also generates numbers with extended precision (53
bits if machine follows IEEE 754 standard) if incPrec = true. nextSeed
will be the seed of the next declared RngStream. */
for (int i 06; ++i) {
Bg[i] = Cg[i] = Ig[i] = nextSeed[i];
}
MatVecModM (A1p127nextSeednextSeedm1);
MatVecModM (A2p127, &nextSeed[3], &nextSeed[3], m2);
}
//-------------------------------------------------------------------------
// Reset Stream to beginning of Stream.
//
void RngStream::ResetStartStream ()
{
for (
int i 06; ++i)
Cg[i] = Bg[i] = Ig[i];
}
//-------------------------------------------------------------------------


// Reset Stream to beginning of SubStream.
//
void RngStream::ResetStartSubstream ()
{
for (
int i 06; ++i)
Cg[i] = Bg[i];
}
//-------------------------------------------------------------------------
// Reset Stream to NextSubStream.
//
void RngStream::ResetNextSubstream ()
{
MatVecModM(A1p76BgBgm1);
MatVecModM(A2p76, &Bg[3], &Bg[3], m2);
for (
int i 06; ++i)
Cg[i] = Bg[i];
}
//-------------------------------------------------------------------------
void RngStream::SetPackageSeed (const unsigned long seed[6])
{
if (
CheckSeed (seed)) exit (EXIT_FAILURE);
for (
int i 06; ++i)
nextSeed[i] = seed[i];
}
//-------------------------------------------------------------------------
void RngStream::SetSeed (const unsigned long seed[6])
{
if (
CheckSeed (seed)) exit (EXIT_FAILURE);
for (
int i 06; ++i)
Cg[i] = Bg[i] = Ig[i] = seed[i];
}
//-------------------------------------------------------------------------
// if e > 0, let n = 2^e + c;
// if e < 0, let n = -2^(-e) + c;
// if e = 0, let n = c.
// Jump n steps forward if n >0, backwards if n <0.
//
void RngStream::AdvanceState (long elong c)
{
double B1[3][3], C1[3][3], B2[3][3], C2[3][3];

if (
0) {
MatTwoPowModM (A1p0B1m1e);
MatTwoPowModM (A2p0B2m2e);
} else if (
0) {
MatTwoPowModM (InvA1B1m1, -e);
MatTwoPowModM (InvA2B2m2, -e);
}
if (
>= 0) {
MatPowModM (A1p0C1m1c);
MatPowModM (A2p0C2m2c);
} else {
MatPowModM (InvA1C1m1, -c);
MatPowModM (InvA2C2m2, -c);
}
if (
e) {
MatMatModM (B1C1C1m1);
MatMatModM (B2C2C2m2);
}
MatVecModM (C1CgCgm1);
MatVecModM (C2, &Cg[3], &Cg[3], m2);
}
//-------------------------------------------------------------------------
void RngStream::GetState (unsigned long seed[6]) const
{
for (
int i 06; ++i)
seed[i] = static_cast<unsigned long> (Cg[i]);
}
//-------------------------------------------------------------------------
void RngStream::WriteState () const
{
cout << "The current state of the Rngstream";
if (
name.size() > 0)
cout <<""<< name;
cout << ":\n Cg ={";
for (
int i 0;<5i++) {
cout << static_cast<unsigned long> (Cg [i]) << ", ";
}
cout << static_cast<unsigned long> (Cg [5]) << " }\n\n";
}

//-------------------------------------------------------------------------
void RngStream::WriteStateFull () const
{
int i;
cout << "The RngStream";
if (
name.size() > 0)
cout <<""<< name;
cout << ":\n anti = " << (anti "true" "false") << "\n";
cout << " incPrec = " << (incPrec "true" "false") << "\n";
cout << " Ig = { ";
for (
=0;<5i++) {
cout << static_cast<unsigned long> (Ig [i]) << ", ";
}
cout << static_cast<unsigned long> (Ig [5]) << " }\n";
cout << " Bg = { ";
for (
=0;<5i++) {
cout << static_cast<unsigned long> (Bg [i]) << ", ";
}
cout << static_cast<unsigned long> (Bg [5]) << " }\n";
cout << " Cg = { ";
for (
=0;<5i++) {
cout << static_cast<unsigned long> (Cg [i]) << ", ";
}
cout << static_cast<unsigned long> (Cg [5]) << " }\n\n";
}
//-------------------------------------------------------------------------
void RngStream::IncreasedPrecis (bool incp)
{
incPrec incp;
}
//-------------------------------------------------------------------------
void RngStream::SetAntithetic (bool a)
{
anti a;
}
//-------------------------------------------------------------------------
// Generate the next random number.
//

double RngStream::RandU01 ()
{
if (
incPrec)
return 
U01d();
else
return 
U01();
}
//-------------------------------------------------------------------------
// Generate the next random integer.
//
long RngStream::RandInt (long lowlong high)
{
return 
low static_cast<long> ((high low 1) * RandU01 ());
};

//--End of Class Implementation-------------------------------

#endif 
dimbulb is offline  
 

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73