random.hpp 20 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497
  1. /// \file
  2. // Range v3 library
  3. //
  4. // Copyright Casey Carter 2016
  5. //
  6. // Use, modification and distribution is subject to the
  7. // Boost Software License, Version 1.0. (See accompanying
  8. // file LICENSE_1_0.txt or copy at
  9. // http://www.boost.org/LICENSE_1_0.txt)
  10. //
  11. // Project home: https://github.com/ericniebler/range-v3
  12. //
  13. /*
  14. * Random-Number Utilities (randutil)
  15. * Addresses common issues with C++11 random number generation.
  16. * Makes good seeding easier, and makes using RNGs easy while retaining
  17. * all the power.
  18. *
  19. * The MIT License (MIT)
  20. *
  21. * Copyright (c) 2015 Melissa E. O'Neill
  22. *
  23. * Permission is hereby granted, free of charge, to any person obtaining a copy
  24. * of this software and associated documentation files (the "Software"), to deal
  25. * in the Software without restriction, including without limitation the rights
  26. * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
  27. * copies of the Software, and to permit persons to whom the Software is
  28. * furnished to do so, subject to the following conditions:
  29. *
  30. * The above copyright notice and this permission notice shall be included in
  31. * all copies or substantial portions of the Software.
  32. *
  33. * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  34. * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  35. * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
  36. * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  37. * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
  38. * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
  39. * SOFTWARE.
  40. */
  41. #ifndef RANGES_V3_UTILITY_RANDOM_HPP
  42. #define RANGES_V3_UTILITY_RANDOM_HPP
  43. #include <array>
  44. #include <cstddef>
  45. #include <cstdint>
  46. #include <initializer_list>
  47. #include <new>
  48. #include <random>
  49. #include <meta/meta.hpp>
  50. #include <concepts/concepts.hpp>
  51. #include <range/v3/range_fwd.hpp>
  52. #include <range/v3/algorithm/copy.hpp>
  53. #include <range/v3/algorithm/generate.hpp>
  54. #include <range/v3/functional/invoke.hpp>
  55. #include <range/v3/functional/reference_wrapper.hpp>
  56. #include <range/v3/iterator/concepts.hpp>
  57. #if !RANGES_CXX_THREAD_LOCAL
  58. #include <mutex>
  59. #endif
  60. #include <range/v3/detail/prologue.hpp>
  61. RANGES_DIAGNOSTIC_PUSH
  62. RANGES_DIAGNOSTIC_IGNORE_CXX17_COMPAT
  63. namespace ranges
  64. {
  65. /// \addtogroup group-numerics
  66. /// @{
  67. // clang-format off
  68. /// \concept uniform_random_bit_generator_
  69. /// \brief The \c uniform_random_bit_generator_ concept
  70. template<typename Gen>
  71. CPP_requires(uniform_random_bit_generator_,
  72. requires() //
  73. (
  74. Gen::min(),
  75. Gen::max()
  76. ));
  77. /// \concept uniform_random_bit_generator_
  78. /// \brief The \c uniform_random_bit_generator_ concept
  79. template(typename Gen)(
  80. concept (uniform_random_bit_generator_)(Gen),
  81. unsigned_integral<invoke_result_t<Gen &>> AND
  82. same_as<invoke_result_t<Gen &>, decltype(Gen::min())> AND
  83. same_as<invoke_result_t<Gen &>, decltype(Gen::max())>);
  84. /// \concept uniform_random_bit_generator
  85. /// \brief The \c uniform_random_bit_generator concept
  86. template<typename Gen>
  87. CPP_concept uniform_random_bit_generator =
  88. invocable<Gen &> &&
  89. CPP_requires_ref(ranges::uniform_random_bit_generator_, Gen) &&
  90. CPP_concept_ref(ranges::uniform_random_bit_generator_, Gen);
  91. // clang-format on
  92. /// @}
  93. /// \cond
  94. namespace detail
  95. {
  96. namespace randutils
  97. {
  98. inline std::array<std::uint32_t, 8> get_entropy()
  99. {
  100. std::array<std::uint32_t, 8> seeds;
  101. // Hopefully high-quality entropy from random_device.
  102. #if defined(__GLIBCXX__) && defined(RANGES_WORKAROUND_VALGRIND_RDRAND)
  103. std::random_device rd{"/dev/urandom"};
  104. #else
  105. std::random_device rd;
  106. #endif
  107. std::uniform_int_distribution<std::uint32_t> dist{};
  108. ranges::generate(seeds, [&] { return dist(rd); });
  109. return seeds;
  110. }
  111. template(typename I)(
  112. requires unsigned_integral<I>)
  113. constexpr I fast_exp(I x, I power, I result = I{1})
  114. {
  115. return power == I{0}
  116. ? result
  117. : randutils::fast_exp(
  118. x * x, power >> 1, result * (power & I{1} ? x : 1));
  119. }
  120. //////////////////////////////////////////////////////////////////////////////
  121. //
  122. // seed_seq_fe
  123. //
  124. //////////////////////////////////////////////////////////////////////////////
  125. /*
  126. * seed_seq_fe implements a fixed-entropy seed sequence; it conforms to all
  127. * the requirements of a Seed Sequence concept.
  128. *
  129. * seed_seq_fe<N> implements a seed sequence which seeds based on a store of
  130. * N * 32 bits of entropy. Typically, it would be initialized with N or more
  131. * integers.
  132. *
  133. * seed_seq_fe128 and seed_seq_fe256 are provided as convenience typedefs for
  134. * 128- and 256-bit entropy stores respectively. These variants outperform
  135. * std::seed_seq, while being better mixing the bits it is provided as
  136. * entropy. In almost all common use cases, they serve as better drop-in
  137. * replacements for seed_seq.
  138. *
  139. * Technical details
  140. *
  141. * Assuming it constructed with M seed integers as input, it exhibits the
  142. * following properties
  143. *
  144. * * Diffusion/Avalanche: A single-bit change in any of the M inputs has a
  145. * 50% chance of flipping every bit in the bitstream produced by generate.
  146. * Initializing the N-word entropy store with M words requires O(N * M)
  147. * time precisely because of the avalanche requirements. Once constructed,
  148. * calls to generate are linear in the number of words generated.
  149. *
  150. * * Bias freedom/Bijection: If M == N, the state of the entropy store is a
  151. * bijection from the M inputs (i.e., no states occur twice, none are
  152. * omitted). If M > N the number of times each state can occur is the same
  153. * (each state occurs 2**(32*(M-N)) times, where ** is the power function).
  154. * If M < N, some states cannot occur (bias) but no state occurs more
  155. * than once (it's impossible to avoid bias if M < N; ideally N should not
  156. * be chosen so that it is more than M).
  157. *
  158. * Likewise, the generate function has similar properties (with the entropy
  159. * store as the input data). If more outputs are requested than there is
  160. * entropy, some outputs cannot occur. For example, the Mersenne Twister
  161. * will request 624 outputs, to initialize its 19937-bit state, which is
  162. * much larger than a 128-bit or 256-bit entropy pool. But in practice,
  163. * limiting the Mersenne Twister to 2**128 possible initializations gives
  164. * us enough initializations to give a unique initialization to trillions
  165. * of computers for billions of years. If you really have 624 words of
  166. * *real* high-quality entropy you want to use, you probably don't need
  167. * an entropy mixer like this class at all. But if you *really* want to,
  168. * nothing is stopping you from creating a randutils::seed_seq_fe<624>.
  169. *
  170. * * As a consequence of the above properties, if all parts of the provided
  171. * seed data are kept constant except one, and the remaining part is varied
  172. * through K different states, K different output sequences will be
  173. * produced.
  174. *
  175. * * Also, because the amount of entropy stored is fixed, this class never
  176. * performs dynamic allocation and is free of the possibility of generating
  177. * an exception.
  178. *
  179. * Ideas used to implement this code include hashing, a simple PCG generator
  180. * based on an MCG base with an XorShift output function and permutation
  181. * functions on tuples.
  182. *
  183. * More detail at
  184. * http://www.pcg-random.org/posts/developing-a-seed_seq-alternative.html
  185. */
  186. template<std::size_t count, typename IntRep = std::uint32_t>
  187. struct seed_seq_fe
  188. {
  189. public:
  190. CPP_assert(unsigned_integral<IntRep>);
  191. typedef IntRep result_type;
  192. private:
  193. static constexpr std::size_t mix_rounds = 1 + (count <= 2);
  194. static constexpr std::uint32_t INIT_A = 0x43b0d7e5;
  195. static constexpr std::uint32_t MULT_A = 0x931e8875;
  196. static constexpr std::uint32_t INIT_B = 0x8b51f9dd;
  197. static constexpr std::uint32_t MULT_B = 0x58f38ded;
  198. static constexpr std::uint32_t MIX_MULT_L = 0xca01f9dd;
  199. static constexpr std::uint32_t MIX_MULT_R = 0x4973f715;
  200. static constexpr std::uint32_t XSHIFT = sizeof(IntRep) * 8 / 2;
  201. std::array<IntRep, count> mixer_;
  202. template(typename I, typename S)(
  203. requires input_iterator<I> AND sentinel_for<S, I> AND
  204. convertible_to<iter_reference_t<I>, IntRep>)
  205. void mix_entropy(I first, S last)
  206. {
  207. auto hash_const = INIT_A;
  208. auto hash = [&](IntRep value) RANGES_INTENDED_MODULAR_ARITHMETIC {
  209. value ^= hash_const;
  210. hash_const *= MULT_A;
  211. value *= hash_const;
  212. value ^= value >> XSHIFT;
  213. return value;
  214. };
  215. auto mix = [](IntRep x, IntRep y) RANGES_INTENDED_MODULAR_ARITHMETIC {
  216. IntRep result = MIX_MULT_L * x - MIX_MULT_R * y;
  217. result ^= result >> XSHIFT;
  218. return result;
  219. };
  220. for(auto & elem : mixer_)
  221. {
  222. if(first != last)
  223. {
  224. elem = hash(static_cast<IntRep>(*first));
  225. ++first;
  226. }
  227. else
  228. elem = hash(IntRep{0});
  229. }
  230. for(auto & src : mixer_)
  231. for(auto & dest : mixer_)
  232. if(&src != &dest)
  233. dest = mix(dest, hash(src));
  234. for(; first != last; ++first)
  235. for(auto & dest : mixer_)
  236. dest = mix(dest, hash(static_cast<IntRep>(*first)));
  237. }
  238. public:
  239. seed_seq_fe(const seed_seq_fe &) = delete;
  240. void operator=(const seed_seq_fe &) = delete;
  241. template(typename T)(
  242. requires convertible_to<T const &, IntRep>)
  243. seed_seq_fe(std::initializer_list<T> init)
  244. {
  245. seed(init.begin(), init.end());
  246. }
  247. template(typename I, typename S)(
  248. requires input_iterator<I> AND sentinel_for<S, I> AND
  249. convertible_to<iter_reference_t<I>, IntRep>)
  250. seed_seq_fe(I first, S last)
  251. {
  252. seed(first, last);
  253. }
  254. // generating functions
  255. template(typename I, typename S)(
  256. requires random_access_iterator<I> AND sentinel_for<S, I>)
  257. RANGES_INTENDED_MODULAR_ARITHMETIC //
  258. void generate(I first, S const last) const
  259. {
  260. auto src_begin = mixer_.begin();
  261. auto src_end = mixer_.end();
  262. auto src = src_begin;
  263. auto hash_const = INIT_B;
  264. for(; first != last; ++first)
  265. {
  266. auto dataval = *src;
  267. if(++src == src_end)
  268. src = src_begin;
  269. dataval ^= hash_const;
  270. hash_const *= MULT_B;
  271. dataval *= hash_const;
  272. dataval ^= dataval >> XSHIFT;
  273. *first = dataval;
  274. }
  275. }
  276. constexpr std::size_t size() const
  277. {
  278. return count;
  279. }
  280. template(typename O)(
  281. requires weakly_incrementable<O> AND
  282. indirectly_copyable<decltype(mixer_.begin()), O>)
  283. RANGES_INTENDED_MODULAR_ARITHMETIC void param(O dest) const
  284. {
  285. constexpr IntRep INV_A = randutils::fast_exp(MULT_A, IntRep(-1));
  286. constexpr IntRep MIX_INV_L =
  287. randutils::fast_exp(MIX_MULT_L, IntRep(-1));
  288. auto mixer_copy = mixer_;
  289. for(std::size_t round = 0; round < mix_rounds; ++round)
  290. {
  291. // Advance to the final value. We'll backtrack from that.
  292. auto hash_const =
  293. INIT_A * randutils::fast_exp(MULT_A, IntRep(count * count));
  294. for(auto src = mixer_copy.rbegin(); src != mixer_copy.rend();
  295. ++src)
  296. for(auto rdest = mixer_copy.rbegin();
  297. rdest != mixer_copy.rend();
  298. ++rdest)
  299. if(src != rdest)
  300. {
  301. IntRep revhashed = *src;
  302. auto mult_const = hash_const;
  303. hash_const *= INV_A;
  304. revhashed ^= hash_const;
  305. revhashed *= mult_const;
  306. revhashed ^= revhashed >> XSHIFT;
  307. IntRep unmixed = *rdest;
  308. unmixed ^= unmixed >> XSHIFT;
  309. unmixed += MIX_MULT_R * revhashed;
  310. unmixed *= MIX_INV_L;
  311. *rdest = unmixed;
  312. }
  313. for(auto i = mixer_copy.rbegin(); i != mixer_copy.rend(); ++i)
  314. {
  315. IntRep unhashed = *i;
  316. unhashed ^= unhashed >> XSHIFT;
  317. unhashed *= randutils::fast_exp(hash_const, IntRep(-1));
  318. hash_const *= INV_A;
  319. unhashed ^= hash_const;
  320. *i = unhashed;
  321. }
  322. }
  323. ranges::copy(mixer_copy, dest);
  324. }
  325. template(typename I, typename S)(
  326. requires input_iterator<I> AND sentinel_for<S, I> AND
  327. convertible_to<iter_reference_t<I>, IntRep>)
  328. void seed(I first, S last)
  329. {
  330. mix_entropy(first, last);
  331. // For very small sizes, we do some additional mixing. For normal
  332. // sizes, this loop never performs any iterations.
  333. for(std::size_t i = 1; i < mix_rounds; ++i)
  334. stir();
  335. }
  336. seed_seq_fe & stir()
  337. {
  338. mix_entropy(mixer_.begin(), mixer_.end());
  339. return *this;
  340. }
  341. };
  342. using seed_seq_fe128 = seed_seq_fe<4, std::uint32_t>;
  343. using seed_seq_fe256 = seed_seq_fe<8, std::uint32_t>;
  344. //////////////////////////////////////////////////////////////////////////////
  345. //
  346. // auto_seeded
  347. //
  348. //////////////////////////////////////////////////////////////////////////////
  349. /*
  350. * randutils::auto_seeded
  351. *
  352. * Extends a seed sequence class with a nondeterministic default
  353. * constructor. Uses a variety of local sources of entropy to portably
  354. * initialize any seed sequence to a good default state.
  355. *
  356. * In normal use, it's accessed via one of the following type aliases, which
  357. * use seed_seq_fe128 and seed_seq_fe256 above.
  358. *
  359. * randutils::auto_seed_128
  360. * randutils::auto_seed_256
  361. *
  362. * It's discussed in detail at
  363. * http://www.pcg-random.org/posts/simple-portable-cpp-seed-entropy.html
  364. * and its motivation (why you can't just use std::random_device) here
  365. * http://www.pcg-random.org/posts/cpps-random_device.html
  366. */
  367. template<typename SeedSeq>
  368. struct auto_seeded : public SeedSeq
  369. {
  370. auto_seeded()
  371. : auto_seeded(randutils::get_entropy())
  372. {}
  373. template<std::size_t N>
  374. auto_seeded(std::array<std::uint32_t, N> const & seeds)
  375. : SeedSeq(seeds.begin(), seeds.end())
  376. {}
  377. using SeedSeq::SeedSeq;
  378. const SeedSeq & base() const
  379. {
  380. return *this;
  381. }
  382. SeedSeq & base()
  383. {
  384. return *this;
  385. }
  386. };
  387. using auto_seed_128 = auto_seeded<seed_seq_fe128>;
  388. using auto_seed_256 = auto_seeded<seed_seq_fe256>;
  389. } // namespace randutils
  390. using default_URNG = meta::if_c<(sizeof(void *) >= sizeof(long long)),
  391. std::mt19937_64, std::mt19937>;
  392. #if !RANGES_CXX_THREAD_LOCAL
  393. template<typename URNG>
  394. class sync_URNG : private URNG
  395. {
  396. mutable std::mutex mtx_;
  397. public:
  398. using URNG::URNG;
  399. sync_URNG() = default;
  400. using typename URNG::result_type;
  401. result_type operator()()
  402. {
  403. std::lock_guard<std::mutex> guard{mtx_};
  404. return static_cast<URNG &>(*this)();
  405. }
  406. using URNG::max;
  407. using URNG::min;
  408. };
  409. using default_random_engine = sync_URNG<default_URNG>;
  410. #else
  411. using default_random_engine = default_URNG;
  412. #endif
  413. template<typename T = void>
  414. default_random_engine & get_random_engine()
  415. {
  416. using Seeder = meta::if_c<(sizeof(default_URNG) > 16),
  417. randutils::auto_seed_256,
  418. randutils::auto_seed_128>;
  419. #if RANGES_CXX_THREAD_LOCAL >= RANGES_CXX_THREAD_LOCAL_11
  420. static thread_local default_random_engine engine{Seeder{}.base()};
  421. #elif RANGES_CXX_THREAD_LOCAL
  422. static __thread bool initialized = false;
  423. static __thread meta::_t<std::aligned_storage<sizeof(default_random_engine),
  424. alignof(default_random_engine)>>
  425. storage;
  426. if(!initialized)
  427. {
  428. ::new(static_cast<void *>(&storage))
  429. default_random_engine{Seeder{}.base()};
  430. initialized = true;
  431. }
  432. auto & engine = reinterpret_cast<default_random_engine &>(storage);
  433. #else
  434. static default_random_engine engine{Seeder{}.base()};
  435. #endif // RANGES_CXX_THREAD_LOCAL
  436. return engine;
  437. }
  438. } // namespace detail
  439. /// \endcond
  440. } // namespace ranges
  441. RANGES_DIAGNOSTIC_POP
  442. #include <range/v3/detail/epilogue.hpp>
  443. #endif