429 lines
12 KiB
C++
429 lines
12 KiB
C++
/* Code was taken from here:
|
|
https://prng.di.unimi.it
|
|
And wrapped as a bit generator for C++'s random module' */
|
|
|
|
/* Written in 2019 by David Blackman and Sebastiano Vigna (vigna@acm.org)
|
|
|
|
To the extent possible under law, the author has dedicated all copyright
|
|
and related and neighboring rights to this software to the public domain
|
|
worldwide. This software is distributed without any warranty.
|
|
|
|
See <http://creativecommons.org/publicdomain/zero/1.0/>. */
|
|
#ifndef XOSHIRO_H
|
|
#define XOSHIRO_H
|
|
|
|
#include <iostream>
|
|
#include <cstring>
|
|
#include <cstdint>
|
|
|
|
namespace Xoshiro {
|
|
|
|
|
|
#if __cplusplus >= 202001L
|
|
# define HAS_CPP20
|
|
#endif
|
|
|
|
static inline std::uint64_t rotl64(const std::uint64_t x, const int k) {
|
|
return (x << k) | (x >> (64 - k));
|
|
}
|
|
|
|
static inline std::uint32_t rotl32(const std::uint32_t x, const int k) {
|
|
return (x << k) | (x >> (32 - k));
|
|
}
|
|
|
|
/* these are in order to avoid gcc warnings about 'strict aliasing rules' */
|
|
static inline std::uint32_t extract_32bits_from64_left(const std::uint64_t x)
|
|
{
|
|
std::uint32_t out;
|
|
std::memcpy(&out, reinterpret_cast<const std::uint32_t*>(&x), sizeof(std::uint32_t));
|
|
return out;
|
|
}
|
|
|
|
static inline std::uint32_t extract_32bits_from64_right(const std::uint64_t x)
|
|
{
|
|
std::uint32_t out;
|
|
std::memcpy(&out, reinterpret_cast<const std::uint32_t*>(&x) + 1, sizeof(std::uint32_t));
|
|
return out;
|
|
}
|
|
|
|
static inline void assign_32bits_to64_left(std::uint64_t assign_to, const std::uint32_t take_from)
|
|
{
|
|
std::memcpy(reinterpret_cast<std::uint32_t*>(&assign_to), &take_from, sizeof(std::uint32_t));
|
|
}
|
|
|
|
static inline void assign_32bits_to64_right(std::uint64_t assign_to, const std::uint32_t take_from)
|
|
{
|
|
std::memcpy(reinterpret_cast<std::uint32_t*>(&assign_to) + 1, &take_from, sizeof(std::uint32_t));
|
|
}
|
|
|
|
constexpr static const uint64_t JUMP_X256PP[] = { 0x180ec6d33cfd0aba, 0xd5a61266f0c9392c, 0xa9582618e03fc9aa, 0x39abdc4529b1661c };
|
|
|
|
constexpr static const uint64_t LONG_JUMP_X256PP[] = { 0x76e15d3efefdcbbf, 0xc5004e441c522fb3, 0x77710069854ee241, 0x39109bb02acbe635 };
|
|
|
|
constexpr static const uint32_t JUMP_X128PP[] = { 0x8764000b, 0xf542d2d3, 0x6fa035c3, 0x77f2db5b };
|
|
|
|
constexpr static const uint32_t LONG_JUMP_X128PP[] = { 0xb523952e, 0x0b6f099f, 0xccf5a0ef, 0x1c580662 };
|
|
|
|
template <class int_t, class rng_t>
|
|
static inline void jump_state(const int_t jump_table[4], rng_t &rng)
|
|
{
|
|
int_t s0 = 0;
|
|
int_t s1 = 0;
|
|
int_t s2 = 0;
|
|
int_t s3 = 0;
|
|
for (int i = 0; i < 4; i++)
|
|
{
|
|
for (int b = 0; b < static_cast<int>(8*sizeof(int_t)); b++) // static_cast has been added by Tourmaline Engine!
|
|
{
|
|
if (jump_table[i] & ((int_t)1) << b)
|
|
{
|
|
s0 ^= rng.state[0];
|
|
s1 ^= rng.state[1];
|
|
s2 ^= rng.state[2];
|
|
s3 ^= rng.state[3];
|
|
}
|
|
rng();
|
|
}
|
|
}
|
|
|
|
rng.state[0] = s0;
|
|
rng.state[1] = s1;
|
|
rng.state[2] = s2;
|
|
rng.state[3] = s3;
|
|
}
|
|
|
|
/* This is a fixed-increment version of Java 8's SplittableRandom generator
|
|
See http://dx.doi.org/10.1145/2714064.2660195 and
|
|
http://docs.oracle.com/javase/8/docs/api/java/util/SplittableRandom.html
|
|
|
|
It is a very fast generator passing BigCrush, and it can be useful if
|
|
for some reason you absolutely want 64 bits of state. */
|
|
static inline std::uint64_t splitmix64(const std::uint64_t seed)
|
|
{
|
|
std::uint64_t z = (seed + 0x9e3779b97f4a7c15);
|
|
z = (z ^ (z >> 30)) * 0xbf58476d1ce4e5b9;
|
|
z = (z ^ (z >> 27)) * 0x94d049bb133111eb;
|
|
return z ^ (z >> 31);
|
|
}
|
|
|
|
/* This is xoshiro256++ 1.0, one of our all-purpose, rock-solid generators.
|
|
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.
|
|
|
|
For generating just floating-point numbers, xoshiro256+ is even faster.
|
|
|
|
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 s. */
|
|
class Xoshiro256PP
|
|
{
|
|
public:
|
|
using result_type = std::uint64_t;
|
|
std::uint64_t state[4] = {0x3d23dce41c588f8c, 0x10c770bb8da027b0, 0xc7a4c5e87c63ba25, 0xa830f83239465a2e};
|
|
|
|
constexpr static result_type min()
|
|
{
|
|
return 0;
|
|
}
|
|
|
|
constexpr static result_type max()
|
|
{
|
|
return UINT64_MAX;
|
|
}
|
|
|
|
Xoshiro256PP() = default;
|
|
|
|
~Xoshiro256PP() noexcept = default;
|
|
|
|
Xoshiro256PP(Xoshiro256PP &other) = default;
|
|
|
|
Xoshiro256PP& operator=(const Xoshiro256PP &other) = default;
|
|
|
|
Xoshiro256PP(Xoshiro256PP &&) noexcept = default;
|
|
|
|
Xoshiro256PP& operator=(Xoshiro256PP &&) noexcept = default;
|
|
|
|
void seed(const std::uint64_t seed)
|
|
{
|
|
this->state[0] = splitmix64(splitmix64(seed));
|
|
this->state[1] = splitmix64(this->state[0]);
|
|
this->state[2] = splitmix64(this->state[1]);
|
|
this->state[3] = splitmix64(this->state[2]);
|
|
}
|
|
|
|
void seed(const std::uint64_t seed[4])
|
|
{
|
|
std::memcpy(this->state, seed, 4*sizeof(std::uint64_t));
|
|
}
|
|
|
|
template<class Sseq>
|
|
void seed(Sseq& seq)
|
|
{
|
|
seq.generate(reinterpret_cast<std::uint32_t*>(&this->state[0]),
|
|
reinterpret_cast<std::uint32_t*>(&this->state[0] + 4));
|
|
}
|
|
|
|
Xoshiro256PP(const std::uint64_t seed)
|
|
{
|
|
this->seed(seed);
|
|
}
|
|
|
|
Xoshiro256PP(const std::uint64_t seed[4])
|
|
{
|
|
this->seed(seed);
|
|
}
|
|
|
|
template<class Sseq>
|
|
Xoshiro256PP(Sseq& seq)
|
|
{
|
|
this->seed(seq);
|
|
}
|
|
|
|
result_type operator()()
|
|
{
|
|
const std::uint64_t result = rotl64(this->state[0] + this->state[3], 23) + this->state[0];
|
|
const std::uint64_t t = this->state[1] << 17;
|
|
this->state[2] ^= this->state[0];
|
|
this->state[3] ^= this->state[1];
|
|
this->state[1] ^= this->state[2];
|
|
this->state[0] ^= this->state[3];
|
|
this->state[2] ^= t;
|
|
this->state[3] = rotl64(this->state[3], 45);
|
|
return result;
|
|
}
|
|
|
|
void discard(unsigned long long z)
|
|
{
|
|
for (unsigned long long ix = 0; ix < z; ix++)
|
|
this->operator()();
|
|
}
|
|
|
|
Xoshiro256PP jump()
|
|
{
|
|
Xoshiro256PP new_gen = *this;
|
|
jump_state(JUMP_X256PP, new_gen);
|
|
return new_gen;
|
|
}
|
|
|
|
Xoshiro256PP long_jump()
|
|
{
|
|
Xoshiro256PP new_gen = *this;
|
|
jump_state(LONG_JUMP_X256PP, new_gen);
|
|
return new_gen;
|
|
}
|
|
|
|
bool operator==(const Xoshiro256PP &rhs)
|
|
{
|
|
return std::memcmp(this->state, rhs.state, 4*sizeof(std::uint64_t));
|
|
}
|
|
|
|
#ifndef HAS_CPP20
|
|
bool operator!=(const Xoshiro256PP &rhs)
|
|
{
|
|
return !std::memcmp(this->state, rhs.state, 4*sizeof(std::uint64_t));
|
|
}
|
|
#endif
|
|
|
|
template< class CharT, class Traits >
|
|
friend std::basic_ostream<CharT,Traits>&
|
|
operator<<(std::basic_ostream<CharT,Traits>& ost, const Xoshiro256PP& e)
|
|
{
|
|
ost.write(reinterpret_cast<const char*>(&e.state[0]), sizeof(std::uint64_t));
|
|
ost.put(' ');
|
|
ost.write(reinterpret_cast<const char*>(&e.state[1]), sizeof(std::uint64_t));
|
|
ost.put(' ');
|
|
ost.write(reinterpret_cast<const char*>(&e.state[2]), sizeof(std::uint64_t));
|
|
ost.put(' ');
|
|
ost.write(reinterpret_cast<const char*>(&e.state[3]), sizeof(std::uint64_t));
|
|
return ost;
|
|
}
|
|
template< class CharT, class Traits >
|
|
friend std::basic_istream<CharT,Traits>&
|
|
operator>>(std::basic_istream<CharT,Traits>& ist, Xoshiro256PP& e)
|
|
{
|
|
ist.read(reinterpret_cast<char*>(&e.state[0]), sizeof(std::uint64_t));
|
|
ist.get();
|
|
ist.read(reinterpret_cast<char*>(&e.state[1]), sizeof(std::uint64_t));
|
|
ist.get();
|
|
ist.read(reinterpret_cast<char*>(&e.state[2]), sizeof(std::uint64_t));
|
|
ist.get();
|
|
ist.read(reinterpret_cast<char*>(&e.state[3]), sizeof(std::uint64_t));
|
|
return ist;
|
|
}
|
|
};
|
|
|
|
/* This is xoshiro128++ 1.0, one of our 32-bit all-purpose, rock-solid
|
|
generators. It has excellent speed, a state size (128 bits) that is
|
|
large enough for mild parallelism, and it passes all tests we are aware
|
|
of.
|
|
|
|
For generating just single-precision (i.e., 32-bit) floating-point
|
|
numbers, xoshiro128+ is even faster.
|
|
|
|
The state must be seeded so that it is not everywhere zero. */
|
|
class Xoshiro128PP
|
|
{
|
|
public:
|
|
using result_type = std::uint32_t;
|
|
std::uint32_t state[4] = {0x1c588f8c, 0x3d23dce4, 0x8da027b0, 0x10c770bb};
|
|
|
|
constexpr static result_type min()
|
|
{
|
|
return 0;
|
|
}
|
|
|
|
constexpr static result_type max()
|
|
{
|
|
return UINT32_MAX;
|
|
}
|
|
|
|
Xoshiro128PP() = default;
|
|
|
|
~Xoshiro128PP() noexcept = default;
|
|
|
|
Xoshiro128PP(Xoshiro128PP &other) = default;
|
|
|
|
Xoshiro128PP& operator=(const Xoshiro128PP &other) = default;
|
|
|
|
Xoshiro128PP(Xoshiro128PP &&) noexcept = default;
|
|
|
|
Xoshiro128PP& operator=(Xoshiro128PP &&) noexcept = default;
|
|
|
|
void seed(const std::uint64_t seed)
|
|
{
|
|
const auto t1 = splitmix64(seed);
|
|
const auto t2 = splitmix64(t1);
|
|
this->state[0] = splitmix64(extract_32bits_from64_left(t1));
|
|
this->state[1] = splitmix64(extract_32bits_from64_right(t1));
|
|
this->state[2] = splitmix64(extract_32bits_from64_left(t2));
|
|
this->state[3] = splitmix64(extract_32bits_from64_right(t2));
|
|
}
|
|
|
|
void seed(const std::uint32_t seed)
|
|
{
|
|
std::uint64_t temp;
|
|
assign_32bits_to64_left(temp, seed);
|
|
assign_32bits_to64_right(temp, seed);
|
|
this->seed(temp);
|
|
}
|
|
|
|
void seed(const std::uint64_t seed[2])
|
|
{
|
|
std::memcpy(this->state, seed, 4*sizeof(std::uint32_t));
|
|
}
|
|
|
|
void seed(const std::uint32_t seed[4])
|
|
{
|
|
std::memcpy(this->state, seed, 4*sizeof(std::uint32_t));
|
|
}
|
|
|
|
template<class Sseq>
|
|
void seed(Sseq& seq)
|
|
{
|
|
seq.generate(&this->state[0], &this->state[0] + 4);
|
|
}
|
|
|
|
Xoshiro128PP(const std::uint32_t seed)
|
|
{
|
|
this->seed(seed);
|
|
}
|
|
|
|
Xoshiro128PP(const std::uint64_t seed)
|
|
{
|
|
this->seed(seed);
|
|
}
|
|
|
|
Xoshiro128PP(const std::uint32_t seed[4])
|
|
{
|
|
this->seed(seed);
|
|
}
|
|
|
|
Xoshiro128PP(const std::uint64_t seed[2])
|
|
{
|
|
this->seed(seed);
|
|
}
|
|
|
|
template<class Sseq>
|
|
Xoshiro128PP(Sseq& seq)
|
|
{
|
|
this->seed(seq);
|
|
}
|
|
|
|
result_type operator()()
|
|
{
|
|
const std::uint32_t result = rotl32(this->state[0] + this->state[3], 7) + this->state[0];
|
|
const std::uint32_t t = this->state[1] << 9;
|
|
this->state[2] ^= this->state[0];
|
|
this->state[3] ^= this->state[1];
|
|
this->state[1] ^= this->state[2];
|
|
this->state[0] ^= this->state[3];
|
|
this->state[2] ^= t;
|
|
this->state[3] = rotl32(this->state[3], 11);
|
|
return result;
|
|
}
|
|
|
|
Xoshiro128PP jump()
|
|
{
|
|
Xoshiro128PP new_gen = *this;
|
|
jump_state(JUMP_X128PP, new_gen);
|
|
return new_gen;
|
|
}
|
|
|
|
Xoshiro128PP long_jump()
|
|
{
|
|
Xoshiro128PP new_gen = *this;
|
|
jump_state(LONG_JUMP_X128PP, new_gen);
|
|
return new_gen;
|
|
}
|
|
|
|
void discard(unsigned long long z)
|
|
{
|
|
for (unsigned long long ix = 0; ix < z; ix++)
|
|
this->operator()();
|
|
}
|
|
|
|
bool operator==(const Xoshiro128PP &rhs)
|
|
{
|
|
return std::memcmp(this->state, rhs.state, 4*sizeof(std::uint32_t));
|
|
}
|
|
|
|
#ifndef HAS_CPP20
|
|
bool operator!=(const Xoshiro128PP &rhs)
|
|
{
|
|
return !std::memcmp(this->state, rhs.state, 4*sizeof(std::uint32_t));
|
|
}
|
|
#endif
|
|
|
|
template< class CharT, class Traits >
|
|
friend std::basic_ostream<CharT,Traits>&
|
|
operator<<(std::basic_ostream<CharT,Traits>& ost, const Xoshiro128PP& e)
|
|
{
|
|
ost.write(reinterpret_cast<const char*>(&e.state[0]), sizeof(std::uint32_t));
|
|
ost.put(' ');
|
|
ost.write(reinterpret_cast<const char*>(&e.state[1]), sizeof(std::uint32_t));
|
|
ost.put(' ');
|
|
ost.write(reinterpret_cast<const char*>(&e.state[2]), sizeof(std::uint32_t));
|
|
ost.put(' ');
|
|
ost.write(reinterpret_cast<const char*>(&e.state[3]), sizeof(std::uint32_t));
|
|
return ost;
|
|
}
|
|
template< class CharT, class Traits >
|
|
friend std::basic_istream<CharT,Traits>&
|
|
operator>>(std::basic_istream<CharT,Traits>& ist, Xoshiro128PP& e)
|
|
{
|
|
ist.read(reinterpret_cast<char*>(&e.state[0]), sizeof(std::uint32_t));
|
|
ist.get();
|
|
ist.read(reinterpret_cast<char*>(&e.state[1]), sizeof(std::uint32_t));
|
|
ist.get();
|
|
ist.read(reinterpret_cast<char*>(&e.state[2]), sizeof(std::uint32_t));
|
|
ist.get();
|
|
ist.read(reinterpret_cast<char*>(&e.state[3]), sizeof(std::uint32_t));
|
|
return ist;
|
|
}
|
|
};
|
|
|
|
}
|
|
|
|
#endif
|