rand/pcg: add pcg_xsh_rr generator

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

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 <danny@nerdcruft.net>
+ */
+
+#include "pcg.hpp"
+
+template class cruft::rand::pcg_xsh_rr<u64,u32>;
+

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 <danny@nerdcruft.net>
+ */
+
+#pragma once
+
+#include "../bitwise.hpp"
+#include "../std.hpp"
+
+#include <limits>
+#include <random>
+
+
+// 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<result_type>::max (); }
+        static constexpr auto min () noexcept { return std::numeric_limits<result_type>::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;
+    };
+}

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<cruft::rand::lcg_t> (void) { return "lcg_
 
 template <> std::string type_to_string<cruft::rand::mwc64x> (void) { return "mwc64x"; }
 
+template <> std::string type_to_string<cruft::rand::pcg_xsh_rr<u64,u32>> (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<cruft::rand::xorshift<uint64_t>> (tap, 0x1234u);
     test_buckets<cruft::rand::lcg_t> (tap, 0x1234u);
     test_buckets<cruft::rand::mwc64x> (tap, 0x1234u);
+    test_buckets<cruft::rand::pcg_xsh_rr<>> (tap, 0x1234u);
 
     return tap.status ();
 }