xoshiro_base.hpp 14 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473
  1. /*
  2. * Copyright Matt Borland 2022 - 2025.
  3. * Distributed under the Boost Software License, Version 1.0. (See
  4. * accompanying file LICENSE_1_0.txt or copy at
  5. * http://www.boost.org/LICENSE_1_0.txt)
  6. *
  7. * See http://www.boost.org for most recent version including documentation.
  8. *
  9. * $Id$
  10. */
  11. #ifndef BOOST_RANDOM_DETAIL_XOSHIRO_BASE
  12. #define BOOST_RANDOM_DETAIL_XOSHIRO_BASE
  13. #include <boost/random/splitmix64.hpp>
  14. #include <boost/random/detail/seed.hpp>
  15. #include <boost/throw_exception.hpp>
  16. #include <boost/config.hpp>
  17. #include <array>
  18. #include <utility>
  19. #include <stdexcept>
  20. #include <limits>
  21. #include <initializer_list>
  22. #include <cstdint>
  23. #include <cstdlib>
  24. #include <string>
  25. #include <ios>
  26. #include <istream>
  27. #include <type_traits>
  28. #include <iterator>
  29. namespace boost {
  30. namespace random {
  31. namespace detail {
  32. // N is the number of words (e.g. for xoshiro 256 N=4)
  33. template <typename Derived, std::size_t N, typename OutputType = std::uint64_t, typename BlockType = std::uint64_t>
  34. class xoshiro_base
  35. {
  36. protected:
  37. std::array<BlockType, N> state_;
  38. private:
  39. using xoshiro_type = std::integral_constant<BlockType, N>;
  40. inline std::uint64_t concatenate(std::uint32_t word1, std::uint32_t word2) noexcept
  41. {
  42. return static_cast<std::uint64_t>(word1) << 32U | word2;
  43. }
  44. template <typename Sseq>
  45. inline void sseq_seed_64(Sseq& seq)
  46. {
  47. std::array<std::uint32_t, N * 2> seeds;
  48. seq.generate(seeds.begin(), seeds.end());
  49. for (std::size_t i = 0; i < state_.size(); ++i)
  50. {
  51. state_[i] = concatenate(seeds[2*i], seeds[2*i + 1]);
  52. }
  53. }
  54. template <typename Sseq>
  55. inline void sseq_seed_32(Sseq& seq)
  56. {
  57. seq.generate(state_.begin(), state_.end());
  58. }
  59. inline void jump_impl(const std::integral_constant<std::uint64_t, 4>&) noexcept
  60. {
  61. constexpr std::array<std::uint64_t, 4U> jump = {{ UINT64_C(0x180ec6d33cfd0aba), UINT64_C(0xd5a61266f0c9392c),
  62. UINT64_C(0xa9582618e03fc9aa), UINT64_C(0x39abdc4529b1661c) }};
  63. std::uint64_t s0 = 0;
  64. std::uint64_t s1 = 0;
  65. std::uint64_t s2 = 0;
  66. std::uint64_t s3 = 0;
  67. for (std::uint64_t i = 0; i < jump.size(); i++)
  68. {
  69. for (std::size_t b = 0; b < 64U; b++)
  70. {
  71. if (jump[i] & UINT64_C(1) << b) {
  72. s0 ^= state_[0];
  73. s1 ^= state_[1];
  74. s2 ^= state_[2];
  75. s3 ^= state_[3];
  76. }
  77. next();
  78. }
  79. }
  80. state_[0] = s0;
  81. state_[1] = s1;
  82. state_[2] = s2;
  83. state_[3] = s3;
  84. }
  85. inline void jump_impl(const std::integral_constant<std::uint64_t, 8>&) noexcept
  86. {
  87. constexpr std::array<std::uint64_t, 8U> jump = {{ UINT64_C(0x33ed89b6e7a353f9), UINT64_C(0x760083d7955323be),
  88. UINT64_C(0x2837f2fbb5f22fae), UINT64_C(0x4b8c5674d309511c),
  89. UINT64_C(0xb11ac47a7ba28c25), UINT64_C(0xf1be7667092bcc1c),
  90. UINT64_C(0x53851efdb6df0aaf), UINT64_C(0x1ebbc8b23eaf25db) }};
  91. std::array<std::uint64_t, 8U> t = {{ 0, 0, 0, 0, 0, 0, 0, 0 }};
  92. for (std::size_t i = 0; i < jump.size(); ++i)
  93. {
  94. for (std::size_t b = 0; b < 64U; ++b)
  95. {
  96. if (jump[i] & UINT64_C(1) << b)
  97. {
  98. for (std::size_t w = 0; w < state_.size(); ++w)
  99. {
  100. t[w] ^= state_[w];
  101. }
  102. }
  103. next();
  104. }
  105. }
  106. state_ = t;
  107. }
  108. inline void long_jump_impl(const std::integral_constant<std::uint64_t, 4>&) noexcept
  109. {
  110. constexpr std::array<std::uint64_t, 4> long_jump = {{ UINT64_C(0x76e15d3efefdcbbf), UINT64_C(0xc5004e441c522fb3),
  111. UINT64_C(0x77710069854ee241), UINT64_C(0x39109bb02acbe635) }};
  112. std::uint64_t s0 = 0;
  113. std::uint64_t s1 = 0;
  114. std::uint64_t s2 = 0;
  115. std::uint64_t s3 = 0;
  116. for (std::size_t i = 0; i < long_jump.size(); i++)
  117. {
  118. for (std::size_t b = 0; b < 64; b++)
  119. {
  120. if (long_jump[i] & UINT64_C(1) << b)
  121. {
  122. s0 ^= state_[0];
  123. s1 ^= state_[1];
  124. s2 ^= state_[2];
  125. s3 ^= state_[3];
  126. }
  127. next();
  128. }
  129. }
  130. state_[0] = s0;
  131. state_[1] = s1;
  132. state_[2] = s2;
  133. state_[3] = s3;
  134. }
  135. inline void long_jump_impl(const std::integral_constant<std::uint64_t, 8>&) noexcept
  136. {
  137. constexpr std::array<std::uint64_t, 8U> long_jump = {{ UINT64_C(0x11467fef8f921d28), UINT64_C(0xa2a819f2e79c8ea8),
  138. UINT64_C(0xa8299fc284b3959a), UINT64_C(0xb4d347340ca63ee1),
  139. UINT64_C(0x1cb0940bedbff6ce), UINT64_C(0xd956c5c4fa1f8e17),
  140. UINT64_C(0x915e38fd4eda93bc), UINT64_C(0x5b3ccdfa5d7daca5) }};
  141. std::array<std::uint64_t, 8U> t = {{ 0, 0, 0, 0, 0, 0, 0, 0 }};
  142. for (std::size_t i = 0; i < long_jump.size(); ++i)
  143. {
  144. for (std::size_t b = 0; b < 64U; ++b)
  145. {
  146. if (long_jump[i] & UINT64_C(1) << b)
  147. {
  148. for (std::size_t w = 0; w < state_.size(); ++w)
  149. {
  150. t[w] ^= state_[w];
  151. }
  152. }
  153. next();
  154. }
  155. }
  156. state_ = t;
  157. }
  158. inline void jump_impl(const std::integral_constant<std::uint32_t, 4>&) noexcept
  159. {
  160. constexpr std::array<std::uint32_t, 4> jump = {{ UINT32_C(0x8764000b), UINT32_C(0xf542d2d3),
  161. UINT32_C(0x6fa035c3), UINT32_C(0x77f2db5b) }};
  162. std::uint32_t s0 = 0;
  163. std::uint32_t s1 = 0;
  164. std::uint32_t s2 = 0;
  165. std::uint32_t s3 = 0;
  166. for (std::size_t i = 0; i < jump.size(); i++)
  167. {
  168. for (std::size_t b = 0; b < 32U; b++)
  169. {
  170. if (jump[i] & UINT32_C(1) << b)
  171. {
  172. s0 ^= state_[0];
  173. s1 ^= state_[1];
  174. s2 ^= state_[2];
  175. s3 ^= state_[3];
  176. }
  177. next();
  178. }
  179. }
  180. state_[0] = s0;
  181. state_[1] = s1;
  182. state_[2] = s2;
  183. state_[3] = s3;
  184. }
  185. inline void long_jump_impl(const std::integral_constant<std::uint32_t, 4>&) noexcept
  186. {
  187. constexpr std::array<std::uint32_t, 4> jump = {{ UINT32_C(0xb523952e), UINT32_C(0x0b6f099f),
  188. UINT32_C(0xccf5a0ef), UINT32_C(0x1c580662) }};
  189. std::uint32_t s0 = 0;
  190. std::uint32_t s1 = 0;
  191. std::uint32_t s2 = 0;
  192. std::uint32_t s3 = 0;
  193. for (std::size_t i = 0; i < jump.size(); i++)
  194. {
  195. for (std::size_t b = 0; b < 32; b++)
  196. {
  197. if (jump[i] & UINT32_C(1) << b)
  198. {
  199. s0 ^= state_[0];
  200. s1 ^= state_[1];
  201. s2 ^= state_[2];
  202. s3 ^= state_[3];
  203. }
  204. next();
  205. }
  206. }
  207. state_[0] = s0;
  208. state_[1] = s1;
  209. state_[2] = s2;
  210. state_[3] = s3;
  211. }
  212. public:
  213. using result_type = OutputType;
  214. using seed_type = BlockType;
  215. static constexpr bool has_fixed_range {false};
  216. /** Seeds the generator using the default seed of boost::random::splitmix64 */
  217. void seed()
  218. {
  219. splitmix64 gen;
  220. for (auto& i : state_)
  221. {
  222. i = static_cast<seed_type>(gen());
  223. }
  224. }
  225. /** Seeds the generator with a user provided seed. */
  226. void seed(const seed_type value)
  227. {
  228. splitmix64 gen(value);
  229. for (auto& i : state_)
  230. {
  231. i = static_cast<seed_type>(gen());
  232. }
  233. }
  234. /**
  235. * Seeds the generator with 32-bit values produced by @c seq.generate().
  236. */
  237. template <typename Sseq, typename std::enable_if<!std::is_convertible<Sseq, seed_type>::value, bool>::type = true>
  238. void seed(Sseq& seq)
  239. {
  240. BOOST_IF_CONSTEXPR (std::is_same<BlockType, std::uint64_t>::value)
  241. {
  242. sseq_seed_64(seq);
  243. }
  244. else
  245. {
  246. sseq_seed_32(seq);
  247. }
  248. }
  249. /** Sets the state of the generator using values from an iterator range. */
  250. template <typename FIter>
  251. void seed(FIter first, FIter last)
  252. {
  253. static_assert(std::is_integral<typename std::iterator_traits<FIter>::value_type>::value,
  254. "Value type must be a built-in integer type" );
  255. std::size_t offset = 0;
  256. while (first != last && offset < state_.size())
  257. {
  258. state_[offset++] = static_cast<seed_type>(*first++);
  259. }
  260. if (offset != state_.size())
  261. {
  262. boost::throw_exception(std::invalid_argument("Not enough elements in call to seed."));
  263. }
  264. }
  265. /**
  266. * Constructs a @c xoshiro and calls @c seed().
  267. */
  268. xoshiro_base() { seed(); }
  269. /** Seeds the generator with a user provided seed. */
  270. explicit xoshiro_base(const seed_type value)
  271. {
  272. seed(value);
  273. }
  274. template <typename FIter>
  275. xoshiro_base(FIter& first, FIter last) { seed(first, last); }
  276. /**
  277. * Seeds the generator with 64-bit values produced by @c seq.generate().
  278. *
  279. * @xmlnote
  280. * The copy constructor will always be preferred over
  281. * the templated constructor.
  282. * @endxmlnote
  283. */
  284. template <typename Sseq, typename std::enable_if<!std::is_convertible<Sseq, xoshiro_base>::value, bool>::type = true>
  285. explicit xoshiro_base(Sseq& seq)
  286. {
  287. seed(seq);
  288. }
  289. // Hit all of our rule of 5 explicitly to ensure old platforms work correctly
  290. ~xoshiro_base() = default;
  291. xoshiro_base(const xoshiro_base& other) noexcept { state_ = other.state(); }
  292. xoshiro_base& operator=(const xoshiro_base& other) noexcept { state_ = other.state(); return *this; }
  293. xoshiro_base(xoshiro_base&& other) noexcept { state_ = other.state(); }
  294. xoshiro_base& operator=(xoshiro_base&& other) noexcept { state_ = other.state(); return *this; }
  295. inline result_type next() noexcept
  296. {
  297. return static_cast<Derived*>(this)->next();
  298. }
  299. /** This is the jump function for the generator. It is equivalent
  300. * to 2^128 calls to next(); it can be used to generate 2^128
  301. * non-overlapping subsequences for parallel computations. */
  302. inline void jump() noexcept
  303. {
  304. jump_impl(xoshiro_type());
  305. }
  306. /** This is the long-jump function for the generator. It is equivalent to
  307. * 2^192 calls to next(); it can be used to generate 2^64 starting points,
  308. * from each of which jump() will generate 2^64 non-overlapping
  309. * subsequences for parallel distributed computations. */
  310. inline void long_jump() noexcept
  311. {
  312. long_jump_impl(xoshiro_type());
  313. }
  314. /** Returns the next value of the generator. */
  315. inline result_type operator()() noexcept
  316. {
  317. return next();
  318. }
  319. /** Advances the state of the generator by @c z. */
  320. inline void discard(const std::uint64_t z) noexcept
  321. {
  322. for (std::uint64_t i {}; i < z; ++i)
  323. {
  324. next();
  325. }
  326. }
  327. /**
  328. * Returns true if the two generators will produce identical
  329. * sequences of values.
  330. */
  331. inline friend bool operator==(const xoshiro_base& lhs, const xoshiro_base& rhs) noexcept
  332. {
  333. return lhs.state_ == rhs.state_;
  334. }
  335. /**
  336. * Returns true if the two generators will produce different
  337. * sequences of values.
  338. */
  339. inline friend bool operator!=(const xoshiro_base& lhs, const xoshiro_base& rhs) noexcept
  340. {
  341. return lhs.state_ != rhs.state_;
  342. }
  343. /** Writes a @c xorshiro to a @c std::ostream. */
  344. template <typename CharT, typename Traits>
  345. inline friend std::basic_ostream<CharT, Traits>& operator<<(std::basic_ostream<CharT, Traits>& ost,
  346. const xoshiro_base& e)
  347. {
  348. for (std::size_t i {}; i < e.state_.size(); ++i)
  349. {
  350. ost << e.state_[i] << ' ';
  351. }
  352. return ost;
  353. }
  354. /** Reads a @c xorshiro from a @c std::istream. */
  355. template <typename CharT, typename Traits>
  356. inline friend std::basic_istream<CharT, Traits>& operator>>(std::basic_istream<CharT, Traits>& ist,
  357. xoshiro_base& e)
  358. {
  359. for (std::size_t i {}; i < e.state_.size(); ++i)
  360. {
  361. ist >> e.state_[i] >> std::ws;
  362. }
  363. return ist;
  364. }
  365. /** Fills a range with random values */
  366. template <typename FIter>
  367. inline void generate(FIter first, FIter last) noexcept
  368. {
  369. using iter_type = typename std::iterator_traits<FIter>::value_type;
  370. while (first != last)
  371. {
  372. *first++ = static_cast<iter_type>(next());
  373. }
  374. }
  375. /**
  376. * Returns the largest value that the @c xorshiro
  377. * can produce.
  378. */
  379. static constexpr result_type (max)() noexcept
  380. {
  381. return (std::numeric_limits<result_type>::max)();
  382. }
  383. /**
  384. * Returns the smallest value that the @c xorshiro
  385. * can produce.
  386. */
  387. static constexpr result_type (min)() noexcept
  388. {
  389. return (std::numeric_limits<result_type>::min)();
  390. }
  391. inline std::array<BlockType, N> state() const noexcept
  392. {
  393. return state_;
  394. }
  395. };
  396. } // namespace detail
  397. } // namespace random
  398. } // namespace boost
  399. #endif