random/detail/xoshiro_base.hpp

0001 /*
0002  * Copyright Matt Borland 2022 - 2025.
0003  * Distributed under the Boost Software License, Version 1.0. (See
0004         * accompanying file LICENSE_1_0.txt or copy at
0005  * http://www.boost.org/LICENSE_1_0.txt)
0006  *
0007  * See http://www.boost.org for most recent version including documentation.
0008  *
0009  * $Id$
0010  */
0011
0012 #ifndef BOOST_RANDOM_DETAIL_XOSHIRO_BASE
0013 #define BOOST_RANDOM_DETAIL_XOSHIRO_BASE
0014
0015 #include <boost/random/splitmix64.hpp>
0016 #include <boost/random/detail/seed.hpp>
0017 #include <boost/throw_exception.hpp>
0018 #include <boost/config.hpp>
0019 #include <array>
0020 #include <utility>
0021 #include <stdexcept>
0022 #include <limits>
0023 #include <initializer_list>
0024 #include <cstdint>
0025 #include <cstdlib>
0026 #include <string>
0027 #include <ios>
0028 #include <istream>
0029 #include <type_traits>
0030 #include <iterator>
0031
0032 namespace boost {
0033 namespace random {
0034 namespace detail {
0035
0036 // N is the number of words (e.g. for xoshiro 256 N=4)
0037 template <typename Derived, std::size_t N, typename OutputType = std::uint64_t, typename BlockType = std::uint64_t>
0038 class xoshiro_base
0039 {
0040 protected:
0041
0042     std::array<BlockType, N> state_;
0043
0044 private:
0045
0046     using xoshiro_type = std::integral_constant<BlockType, N>;
0047
0048     inline std::uint64_t concatenate(std::uint32_t word1, std::uint32_t word2) noexcept
0049     {
0050         return static_cast<std::uint64_t>(word1) << 32U | word2;
0051     }
0052
0053     template <typename Sseq>
0054     inline void sseq_seed_64(Sseq& seq)
0055     {
0056         for (auto& i : state_)
0057         {
0058             std::array<std::uint32_t, 2> seeds;
0059             seq.generate(seeds.begin(), seeds.end());
0060
0061             i = concatenate(seeds[0], seeds[1]);
0062         }
0063     }
0064
0065     template <typename Sseq>
0066     inline void sseq_seed_32(Sseq& seq)
0067     {
0068         seq.generate(state_.begin(), state_.end());
0069     }
0070
0071     inline void jump_impl(const std::integral_constant<std::uint64_t, 4>&) noexcept
0072     {
0073         constexpr std::array<std::uint64_t, 4U> jump = {{ UINT64_C(0x180ec6d33cfd0aba), UINT64_C(0xd5a61266f0c9392c),
0074                                                           UINT64_C(0xa9582618e03fc9aa), UINT64_C(0x39abdc4529b1661c) }};
0075
0076         std::uint64_t s0 = 0;
0077         std::uint64_t s1 = 0;
0078         std::uint64_t s2 = 0;
0079         std::uint64_t s3 = 0;
0080
0081         for (std::uint64_t i = 0; i < jump.size(); i++)
0082         {
0083             for (std::size_t b = 0; b < 64U; b++)
0084             {
0085                 if (jump[i] & UINT64_C(1) << b) {
0086                     s0 ^= state_[0];
0087                     s1 ^= state_[1];
0088                     s2 ^= state_[2];
0089                     s3 ^= state_[3];
0090                 }
0091
0092                 next();
0093             }
0094         }
0095
0096         state_[0] = s0;
0097         state_[1] = s1;
0098         state_[2] = s2;
0099         state_[3] = s3;
0100     }
0101
0102     inline void jump_impl(const std::integral_constant<std::uint64_t, 8>&) noexcept
0103     {
0104         constexpr std::array<std::uint64_t, 8U> jump = {{ UINT64_C(0x33ed89b6e7a353f9), UINT64_C(0x760083d7955323be),
0105                                                           UINT64_C(0x2837f2fbb5f22fae), UINT64_C(0x4b8c5674d309511c),
0106                                                           UINT64_C(0xb11ac47a7ba28c25), UINT64_C(0xf1be7667092bcc1c),
0107                                                           UINT64_C(0x53851efdb6df0aaf), UINT64_C(0x1ebbc8b23eaf25db) }};
0108
0109         std::array<std::uint64_t, 8U> t = {{ 0, 0, 0, 0, 0, 0, 0, 0 }};
0110
0111         for (std::size_t i = 0; i < jump.size(); ++i)
0112         {
0113             for (std::size_t b = 0; b < 64U; ++b)
0114             {
0115                 if (jump[i] & UINT64_C(1) << b)
0116                 {
0117                     for (std::size_t w = 0; w < state_.size(); ++w)
0118                     {
0119                         t[w] ^= state_[w];
0120                     }
0121                 }
0122
0123                 next();
0124             }
0125         }
0126
0127         state_ = t;
0128     }
0129
0130     inline void long_jump_impl(const std::integral_constant<std::uint64_t, 4>&) noexcept
0131     {
0132         constexpr std::array<std::uint64_t, 4> long_jump = {{ UINT64_C(0x76e15d3efefdcbbf), UINT64_C(0xc5004e441c522fb3),
0133                                                               UINT64_C(0x77710069854ee241), UINT64_C(0x39109bb02acbe635) }};
0134
0135         std::uint64_t s0 = 0;
0136         std::uint64_t s1 = 0;
0137         std::uint64_t s2 = 0;
0138         std::uint64_t s3 = 0;
0139
0140         for (std::size_t i = 0; i < long_jump.size(); i++)
0141         {
0142             for (std::size_t b = 0; b < 64; b++)
0143             {
0144                 if (long_jump[i] & UINT64_C(1) << b)
0145                 {
0146                     s0 ^= state_[0];
0147                     s1 ^= state_[1];
0148                     s2 ^= state_[2];
0149                     s3 ^= state_[3];
0150                 }
0151
0152                 next();
0153             }
0154         }
0155
0156         state_[0] = s0;
0157         state_[1] = s1;
0158         state_[2] = s2;
0159         state_[3] = s3;
0160     }
0161
0162     inline void long_jump_impl(const std::integral_constant<std::uint64_t, 8>&) noexcept
0163     {
0164         constexpr std::array<std::uint64_t, 8U> long_jump = {{ UINT64_C(0x11467fef8f921d28), UINT64_C(0xa2a819f2e79c8ea8),
0165                                                                UINT64_C(0xa8299fc284b3959a), UINT64_C(0xb4d347340ca63ee1),
0166                                                                UINT64_C(0x1cb0940bedbff6ce), UINT64_C(0xd956c5c4fa1f8e17),
0167                                                                UINT64_C(0x915e38fd4eda93bc), UINT64_C(0x5b3ccdfa5d7daca5) }};
0168
0169         std::array<std::uint64_t, 8U> t = {{ 0, 0, 0, 0, 0, 0, 0, 0 }};
0170
0171         for (std::size_t i = 0; i < long_jump.size(); ++i)
0172         {
0173             for (std::size_t b = 0; b < 64U; ++b)
0174             {
0175                 if (long_jump[i] & UINT64_C(1) << b)
0176                 {
0177                     for (std::size_t w = 0; w < state_.size(); ++w)
0178                     {
0179                         t[w] ^= state_[w];
0180                     }
0181                 }
0182
0183                 next();
0184             }
0185         }
0186
0187         state_ = t;
0188     }
0189
0190     inline void jump_impl(const std::integral_constant<std::uint32_t, 4>&) noexcept
0191     {
0192         constexpr std::array<std::uint32_t, 4> jump = {{ UINT32_C(0x8764000b), UINT32_C(0xf542d2d3),
0193                                                          UINT32_C(0x6fa035c3), UINT32_C(0x77f2db5b) }};
0194
0195         std::uint32_t s0 = 0;
0196         std::uint32_t s1 = 0;
0197         std::uint32_t s2 = 0;
0198         std::uint32_t s3 = 0;
0199
0200         for (std::size_t i = 0; i < jump.size(); i++)
0201         {
0202             for (std::size_t b = 0; b < 32U; b++)
0203             {
0204                 if (jump[i] & UINT32_C(1) << b)
0205                 {
0206                     s0 ^= state_[0];
0207                     s1 ^= state_[1];
0208                     s2 ^= state_[2];
0209                     s3 ^= state_[3];
0210                 }
0211
0212                 next();
0213             }
0214         }
0215
0216         state_[0] = s0;
0217         state_[1] = s1;
0218         state_[2] = s2;
0219         state_[3] = s3;
0220     }
0221
0222     inline void long_jump_impl(const std::integral_constant<std::uint32_t, 4>&) noexcept
0223     {
0224         constexpr std::array<std::uint32_t, 4> jump = {{ UINT32_C(0xb523952e), UINT32_C(0x0b6f099f),
0225                                                          UINT32_C(0xccf5a0ef), UINT32_C(0x1c580662) }};
0226
0227         std::uint32_t s0 = 0;
0228         std::uint32_t s1 = 0;
0229         std::uint32_t s2 = 0;
0230         std::uint32_t s3 = 0;
0231
0232         for (std::size_t i = 0; i < jump.size(); i++)
0233         {
0234             for (std::size_t b = 0; b < 32; b++)
0235             {
0236                 if (jump[i] & UINT32_C(1) << b)
0237                 {
0238                     s0 ^= state_[0];
0239                     s1 ^= state_[1];
0240                     s2 ^= state_[2];
0241                     s3 ^= state_[3];
0242                 }
0243
0244                 next();
0245             }
0246         }
0247
0248         state_[0] = s0;
0249         state_[1] = s1;
0250         state_[2] = s2;
0251         state_[3] = s3;
0252     }
0253
0254 public:
0255
0256     using result_type = OutputType;
0257     using seed_type = BlockType;
0258
0259     static constexpr bool has_fixed_range {false};
0260
0261     /** Seeds the generator using the default seed of boost::random::splitmix64 */
0262     void seed()
0263     {
0264         splitmix64 gen;
0265         for (auto& i : state_)
0266         {
0267             i = static_cast<seed_type>(gen());
0268         }
0269     }
0270
0271     /** Seeds the generator with a user provided seed. */
0272     void seed(const seed_type value)
0273     {
0274         splitmix64 gen(value);
0275         for (auto& i : state_)
0276         {
0277             i = static_cast<seed_type>(gen());
0278         }
0279     }
0280
0281     /**
0282      * Seeds the generator with 32-bit values produced by @c seq.generate().
0283      */
0284     template <typename Sseq, typename std::enable_if<!std::is_convertible<Sseq, seed_type>::value, bool>::type = true>
0285     void seed(Sseq& seq)
0286     {
0287         BOOST_IF_CONSTEXPR (std::is_same<BlockType, std::uint64_t>::value)
0288         {
0289             sseq_seed_64(seq);
0290         }
0291         else
0292         {
0293             sseq_seed_32(seq);
0294         }
0295     }
0296
0297     /** Sets the state of the generator using values from an iterator range. */
0298     template <typename FIter>
0299     void seed(FIter first, FIter last)
0300     {
0301         static_assert(std::is_integral<typename std::iterator_traits<FIter>::value_type>::value,
0302                       "Value type must be a built-in integer type" );
0303
0304         std::size_t offset = 0;
0305         while (first != last && offset < state_.size())
0306         {
0307             state_[offset++] = static_cast<seed_type>(*first++);
0308         }
0309
0310         if (offset != state_.size())
0311         {
0312             boost::throw_exception(std::invalid_argument("Not enough elements in call to seed."));
0313         }
0314     }
0315
0316     /**
0317      * Constructs a @c xoshiro and calls @c seed().
0318      */
0319     xoshiro_base() { seed(); }
0320
0321     /** Seeds the generator with a user provided seed. */
0322     explicit xoshiro_base(const seed_type value)
0323     {
0324         seed(value);
0325     }
0326
0327     template <typename FIter>
0328     xoshiro_base(FIter& first, FIter last) { seed(first, last); }
0329
0330     /**
0331      * Seeds the generator with 64-bit values produced by @c seq.generate().
0332      *
0333      * @xmlnote
0334      * The copy constructor will always be preferred over
0335      * the templated constructor.
0336      * @endxmlnote
0337      */
0338     template <typename Sseq, typename std::enable_if<!std::is_convertible<Sseq, xoshiro_base>::value, bool>::type = true>
0339     explicit xoshiro_base(Sseq& seq)
0340     {
0341         seed(seq);
0342     }
0343
0344     // Hit all of our rule of 5 explicitly to ensure old platforms work correctly
0345     ~xoshiro_base() = default;
0346     xoshiro_base(const xoshiro_base& other) noexcept { state_ = other.state(); }
0347     xoshiro_base& operator=(const xoshiro_base& other) noexcept { state_ = other.state(); return *this; }
0348     xoshiro_base(xoshiro_base&& other) noexcept { state_ = other.state(); }
0349     xoshiro_base& operator=(xoshiro_base&& other) noexcept { state_ = other.state(); return *this; }
0350
0351     inline result_type next() noexcept
0352     {
0353         return static_cast<Derived*>(this)->next();
0354     }
0355
0356     /** This is the jump function for the generator. It is equivalent
0357      *  to 2^128 calls to next(); it can be used to generate 2^128
0358      *  non-overlapping subsequences for parallel computations. */
0359     inline void jump() noexcept
0360     {
0361         jump_impl(xoshiro_type());
0362     }
0363
0364     /** This is the long-jump function for the generator. It is equivalent to
0365      *  2^192 calls to next(); it can be used to generate 2^64 starting points,
0366      *  from each of which jump() will generate 2^64 non-overlapping
0367      *  subsequences for parallel distributed computations. */
0368     inline void long_jump() noexcept
0369     {
0370         long_jump_impl(xoshiro_type());
0371     }
0372
0373     /**  Returns the next value of the generator. */
0374     inline result_type operator()() noexcept
0375     {
0376         return next();
0377     }
0378
0379     /** Advances the state of the generator by @c z. */
0380     inline void discard(const std::uint64_t z) noexcept
0381     {
0382         for (std::uint64_t i {}; i < z; ++i)
0383         {
0384             next();
0385         }
0386     }
0387
0388     /**
0389      * Returns true if the two generators will produce identical
0390      * sequences of values.
0391      */
0392     inline friend bool operator==(const xoshiro_base& lhs, const xoshiro_base& rhs) noexcept
0393     {
0394         return lhs.state_ == rhs.state_;
0395     }
0396
0397     /**
0398      * Returns true if the two generators will produce different
0399      * sequences of values.
0400      */
0401     inline friend bool operator!=(const xoshiro_base& lhs, const xoshiro_base& rhs) noexcept
0402     {
0403         return lhs.state_ != rhs.state_;
0404     }
0405
0406     /**  Writes a @c xorshiro to a @c std::ostream. */
0407     template <typename CharT, typename Traits>
0408     inline friend std::basic_ostream<CharT, Traits>& operator<<(std::basic_ostream<CharT, Traits>& ost,
0409                                                                 const xoshiro_base& e)
0410     {
0411         for (std::size_t i {}; i < e.state_.size(); ++i)
0412         {
0413             ost << e.state_[i] << ' ';
0414         }
0415
0416         return ost;
0417     }
0418
0419     /**  Reads a @c xorshiro from a @c std::istream. */
0420     template <typename CharT, typename Traits>
0421     inline friend std::basic_istream<CharT, Traits>& operator>>(std::basic_istream<CharT, Traits>& ist,
0422                                                                 xoshiro_base& e)
0423     {
0424         for (std::size_t i {}; i < e.state_.size(); ++i)
0425         {
0426             ist >> e.state_[i] >> std::ws;
0427         }
0428
0429         return ist;
0430     }
0431
0432     /** Fills a range with random values */
0433     template <typename FIter>
0434     inline void generate(FIter first, FIter last) noexcept
0435     {
0436         using iter_type = typename std::iterator_traits<FIter>::value_type;
0437
0438         while (first != last)
0439         {
0440             *first++ = static_cast<iter_type>(next());
0441         }
0442     }
0443
0444     /**
0445      * Returns the largest value that the @c xorshiro
0446      * can produce.
0447      */
0448     static constexpr result_type (max)() noexcept
0449     {
0450         return (std::numeric_limits<result_type>::max)();
0451     }
0452
0453     /**
0454      * Returns the smallest value that the @c xorshiro
0455      * can produce.
0456      */
0457     static constexpr result_type (min)() noexcept
0458     {
0459         return (std::numeric_limits<result_type>::min)();
0460     }
0461
0462     inline std::array<BlockType, N> state() const noexcept
0463     {
0464         return state_;
0465     }
0466 };
0467
0468 } // namespace detail
0469 } // namespace random
0470 } // namespace boost
0471
0472 #endif