Coverage Report

Created: 2026-03-31 07:01

next uncovered line (L), next uncovered region (R), next uncovered branch (B)
/src/quantlib/ql/math/randomnumbers/burley2020sobolrsg.cpp
Line
Count
Source
1
/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3
/*
4
 Copyright (C) 2023 Peter Caspers
5
6
 This file is part of QuantLib, a free-software/open-source library
7
 for financial quantitative analysts and developers - http://quantlib.org/
8
9
 QuantLib is free software: you can redistribute it and/or modify it
10
 under the terms of the QuantLib license.  You should have received a
11
 copy of the license along with this program; if not, please email
12
 <quantlib-dev@lists.sf.net>. The license is also available online at
13
 <https://www.quantlib.org/license.shtml>.
14
15
 This program is distributed in the hope that it will be useful, but WITHOUT
16
 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17
 FOR A PARTICULAR PURPOSE.  See the license for more details.
18
*/
19
20
#include <ql/math/randomnumbers/burley2020sobolrsg.hpp>
21
#include <ql/math/randomnumbers/mt19937uniformrng.hpp>
22
#include <ql/errors.hpp>
23
24
namespace QuantLib {
25
26
    Burley2020SobolRsg::Burley2020SobolRsg(Size dimensionality,
27
                                           unsigned long seed,
28
                                           SobolRsg::DirectionIntegers directionIntegers,
29
                                           unsigned long scrambleSeed)
30
0
    : dimensionality_(dimensionality), seed_(seed), directionIntegers_(directionIntegers),
31
0
      integerSequence_(dimensionality), sequence_(std::vector<Real>(dimensionality), 1.0) {
32
0
        reset();
33
0
        group4Seeds_.resize((dimensionality_ - 1) / 4 + 1);
34
0
        MersenneTwisterUniformRng mt(scrambleSeed);
35
0
        for (auto& s : group4Seeds_) {
36
0
            s = static_cast<std::uint32_t>(mt.nextInt32());
37
0
        }
38
0
    }
39
40
0
    void Burley2020SobolRsg::reset() const {
41
0
        sobolRsg_ = ext::make_shared<SobolRsg>(dimensionality_, seed_, directionIntegers_, false);
42
0
        nextSequenceCounter_ = 0;
43
0
    }
44
45
0
    const std::vector<std::uint32_t>& Burley2020SobolRsg::skipTo(std::uint32_t n) const {
46
0
        nextSequenceCounter_ = n;
47
0
        nextInt32Sequence();
48
0
        --nextSequenceCounter_;
49
0
        return integerSequence_;
50
0
    }
51
52
    namespace {
53
54
        // for reverseBits() see http://graphics.stanford.edu/~seander/bithacks.html#BitReverseTable
55
56
        const std::uint8_t bitReverseTable[] = {
57
            0U,   128U, 64U,  192U, 32U,  160U, 96U,  224U, 16U,  144U, 80U,  208U, 48U,  176U,
58
            112U, 240U, 8U,   136U, 72U,  200U, 40U,  168U, 104U, 232U, 24U,  152U, 88U,  216U,
59
            56U,  184U, 120U, 248U, 4U,   132U, 68U,  196U, 36U,  164U, 100U, 228U, 20U,  148U,
60
            84U,  212U, 52U,  180U, 116U, 244U, 12U,  140U, 76U,  204U, 44U,  172U, 108U, 236U,
61
            28U,  156U, 92U,  220U, 60U,  188U, 124U, 252U, 2U,   130U, 66U,  194U, 34U,  162U,
62
            98U,  226U, 18U,  146U, 82U,  210U, 50U,  178U, 114U, 242U, 10U,  138U, 74U,  202U,
63
            42U,  170U, 106U, 234U, 26U,  154U, 90U,  218U, 58U,  186U, 122U, 250U, 6U,   134U,
64
            70U,  198U, 38U,  166U, 102U, 230U, 22U,  150U, 86U,  214U, 54U,  182U, 118U, 246U,
65
            14U,  142U, 78U,  206U, 46U,  174U, 110U, 238U, 30U,  158U, 94U,  222U, 62U,  190U,
66
            126U, 254U, 1U,   129U, 65U,  193U, 33U,  161U, 97U,  225U, 17U,  145U, 81U,  209U,
67
            49U,  177U, 113U, 241U, 9U,   137U, 73U,  201U, 41U,  169U, 105U, 233U, 25U,  153U,
68
            89U,  217U, 57U,  185U, 121U, 249U, 5U,   133U, 69U,  197U, 37U,  165U, 101U, 229U,
69
            21U,  149U, 85U,  213U, 53U,  181U, 117U, 245U, 13U,  141U, 77U,  205U, 45U,  173U,
70
            109U, 237U, 29U,  157U, 93U,  221U, 61U,  189U, 125U, 253U, 3U,   131U, 67U,  195U,
71
            35U,  163U, 99U,  227U, 19U,  147U, 83U,  211U, 51U,  179U, 115U, 243U, 11U,  139U,
72
            75U,  203U, 43U,  171U, 107U, 235U, 27U,  155U, 91U,  219U, 59U,  187U, 123U, 251U,
73
            7U,   135U, 71U,  199U, 39U,  167U, 103U, 231U, 23U,  151U, 87U,  215U, 55U,  183U,
74
            119U, 247U, 15U,  143U, 79U,  207U, 47U,  175U, 111U, 239U, 31U,  159U, 95U,  223U,
75
            63U,  191U, 127U, 255U};
76
77
0
        inline std::uint32_t reverseBits(std::uint32_t x) {
78
0
            return (bitReverseTable[x & 0xff] << 24) | (bitReverseTable[(x >> 8) & 0xff] << 16) |
79
0
                   (bitReverseTable[(x >> 16) & 0xff] << 8) | (bitReverseTable[(x >> 24) & 0xff]);
80
0
        }
81
82
0
        inline std::uint32_t laine_karras_permutation(std::uint32_t x, std::uint32_t seed) {
83
0
            x += seed;
84
0
            x ^= x * 0x6c50b47cU;
85
0
            x ^= x * 0xb82f1e52U;
86
0
            x ^= x * 0xc7afe638U;
87
0
            x ^= x * 0x8d22f6e6U;
88
0
            return x;
89
0
        }
90
91
0
        inline std::uint32_t nested_uniform_scramble(std::uint32_t x, std::uint32_t seed) {
92
0
            x = reverseBits(x);
93
0
            x = laine_karras_permutation(x, seed);
94
0
            x = reverseBits(x);
95
0
            return x;
96
0
        }
97
98
        // the results depend a lot on the details of the hash_combine() function that is used
99
        // we use hash_combine() calling hash(), hash_mix() as implemented here:
100
        // https://github.com/boostorg/container_hash/blob/boost-1.83.0/include/boost/container_hash/hash.hpp#L560
101
        // https://github.com/boostorg/container_hash/blob/boost-1.83.0/include/boost/container_hash/hash.hpp#L115
102
        // https://github.com/boostorg/container_hash/blob/boost-1.83.0/include/boost/container_hash/detail/hash_mix.hpp#L67
103
104
0
        inline std::uint64_t local_hash_mix(std::uint64_t x) {
105
0
            const std::uint64_t m = 0xe9846af9b1a615d;
106
0
            x ^= x >> 32;
107
0
            x *= m;
108
0
            x ^= x >> 32;
109
0
            x *= m;
110
0
            x ^= x >> 28;
111
0
            return x;
112
0
        }
113
114
0
        inline std::uint64_t local_hash(const std::uint64_t v) {
115
0
            std::uint64_t seed = 0;
116
0
            seed = (v >> 32) + local_hash_mix(seed);
117
0
            seed = (v & 0xFFFFFFFF) + local_hash_mix(seed);
118
0
            return seed;
119
0
        }
120
121
0
        inline std::uint64_t local_hash_combine(std::uint64_t x, const uint64_t v) {
122
0
            return local_hash_mix(x + 0x9e3779b9 + local_hash(v));
123
0
        }
124
    }
125
126
0
    const std::vector<std::uint32_t>& Burley2020SobolRsg::nextInt32Sequence() const {
127
0
        auto n = nested_uniform_scramble(nextSequenceCounter_, group4Seeds_[0]);
128
0
        const auto& seq = sobolRsg_->skipTo(n);
129
0
        std::copy(seq.begin(), seq.end(), integerSequence_.begin());
130
0
        Size i = 0, group = 0;
131
0
        do {
132
0
            std::uint64_t seed = group4Seeds_[group++];
133
0
            for (Size g = 0; g < 4 && i < dimensionality_; ++g, ++i) {
134
0
                seed = local_hash_combine(seed, g);
135
0
                integerSequence_[i] =
136
0
                    nested_uniform_scramble(integerSequence_[i], static_cast<std::uint32_t>(seed));
137
0
            }
138
0
        } while (i < dimensionality_);
139
0
        QL_REQUIRE(++nextSequenceCounter_ != 0,
140
0
                   "Burley2020SobolRsg::nextInt32Sequence(): period exceeded");
141
0
        return integerSequence_;
142
0
    }
143
144
0
    const SobolRsg::sample_type& Burley2020SobolRsg::nextSequence() const {
145
0
        const std::vector<std::uint32_t>& v = nextInt32Sequence();
146
        // normalize to get a double in (0,1)
147
0
        for (Size k = 0; k < dimensionality_; ++k) {
148
0
            sequence_.value[k] = (static_cast<double>(v[k]) + 0.5) / 4294967296.0;
149
0
        }
150
0
        return sequence_;
151
0
    }
152
}