From 6811a197325a094db06a051860c85daf438b2435 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Fri, 1 Dec 2017 17:32:27 -0500 Subject: [PATCH] Use unsigned integers in the RNG --- src/random_lcg.C | 49 +++++++++++++++++++++++------------------------- src/random_lcg.h | 10 +++++----- 2 files changed, 28 insertions(+), 31 deletions(-) diff --git a/src/random_lcg.C b/src/random_lcg.C index c78c40748d..8383401f9b 100644 --- a/src/random_lcg.C +++ b/src/random_lcg.C @@ -6,12 +6,12 @@ int64_t seed = 1; // LCG parameters -const int64_t prn_mult = 2806196910506780709LL; // multiplication factor, g -const int64_t prn_add = 1; // additive factor, c -const int64_t prn_mod = -0x8000000000000000; // -2^63 -const int64_t prn_mask = 0x7fffffffffffffff; // 2^63 - 1 -const int64_t prn_stride = 152917LL; // stride between particles -const double prn_norm = pow(2, -63); // 2^-63 +const uint64_t prn_mult = 2806196910506780709LL; // multiplication factor, g +const uint64_t prn_add = 1; // additive factor, c +const uint64_t prn_mod = 0x8000000000000000; // 2^63 +const uint64_t prn_mask = 0x7fffffffffffffff; // 2^63 - 1 +const uint64_t prn_stride = 152917LL; // stride between particles +const double prn_norm = pow(2, -63); // 2^-63 // Module constants const int N_STREAMS = 5; @@ -22,8 +22,8 @@ const int STREAM_URR_PTABLE = 3; const int STREAM_VOLUME = 4; // Current PRNG state -int64_t prn_seed[N_STREAMS]; // current seed -int stream; // current RNG stream +uint64_t prn_seed[N_STREAMS]; // current seed +int stream; // current RNG stream #pragma omp threadprivate(prn_seed, stream) @@ -40,8 +40,7 @@ prn() // Once the integer is calculated, we just need to divide by 2**m, // represented here as multiplying by a pre-calculated factor - double pseudo_rn = prn_seed[stream] * prn_norm; - return pseudo_rn; + return prn_seed[stream] * prn_norm; } //============================================================================== @@ -49,10 +48,9 @@ prn() //============================================================================== extern "C" double -future_prn(int64_t n) +future_prn(uint64_t n) { - double pseudo_rn = future_seed(n, prn_seed[stream]) * prn_norm; - return pseudo_rn; + return future_seed(n, prn_seed[stream]) * prn_norm; } //============================================================================== @@ -60,7 +58,7 @@ future_prn(int64_t n) //============================================================================== extern "C" void -set_particle_seed(int64_t id) +set_particle_seed(uint64_t id) { for (int i = 0; i < N_STREAMS; i++) { prn_seed[i] = future_seed(id * prn_stride, seed + i); @@ -72,7 +70,7 @@ set_particle_seed(int64_t id) //============================================================================== extern "C" void -advance_prn_seed(int64_t n) +advance_prn_seed(uint64_t n) { prn_seed[stream] = future_seed(n, prn_seed[stream]); } @@ -81,16 +79,16 @@ advance_prn_seed(int64_t n) // FUTURE_SEED //============================================================================== -extern "C" int64_t -future_seed(int64_t n, int64_t seed) +extern "C" uint64_t +future_seed(uint64_t n, uint64_t seed) { // In cases where we want to skip backwards, we add the period of the random // number generator until the number of PRNs to skip is positive since // skipping ahead that much is the same as skipping backwards by the original // amount. - int64_t nskip = n; - while (nskip < 0) nskip += prn_mod; + uint64_t nskip = n; + while (nskip > prn_mod) nskip += prn_mod; // Make sure nskip is less than 2^M. nskip &= prn_mask; @@ -102,10 +100,10 @@ future_seed(int64_t n, int64_t seed) // and C which can then be used to find x_N = G*x_0 + C mod 2^M. // Initialize constants - int64_t g = prn_mult; - int64_t c = prn_add; - int64_t g_new = 1; - int64_t c_new = 0; + uint64_t g = prn_mult; + uint64_t c = prn_add; + uint64_t g_new = 1; + uint64_t c_new = 0; while (nskip > 0) { // Check if the least significant bit is 1. @@ -121,8 +119,7 @@ future_seed(int64_t n, int64_t seed) } // With G and C, we can now find the new seed. - int64_t new_seed = (g_new * seed + c_new) & prn_mask; - return new_seed; + return (g_new * seed + c_new) & prn_mask; } //============================================================================== @@ -140,7 +137,7 @@ prn_set_stream(int i) //============================================================================== extern "C" int -openmc_set_seed(int64_t new_seed) +openmc_set_seed(uint64_t new_seed) { seed = new_seed; #pragma omp parallel diff --git a/src/random_lcg.h b/src/random_lcg.h index 90d1c253ef..889195f109 100644 --- a/src/random_lcg.h +++ b/src/random_lcg.h @@ -14,21 +14,21 @@ extern "C" double prn(); // current seed. //============================================================================== -extern "C" double future_prn(int64_t n); +extern "C" double future_prn(uint64_t n); //============================================================================== // SET_PARTICLE_SEED sets the seed to a unique value based on the ID of the // particle. //============================================================================== -extern "C" void set_particle_seed(int64_t id); +extern "C" void set_particle_seed(uint64_t id); //============================================================================== // ADVANCE_PRN_SEED advances the random number seed 'n' times from the current // seed. //============================================================================== -extern "C" void advance_prn_seed(int64_t n); +extern "C" void advance_prn_seed(uint64_t n); //============================================================================== // FUTURE_SEED advances the random number seed 'skip' times. This is usually @@ -37,7 +37,7 @@ extern "C" void advance_prn_seed(int64_t n); // are used. //============================================================================== -extern "C" int64_t future_seed(int64_t n, int64_t seed); +extern "C" uint64_t future_seed(uint64_t n, uint64_t seed); //============================================================================== // PRN_SET_STREAM changes the random number stream. If random numbers are needed @@ -51,6 +51,6 @@ extern "C" void prn_set_stream(int n); // API FUNCTIONS //============================================================================== -extern "C" int openmc_set_seed(int64_t new_seed); +extern "C" int openmc_set_seed(uint64_t new_seed); #endif // RANDOM_LCG_H