123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117 |
- /*
- * 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(initial_state, initseq)
- // pcg32_srandom_r(rng, initial_state, 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 initial_state, uint64_t initseq)
- {
- rng->state = 0U;
- rng->inc = (initseq << 1u) | 1u;
- pcg32_random_r(rng);
- rng->state += initial_state;
- 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)
- {
- uint64_t oldstate = rng->state;
- rng->state = oldstate * 6364136223846793005ULL + rng->inc;
- uint32_t xorshifted = (uint32_t)(((oldstate >> 18u) ^ oldstate) >> 27u);
- uint32_t rot = (uint32_t)(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);
- }
|