Add 512 bit xoshiro family

doc/generator_defs.qbk

@@ -36,6 +36,9 @@
 [template ranlux64_4_01_speed[] 2%]
 [template mt19937ar_c_speed[] 97%]
 [template splitmix64_speed[] 155%]
-[template xoshiro256pp_speed[] 90%]
-[template xoshiro256d_speed[] 101%]
-[template xoshiro256mm_speed[] 82%]
+[template xoshiro256pp_speed[] 132%]
+[template xoshiro256d_speed[] 139%]
+[template xoshiro256mm_speed[] 116%]
+[template xoshiro512pp_speed[] 133%]
+[template xoshiro512d_speed[] 144%]
+[template xoshiro512mm_speed[] 113%]

doc/generators.qbk

@@ -87,6 +87,9 @@ numbers mean faster random number generation.
   [[__xoshiro256pp] [2[sup 256] [`4*sizeof(uint64_t)`] [[xoshiro256pp_speed]] [[xoshiro256++ from https://prng.di.unimi.it]]
   [[__xoshiro256d] [2[sup 256] [`4*sizeof(uint64_t)`] [[xoshiro256d_speed]] [[This generator returns doubles instead of uint64_t. It is modified xoshiro256+ from https://prng.di.unimi.it]]
   [[__xoshiro256mm] [2[sup 256] [`4*sizeof(uint64_t)`] [[xoshiro256mm_speed]] [[xoshiro256** from https://prng.di.unimi.it]]
+  [[__xoshiro512pp] [2[sup 512] [`8*sizeof(uint64_t)`] [[xoshiro512pp_speed]] [[xoshiro512++ from https://prng.di.unimi.it]]
+  [[__xoshiro512d] [2[sup 512] [`8*sizeof(uint64_t)`] [[xoshiro512d_speed]] [[This generator returns doubles instead of uint64_t. It is modified xoshiro512+ from https://prng.di.unimi.it]]
+  [[__xoshiro512mm] [2[sup 512] [`8*sizeof(uint64_t)`] [[xoshiro512mm_speed]] [[xoshiro512** from https://prng.di.unimi.it]]
 ]
 
 As observable from the table, there is generally a quality/performance/memory
@@ -98,8 +101,8 @@ generators for a given application of Monte Carlo simulation will improve
 the confidence in the results.
 
 If the names of the generators don't ring any bell and you have no idea
-which generator to use, it is reasonable to employ __mt19937 for a start: It
-is fast and has acceptable quality.
+which generator to use, it is reasonable to employ __xoshiro256pp for a start: It
+is fast and of high quality.
 
 [note These random number generators are not intended for use in applications
 where non-deterministic random numbers are required. See __random_device

include/boost/random/detail/xoshiro_base.hpp

@@ -27,7 +27,6 @@
 #include <istream>
 #include <type_traits>
 #include <iterator>
-#include <iostream>
 
 namespace boost {
 namespace random {
@@ -79,6 +78,34 @@ class xoshiro_base
         state_[3] = s3;
     }
 
+    inline void jump_impl(const std::integral_constant<std::size_t, 8>&) noexcept
+    {
+        constexpr std::array<std::uint64_t, 8U> jump = {{ UINT64_C(0x33ed89b6e7a353f9), UINT64_C(0x760083d7955323be),
+                                                          UINT64_C(0x2837f2fbb5f22fae), UINT64_C(0x4b8c5674d309511c),
+                                                          UINT64_C(0xb11ac47a7ba28c25), UINT64_C(0xf1be7667092bcc1c),
+                                                          UINT64_C(0x53851efdb6df0aaf), UINT64_C(0x1ebbc8b23eaf25db) }};
+
+        std::array<std::uint64_t, 8U> t = {{ 0, 0, 0, 0, 0, 0, 0, 0 }};
+
+        for (std::size_t i = 0; i < jump.size(); ++i)
+        {
+            for (std::size_t b = 0; b < 64U; ++b)
+            {
+                if (jump[i] & UINT64_C(1) << b)
+                {
+                    for (std::size_t w = 0; w < state_.size(); ++w)
+                    {
+                        t[w] ^= state_[w];
+                    }
+                }
+
+                next();
+            }
+        }
+
+        state_ = t;
+    }
+
     inline void long_jump_impl(const std::integral_constant<std::size_t, 4>&) noexcept
     {
         constexpr std::array<std::uint64_t, 4> long_jump = {{ UINT64_C(0x76e15d3efefdcbbf), UINT64_C(0xc5004e441c522fb3),
@@ -110,6 +137,34 @@ class xoshiro_base
         state_[2] = s2;
         state_[3] = s3;
     }
+
+    inline void long_jump_impl(const std::integral_constant<std::size_t, 8>&) noexcept
+    {
+        constexpr std::array<std::uint64_t, 8U> long_jump = {{ UINT64_C(0x11467fef8f921d28), UINT64_C(0xa2a819f2e79c8ea8),
+                                                               UINT64_C(0xa8299fc284b3959a), UINT64_C(0xb4d347340ca63ee1),
+                                                               UINT64_C(0x1cb0940bedbff6ce), UINT64_C(0xd956c5c4fa1f8e17),
+                                                               UINT64_C(0x915e38fd4eda93bc), UINT64_C(0x5b3ccdfa5d7daca5) }};
+
+        std::array<std::uint64_t, 8U> t = {{ 0, 0, 0, 0, 0, 0, 0, 0 }};
+
+        for (std::size_t i = 0; i < long_jump.size(); ++i)
+        {
+            for (std::size_t b = 0; b < 64U; ++b)
+            {
+                if (long_jump[i] & UINT64_C(1) << b)
+                {
+                    for (std::size_t w = 0; w < state_.size(); ++w)
+                    {
+                        t[w] ^= state_[w];
+                    }
+                }
+
+                next();
+            }
+        }
+
+        state_ = t;
+    }
 
 public:
 

include/boost/random/xoshiro.hpp

@@ -93,14 +93,32 @@ class xoshiro256d final : public detail::xoshiro_base<xoshiro256d, 4, double>
         return result;
     }
 
-    inline double next() noexcept
+    inline result_type next() noexcept
     {
         #if (__cplusplus >= 201703L || _MSVC_LANG >= 201703L) && defined(__cpp_hex_float) && __cpp_hex_float >= 201603L
         return static_cast<double>((next_int() >> 11)) * 0x1.0p-53;
         #else
         return static_cast<double>((next_int() >> 11)) * 1.11022302462515654e-16;
         #endif
     }
+
+    static constexpr result_type (min)() noexcept
+    {
+        #if (__cplusplus >= 201703L || _MSVC_LANG >= 201703L) && defined(__cpp_hex_float) && __cpp_hex_float >= 201603L
+        return static_cast<double>((std::numeric_limits<std::uint64_t>::min)()) * 0x1.0p-53;
+        #else
+        return static_cast<double>((std::numeric_limits<std::uint64_t>::min)()) * 1.11022302462515654e-16;
+        #endif
+    }
+
+    static constexpr result_type (max)() noexcept
+    {
+        #if (__cplusplus >= 201703L || _MSVC_LANG >= 201703L) && defined(__cpp_hex_float) && __cpp_hex_float >= 201603L
+        return static_cast<double>((std::numeric_limits<std::uint64_t>::max)()) * 0x1.0p-53;
+        #else
+        return static_cast<double>((std::numeric_limits<std::uint64_t>::max)()) * 1.11022302462515654e-16;
+        #endif
+    }
 };
 
 /**
@@ -140,6 +158,173 @@ class xoshiro256mm final : public detail::xoshiro_base<xoshiro256mm, 4>
     }
 };
 
+/**
+ * This is xoshiro512++ 1.0, one of our all-purpose, rock-solid
+ * generators. It has excellent (about 1ns) speed, a state (512 bits) that
+ * is large enough for any parallel application, and it passes all tests
+ * we are aware of.
+ *
+ * For generating just floating-point numbers, xoshiro512+ is even faster.
+ *
+ * The state must be seeded so that it is not everywhere zero. If you have
+ * a 64-bit seed, we suggest to seed a splitmix64 generator and use its
+ * output to fill s.
+ */
+
+class xoshiro512pp final : public detail::xoshiro_base<xoshiro512pp, 8>
+{
+private:
+
+    using Base = detail::xoshiro_base<xoshiro512pp, 8>;
+
+public:
+
+    using Base::Base;
+
+    inline result_type next() noexcept
+    {
+        const std::uint64_t result = boost::core::rotl(state_[0] + state_[2], 17) + state_[2];
+
+        const std::uint64_t t = state_[1] << 11;
+
+        state_[2] ^= state_[0];
+        state_[5] ^= state_[1];
+        state_[1] ^= state_[2];
+        state_[7] ^= state_[3];
+        state_[3] ^= state_[4];
+        state_[4] ^= state_[5];
+        state_[0] ^= state_[6];
+        state_[6] ^= state_[7];
+
+        state_[6] ^= t;
+
+        state_[7] = boost::core::rotl(state_[7], 21);
+
+        return result;
+    }
+};
+
+/* This is xoshiro512** 1.0, one of our all-purpose, rock-solid generators
+ * with increased state size. It has excellent (about 1ns) speed, a state
+ * (512 bits) that is large enough for any parallel application, and it
+ * passes all tests we are aware of.
+ *
+ * For generating just floating-point numbers, xoshiro512+ is even faster.
+ *
+ * The state must be seeded so that it is not everywhere zero. If you have
+ * a 64-bit seed, we suggest to seed a splitmix64 generator and use its
+ *  output to fill s.
+ */
+
+class xoshiro512mm final : public detail::xoshiro_base<xoshiro512mm, 8>
+{
+private:
+
+    using Base = detail::xoshiro_base<xoshiro512mm, 8>;
+
+public:
+
+    using Base::Base;
+
+    inline result_type next() noexcept
+    {
+        const std::uint64_t result = boost::core::rotl(state_[1] * 5, 7) * 9;
+
+        const std::uint64_t t = state_[1] << 11;
+
+        state_[2] ^= state_[0];
+        state_[5] ^= state_[1];
+        state_[1] ^= state_[2];
+        state_[7] ^= state_[3];
+        state_[3] ^= state_[4];
+        state_[4] ^= state_[5];
+        state_[0] ^= state_[6];
+        state_[6] ^= state_[7];
+
+        state_[6] ^= t;
+
+        state_[7] = boost::core::rotl(state_[7], 21);
+
+        return result;
+    }
+};
+
+/* This is xoshiro512+ 1.0, our generator for floating-point numbers with
+ * increased state size. We suggest to use its upper bits for
+ * floating-point generation, as it is slightly faster than xoshiro512**.
+ * It passes all tests we are aware of except for the lowest three bits,
+ * which might fail linearity tests (and just those), so if low linear
+ * complexity is not considered an issue (as it is usually the case) it
+ * can be used to generate 64-bit outputs, too.
+ *
+ * We suggest to use a sign test to extract a random Boolean value, and
+ * right shifts to extract subsets of bits.
+ *
+ * The state must be seeded so that it is not everywhere zero. If you have
+ * a 64-bit seed, we suggest to seed a splitmix64 generator and use its
+ *  output to fill s.
+ */
+
+class xoshiro512d final : public detail::xoshiro_base<xoshiro512d, 8, double>
+{
+private:
+
+    using Base = detail::xoshiro_base<xoshiro512d, 8, double>;
+
+public:
+
+    using Base::Base;
+
+    inline std::uint64_t next_int() noexcept
+    {
+        const std::uint64_t result = state_[0] + state_[2];
+
+        const std::uint64_t t = state_[1] << 11;
+
+        state_[2] ^= state_[0];
+        state_[5] ^= state_[1];
+        state_[1] ^= state_[2];
+        state_[7] ^= state_[3];
+        state_[3] ^= state_[4];
+        state_[4] ^= state_[5];
+        state_[0] ^= state_[6];
+        state_[6] ^= state_[7];
+
+        state_[6] ^= t;
+
+        state_[7] = boost::core::rotl(state_[7], 21);
+
+        return result;
+    }
+
+    inline result_type next() noexcept
+    {
+        #if (__cplusplus >= 201703L || _MSVC_LANG >= 201703L) && defined(__cpp_hex_float) && __cpp_hex_float >= 201603L
+        return static_cast<double>((next_int() >> 11)) * 0x1.0p-53;
+        #else
+        return static_cast<double>((next_int() >> 11)) * 1.11022302462515654e-16;
+        #endif
+    }
+
+    static constexpr result_type (min)() noexcept
+    {
+        #if (__cplusplus >= 201703L || _MSVC_LANG >= 201703L) && defined(__cpp_hex_float) && __cpp_hex_float >= 201603L
+        return static_cast<double>((std::numeric_limits<std::uint64_t>::min)()) * 0x1.0p-53;
+        #else
+        return static_cast<double>((std::numeric_limits<std::uint64_t>::min)()) * 1.11022302462515654e-16;
+        #endif
+    }
+
+    static constexpr result_type (max)() noexcept
+    {
+        #if (__cplusplus >= 201703L || _MSVC_LANG >= 201703L) && defined(__cpp_hex_float) && __cpp_hex_float >= 201603L
+        return static_cast<double>((std::numeric_limits<std::uint64_t>::max)()) * 0x1.0p-53;
+        #else
+        return static_cast<double>((std::numeric_limits<std::uint64_t>::max)()) * 1.11022302462515654e-16;
+        #endif
+    }
+};
+
 } // namespace random
 } // namespace boost
 

performance/random_speed.cpp

@@ -385,6 +385,9 @@ int main(int argc, char*argv[])
   run(iter, "xoshiro256pp", boost::random::xoshiro256pp());
   run(iter, "xoshiro256d", boost::random::xoshiro256d());
   run(iter, "xoshiro256mm", boost::random::xoshiro256mm());
+  run(iter, "xoshiro512pp", boost::random::xoshiro512pp());
+  run(iter, "xoshiro512d", boost::random::xoshiro512d());
+  run(iter, "xoshiro512mm", boost::random::xoshiro512mm());
 
 #ifdef HAVE_MT19937INT_C
   // requires the original mt19937int.c

test/Jamfile.v2

@@ -71,6 +71,12 @@ run test_xoshiro256d.cpp /boost/test//boost_unit_test_framework ;
 run test_comp_xoshiro256d.cpp ;
 run test_xoshiro256mm.cpp /boost/test//boost_unit_test_framework ;
 run test_comp_xoshiro256mm.cpp ;
+run test_xoshiro512pp.cpp /boost/test//boost_unit_test_framework ;
+run test_comp_xoshiro512pp.cpp ;
+run test_xoshiro512mm.cpp /boost/test//boost_unit_test_framework ;
+run test_comp_xoshiro512mm.cpp ;
+run test_xoshiro512d.cpp /boost/test//boost_unit_test_framework ;
+run test_comp_xoshiro512d.cpp ;
 
 run niederreiter_base2_validate.cpp /boost/test//boost_unit_test_framework ;
 run sobol_validate.cpp /boost/test//boost_unit_test_framework ;