Appendix B: CLCG4 Random Number Generator (C# Language Version)
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
namespace RandomGenerator
{
public class RandomGenerator
{
public static long Maxgen = 100;
public static long i;
public static long j;
public static long [] aw = new long[4];
public static long [] avw = new long[4];
public static long [] a = new long[4] {45991, 207707, 138556, 49689};
public static long [] m = new long[4] {2147483647, 2147483543, 2147483423, 2147483323);
public static long [,] Ig = new long[4, 101];
public static long [,] Lg = new long[4, 101];
public static long [,] Cg = new long[4, 101];
public static long H = 32768; // - 2^15: use 1n MultModM.
public enum SeedType
{
InitialSeed,
LastSeed,
NewSeed,
}
public static int gSeedl;
public static int gSeed2;
public static int gSeed3;
public static int gSeed4;
/************************************************************************
This is the Random Number Generator Based on the Combination of Four LCGs (Linear Congruential Generators) by Pierre L'euyer and Terry H. Andres.
*
* This is a portable package for uniform random number generation based on a backbone generator with period length near 2^121, which is a combination
* of four linear congruential generators. The package provides for multiple
* (virtual) generators evolving in parallel. Each generator also has many
* disjoint subsequences, and software tools are provided to reset the state
* of any generator to the beginning of its first, previous, or current
* subsequence. Such facilities are helpful to maintain synchronization for
* implementing variance reduction methods in simulation.
* For those who are interested in the details or the theories behind the
* implementation, you ran refer to the following papers by Prof. L'Ecuyer:
" 1. A Random Number Generator Based on the Combination of Four LCGs.
(This can be downloaded from the "papers" page of Prof. L'Ecuyer's web page: http://www.iro.umontreal.ca/~lecuyer/myftp/papers/clcg4.ps)
• 2. Implementing a Random Number Package with Splitting Facilities
(ACM Transactions on Math. Software, Vol. 17, No. 1, March 91, p. 98-111)
• 3. Efficient and portable combined random number generator5.
(Communications of the ACM, 31(6): 742-749 and 774, 1988.)
* and the papers referred to therein.
public static float [] getBlockOfRandomNumbers(int seedl, int seed2, int seed3, int seed4, int iCount, int iTestSwitch)
{
/ ************************************************************************
*
* The Seeds MUST satisfy the following constraints:
long[] s = new long[4];
if (iTestSwitch == 1)
{
// TEST Random Number Generator
//1) Initialize the CLCG4 generator using the defaults.
//You will have access to 100 sub-generators in an array indexed from 0 to 99.
//2) Advance the generator 20,000 numbers for generators 0, 1, 2, 3
//3) The next five numbers from each generator segment and they should match the numbers below.
// Skipping 20,000 numbers for generator 0 - 3...
Gen 0
Gen 1
Gen 2
Gen 3
// 1
0.5466955
0.8833518
0.8046209
0.7498480
// 2
0.3769046
0.9406571
0.7508963
0.0273824
// 3
0.8909396
0.5905970
0.6122372
0.0797647
// 4
0.8216189
0.1320690
0.2008421
0.4063177
// 5
0.0613126
0.9316329
0.2755311
0.1423134
gSeed1 = 11111111; //<-- Initialize using the defaults of the TEST Requirements.
gSeed2 = 22222222; //<-- Initialize using the defaults of the TEST Requirements.
gSeed3 = 33333333; //<-- Initialize using the defaults of the TEST Requirements.
gSeed4 = 44444444; //<-- Initialize using the defaults of the TEST Requirements.
iCount = 5; //<-- This required by the Rules of the TEST Requirements
float[] fReturnBlock2 = new float[20]; //<-- 20 to return 4 Blocks of 5 from the 4 Generator Segments for TESTING.
InitDefault();
// Advance the generator 20,000 numbers for generators 0, 1, 2, 3
for (i = 1; i <= 20000; i++) GenVal(0);
for (i = 1; i <= 20000; i++) GenVal(1);
for (i = 1; i <= 20000; i++) GenVal(2);
for (i = 1; i <= 20000; i++) GenVal(3);
//IF want to see the next 5 Test RAW Result Values for each of the 4 Generator Segments output in a Console Applicaton in which you have reference the RandomGenertor.dll then un-remark the 11 lines below.
//Console.WriteLine("Gen 0");
//for (i = 1; i <= iCount; i++) Console.WriteLine("{0:F2}", i + ": " + GenVal(0));
//Console.WriteLine("----------------");
//Console.WriteLine("Gen 1");
//for (i = 1; i <= iCount; i++) Console.WriteLine("{0:F2}", i + ": " + GenVal(1));
//Console.WriteLine("----------------");
//Console.WriteLine("Gen 2");
//for (i = 1; i <= iCount; i++) Console.WriteLine("{0:F2}", i + ": " + GenVal(2));
//Console.WriteLine("----------------");
//Console.WriteLine("Gen 3");
//for (i = 1; i <= iCount; i++) Console.WriteLine("{0:F2}", i + ": " + GenVal(3));
//Return the next 5 Test Result Values for each of the 4 Generator Segments
for (i = 0; i < iCount; i++)
{
fReturnBlock2[i] = (float)GenVal(0);
}
for (i = 0; i < iCount; i++)
{
fReturnBlock2[i + 5] = (float)GenVal(1);
}
for (i = 0; i < iCount; i++)
{
fReturnBlock2[i + 10] = (float)GenVal(2);
}
for (i = 0; i < iCount; i++)
{
fReturnBlock2[i + 15] = (float)GenVal(3);
}
for (i = 0; i < iCount; i++)
{
fReturnBlock2[i] = (float)GenVal(0);
}
return fReturnBlock2;
}
static long MultModM(long s, long t, long M)
/* Returns (s*t) MOD M. Assumes that -M < s < M and -M < t < M.
/* See L'Ecuyer and Cote (1991).
{
long R, S0, S1, q, qh, rh, k;
if (s < 0) s += M;
if (t < 0) t += M;
if (s < H) {S0 = s; R = 0;}
else
{
S1 = s / H; S0 = s - H * S1;
qh = M / H; rh = M - H * qh;
if (S1 >= H)
{
S1 -= H; k = t / qh; R = H * (t - k * qh) - k * rh;
while (R < 0) R += M;
}
else R = 0;
if (S1 != 0)
{
q = M / S1; k = t / q; R -= k * (M - S1 * q);
if (R > 0) R -= M;
R += S1 * (t - k * q);
while (R < 0) R += M;
}
k = R / qh; R = H * (R - k * qh) - k * rh;
while (R < 0) R += M;
}
if (S0 != 0)
{
q = M / S0; k = t / q; R -= k * (M - S0 * q);
if (R > 0) R -= M;
R += S0 * (t - k * q);
while (R < 0) R += M;
}
public static void SetSeed(short g, long[] s)
{
if (g > Maxgen) Console.WriteLine("ERROR: SetSeed with g > Maxgen \n");
for (j = 0; j < 4; j++) Ig[j, g] = s[j];
InitGenerator(g, SeedType.InitialSeed);
}
public static void WriteState(short g)
{
Console.WriteLine("\n State of generator g = %u :", g);
for (j = 0; j < 4; j++)
{
Console.WriteLine("\n Cg[%u] = %lu", j, Cg[j, g]);
}
Console.WriteLine("\n");
}
public static void GetState(short g, long[] s)
{
for (j = 0; j < 4; j++) s[j] = Cg[j, g];
}
public static void InitGenerator(long g, SeedType Where)
{
if (g > Maxgen) Console.WriteLine("ERROR: InitGenerator with g > Maxgen \n");
for (j = 0; j < 4; j++)
{
switch (Where)
{
case SeedType.InitialSeed: Lg[j, g] = Ig[j, g]; break; case SeedType.NewSeed: Lg[j, g] = MultModM(aw[j], Lg[j, g], m[j]); break; case SeedType.LastSeed: break;
}
Cg[j, g] = Lg[j, g];
}
}
public static void SetInitialSeed(long[] s)
{
long g;
for (j = 0; j < 4; j++) Ig[j, 0] = s[j];
InitGenerator(0, SeedType.InitialSeed);
for (g = 1; g <= Maxgen; g++)
{
for (j = 0; j < 4; j++) Ig[j, g] = MultModM(avw[j], Ig[j, (g - 1)], m[j]);
InitGenerator(g, SeedType.InitialSeed);
}
}
public static void Init(long v, long w)
{
long[] sd = new long[4] {gSeed1, gSeed2, gSeed3, gSeed4};
for (j = 0; j < 4; j++)
{
aw[j] = a[j];
for (i = 1; i <= w; i++) aw[j] = MultModM(aw[j], aw[j], m[j]);
avw[j] = aw[j];
for (i = 1; i <= v; i++) avw[j] = MultModM(avw[j], avw[j], m[j]);
}
SetInitialSeed(sd);
}
public static double GenVal(short g)
{
long k, s;
double u;
u = 0.0;
if (g > Maxgen) Console.WriteLine("ERROR: Genval with g > Maxgen \n");
s = Cg[0, g]; k = s / 46693;
s = 45991 * (s - k * 46693) - k * 25884;
if (s < 0) s = s + 2147483647; Cg[0, g] = s;
u = u + 4.65661287524579692e-10 * s;
s = Cg[1, g]; k = s / 10339;
s = 207707 * (s - k * 10339) - k * 870;
if (s < 0) s = s + 2147483543; Cg[1, g] = s;
u = u - 4.65661310075985993e-10 * s;
if (u < 0) u = u + 1.0;
s = Cg[2, g]; k = s / 15499;
s = 138556 * (s - k * 15499) - k * 3979;
if (s < 0.0) s = s + 2147483423; Cg[2, g] = s;
u = u + 4.65661336096842131e-10 * s;
if (u >= 1.0) u = u - 1.0;
s = Cg[3, g]; k = s / 43218;
s = 49689 * (s - k * 43218) - k * 24121;
if (s < 0) s = s + 2147483323; Cg[3, g] = s;
u = u - 4.65661357780891134e-10 * s;
if (u < 0) u = u + 1.0;
return (u);
}
public static void InitDefault()
{
Init(31, 41);
}
}
}