diff --git a/doc/generators.qbk b/doc/generators.qbk index 863a6741e..90391494b 100644 --- a/doc/generators.qbk +++ b/doc/generators.qbk @@ -116,6 +116,13 @@ If you need a 32-bit return instead of 64-bit prefer __xoshiro128pp. If you need either double or float as the return type prefer __xoshiro256d or __xoshiro128f respectively. +The xoshiro generators compute `discard(z)` in time logarithmic in `z`. +Besides `jump()` and `long_jump()`, they provide `jump(c, e)`, which advances +the generator by c*2[sup e] steps in time logarithmic in c plus linear in e, +to obtain non-overlapping subsequences at arbitrary distances: for example, +for the xoshiro256 generators `jump(1, 128)` and `jump(1, 192)` are equivalent +to `jump()` and `long_jump()`. + [note These random number generators are not intended for use in applications where non-deterministic random numbers are required. See __random_device for a choice of (hopefully) non-deterministic random number generators.] diff --git a/include/boost/random/detail/f2x_arithmetic.hpp b/include/boost/random/detail/f2x_arithmetic.hpp new file mode 100644 index 000000000..36d2f3250 --- /dev/null +++ b/include/boost/random/detail/f2x_arithmetic.hpp @@ -0,0 +1,115 @@ +/* + * Copyright Sebastiano Vigna 2026. + * Distributed under the Boost Software License, Version 1.0. (See + * accompanying file LICENSE_1_0.txt or copy at + * http://www.boost.org/LICENSE_1_0.txt) + * + * See http://www.boost.org for most recent version including documentation. + */ + +#ifndef BOOST_RANDOM_DETAIL_F2X_ARITHMETIC_HPP +#define BOOST_RANDOM_DETAIL_F2X_ARITHMETIC_HPP + +#include +#include +#include + +namespace boost { +namespace random { +namespace detail { + +// Arithmetic in the ring F_2[x] / (p), where p is the characteristic +// polynomial of the transition map of a linear generator over F_2 with +// 64 * W bits of state. Calling next() once is equivalent to multiplying the +// state by x in this ring; thus, advancing the generator by n steps is +// equivalent to multiplying by x^n mod p, and the polynomial x^n mod p can be +// applied to the state with the same accumulate-and-step loop used by the +// predefined jump functions. This is a translation of f2x.c from +// https://prng.di.unimi.it. +// +// Polynomials are arrays of W words; coefficient k is bit k % 64 of word +// k / 64. The characteristic polynomial is monic of degree 64 * W, and only +// its 64 * W low-order coefficients are stored. + +// a <- a * x mod p +template +inline void f2x_mul_x(std::array& a, const std::array& p) noexcept +{ + std::uint64_t carry = 0; + for (std::size_t i = 0; i < W; ++i) + { + const std::uint64_t next_carry = a[i] >> 63; + a[i] = a[i] << 1 | carry; + carry = next_carry; + } + + // carry is now the coefficient of x^(64 * W) + if (carry != 0) + { + for (std::size_t i = 0; i < W; ++i) + { + a[i] ^= p[i]; + } + } +} + +// Returns a * b mod p, by Horner's rule on the coefficients of a. +template +inline std::array f2x_mul_mod(const std::array& a, + const std::array& b, + const std::array& p) noexcept +{ + std::array r {}; + for (std::size_t k = 64 * W; k-- > 0;) + { + f2x_mul_x(r, p); + if ((a[k / 64] >> (k % 64)) & 1) + { + for (std::size_t i = 0; i < W; ++i) + { + r[i] ^= b[i]; + } + } + } + + return r; +} + +// Returns x^(c * 2^e) mod p, by square-and-multiply on the bits of c, +// followed by e squarings. +template +inline std::array f2x_jump_poly(const std::uint64_t c, const unsigned int e, + const std::array& p) noexcept +{ + std::array r {}; + r[0] = 1; + + // Leading zeros of c can be skipped, as they would just square 1 + int k = 63; + while (k >= 0 && ((c >> k) & 1) == 0) + { + --k; + } + + for (; k >= 0; --k) + { + r = f2x_mul_mod(r, r, p); + if ((c >> k) & 1) + { + f2x_mul_x(r, p); + } + } + + for (unsigned int i = 0; i < e; ++i) + { + r = f2x_mul_mod(r, r, p); + } + + return r; +} + +} // namespace detail +} // namespace random +} // namespace boost + +#endif // BOOST_RANDOM_DETAIL_F2X_ARITHMETIC_HPP diff --git a/include/boost/random/detail/xoshiro_base.hpp b/include/boost/random/detail/xoshiro_base.hpp index f9a794e7e..d8743da96 100644 --- a/include/boost/random/detail/xoshiro_base.hpp +++ b/include/boost/random/detail/xoshiro_base.hpp @@ -13,6 +13,7 @@ #define BOOST_RANDOM_DETAIL_XOSHIRO_BASE #include +#include #include #include #include @@ -251,6 +252,59 @@ class xoshiro_base state_[3] = s3; } + // The characteristic polynomials of the linear engines, without the + // leading term (see f2x_arithmetic.hpp) + static std::array charpoly(const std::integral_constant&) noexcept + { + return {{ UINT64_C(0x9d116f2bb0f0f001), UINT64_C(0x0280002bcefd1a5e), + UINT64_C(0x04b4edcf26259f85), UINT64_C(0x0003c03c3f3ecb19) }}; + } + + static std::array charpoly(const std::integral_constant&) noexcept + { + return {{ UINT64_C(0xcf3cff0c00000001), UINT64_C(0x7fdc78d886f00c63), + UINT64_C(0xf05e63fca6d7b781), UINT64_C(0x7a67058e7bbab6f0), + UINT64_C(0xf11eef832e32518f), UINT64_C(0x51ba7c47edc758ad), + UINT64_C(0x8f2d27268ce4b20b), UINT64_C(0x0000500055d8b77f) }}; + } + + static std::array charpoly(const std::integral_constant&) noexcept + { + return {{ UINT64_C(0x1b489db6de18fc01), UINT64_C(0x00fc65a2006254b1) }}; + } + + // Applies the jump polynomial poly (x^n mod the characteristic + // polynomial, for the desired distance n) to the state, using the same + // accumulate-and-step loop as jump(). + template + inline void jump_apply(const std::array& poly) noexcept + { + std::array t {}; + + for (std::size_t i = 0; i < W; ++i) + { + for (std::size_t b = 0; b < 64U; ++b) + { + if (poly[i] & UINT64_C(1) << b) + { + for (std::size_t w = 0; w < state_.size(); ++w) + { + t[w] ^= state_[w]; + } + } + + next(); + } + } + + state_ = t; + } + + // Below this distance, discard() steps the generator instead of jumping. + // The threshold is proportional to the size of the state, and it is close + // to the break-even point on current x86-64 processors. + static constexpr std::uint64_t discard_threshold = UINT64_C(2048) * N * sizeof(BlockType) / sizeof(std::uint64_t); + public: using result_type = OutputType; @@ -354,34 +408,66 @@ class xoshiro_base } /** This is the jump function for the generator. It is equivalent - * to 2^128 calls to next(); it can be used to generate 2^128 - * non-overlapping subsequences for parallel computations. */ + * to 2^(n/2) calls to next(), where n is the number of bits of state + * (2^64 for xoshiro128, 2^128 for xoshiro256, and 2^256 for + * xoshiro512); it can be used to generate 2^(n/2) non-overlapping + * subsequences for parallel computations. */ inline void jump() noexcept { jump_impl(xoshiro_type()); } /** This is the long-jump function for the generator. It is equivalent to - * 2^192 calls to next(); it can be used to generate 2^64 starting points, - * from each of which jump() will generate 2^64 non-overlapping + * 2^(3n/4) calls to next(), where n is the number of bits of state + * (2^96 for xoshiro128, 2^192 for xoshiro256, and 2^384 for + * xoshiro512); it can be used to generate 2^(n/4) starting points, + * from each of which jump() will generate 2^(n/4) non-overlapping * subsequences for parallel distributed computations. */ inline void long_jump() noexcept { long_jump_impl(xoshiro_type()); } + /** + * Advances the state of the generator by c * 2^e steps, in time + * logarithmic in c plus linear in e, by computing x^(c * 2^e) modulo + * the characteristic polynomial of the generator. This makes it possible + * to obtain non-overlapping subsequences at distances larger than + * discard() can express: for example, for xoshiro256 generators + * jump(1, 128) is equivalent to jump(), and jump(1, 192) is equivalent + * to long_jump(). For the jump to be meaningful, c * 2^e should be + * smaller than the period of the generator. + */ + inline void jump(const std::uint64_t c, const unsigned int e) noexcept + { + if (c != 0) + { + jump_apply(f2x_jump_poly(c, e, charpoly(xoshiro_type()))); + } + } + /** Returns the next value of the generator. */ inline result_type operator()() noexcept { return next(); } - /** Advances the state of the generator by @c z. */ + /** + * Advances the state of the generator by @c z steps. Large distances are + * computed in time logarithmic in @c z, using jump(z, 0). + */ inline void discard(const std::uint64_t z) noexcept { - for (std::uint64_t i {}; i < z; ++i) + if (z < discard_threshold) + { + for (std::uint64_t i {}; i < z; ++i) + { + next(); + } + } + else { - next(); + jump(z, 0); } } diff --git a/test/Jamfile.v2 b/test/Jamfile.v2 index b8fcc5799..d511a85e5 100644 --- a/test/Jamfile.v2 +++ b/test/Jamfile.v2 @@ -85,6 +85,7 @@ run test_xoshiro128mm.cpp /boost/test//boost_unit_test_framework ; run test_comp_xoshiro128mm.cpp ; run test_xoshiro128f.cpp /boost/test//boost_unit_test_framework ; run test_comp_xoshiro128f.cpp ; +run test_xoshiro_jump.cpp ; run test_mwc192x.cpp /boost/test//boost_unit_test_framework ; run test_comp_mwc192x.cpp ; run test_comp_mwc192x_portable.cpp ; diff --git a/test/test_xoshiro_jump.cpp b/test/test_xoshiro_jump.cpp new file mode 100644 index 000000000..5581e5589 --- /dev/null +++ b/test/test_xoshiro_jump.cpp @@ -0,0 +1,259 @@ +/* + * Copyright Sebastiano Vigna 2026. + * Distributed under the Boost Software License, Version 1.0. (See + * accompanying file LICENSE_1_0.txt or copy at + * http://www.boost.org/LICENSE_1_0.txt) + * + * Checks of discard() and jump(c, e) of the xoshiro generators. + */ + +#include +#include +#include +#include +#include +#include + +using namespace boost::random; + +// Builds a generator from its state. +template +Engine make_engine(std::initializer_list state) +{ + const std::vector v(state); + Engine e; + e.seed(v.begin(), v.end()); + return e; +} + +// Checks that discard(d) and jump(d, 0) are equivalent to d calls to +// operator()(), across the thresholds below which discard() steps the +// generator. +template +void check_discard(const Engine& start) +{ + static const std::uint64_t distances[] = { 0, 1, 2, 3, 7, 64, 1000, 4095, 4096, 4097, 8191, 8192, 8193, 9307, + 16383, 16384, 16385, 65537, 100000 }; + for (const std::uint64_t d : distances) + { + Engine a = start; + for (std::uint64_t i = 0; i < d; ++i) + { + a(); + } + + Engine b = start; + b.discard(d); + BOOST_TEST(a == b); + + Engine c = start; + c.jump(d, 0); + BOOST_TEST(a == c); + } +} + +// Checks that discard(u) followed by discard(v) is equivalent to discard(u + v). +template +void check_discard_composition(const Engine& start) +{ + static const std::uint64_t pairs[][2] = { + { 1, UINT64_C(0xfffffffffffffffe) }, + { 100000, UINT64_C(0x123456789abcdef) }, + { UINT64_C(0x7fffffffffffffff), UINT64_C(0x8000000000000000) }, + { UINT64_C(0x9e3779b97f4a7c15), UINT64_C(0x0123456789abcdef) } + }; + for (const auto& p : pairs) + { + Engine a = start; + a.discard(p[0]); + a.discard(p[1]); + Engine b = start; + b.discard(p[0] + p[1]); + BOOST_TEST(a == b); + } +} + +// Checks jump(c, e) for a generator with 'bits' bits of state, whose jump() +// and long_jump() advance it by 2^jump_log and 2^long_jump_log steps. +template +void check_jump(const Engine& start, const unsigned int bits, const unsigned int jump_log, const unsigned int long_jump_log) +{ + // jump(1, jump_log) and jump(1, long_jump_log) are jump() and long_jump() + { + Engine a = start; + a.jump(1, jump_log); + Engine b = start; + b.jump(); + BOOST_TEST(a == b); + + a = start; + a.jump(1, long_jump_log); + b = start; + b.long_jump(); + BOOST_TEST(a == b); + } + + // jump(0, e) does nothing + { + Engine a = start; + a.jump(0, 0); + BOOST_TEST(a == start); + a.jump(0, jump_log); + BOOST_TEST(a == start); + } + + // As the period is 2^bits - 1, a jump of 2^bits steps is a single step + { + Engine a = start; + a.jump(1, bits); + Engine b = start; + b(); + BOOST_TEST(a == b); + } + + // jump(c, e) is equivalent to jump(c * 2^e, 0) when c * 2^e fits in 64 bits + static const std::uint64_t ce[][2] = { { 3, 8 }, { 5, 20 }, { UINT64_C(0x1234567), 30 }, { 1, 63 }, { UINT64_C(0x9e3779b97f4a7c15), 0 } }; + for (const auto& p : ce) + { + Engine a = start; + a.jump(p[0], static_cast(p[1])); + Engine b = start; + b.jump(p[0] << p[1], 0); + BOOST_TEST(a == b); + } + + // jump(c1, e) followed by jump(c2, e) is equivalent to jump(c1 + c2, e) + const unsigned int es[] = { 0, 1, 64, jump_log, long_jump_log }; + for (const unsigned int e : es) + { + Engine a = start; + a.jump(3, e); + a.jump(UINT64_C(0x123456789abcdef), e); + Engine b = start; + b.jump(3 + UINT64_C(0x123456789abcdef), e); + BOOST_TEST(a == b); + } +} + +template +void check_all(const Engine& start, const unsigned int bits, const unsigned int jump_log, const unsigned int long_jump_log) +{ + check_discard(start); + check_discard(Engine(42)); + check_discard_composition(start); + check_jump(start, bits, jump_log, long_jump_log); +} + +// The states below were computed with jump_ce() of the reference +// implementations at https://prng.di.unimi.it. They depend only on the +// linear engine, so they are the same for all generators with the same state. + +template +void check_reference_128() +{ + const Engine start = make_engine({ UINT32_C(0x89abcdef), UINT32_C(0x92a217f3), UINT32_C(0x9b9861f7), UINT32_C(0xa48eabfb) }); + Engine e = start; + e.discard(12345); + BOOST_TEST(e == make_engine({ UINT32_C(0xfde876d9), UINT32_C(0x67184059), UINT32_C(0xa45d92c9), UINT32_C(0xb9a7a002) })); + e = start; + e.discard(UINT64_C(0x9e3779b97f4a7c15)); + BOOST_TEST(e == make_engine({ UINT32_C(0x791642c5), UINT32_C(0x575b3521), UINT32_C(0x1391e379), UINT32_C(0x9f89430e) })); + e = start; + e.discard(UINT64_C(0xffffffffffffffff)); + BOOST_TEST(e == make_engine({ UINT32_C(0x75b734e2), UINT32_C(0x79451d39), UINT32_C(0x7e7614d4), UINT32_C(0x0a70b3f2) })); + e = start; + e.jump(3, 70); + BOOST_TEST(e == make_engine({ UINT32_C(0x6a17ac7b), UINT32_C(0xae4fb410), UINT32_C(0x86f3bf52), UINT32_C(0xde12fa78) })); + e = start; + e.jump(UINT64_C(0x123456789abcdef), 71); + BOOST_TEST(e == make_engine({ UINT32_C(0xd9ef01fc), UINT32_C(0xbcf9da3f), UINT32_C(0xd7172609), UINT32_C(0x64c3a091) })); + + check_all(start, 128, 64, 96); +} + +template +void check_reference_256() +{ + const Engine start = make_engine({ UINT64_C(0x0123456789abcdef), UINT64_C(0xa07e048892a217f3), + UINT64_C(0x3fd8c3a99b9861f7), UINT64_C(0xdf3382caa48eabfb) }); + Engine e = start; + e.discard(12345); + BOOST_TEST(e == make_engine({ UINT64_C(0x69cf32f8d724cd8b), UINT64_C(0x7c331e0a4fec602e), + UINT64_C(0xd7d968ad7dcb2779), UINT64_C(0xf1e6137266432579) })); + e = start; + e.discard(UINT64_C(0x9e3779b97f4a7c15)); + BOOST_TEST(e == make_engine({ UINT64_C(0x70281a74ac355027), UINT64_C(0xa42ca095899e9167), + UINT64_C(0x9d1835a0865afdd1), UINT64_C(0xdf4db42b007e379d) })); + e = start; + e.discard(UINT64_C(0xffffffffffffffff)); + BOOST_TEST(e == make_engine({ UINT64_C(0x91e3122b2590ba02), UINT64_C(0xa608df3af94e5571), + UINT64_C(0x83d5db190ae2b3c3), UINT64_C(0xc2a8e84c4217a2b6) })); + e = start; + e.jump(3, 70); + BOOST_TEST(e == make_engine({ UINT64_C(0xdd829b8601731221), UINT64_C(0x8ad107e028c5a37b), + UINT64_C(0xbfc8df29986dd8b4), UINT64_C(0xff9b643b5cb48331) })); + e = start; + e.jump(UINT64_C(0x123456789abcdef), 135); + BOOST_TEST(e == make_engine({ UINT64_C(0x13da6f31fe6072de), UINT64_C(0x23d115d50961ce4c), + UINT64_C(0x990383bc8f71e2f3), UINT64_C(0x29e09f9e78b1d2b4) })); + + check_all(start, 256, 128, 192); +} + +template +void check_reference_512() +{ + const Engine start = make_engine({ UINT64_C(0x0123456789abcdef), UINT64_C(0xa07e048892a217f3), + UINT64_C(0x3fd8c3a99b9861f7), UINT64_C(0xdf3382caa48eabfb), + UINT64_C(0x7e8e41ebad84f5ff), UINT64_C(0x1de9010cb67b4003), + UINT64_C(0xbd43c02dbf718a07), UINT64_C(0x5c9e7f4ec867d40b) }); + Engine e = start; + e.discard(12345); + BOOST_TEST(e == make_engine({ UINT64_C(0x670dfab1a4dbd932), UINT64_C(0xde097f4981d3f4d1), + UINT64_C(0xec2e66f685013683), UINT64_C(0xca4df29dc478e3d5), + UINT64_C(0x7d0ffb16c7d5db22), UINT64_C(0x9cb418f689f884dc), + UINT64_C(0x3269b7eed6746a30), UINT64_C(0x6dd164f43c3c3299) })); + e = start; + e.discard(UINT64_C(0x9e3779b97f4a7c15)); + BOOST_TEST(e == make_engine({ UINT64_C(0x63d256ba236fe699), UINT64_C(0x9eb057bedbf8a9e8), + UINT64_C(0xd55b8f541d8df6bb), UINT64_C(0xcb3cda5604a5c200), + UINT64_C(0x6484c0edfd9c1574), UINT64_C(0x0fb55f7648bc7ce1), + UINT64_C(0x44e19ae84b723000), UINT64_C(0xdb56ea6746568ebb) })); + e = start; + e.discard(UINT64_C(0xffffffffffffffff)); + BOOST_TEST(e == make_engine({ UINT64_C(0x03d11e951e72419a), UINT64_C(0x775e3fd2797ae991), + UINT64_C(0x6e525c9f87143faf), UINT64_C(0x40f37b28d4a59ecb), + UINT64_C(0xcc5ffade71259ecb), UINT64_C(0x9f8fc85247f83135), + UINT64_C(0xbe8a071b77f1aa3f), UINT64_C(0x4c23cb4260054391) })); + e = start; + e.jump(3, 70); + BOOST_TEST(e == make_engine({ UINT64_C(0x15710a2cbcebaddd), UINT64_C(0xca329d2b933fa0e3), + UINT64_C(0x86192c2a1622632f), UINT64_C(0xd7dab9c5a634aff7), + UINT64_C(0x5bc1fd588730cf0f), UINT64_C(0x2da6036b5f331de7), + UINT64_C(0x3630f4bd57b7278c), UINT64_C(0x85e80425a2c59fd8) })); + e = start; + e.jump(UINT64_C(0x123456789abcdef), 263); + BOOST_TEST(e == make_engine({ UINT64_C(0x1716bf48194e0779), UINT64_C(0xeeee48e6c78df93d), + UINT64_C(0x29f9411264d13131), UINT64_C(0x54ccce075543e94f), + UINT64_C(0x960e7688bcef2b39), UINT64_C(0x72ff68e6a85adbe6), + UINT64_C(0xd2406b47887bc7ab), UINT64_C(0xf71e945d33fd21f2) })); + + check_all(start, 512, 256, 384); +} + +int main() +{ + check_reference_128(); + check_reference_128(); + check_reference_128(); + + check_reference_256(); + check_reference_256(); + check_reference_256(); + + check_reference_512(); + check_reference_512(); + check_reference_512(); + + return boost::report_errors(); +}