diff options
Diffstat (limited to 'Source/Utility/pcg_basic.cpp')
-rw-r--r-- | Source/Utility/pcg_basic.cpp | 121 |
1 files changed, 121 insertions, 0 deletions
diff --git a/Source/Utility/pcg_basic.cpp b/Source/Utility/pcg_basic.cpp new file mode 100644 index 00000000..27147e2c --- /dev/null +++ b/Source/Utility/pcg_basic.cpp @@ -0,0 +1,121 @@ +//----------------------------------------------------------------------------- +// Name: pcg_basic.cpp +// Description: +// License: Read below +//----------------------------------------------------------------------------- +/* + * PCG Random Number Generation for C. + * + * Copyright 2014 Melissa O'Neill <oneill@pcg-random.org> + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + * For additional information about the PCG random number generation scheme, + * including its license and other licensing options, visit + * + * http://www.pcg-random.org + */ + +/* + * This code is derived from the full C implementation, which is in turn + * derived from the canonical C++ PCG implementation. The C++ version + * has many additional features and is preferable if you can use C++ in + * your project. + */ + +#include "pcg_basic.h" + +// state for global RNGs + +static pcg32_random_t pcg32_global = PCG32_INITIALIZER; + +// pcg32_srandom(initstate, initseq) +// pcg32_srandom_r(rng, initstate, initseq): +// Seed the rng. Specified in two parts, state initializer and a +// sequence selection constant (a.k.a. stream id) + +void pcg32_srandom_r(pcg32_random_t* rng, uint64_t initstate, uint64_t initseq) +{ + rng->state = 0U; + rng->inc = (initseq << 1u) | 1u; + pcg32_random_r(rng); + rng->state += initstate; + pcg32_random_r(rng); +} + +void pcg32_srandom(uint64_t seed, uint64_t seq) +{ + pcg32_srandom_r(&pcg32_global, seed, seq); +} + +// pcg32_random() +// pcg32_random_r(rng) +// Generate a uniformly distributed 32-bit random number + +uint32_t pcg32_random_r(pcg32_random_t* rng) +{ + uint32_t oldstate = (uint32_t) rng->state; + rng->state = oldstate * 6364136223846793005ULL + rng->inc; + uint32_t xorshifted = ((oldstate >> 18u) ^ oldstate) >> 27u; + uint32_t rot = oldstate >> 59u; + return (xorshifted >> rot) | (xorshifted << ((-rot) & 31)); +} + +uint32_t pcg32_random() +{ + return pcg32_random_r(&pcg32_global); +} + + +// pcg32_boundedrand(bound): +// pcg32_boundedrand_r(rng, bound): +// Generate a uniformly distributed number, r, where 0 <= r < bound + +uint32_t pcg32_boundedrand_r(pcg32_random_t* rng, uint32_t bound) +{ + // To avoid bias, we need to make the range of the RNG a multiple of + // bound, which we do by dropping output less than a threshold. + // A naive scheme to calculate the threshold would be to do + // + // uint32_t threshold = 0x100000000ull % bound; + // + // but 64-bit div/mod is slower than 32-bit div/mod (especially on + // 32-bit platforms). In essence, we do + // + // uint32_t threshold = (0x100000000ull-bound) % bound; + // + // because this version will calculate the same modulus, but the LHS + // value is less than 2^32. + + uint32_t threshold = -bound % bound; + + // Uniformity guarantees that this loop will terminate. In practice, it + // should usually terminate quickly; on average (assuming all bounds are + // equally likely), 82.25% of the time, we can expect it to require just + // one iteration. In the worst case, someone passes a bound of 2^31 + 1 + // (i.e., 2147483649), which invalidates almost 50% of the range. In + // practice, bounds are typically small and only a tiny amount of the range + // is eliminated. + for (;;) { + uint32_t r = pcg32_random_r(rng); + if (r >= threshold) + return r % bound; + } +} + + +uint32_t pcg32_boundedrand(uint32_t bound) +{ + return pcg32_boundedrand_r(&pcg32_global, bound); +} + |