Advertisement
ulfben

xoshiro256 with public interface

Nov 19th, 2023 (edited)
1,385
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 7.34 KB | None | 0 0
  1. #pragma once
  2. #include <array>
  3. #include <cassert>
  4. #include <bitset>
  5. #include <chrono> //for entropy in createSeeds
  6. #include <concepts>
  7. #include <cstdint>
  8. #include <limits>
  9. #include <span>
  10. #include <thread> //for entropy in createSeeds
  11. #include <ctime> //for entropy in createSeeds
  12. #include <type_traits>
  13.  
  14. // The "xoshiro256** 1.0" generator.
  15. // public interface, rejection sampling and seeding utilities by Ulf Benjaminsson (2023)
  16. // https://ulfbenjaminsson.com/
  17. // Based on C++ port by Arthur O'Dwyer (2021).
  18. // https://quuxplusone.github.io/blog/2021/11/23/xoshiro/
  19. // of the C version by David Blackman and Sebastiano Vigna (2018),
  20. // https://prng.di.unimi.it/xoshiro256starstar.c
  21.  
  22. class RNG{
  23. public:
  24.     using u64 = std::uint_fast64_t;
  25.     static constexpr std::size_t SEED_COUNT = 4;
  26.     using State = std::array<u64, SEED_COUNT>;
  27.     using Span = std::span<const u64, SEED_COUNT>;
  28.     static constexpr auto USE_REJECTION_SAMPLING = false;
  29.     // When enabled, all ranged functions will use rejection sampling to ensure a more uniform
  30.     // distribution of random numbers across large ranges. The concern with large ranges is that
  31.     // methods like modulo reduction (randNum % range) might not evenly distribute numbers
  32.     // across the range, leading to bias. As a rule of thumb, if the range is more than half of
  33.     // the maximum output (e.g., more than 2^63), you might start to see the benefits of using
  34.     // rejection sampling to ensure uniformity.
  35.  
  36.     constexpr explicit RNG(u64 seed) noexcept{
  37.         s[0] = splitmix64(seed);
  38.         seed += 0x9E3779B97F4A7C15uLL;
  39.         s[1] = splitmix64(seed);
  40.         seed += 0x7F4A7C15uLL;
  41.         s[2] = splitmix64(seed);
  42.         seed += 0x9E3779B9uLL;
  43.         s[3] = splitmix64(seed);
  44.     }
  45.     constexpr explicit RNG(double seed) noexcept
  46.         : RNG(static_cast<u64>(seed)){}
  47.  
  48.     constexpr explicit RNG(Span seeds) noexcept{
  49.         std::ranges::copy(seeds, s.begin());
  50.     }
  51.  
  52.     static constexpr u64 max() noexcept{
  53.         return std::numeric_limits<u64>::max();
  54.     }
  55.  
  56.     constexpr u64 next() noexcept{
  57.         return nextU64();
  58.     }
  59.  
  60.     constexpr bool coinToss() noexcept{
  61.         return next() & 1; //checks the least significant bit
  62.     }
  63.  
  64.     template<std::floating_point Real>
  65.     constexpr Real normalized() noexcept{
  66.         return static_cast<Real>(next()) / static_cast<Real>(max());
  67.     }
  68.  
  69.     template<std::floating_point Real>
  70.     constexpr Real inRange(Real range) noexcept{
  71.         return normalized<Real>() * range;
  72.     }
  73.  
  74.     template<std::floating_point Real>
  75.     constexpr Real inRange(Real from, Real to) noexcept{
  76.         assert(from < to && "RNG: inverted range.");
  77.         return from + normalized<Real>() * (to - from);
  78.     }
  79.  
  80.     template<std::integral T>
  81.     constexpr T inRange(T from, T to) noexcept{
  82.         assert(from < to && "RNG: inverted range.");
  83.         using UT = std::make_unsigned_t<T>;
  84.         UT range = static_cast<UT>(to - from);
  85.         return static_cast<T>(inRange(range)) + from;
  86.     }
  87.  
  88.     template<std::integral T>
  89.     constexpr T inRange(T range) noexcept{
  90.         using UT = std::make_unsigned_t<T>;
  91.         UT num = static_cast<UT>(inRange(static_cast<u64>(std::abs(range))));
  92.         return (range < 0) ? -static_cast<T>(num) : static_cast<T>(num);
  93.     }
  94.  
  95.     constexpr u64 inRange(u64 range) noexcept{
  96.         if(range == 0){
  97.             assert(false && "RNG::inRange called with empty range!");
  98.             return 0;
  99.         }
  100.         if constexpr(USE_REJECTION_SAMPLING){
  101.             return uniformRandom(range);
  102.         }
  103.         return next() / (max() / range);
  104.     }
  105.  
  106.     constexpr Span state() const{
  107.         return {s};
  108.     }
  109.  
  110.     /* the jump() function is equivalent to 2^128 calls to next();
  111.     it can be used to generate 2^128 non-overlapping subsequences
  112.     for parallel computations. */
  113.     constexpr void jump() noexcept{
  114.         constexpr std::array<std::bitset<64>, SEED_COUNT> JUMP{
  115.             0x180ec6d33cfd0abaULL, 0xd5a61266f0c9392cULL, 0xa9582618e03fc9aaULL, 0x39abdc4529b1661cULL
  116.         };
  117.         State temp{0};
  118.         for(const auto& bits : JUMP){
  119.             for(std::size_t b = 0; b < 64; ++b){
  120.                 if(bits.test(b)){
  121.                     temp[0] ^= s[0];
  122.                     temp[1] ^= s[1];
  123.                     temp[2] ^= s[2];
  124.                     temp[3] ^= s[3];
  125.                 }
  126.                 next();
  127.             }
  128.         }
  129.         s = temp;
  130.     }
  131.  
  132.     //convenience function to allow hashing from const values.
  133.     static constexpr u64 splitmix64_hash(u64 x) noexcept{        
  134.         return splitmix64(x);
  135.     }
  136.    
  137. private:
  138.     State s{};
  139.  
  140.     constexpr u64 nextU64() noexcept{
  141.         const u64 result = rotl(s[1] * 5, 7) * 9;
  142.         const u64 t = s[1] << 17;
  143.         s[2] ^= s[0];
  144.         s[3] ^= s[1];
  145.         s[1] ^= s[2];
  146.         s[0] ^= s[3];
  147.         s[2] ^= t;
  148.         s[3] = rotl(s[3], 45);
  149.         return result;
  150.     }
  151.  
  152.     static constexpr u64 rotl(u64 x, int k) noexcept{
  153.         return (x << k) | (x >> (64 - k));
  154.     }
  155.  
  156.     static constexpr u64 splitmix64(u64& x) noexcept{
  157.         u64 z = (x += 0x9e3779b97f4a7c15uLL);
  158.         z = (z ^ (z >> 30)) * 0xbf58476d1ce4e5b9uLL;
  159.         z = (z ^ (z >> 27)) * 0x94d049bb133111ebuLL;
  160.         return z ^ (z >> 31);
  161.     }
  162.  
  163.     //uses rejection sampling to ensure fair scaling. Costlier than inRange(u64)
  164.     constexpr u64 uniformRandom(u64 range) noexcept{
  165.         if(range == 0){
  166.             assert(false && "RNG::uniformRandom called with empty range!");
  167.             return 0;
  168.         }
  169.         const u64 rangeLimit = max() - range + 1;
  170.         u64 n = next();
  171.         while((n - (n % range)) >= rangeLimit){
  172.             n = next();
  173.         }
  174.         return n % range;
  175.     }
  176. };
  177.  
  178. //Some strategies for seeding the full 256-bit state of xoshiro256.
  179. // createSeeds uses date, time-since-launch, CPU time, thread ID, and a memory address as sources of entropy.
  180. // each value is hashed using splitmix64.
  181. static typename RNG::State createSeeds() noexcept{
  182.     using u64 = RNG::u64;
  183.     using namespace std::chrono;
  184.     const auto current_date = static_cast<u64>(system_clock::now().time_since_epoch().count());
  185.     const auto uptime = static_cast<u64>(high_resolution_clock::now().time_since_epoch().count());
  186.     const auto cpu_time = static_cast<u64>(std::clock());    
  187.     const auto mixed_time = static_cast<u64>((uptime << 1) ^ current_date);
  188.     const auto thread_id = static_cast<u64>(std::hash<std::thread::id>{}(std::this_thread::get_id()));    
  189.     const int local{};
  190.     const auto local_addr = static_cast<u64>(reinterpret_cast<std::uintptr_t>(&local));
  191.     return {
  192.         RNG::splitmix64_hash(mixed_time),
  193.         RNG::splitmix64_hash(thread_id),
  194.         RNG::splitmix64_hash(cpu_time),
  195.         RNG::splitmix64_hash(local_addr),
  196.     };
  197. }
  198.  
  199. /* createSeedsB *additionally* mix those sources with entropy from std::random_device. This is likely significantly slow(er).
  200. #include <random>
  201. static typename RNG::State createSeedsB() {
  202.     using u64 = RNG::u64;
  203.     std::random_device rd;
  204.     auto seeds = createSeeds();
  205.     seeds[0] = static_cast<u64>(rd()) ^ seeds[0];
  206.     seeds[1] = static_cast<u64>(rd()) ^ (seeds[0] << 1);
  207.     seeds[2] = static_cast<u64>(rd()) ^ (seeds[1] << 1);
  208.     seeds[3] = static_cast<u64>(rd()) ^ (seeds[2] << 1);
  209.     return seeds;
  210. }*/
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement