From 508c4ab81a7068986ae818768b1444bbbd83b9c7 Mon Sep 17 00:00:00 2001 From: rtk0c Date: Sat, 17 May 2025 14:13:10 -0700 Subject: Wrap random source Havne't dont too much investigation if PCG is a good choice; this give some freedom to move later --- src/pcg.hpp | 207 ------------------------------------------------------------ 1 file changed, 207 deletions(-) delete mode 100644 src/pcg.hpp (limited to 'src/pcg.hpp') diff --git a/src/pcg.hpp b/src/pcg.hpp deleted file mode 100644 index 8e74d11..0000000 --- a/src/pcg.hpp +++ /dev/null @@ -1,207 +0,0 @@ -/* - * Tiny self-contained version of the PCG Random Number Generation for C++ - * put together from pieces of the much larger C/C++ codebase. - * Wenzel Jakob, February 2015 - * - * The PCG random number generator was developed by Melissa O'Neill - * - * - * 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 - */ - -#pragma once - -#define PCG32_DEFAULT_STATE 0x853c49e6748fea9bULL -#define PCG32_DEFAULT_STREAM 0xda3e39cb94b95bdbULL -#define PCG32_MULT 0x5851f42d4c957f2dULL - -#include -#include -#include -#include - -/// PCG32 Pseudorandom number generator -struct Pcg32 { - /// Initialize the pseudorandom number generator with default seed - Pcg32() : state(PCG32_DEFAULT_STATE), inc(PCG32_DEFAULT_STREAM) {} - - /// Initialize the pseudorandom number generator with the \ref seed() function - Pcg32(uint64_t initstate, uint64_t initseq = 1u) { seed(initstate, initseq); } - - /** - * \brief Seed the pseudorandom number generator - * - * Specified in two parts: a state initializer and a sequence selection - * constant (a.k.a. stream id) - */ - void seed(uint64_t initstate, uint64_t initseq = 1) { - state = 0U; - inc = (initseq << 1u) | 1u; - next_u32(); - state += initstate; - next_u32(); - } - - /// Generate a uniformly distributed unsigned 32-bit random number - uint32_t next_u32() { - uint64_t oldstate = state; - state = oldstate * PCG32_MULT + inc; - uint32_t xorshifted = (uint32_t)(((oldstate >> 18u) ^ oldstate) >> 27u); - uint32_t rot = (uint32_t)(oldstate >> 59u); - return (xorshifted >> rot) | (xorshifted << ((~rot + 1u) & 31)); - } - - /// Generate a uniformly distributed number, r, where 0 <= r < bound - uint32_t next_u32(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 + 1u) % 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 = next_u32(); - if (r >= threshold) - return r % bound; - } - } - - /// Generate a single precision floating point value on the interval [0, 1) - float next_f32() { - /* Trick from MTGP: generate an uniformly distributed - single precision number in [1,2) and subtract 1. */ - union { - uint32_t u; - float f; - } x; - x.u = (next_u32() >> 9) | 0x3f800000u; - return x.f - 1.0f; - } - - /** - * \brief Generate a double precision floating point value on the interval [0, 1) - * - * \remark Since the underlying random number generator produces 32 bit output, - * only the first 32 mantissa bits will be filled (however, the resolution is still - * finer than in \ref nextFloat(), which only uses 23 mantissa bits) - */ - double next_f64() { - /* Trick from MTGP: generate an uniformly distributed - double precision number in [1,2) and subtract 1. */ - union { - uint64_t u; - double d; - } x; - x.u = ((uint64_t)next_u32() << 20) | 0x3ff0000000000000ULL; - return x.d - 1.0; - } - - /** - * \brief Multi-step advance function (jump-ahead, jump-back) - * - * The method used here is based on Brown, "Random Number Generation - * with Arbitrary Stride", Transactions of the American Nuclear - * Society (Nov. 1994). The algorithm is very similar to fast - * exponentiation. - */ - void advance(int64_t delta_) { - uint64_t - cur_mult = PCG32_MULT, - cur_plus = inc, - acc_mult = 1u, - acc_plus = 0u; - - /* Even though delta is an unsigned integer, we can pass a signed - integer to go backwards, it just goes "the long way round". */ - uint64_t delta = (uint64_t)delta_; - - while (delta > 0) { - if (delta & 1) { - acc_mult *= cur_mult; - acc_plus = acc_plus * cur_mult + cur_plus; - } - cur_plus = (cur_mult + 1) * cur_plus; - cur_mult *= cur_mult; - delta /= 2; - } - state = acc_mult * state + acc_plus; - } - - /** - * \brief Draw uniformly distributed permutation and permute the - * given STL container - * - * From: Knuth, TAoCP Vol. 2 (3rd 3d), Section 3.4.2 - */ - template - void shuffle(Iterator begin, Iterator end) { - for (Iterator it = end - 1; it > begin; --it) - std::iter_swap(it, begin + next_u32((uint32_t)(it - begin + 1))); - } - - /// Compute the distance between two PCG32 pseudorandom number generators - int64_t operator-(const Pcg32& other) const { - assert(inc == other.inc); - - uint64_t - cur_mult = PCG32_MULT, - cur_plus = inc, - cur_state = other.state, - the_bit = 1u, - distance = 0u; - - while (state != cur_state) { - if ((state & the_bit) != (cur_state & the_bit)) { - cur_state = cur_state * cur_mult + cur_plus; - distance |= the_bit; - } - assert((state & the_bit) == (cur_state & the_bit)); - the_bit <<= 1; - cur_plus = (cur_mult + 1ULL) * cur_plus; - cur_mult *= cur_mult; - } - - return (int64_t)distance; - } - - /// Equality operator - bool operator==(const Pcg32& other) const { return state == other.state && inc == other.inc; } - - /// Inequality operator - bool operator!=(const Pcg32& other) const { return state != other.state || inc != other.inc; } - - uint64_t state; // RNG state. All values are possible. - uint64_t inc; // Controls which RNG sequence (stream) is selected. Must *always* be odd. -}; -- cgit v1.2.3-70-g09d2