diff --git a/CMakeLists.txt b/CMakeLists.txt index 7b27b5f4..e4e0ccc4 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -395,6 +395,8 @@ list ( rand/xorshift.hpp rand/mwc64x.cpp rand/mwc64x.hpp + rand/pcg.cpp + rand/pcg.hpp random.cpp random.hpp range.cpp diff --git a/rand/pcg.cpp b/rand/pcg.cpp new file mode 100644 index 00000000..865d7782 --- /dev/null +++ b/rand/pcg.cpp @@ -0,0 +1,12 @@ +/* + * This Source Code Form is subject to the terms of the Mozilla Public + * License, v. 2.0. If a copy of the MPL was not distributed with this + * file, You can obtain one at http://mozilla.org/MPL/2.0/. + * + * Copyright 2019 Danny Robson + */ + +#include "pcg.hpp" + +template class cruft::rand::pcg_xsh_rr; + diff --git a/rand/pcg.hpp b/rand/pcg.hpp new file mode 100644 index 00000000..76c3fbc9 --- /dev/null +++ b/rand/pcg.hpp @@ -0,0 +1,86 @@ +/* + * This Source Code Form is subject to the terms of the Mozilla Public + * License, v. 2.0. If a copy of the MPL was not distributed with this + * file, You can obtain one at http://mozilla.org/MPL/2.0/. + * + * Copyright 2019 Danny Robson + */ + +#pragma once + +#include "../bitwise.hpp" +#include "../std.hpp" + +#include +#include + + +// Implementations of algorithms from the paper: "PCG: A Family of Simple +// Fast Space-Efficient Statistically Good Algorithms for Random Number +// Generation" +namespace cruft::rand { + // Implements the PCG-XSH-RR algorithm; xorshift-high, random rotation. + template < + typename StateT = u64, + typename OutputT = u32 + > + class pcg_xsh_rr { + public: + using result_type = OutputT; + + explicit pcg_xsh_rr (std::seed_seq seed) + : state (increment) + { + // Overlap a series of unsigned words with state words. We might + // over-provision but that shouldn't be an issue with the sizes we + // support. + union { + unsigned u[sizeof (StateT) / sizeof (unsigned) + 1]; + StateT s; + } words; + + seed.generate (std::begin (words.u), std::end (words.u)); + state += words.s; + } + + explicit pcg_xsh_rr (StateT seed) noexcept + : state (increment + seed) + { (*this)(); } + + + static constexpr auto ShiftV = 5u; + + static constexpr auto state_bits = sizeof (StateT ) * 8u; + static constexpr auto state_shift = sizeof (StateT ) * 8u - ShiftV; + static constexpr auto output_shift = sizeof (OutputT) * 8u - ShiftV; + static constexpr auto xshift = (state_bits - output_shift) / 2u; + + + result_type operator() (void) noexcept + { + StateT x = state; + StateT const rrotate = x >> state_shift; + + state = x * multiplier + increment; + x ^= x >> xshift; + + // Be very careful to cast to result_type before doing the + // rotation otherwise we'll rotate the zeroed upper bits into the + // returned value. + return rotater (result_type (x >> output_shift), rrotate); + } + + + static constexpr auto max () noexcept { return std::numeric_limits::max (); } + static constexpr auto min () noexcept { return std::numeric_limits::min (); } + + private: + // Initialiser from the PCG example code. + StateT state = 0x4d595df4d0f33173; + + // Using multiplier and increment from the PCG paper, which appear to + // also correspond to constants chosen by Knuth. + static constexpr StateT const multiplier = 6364136223846793005u; + static constexpr StateT const increment = 1442695040888963407u; + }; +} diff --git a/test/rand/buckets.cpp b/test/rand/buckets.cpp index 7ba6a34b..6275f46f 100644 --- a/test/rand/buckets.cpp +++ b/test/rand/buckets.cpp @@ -1,6 +1,7 @@ #include "rand/xorshift.hpp" #include "rand/lcg.hpp" #include "rand/mwc64x.hpp" +#include "rand/pcg.hpp" #include "tap.hpp" #include "maths.hpp" @@ -20,6 +21,8 @@ template <> std::string type_to_string (void) { return "lcg_ template <> std::string type_to_string (void) { return "mwc64x"; } +template <> std::string type_to_string> (void) { return "pcg_xsh_rr<64,32>"; } + /////////////////////////////////////////////////////////////////////////////// /// check random outputs are roughly divisible between a number of fixed width @@ -62,6 +65,7 @@ main (int,char**) test_buckets> (tap, 0x1234u); test_buckets (tap, 0x1234u); test_buckets (tap, 0x1234u); + test_buckets> (tap, 0x1234u); return tap.status (); }