xoshiro265** uniform random number generator (#1769)

LICENSE.TXT

@@ -156,7 +156,7 @@ QuantLib is
 
     Copyright (C) 2021 Anubhav Pandey
     Copyright (C) 2021 Magnus Mencke
-    Copyright (C) 2021, 2022 Ralf Konrad Eckel
+    Copyright (C) 2021, 2022, 2023 Ralf Konrad Eckel
     Copyright (C) 2021, 2022, 2023 Skandinaviska Enskilda Banken AB (publ)
 
     Copyright (C) 2022, 2023 Ignacio Anguita

QuantLib.vcxproj

@@ -1130,6 +1130,7 @@
     <ClInclude Include="ql\math\randomnumbers\sobolbrownianbridgersg.hpp" />
     <ClInclude Include="ql\math\randomnumbers\sobolrsg.hpp" />
     <ClInclude Include="ql\math\randomnumbers\stochasticcollocationinvcdf.hpp" />
+    <ClInclude Include="ql\math\randomnumbers\xoshiro256starstaruniformrng.hpp" />
     <ClInclude Include="ql\math\richardsonextrapolation.hpp" />
     <ClInclude Include="ql\math\rounding.hpp" />
     <ClInclude Include="ql\math\sampledcurve.hpp" />
@@ -2323,6 +2324,7 @@
     <ClCompile Include="ql\math\randomnumbers\sobolbrownianbridgersg.cpp" />
     <ClCompile Include="ql\math\randomnumbers\sobolrsg.cpp" />
     <ClCompile Include="ql\math\randomnumbers\stochasticcollocationinvcdf.cpp" />
+    <ClCompile Include="ql\math\randomnumbers\xoshiro256starstaruniformrng.cpp" />
     <ClCompile Include="ql\math\richardsonextrapolation.cpp" />
     <ClCompile Include="ql\math\rounding.cpp" />
     <ClCompile Include="ql\math\sampledcurve.cpp" />

QuantLib.vcxproj.filters

@@ -1293,6 +1293,9 @@
     <ClInclude Include="ql\math\randomnumbers\stochasticcollocationinvcdf.hpp">
       <Filter>math\randomnumbers</Filter>
     </ClInclude>
+    <ClInclude Include="ql\math\randomnumbers\xoshiro256starstaruniformrng.hpp">
+      <Filter>math\randomnumbers</Filter>
+    </ClInclude>
     <ClInclude Include="ql\math\solvers1d\all.hpp">
       <Filter>math\solvers1D</Filter>
     </ClInclude>
@@ -5003,6 +5006,9 @@
     <ClCompile Include="ql\math\randomnumbers\stochasticcollocationinvcdf.cpp">
       <Filter>math\randomnumbers</Filter>
     </ClCompile>
+    <ClCompile Include="ql\math\randomnumbers\xoshiro256starstaruniformrng.cpp">
+      <Filter>math\randomnumbers</Filter>
+    </ClCompile>
     <ClCompile Include="ql\math\optimization\armijo.cpp">
       <Filter>math\optimization</Filter>
     </ClCompile>

ql/CMakeLists.txt

@@ -425,6 +425,7 @@ set(QL_SOURCES
     math/randomnumbers/sobolbrownianbridgersg.cpp
     math/randomnumbers/sobolrsg.cpp
     math/randomnumbers/stochasticcollocationinvcdf.cpp
+    math/randomnumbers/xoshiro256starstaruniformrng.cpp
     math/richardsonextrapolation.cpp
     math/rounding.cpp
     math/sampledcurve.cpp
@@ -1536,6 +1537,7 @@ set(QL_HEADERS
     math/randomnumbers/sobolbrownianbridgersg.hpp
     math/randomnumbers/sobolrsg.hpp
     math/randomnumbers/stochasticcollocationinvcdf.hpp
+    math/randomnumbers/xoshiro256starstaruniformrng.hpp
     math/richardsonextrapolation.hpp
     math/rounding.hpp
     math/sampledcurve.hpp

ql/math/randomnumbers/Makefile.am

@@ -23,7 +23,8 @@ this_include_HEADERS = \
 	seedgenerator.hpp \
 	sobolbrownianbridgersg.hpp \
 	sobolrsg.hpp \
-	stochasticcollocationinvcdf.hpp
+	stochasticcollocationinvcdf.hpp \
+	xoshiro256starstaruniformrng.hpp
 
 cpp_files = \
     faurersg.cpp \
@@ -37,7 +38,8 @@ cpp_files = \
 	seedgenerator.cpp \
 	sobolbrownianbridgersg.cpp \
 	sobolrsg.cpp \
-	stochasticcollocationinvcdf.cpp
+	stochasticcollocationinvcdf.cpp \
+	xoshiro256starstaruniformrng.cpp
 
 if UNITY_BUILD
 

ql/math/randomnumbers/all.hpp

@@ -21,4 +21,5 @@
 #include <ql/math/randomnumbers/sobolbrownianbridgersg.hpp>
 #include <ql/math/randomnumbers/sobolrsg.hpp>
 #include <ql/math/randomnumbers/stochasticcollocationinvcdf.hpp>
+#include <ql/math/randomnumbers/xoshiro256starstaruniformrng.hpp>
 

ql/math/randomnumbers/xoshiro256starstaruniformrng.cpp

@@ -0,0 +1,68 @@
+/* -*- 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/xoshiro256starstaruniformrng.hpp>
+
+namespace QuantLib {
+
+    namespace {
+
+        // 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/>.
+        */
+        class SplitMix64 {
+          public:
+            explicit SplitMix64(std::uint64_t x) : x_(x) {}
+            std::uint64_t next() const {
+                auto z = (x_ += 0x9e3779b97f4a7c15);
+                z = (z ^ (z >> 30)) * 0xbf58476d1ce4e5b9;
+                z = (z ^ (z >> 27)) * 0x94d049bb133111eb;
+                return z ^ (z >> 31);
+            };
+
+          private:
+            mutable std::uint64_t x_;
+        };
+    }
+
+    Xoshiro256StarStarUniformRng::Xoshiro256StarStarUniformRng(std::uint64_t seed) {
+        SplitMix64 splitMix64(seed != 0 ? seed : SeedGenerator::instance().get());
+        s0_ = splitMix64.next();
+        s1_ = splitMix64.next();
+        s2_ = splitMix64.next();
+        s3_ = splitMix64.next();
+    }
+
+    Xoshiro256StarStarUniformRng::Xoshiro256StarStarUniformRng(std::uint64_t s0,
+                                                               std::uint64_t s1,
+                                                               std::uint64_t s2,
+                                                               std::uint64_t s3)
+    : s0_(s0), s1_(s1), s2_(s2), s3_(s3) {}
+
+}

ql/math/randomnumbers/xoshiro256starstaruniformrng.hpp

@@ -0,0 +1,100 @@
+/* -*- 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/xoshiro256starstar.c
+// that has been used for this class.
+
+/*  Written in 2018 by David Blackman and 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 xoshiro256starstaruniformrng.hpp
+    \brief xoshiro256** uniform random number generator
+*/
+
+#ifndef quantlib_xoshiro256starstar_uniform_rng_hpp
+#define quantlib_xoshiro256starstar_uniform_rng_hpp
+
+#include <ql/methods/montecarlo/sample.hpp>
+#include <ql/types.hpp>
+#include <cstdint>
+
+namespace QuantLib {
+
+    //! Uniform random number generator
+    /*! xoshiro256** random number generator of period 2**256-1
+
+        For more details see
+            https://prng.di.unimi.it/
+        and its reference implementation
+            https://prng.di.unimi.it/xoshiro256starstar.c
+
+        \test the correctness of the returned values is tested by checking them
+               against the reference implementation in c.
+    */
+    class Xoshiro256StarStarUniformRng {
+      public:
+        typedef Sample<Real> sample_type;
+
+        /*! If the given seed is 0, a random seed will be chosen based on clock(). */
+        explicit Xoshiro256StarStarUniformRng(std::uint64_t seed = 0);
+
+        /*! 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(std::uint64_t s0, std::uint64_t s1, std::uint64_t s2, std::uint64_t s3);
+
+        /*! returns a sample with weight 1.0 containing a random number
+         * 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
+        Real nextReal() const { return (Real(nextInt64() >> 11) + 0.5) * (1.0 / Real(1ULL << 53)); }
+
+        //! return a random integer in the [0,0xffffffffffffffffULL]-interval
+        std::uint64_t nextInt64() const {
+            const auto result = rotl(s1_ * 5, 7) * 9;
+
+            const auto t = s1_ << 17;
+
+            s2_ ^= s0_;
+            s3_ ^= s1_;
+            s1_ ^= s2_;
+            s0_ ^= s3_;
+
+            s2_ ^= t;
+
+            s3_ = rotl(s3_, 45);
+
+            return result;
+        }
+
+      private:
+        static std::uint64_t rotl(std::uint64_t x, std::int32_t k) { return (x << k) | (x >> (64 - k)); }
+        mutable std::uint64_t s0_, s1_, s2_, s3_;
+    };
+
+}
+
+#endif

test-suite/CMakeLists.txt

@@ -166,6 +166,7 @@ set(QL_TEST_SOURCES
     varianceswaps.cpp
     volatilitymodels.cpp
     vpp.cpp
+    xoshiro256starstar.cpp
     zabr.cpp
     zerocouponswap.cpp
 )
@@ -340,6 +341,7 @@ set(QL_TEST_HEADERS
     varianceswaps.hpp
     volatilitymodels.hpp
     vpp.hpp
+    xoshiro256starstar.hpp
     zabr.hpp
     zerocouponswap.hpp
 )

test-suite/Makefile.am

@@ -167,6 +167,7 @@ QL_TEST_SRCS = \
 	varianceswaps.cpp \
 	volatilitymodels.cpp \
 	vpp.cpp \
+	xoshiro256starstar.cpp \
 	zabr.cpp \
 	zerocouponswap.cpp
 
@@ -339,6 +340,7 @@ QL_TEST_HDRS = \
 	varianceswaps.hpp \
 	volatilitymodels.hpp \
 	vpp.hpp \
+	xoshiro256starstar.hpp \
 	zabr.hpp \
 	zerocouponswap.hpp