Added SplitMix64 to initialize xoshiro256**

lballabio · Aug 17, 2023 · 0d490e9 · 0d490e9
1 parent a4a4af7
commit 0d490e9
Show file tree

Hide file tree

Showing 10 changed files with 134 additions and 47 deletions.
diff --git a/QuantLib.vcxproj b/QuantLib.vcxproj
@@ -1124,6 +1124,7 @@
     <ClInclude Include="ql\math\randomnumbers\seedgenerator.hpp" />
     <ClInclude Include="ql\math\randomnumbers\sobolbrownianbridgersg.hpp" />
     <ClInclude Include="ql\math\randomnumbers\sobolrsg.hpp" />
+    <ClInclude Include="ql\math\randomnumbers\splitmix64.cpp" />
     <ClInclude Include="ql\math\randomnumbers\stochasticcollocationinvcdf.hpp" />
     <ClInclude Include="ql\math\randomnumbers\xoshiro256starstaruniformrng.hpp" />
     <ClInclude Include="ql\math\richardsonextrapolation.hpp" />
@@ -2316,6 +2317,7 @@
     <ClCompile Include="ql\math\randomnumbers\seedgenerator.cpp" />
     <ClCompile Include="ql\math\randomnumbers\sobolbrownianbridgersg.cpp" />
     <ClCompile Include="ql\math\randomnumbers\sobolrsg.cpp" />
+    <ClCompile Include="ql\math\randomnumbers\splitmix64.cpp" />
     <ClCompile Include="ql\math\randomnumbers\stochasticcollocationinvcdf.cpp" />
     <ClCompile Include="ql\math\randomnumbers\xoshiro256starstaruniformrng.cpp" />
     <ClCompile Include="ql\math\richardsonextrapolation.cpp" />

diff --git a/QuantLib.vcxproj.filters b/QuantLib.vcxproj.filters
@@ -1272,6 +1272,9 @@
     <ClInclude Include="ql\math\randomnumbers\sobolrsg.hpp">
       <Filter>math\randomnumbers</Filter>
     </ClInclude>
+    <ClInclude Include="ql\math\randomnumbers\splitmix64.hpp">
+      <Filter>math\randomnumbers</Filter>
+    </ClInclude>
     <ClInclude Include="ql\math\randomnumbers\stochasticcollocationinvcdf.hpp">
       <Filter>math\randomnumbers</Filter>
     </ClInclude>
@@ -4952,6 +4955,9 @@
     <ClCompile Include="ql\math\randomnumbers\sobolrsg.cpp">
       <Filter>math\randomnumbers</Filter>
     </ClCompile>
+    <ClCompile Include="ql\math\randomnumbers\splitmix64.cpp">
+      <Filter>math\randomnumbers</Filter>
+    </ClCompile>
     <ClCompile Include="ql\math\randomnumbers\stochasticcollocationinvcdf.cpp">
       <Filter>math\randomnumbers</Filter>
     </ClCompile>

diff --git a/ql/CMakeLists.txt b/ql/CMakeLists.txt
@@ -428,6 +428,7 @@ set(QL_SOURCES
     math/randomnumbers/seedgenerator.cpp
     math/randomnumbers/sobolbrownianbridgersg.cpp
     math/randomnumbers/sobolrsg.cpp
+    math/randomnumbers/splitmix64.cpp
     math/randomnumbers/stochasticcollocationinvcdf.cpp
     math/randomnumbers/xoshiro256starstaruniformrng.cpp
     math/richardsonextrapolation.cpp
@@ -1530,6 +1531,7 @@ set(QL_HEADERS
     math/randomnumbers/seedgenerator.hpp
     math/randomnumbers/sobolbrownianbridgersg.hpp
     math/randomnumbers/sobolrsg.hpp
+    math/randomnumbers/splitmix64.hpp
     math/randomnumbers/stochasticcollocationinvcdf.hpp
     math/randomnumbers/xoshiro256starstaruniformrng.hpp
     math/richardsonextrapolation.hpp

diff --git a/ql/math/randomnumbers/Makefile.am b/ql/math/randomnumbers/Makefile.am
@@ -23,6 +23,7 @@ this_include_HEADERS = \
 	seedgenerator.hpp \
 	sobolbrownianbridgersg.hpp \
 	sobolrsg.hpp \
+	splitmix64.hpp \
 	stochasticcollocationinvcdf.hpp \
 	xoshiro256starstaruniformrng.hpp
 
@@ -38,6 +39,7 @@ cpp_files = \
 	seedgenerator.cpp \
 	sobolbrownianbridgersg.cpp \
 	sobolrsg.cpp \
+	splitmix64.cpp \
 	stochasticcollocationinvcdf.cpp \
 	xoshiro256starstaruniformrng.cpp
 

diff --git a/ql/math/randomnumbers/all.hpp b/ql/math/randomnumbers/all.hpp
@@ -20,6 +20,7 @@
 #include <ql/math/randomnumbers/seedgenerator.hpp>
 #include <ql/math/randomnumbers/sobolbrownianbridgersg.hpp>
 #include <ql/math/randomnumbers/sobolrsg.hpp>
+#include <ql/math/randomnumbers/splitmix64.hpp>
 #include <ql/math/randomnumbers/stochasticcollocationinvcdf.hpp>
 #include <ql/math/randomnumbers/xoshiro256starstaruniformrng.hpp>
 
diff --git a/ql/math/randomnumbers/splitmix64.cpp b/ql/math/randomnumbers/splitmix64.cpp
@@ -0,0 +1,25 @@
+/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
+
+/*
+ Copyright (C) 2023 Ralf Konrad Eckel
+
+ This file is part of QuantLib, a free-software/open-source library
+ for financial quantitative analysts and developers - http://quantlib.org/
+
+ QuantLib is free software: you can redistribute it and/or modify it
+ under the terms of the QuantLib license.  You should have received a
+ copy of the license along with this program; if not, please email
+ <[email protected]>. The license is also available online at
+ <http://quantlib.org/license.shtml>.
+
+ This program is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+ FOR A PARTICULAR PURPOSE.  See the license for more details.
+*/
+
+#include <ql/math/randomnumbers/seedgenerator.hpp>
+#include <ql/math/randomnumbers/splitmix64.hpp>
+
+namespace QuantLib {
+    SplitMix64::SplitMix64(uint64_t x) : x_(x) {}
+}
diff --git a/ql/math/randomnumbers/splitmix64.hpp b/ql/math/randomnumbers/splitmix64.hpp
@@ -0,0 +1,57 @@
+/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
+
+/*
+ Copyright (C) 2023 Ralf Konrad Eckel
+
+ This file is part of QuantLib, a free-software/open-source library
+ for financial quantitative analysts and developers - http://quantlib.org/
+
+ QuantLib is free software: you can redistribute it and/or modify it
+ under the terms of the QuantLib license.  You should have received a
+ copy of the license along with this program; if not, please email
+ <[email protected]>. The license is also available online at
+ <http://quantlib.org/license.shtml>.
+
+ This program is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+ FOR A PARTICULAR PURPOSE.  See the license for more details.
+*/
+
+// NOTE: The following copyright notice applies to the
+// original C implementation https://prng.di.unimi.it/splitmix64.c
+// that has been used for this class.
+
+/*  Written in 2015 by Sebastiano Vigna ([email protected])
+
+To the extent possible under law, the author has dedicated all copyright
+and related and neighboring rights to this software to the public domain
+worldwide. This software is distributed without any warranty.
+
+See <http://creativecommons.org/publicdomain/zero/1.0/>. */
+
+/*! \file splitmix64.hpp
+    \brief
+*/
+
+#ifndef quantlib_splitmix64_hpp
+#define quantlib_splitmix64_hpp
+
+#include <stdint.h>
+
+namespace QuantLib {
+    class SplitMix64 {
+      public:
+        SplitMix64(uint64_t x);
+        uint64_t next() const {
+            auto z = (x_ += 0x9e3779b97f4a7c15);
+            z = (z ^ (z >> 30)) * 0xbf58476d1ce4e5b9;
+            z = (z ^ (z >> 27)) * 0x94d049bb133111eb;
+            return z ^ (z >> 31);
+        };
+
+      private:
+        mutable uint64_t x_;
+    };
+}
+
+#endif // quantlib_splitmix64_hpp
diff --git a/ql/math/randomnumbers/xoshiro256starstaruniformrng.cpp b/ql/math/randomnumbers/xoshiro256starstaruniformrng.cpp
@@ -18,36 +18,24 @@
 */
 
 #include <ql/math/randomnumbers/seedgenerator.hpp>
+#include <ql/math/randomnumbers/splitmix64.hpp>
 #include <ql/math/randomnumbers/xoshiro256starstaruniformrng.hpp>
 
 namespace QuantLib {
 
-    Xoshiro256StarStarUniformRng::Xoshiro256StarStarUniformRng(uint64_t seed)
-    : Xoshiro256StarStarUniformRng(seed, seed, seed, seed) {}
+    Xoshiro256StarStarUniformRng::Xoshiro256StarStarUniformRng(uint64_t seed) {
+        SplitMix64 splitMix64(seed != 0 ? seed : SeedGenerator::instance().get());
+        do {
+            s0_ = splitMix64.next();
+            s1_ = splitMix64.next();
+            s2_ = splitMix64.next();
+            s3_ = splitMix64.next();
+        } while (s0_ == 0 && s1_ == 0 && s2_ == 0 && s3_ == 0);
+    }
 
     Xoshiro256StarStarUniformRng::Xoshiro256StarStarUniformRng(uint64_t s0,
                                                                uint64_t s1,
                                                                uint64_t s2,
                                                                uint64_t s3)
-    : s0_(s0), s1_(s1), s2_(s2), s3_(s3) {
-        seedInitialization();
-
-        // We call nextInt64() 1,000 times as the first random numbers
-        // might be "close" to the given seeds.
-        // See https://github.com/lballabio/QuantLib/pull/1769#discussion_r1296997133
-        for (auto i = 0; i < 1'000; ++i) {
-            nextInt64();
-        }
-    }
-
-    void Xoshiro256StarStarUniformRng::seedInitialization() const {
-        if (s0_ == 0 && s1_ == 0 && s2_ == 0 && s3_ == 0) {
-            do {
-                s0_ = SeedGenerator::instance().get();
-                s1_ = SeedGenerator::instance().get();
-                s2_ = SeedGenerator::instance().get();
-                s3_ = SeedGenerator::instance().get();
-            } while (s0_ == 0 || s1_ == 0 || s2_ == 0 || s3_ == 0);
-        }
-    }
+    : s0_(s0), s1_(s1), s2_(s2), s3_(s3) {}
 }
diff --git a/ql/math/randomnumbers/xoshiro256starstaruniformrng.hpp b/ql/math/randomnumbers/xoshiro256starstaruniformrng.hpp
@@ -58,20 +58,16 @@ namespace QuantLib {
       public:
         typedef Sample<Real> sample_type;
 
-        /*! If the given seed is 0, a random seed will be chosen based on clock().
-         * The constructor itself calls nextInt64() 1,000 times as the first random numbers
-         * might be "close" to the given seed.
-         * See https://github.com/lballabio/QuantLib/pull/1769#discussion_r1296997133 */
+        /*! If the given seed is 0, a random seed will be chosen based on clock(). */
         explicit Xoshiro256StarStarUniformRng(uint64_t seed = 0);
 
-        /*! If the given seeds are all 0, a random seed will be chosen based on clock().
-         * The constructor itself calls nextInt64() 1,000 times as the first random numbers
-         * might be "close" to the given seeds.
-         * See https://github.com/lballabio/QuantLib/pull/1769#discussion_r1296997133 */
+        /*! Make sure that s0, s1, s2 and s3 are chosen randomly.
+         * Otherwise, the results of the first random numbers might not be well distributed.
+         * Especially s0 = s1 = s2 = s3 = 0 does not work and will always return 0. */
         Xoshiro256StarStarUniformRng(uint64_t s0, uint64_t s1, uint64_t s2, uint64_t s3);
 
         /*! returns a sample with weight 1.0 containing a random number
-            in the (0.0, 1.0) interval  */
+         * in the (0.0, 1.0) interval */
         sample_type next() const { return {nextReal(), 1.0}; }
 
         //! return a random number in the (0.0, 1.0)-interval
@@ -97,8 +93,6 @@ namespace QuantLib {
 
       private:
         static uint64_t rotl(uint64_t x, int32_t k) { return (x << k) | (x >> (64 - k)); }
-        void seedInitialization() const;
-
         mutable uint64_t s0_, s1_, s2_, s3_;
     };
 }

diff --git a/test-suite/xoshiro256starstar.cpp b/test-suite/xoshiro256starstar.cpp
@@ -19,6 +19,7 @@
 
 #include "xoshiro256starstar.hpp"
 #include "utilities.hpp"
+#include <ql/math/randomnumbers/splitmix64.hpp>
 #include <ql/math/randomnumbers/xoshiro256starstaruniformrng.hpp>
 #include <numeric>
 
@@ -153,6 +154,7 @@ extern "C" {
 // clang-format on
 
 using QuantLib::Real;
+using QuantLib::SplitMix64;
 using QuantLib::Xoshiro256StarStarUniformRng;
 
 void Xoshiro256StarStarTest::testMeanAndStdDevOfNextReal() {
@@ -190,29 +192,37 @@ void Xoshiro256StarStarTest::testAgainstReferenceImplementationInC() {
         "Testing Xoshiro256StarStarUniformRng::nextInt64() against reference implementation in C");
 
     // some random initial seed
-    static const auto s0 = 10108360646465513120ULL;
-    static const auto s1 = 4416403493985791904ULL;
-    static const auto s2 = 7597776674045431742ULL;
-    static const auto s3 = 6431387443075032236ULL;
+    static const auto seed = 10108360646465513120ULL;
+    SplitMix64 splitMix64(seed);
+
+    static const auto s0 = splitMix64.next();
+    static const auto s1 = splitMix64.next();
+    static const auto s2 = splitMix64.next();
+    static const auto s3 = splitMix64.next();
 
-    // simulate the warmup in our implementation
-    // by burning the first 1,000 random numbers in the reference implementation
     s[0] = s0;
     s[1] = s1;
     s[2] = s2;
     s[3] = s3;
-    for (auto i = 0; i < 1'000; ++i) {
-        next();
-    }
 
-    auto rng = Xoshiro256StarStarUniformRng(s0, s1, s2, s3);
+    auto rngFromSeed = Xoshiro256StarStarUniformRng(seed);
+    auto rngFroms0s1s2s3 = Xoshiro256StarStarUniformRng(s0, s1, s2, s3);
     for (auto i = 0; i < 1'000; i++) {
         auto nextRefImpl = next();
-        auto nextOurs = rng.nextInt64();
-        if (nextRefImpl != nextOurs) {
+        auto nextFromSeed = rngFromSeed.nextInt64();
+        auto nextFroms0s1s2s3 = rngFroms0s1s2s3.nextInt64();
+        if (nextRefImpl != nextFromSeed) {
             BOOST_FAIL("Test failed at index "
                        << i << " (expected from reference implementation: " << nextRefImpl
-                       << "ULL, ours: " << nextOurs << "ULL)");
+                       << "ULL, from Xoshiro256StarStarUniformRng(" << seed
+                       << "ULL): " << nextFromSeed << "ULL)");
+        }
+        if (nextFroms0s1s2s3 != nextFromSeed) {
+            BOOST_FAIL("Test failed at index " << i << " (from Xoshiro256StarStarUniformRng("
+                                               << seed << "): " << nextFroms0s1s2s3
+                                               << "ULL, from Xoshiro256StarStarUniformRng(" << s0
+                                               << "ULL, " << s1 << "ULL, " << s2 << "ULL, " << s3
+                                               << "ULL): " << nextFromSeed << "ULL)");
         }
     }
 }