Skip to content
Snippets Groups Projects
Verified Commit 1cc3748a authored by Dorian Stoll's avatar Dorian Stoll
Browse files

zellularautomat: c: Drop custom RNG implementation

parent dc04e5aa
No related branches found
No related tags found
No related merge requests found
......@@ -16,10 +16,6 @@
#include "openssl/evp.h"
#include "ca_common.h"
#include "random.h"
/* determine random integer between 0 and n-1 */
#define randInt(n) ((int)(nextRandomLEcuyer() * n))
void ca_init(int argc, char** argv, int *lines, int *its)
{
......@@ -36,18 +32,18 @@ void ca_init_config(line_t *buf, int lines, int skip_lines)
{
volatile int scratch;
initRandomLEcuyer(424243);
srand(424243);
/* let the RNG spin for some rounds (used for distributed initialization) */
for (int y = 1; y <= skip_lines; y++) {
for (int x = 1; x <= XSIZE; x++) {
scratch = scratch + randInt(100) >= 50;
scratch = scratch + (rand() % 2);
}
}
for (int y = 1; y <= lines; y++) {
for (int x = 1; x <= XSIZE; x++) {
buf[y][x] = randInt(100) >= 50;
buf[y][x] = rand() % 2;
}
}
}
......
......@@ -3,7 +3,6 @@ cc = meson.get_compiler('c')
sources = [
'ca_common.c',
'ca_openmp.c',
'random.c',
]
dependencies = [
......
#include "random.h"
/* (c) 1996,1997 Thomas Worsch, Peter Sanders */
/* =====================================================================
* The pseudo random number generator functions in this file are
* similar to those in
* Numerical Recipes in C, Second Edition, pages 279-282.
*
* The RNGs have names nextRandom<something> and they return a Float64.
* They are initialized by a call to initRandom<something> which
* take an Int32 as seed.
*/
#define EPS 1.2e-7
#define RNMX (1.0-EPS)
#define IM 2147483647
#define AM ((Float64)1.0/IM)
#define IA 16807
#define IQ 127773
#define IR 2836
static Int32 state = 123456789;
void initRandomParkMiller(Int32 seed)
{
state = seed;
/* but we have to make sure that state never ever is set to zero */
if (state==0) { state = 42; }
}
Float64 nextRandomParkMiller(void)
{
Int32 k;
Float64 result;
k = state/IQ;
state = IA*(state-k*IQ)-k*IR;
if (state < 0) { state += IM; }
result = AM*state;
if (result >= 1.0) { result = RNMX; }
return result;
}
/* =====================================================================
* Now a pseudo RNG with a longer period (roughly 10^18) by lEcuyer
* (Numerical Recipes, page 282).
*/
#define IM1 2147483563
#define IM2 2147483399
#define AM1 ((Float64)1.0/IM1)
#define IMM1 (IM1-1)
#define IA1 40014
#define IA2 40692
#define IQ1 53668
#define IQ2 52774
#define IR1 12211
#define IR2 3791
#define NTAB 32
#define NDIV (1+IMM1/NTAB)
static Int32 state1 = 987654321;
static Int32 state2;
static Int32 y;
static Int32 v[NTAB];
/* ------------------------------------------------------------------ */
static void initRandomSeedLEcuyer(Int32 seed)
{
state1 = seed;
if (state1==0) { state1 = 987654321; }
state2 = state1;
}
/* ------------------------------------------------------------------ */
static void initRandomTabLEcuyer(void)
{
Int32 j, k;
for (j=NTAB+7; j>=0; j--) {
k = state1/IQ1;
state1 = IA1*(state1-k*IQ1)-k*IR1;
if (state1 < 0) { state1 += IM1; }
if (j < NTAB) { v[j] = state1; }
}
y = v[0];
}
/* ------------------------------------------------------------------ */
void initRandomLEcuyer(Int32 seed)
{
initRandomSeedLEcuyer(seed);
initRandomTabLEcuyer();
}
/* ------------------------------------------------------------------ */
static Int32 power(Int32 base, Card64 exp, Int32 modulus)
{
Int64 temp = 1;
Card64 mask;
if (base < 0) { return 0; }
/* note that at on each entry into the following loop body, the
actual value of temp is always positive and fits into an Int32 */
for (mask = ((Card64)1) << 63; mask != 0; mask >>= 1) {
temp = (temp * temp) % modulus;
if (exp & mask) {
temp = (temp * base) % modulus;
}
}
return ((Int32) temp);
}
/* ------------------------------------------------------------------ */
static void forwardRandomLEcuyer(Card64 steps)
{
Int32 a;
a = power(IA1, steps, IM1);
state1 = (Int32) ( (((Int64)a) * state1) % IM1);
a = power(IA2, steps, IM2);
state2 = (Int32) ( (((Int64)a) * state2) % IM2);
}
/* ------------------------------------------------------------------ */
/*
* The following initialization function is intended for use on a
* parallel machine.
* Each PE must call initParallelRandomLEcuyer once at the beginning.
* Calls to nextRandomLEcuyer will then return different (and hopefully
* pseudo unrelated) pseudo random numbers on different PEs.
* The calls to initParallelRandomLEcuyer on different PEs
* *MUST* *SATISFY* *THE* *FOLLOWING* *CONDITIONS*:
* - The parameter seed is the same on all PEs.
* - The parameter total is the same on all PEs and it must be
* the total number of PEs using the RNG.
* - The parameter pe must be different on each PE and for each
* i in the set {0,1,...,total-1} there must be exactly one PE
* calling initParallelRandomLEcuyer with parameter pe set to i.
* Please note:
* - There are *NO* run time checks to see whether the above
* conditions have been satisfied.
* - The numbers generated on different PEs are probably only more or
* less unrelated as long as the number of calls of nextRandomLEcuyer
* on each of the PEs is smaller than the period of the RNG
* (approx. 10^18) divided by the number of PEs (probably below 10^4).
* Since 10^14 calls of nextRandomLEcuyer will take quite some time,
* you are probably safe using this RNG.
* - This function uses nextRandomParkMiller. If you want to have
* reproducible results, make sure that you do not use
* RandomParkMiller yourself before calling initParallelRandomLEcuyer.
*/
void initParallelRandomLEcuyer(Int32 seed, Int32 pe, Int32 total)
{
Card64 steps;
initRandomSeedLEcuyer(seed);
/* The period of the RNG is roughly 2.3e18, i.e. 2^61,
which should be distributed onto the PEs approximately equally;
because we do not know the exact value we are careful and take
one half of the average length of the interval per PE: */
steps = (((Card64)1) << 60)/total;
/* For PE number pe we get the starting point: */
steps = steps * pe;
/* Finally the starting point for each PE is randomly shifted
by an amount which is small compared to the length of
its interval (as long as there are much less than 2^30 PEs :-).
Therefore steps will still be far below the end of its interval: */
steps = steps + (Card64) (nextRandomParkMiller() * (((Card64)1) << 30));
/* Now the RNG is initialized for PE pe as if it had already made steps
many steps from the initial seed. */
forwardRandomLEcuyer(steps);
initRandomTabLEcuyer();
}
/* ------------------------------------------------------------------ */
Float64 nextRandomLEcuyer(void)
{
Int32 k;
Float64 result;
int j;
k = state1/IQ1;
state1 = IA1*(state1-k*IQ1)-k*IR1;
if (state1 < 0) { state1 += IM1; }
k = state2/IQ2;
state2 = IA2*(state2-k*IQ2)-k*IR2;
if (state2 < 0) { state2 += IM2; }
j = y/NDIV;
y = v[j] - state2;
v[j] = state1;
if (y < 1) { y += IMM1; }
result = AM1*y;
if (result >= 1.0) { result = RNMX; }
return result;
}
#include <limits.h>
/* (c) 1996,1997 Thomas Worsch, Peter Sanders */
/* =====================================================================
* The pseudo random number generator functions in this file are
* similar to those in
* Numerical Recipes in C, Second Edition, pages 279-282.
*
* The RNGs have names nextRandom<something> and they return a Float64.
* They are initialized by a call to initRandom<something> which
* take an Int32 as seed.
*/
/* C++ compatibility */
#ifdef __cplusplus
#define CC extern "C"
#else
#define CC
#endif
/* try to find out how 32 bit and 64 bit
* interger types look like in this compiler
* this may fail
* e.g., if the compiler does not support 64 data types...
*/
#if UINT_MAX >> 31 == 1
typedef int Int32;
typedef unsigned int Card32;
#elif USHRT_MAX >> 31 == 1
typedef short Int32;
typedef unsigned short Card32;
#elif ULONG_MAX >> 31 == 1
typedef short Int32;
typedef unsigned short Card32;
#else /* provoke error */
typedef nonexisting Int32;
typedef nonexisting Card32;
#endif
#if UINT_MAX >> 63 == 1
typedef int Int64;
typedef unsigned int Card64;
#elif ULONG_MAX >> 63 == 1
typedef long Int64;
typedef unsigned long Card64;
#else
typedef long long int Int64;
typedef unsigned long long int Card64;
#endif
typedef double Float64;
/* =====================================================================
* The Minimal Standard pseudo RNG (Numerical Recipes, page 279)
* by Park and Miller
* I added an initialization function and could therefore remove
* the MASK mechanism from the original ran0 function.
*/
CC void initRandomParkMiller(Int32 seed);
CC Float64 nextRandomParkMiller (void);
/* =====================================================================
* Now a pseudo RNG with a longer period (roughly 10^18) by lEcuyer
* (Numerical Recipes, page 282).
*/
CC void initRandomLEcuyer(Int32 seed);
CC Float64 nextRandomLEcuyer (void);
/* ------------------------------------------------------------------ */
/*
* The following initialization function is intended for use on a
* parallel machine.
* Each PE must call initParallelRandomLEcuyer once at the beginning.
* Calls to nextRandomLEcuyer will then return different (and hopefully
* pseudo unrelated) pseudo random numbers on different PEs.
* The calls to initParallelRandomLEcuyer on different PEs
* *MUST* *SATISFY* *THE* *FOLLOWING* *CONDITIONS*:
* - The parameter seed is the same on all PEs.
* - The parameter total is the same on all PEs and it must be
* the total number of PEs using the RNG.
* - The parameter pe must be different on each PE and for each
* i in the set {0,1,...,total-1} there must be exactly one PE
* calling initParallelRandomLEcuyer with parameter pe set to i.
* Please note:
* - There are *NO* run time checks to see whether the above
* conditions have been satisfied.
* - The numbers generated on different PEs are probably only more or
* less unrelated as long as the number of calls of nextRandomLEcuyer
* on each of the PEs is smaller than the period of the RNG
* (approx. 10^18) divided by the number of PEs (probably below 10^4).
* Since 10^14 calls of nextRandomLEcuyer will take quite some time,
* you are probably safe using this RNG.
* - This function uses nextRandomParkMiller. If you want to have
* reproducible results, make sure that you do not use
* RandomParkMiller yourself before calling initParallelRandomLEcuyer.
*/
CC void initParallelRandomLEcuyer(Int32 seed, Int32 pe, Int32 total);
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment