diff --git a/src/libraries/util/include/util/xoroshiro.h b/src/libraries/util/include/util/xoroshiro.h new file mode 100644 index 0000000..f37ab36 --- /dev/null +++ b/src/libraries/util/include/util/xoroshiro.h @@ -0,0 +1,55 @@ +#pragma once +/// \file xoroshiro.h +/// xoroshiro256++ random number implementation + +#include + +namespace util +{ + +class xoroshiro256pp +{ +public: + static constexpr size_t state_words = 4; + + /// Constructor. + /// \arg seed A non-zero seed for the generator + xoroshiro256pp(uint64_t seed); + + /// Constructor with explicit state. + xoroshiro256pp(const uint64_t state[state_words]); + + /// Copy constructor + xoroshiro256pp(const xoroshiro256pp &other); + + /// Assignment operator + xoroshiro256pp & operator=(const xoroshiro256pp &other); + + /// Get the next random integer from the generator + /// \returns A random integer between 0-2^64, inclusive + uint64_t next(); + + /// 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. + /// \returns A new generator 2^128 calls ahead of this one + inline xoroshiro256pp jump() const { return do_jump(jump_poly); } + + /// 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 subsequences for + /// parallel distributed computations. + /// \returns A new generator 2^192 calls ahead of this one + xoroshiro256pp long_jump() const { return do_jump(longjump_poly); } + +private: + static const uint64_t jump_poly[state_words]; + static const uint64_t longjump_poly[state_words]; + + xoroshiro256pp do_jump(const uint64_t poly[state_words]) const; + + uint64_t m_state[state_words]; +}; + +} // namespace util \ No newline at end of file diff --git a/src/libraries/util/util.module b/src/libraries/util/util.module index 9a002ab..7c58866 100644 --- a/src/libraries/util/util.module +++ b/src/libraries/util/util.module @@ -8,6 +8,7 @@ module("util", "bip_buffer.cpp", "format.cpp", "spinlock.cpp", + "xoroshiro.cpp", ], public_headers = [ "util/allocator.h", @@ -33,4 +34,5 @@ module("util", "util/spinlock.h", "util/util.h", "util/vector.h", + "util/xoroshiro.h", ]) diff --git a/src/libraries/util/xoroshiro.cpp b/src/libraries/util/xoroshiro.cpp new file mode 100644 index 0000000..db8b799 --- /dev/null +++ b/src/libraries/util/xoroshiro.cpp @@ -0,0 +1,103 @@ +// PRNG based on xoshiro256++ 1.0 by David Blackman and Sebastiano Vigna. +// https://prng.di.unimi.it/xoshiro256plusplus.c + +/* This is xoshiro256++ 1.0, an all-purpose, rock-solid generator by + David Blackman and Sebastiano Vigna. + + It has excellent (sub-ns) speed, a state (256 bits) that is large + enough for any parallel application, and it passes all tests we are + aware of. + + The state must be seeded so that it is not everywhere zero. If you have + a 64-bit seed, we suggest to seed a splitmix64 generator and use its + output to fill the state. */ + + +#include +#include + +namespace util { + +const uint64_t xoroshiro256pp::jump_poly[] = { 0x180ec6d33cfd0aba, 0xd5a61266f0c9392c, 0xa9582618e03fc9aa, 0x39abdc4529b1661c }; +const uint64_t xoroshiro256pp::longjump_poly[] = { 0x76e15d3efefdcbbf, 0xc5004e441c522fb3, 0x77710069854ee241, 0x39109bb02acbe635 }; + +static inline uint64_t rotl(const uint64_t x, int k) { + return (x << k) | (x >> (64 - k)); +} + +xoroshiro256pp::xoroshiro256pp(uint64_t seed) +{ + m_state[0] = seed; + m_state[1] = splitmix64(m_state[0]); + m_state[2] = splitmix64(m_state[1]); + m_state[3] = splitmix64(m_state[2]); +} + +xoroshiro256pp::xoroshiro256pp(const uint64_t state[state_words]) : + m_state {state[0], state[1], state[2], state[3]} +{ +} + +xoroshiro256pp::xoroshiro256pp(const xoroshiro256pp &other) : + m_state {other.m_state[0], other.m_state[1], other.m_state[2], other.m_state[3]} +{ +} + +xoroshiro256pp & +xoroshiro256pp::operator=(const xoroshiro256pp &other) +{ + m_state[0] = other.m_state[0]; + m_state[1] = other.m_state[1]; + m_state[2] = other.m_state[2]; + m_state[3] = other.m_state[3]; + return *this; +} + +uint64_t +xoroshiro256pp::next() +{ + const uint64_t result = rotl(m_state[0] + m_state[3], 23) + m_state[0]; + + const uint64_t t = m_state[1] << 17; + + m_state[2] ^= m_state[0]; + m_state[3] ^= m_state[1]; + m_state[1] ^= m_state[2]; + m_state[0] ^= m_state[3]; + + m_state[2] ^= t; + + m_state[3] = rotl(m_state[3], 45); + + return result; +} + +xoroshiro256pp +xoroshiro256pp::do_jump(const uint64_t poly[state_words]) const +{ + uint64_t s0 = 0; + uint64_t s1 = 0; + uint64_t s2 = 0; + uint64_t s3 = 0; + xoroshiro256pp x = *this; + + for(int i = 0; i < state_words; ++i) { + for(int b = 0; b < 64; ++b) { + if (poly[i] & (1ull << b)) { + s0 ^= x.m_state[0]; + s1 ^= x.m_state[1]; + s2 ^= x.m_state[2]; + s3 ^= x.m_state[3]; + } + x.next(); + } + } + + x.m_state[0] = s0; + x.m_state[1] = s1; + x.m_state[2] = s2; + x.m_state[3] = s3; + return x; +} + +} // namespace util \ No newline at end of file