diff --git a/QuantLib.vcxproj b/QuantLib.vcxproj
index 9c207661ab7..864c34a4530 100644
--- a/QuantLib.vcxproj
+++ b/QuantLib.vcxproj
@@ -1124,6 +1124,7 @@
+
@@ -2316,6 +2317,7 @@
+
diff --git a/QuantLib.vcxproj.filters b/QuantLib.vcxproj.filters
index 2b50e5a0965..0f27ab22fda 100644
--- a/QuantLib.vcxproj.filters
+++ b/QuantLib.vcxproj.filters
@@ -1272,6 +1272,9 @@
math\randomnumbers
+
+ math\randomnumbers
+
math\randomnumbers
@@ -4952,6 +4955,9 @@
math\randomnumbers
+
+ math\randomnumbers
+
math\randomnumbers
diff --git a/ql/CMakeLists.txt b/ql/CMakeLists.txt
index 516484d68b6..73e08391731 100644
--- 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
index 21faedb2ad0..9ab47593043 100644
--- 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
index adcaafb4a07..427d8675bf8 100644
--- a/ql/math/randomnumbers/all.hpp
+++ b/ql/math/randomnumbers/all.hpp
@@ -20,6 +20,7 @@
#include
#include
#include
+#include
#include
#include
diff --git a/ql/math/randomnumbers/splitmix64.cpp b/ql/math/randomnumbers/splitmix64.cpp
new file mode 100644
index 00000000000..2b7d70edc94
--- /dev/null
+++ 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
+ . The license is also available online at
+ .
+
+ 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
+#include
+
+namespace QuantLib {
+ SplitMix64::SplitMix64(uint64_t x) : x_(x) {}
+}
diff --git a/ql/math/randomnumbers/splitmix64.hpp b/ql/math/randomnumbers/splitmix64.hpp
new file mode 100644
index 00000000000..8dacd59cfb4
--- /dev/null
+++ 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
+ . The license is also available online at
+ .
+
+ 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 (vigna@acm.org)
+
+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 . */
+
+/*! \file splitmix64.hpp
+ \brief
+*/
+
+#ifndef quantlib_splitmix64_hpp
+#define quantlib_splitmix64_hpp
+
+#include
+
+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
index 804763f68a7..036985a10ba 100644
--- a/ql/math/randomnumbers/xoshiro256starstaruniformrng.cpp
+++ b/ql/math/randomnumbers/xoshiro256starstaruniformrng.cpp
@@ -18,36 +18,24 @@
*/
#include
+#include
#include
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
index 623b1a87bd2..6bc6f32056b 100644
--- a/ql/math/randomnumbers/xoshiro256starstaruniformrng.hpp
+++ b/ql/math/randomnumbers/xoshiro256starstaruniformrng.hpp
@@ -58,20 +58,16 @@ namespace QuantLib {
public:
typedef Sample 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
index f1c698dea35..e69301dacf0 100644
--- a/test-suite/xoshiro256starstar.cpp
+++ b/test-suite/xoshiro256starstar.cpp
@@ -19,6 +19,7 @@
#include "xoshiro256starstar.hpp"
#include "utilities.hpp"
+#include
#include
#include
@@ -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)");
}
}
}