+ xoshiro256ss (copied from kalman project)

This commit is contained in:
Roland Conybeare 2023-09-23 13:07:45 -04:00
commit fadbcd7b54
12 changed files with 1130 additions and 0 deletions

View file

@ -0,0 +1,36 @@
/* @file engine_concept.hpp */
#pragma once
#include <concepts>
#include <random>
namespace xo {
namespace rng {
/* an engine generates psuedo-random bits.
* given
* RngEngine eng = ...;
*
* RngEngine::result_type x = eng();
*
* puts random bits into x.
*/
template <class RngEngine>
concept engine_concept = requires(RngEngine engine, typename RngEngine::result_type r) {
/* note: the first 4 requirements characterize UniformRandomBitGenerator */
typename RngEngine::result_type;
{ RngEngine(r) };
{ engine.min() } -> std::same_as<typename RngEngine::result_type>;
{ engine.max() } -> std::same_as<typename RngEngine::result_type>;
/* must return value in closed interval [.min(), .max()] */
{ engine() } -> std::same_as<typename RngEngine::result_type>;
{ engine.seed() };
{ engine.seed(r) };
{ engine == engine };
{ engine != engine };
} && std::copyable<RngEngine> && std::uniform_random_bit_generator<RngEngine>;
} /*namespace rng*/
} /*namespace xo*/
/* end engine_concept.hpp */

View file

@ -0,0 +1,70 @@
/* @file random_seed.hpp */
#include "indentlog/print/array.hpp"
#include <iostream>
#include <cstdint>
#include <stdlib.h>
namespace xo {
namespace rng {
/* generate a 64-bit random seed using /dev/urandom or similar source.
* This is relatively expensive; at least cost of a system call
* + may block if host has rebooted recently
*
* Require:
* - T is null-constructible.
*
* return value will contain a T-instance in which representation
* has been populated with random bits. Expecting T to be something
* like int32_t, or std::array<uint64_t, ..>
*/
template<typename T>
void random_seed(T * p_seed) {
/* NOTE: arc4random_buf() works on darwin/nix;
* probably need to do something else on intel linux
*/
arc4random_buf(p_seed, sizeof(*p_seed));
} /*random_seed*/
template<typename T>
T random_seed() {
T retval;
random_seed(&retval);
return retval;
} /*random_seed*/
/* RAII-style random-number seed
*
* Usage:
*
* Seed<xoshiro256ss> seed;
*
* auto eng = xoshiro256ss(seed);
* or
* auto rng = UnitIntervalGen<xoshiro256ss>::make(seed);
*/
template<typename Engine>
struct Seed {
using seed_type = typename Engine::seed_type;
Seed() { random_seed(&seed_); }
operator seed_type const & () const { return seed_; }
seed_type seed_;
}; /*Seed*/
template<typename T>
inline std::ostream &
operator<<(std::ostream & os,
Seed<T> const & x)
{
os << x.seed_;
return os;
} /*operator<<*/
} /*namespace rng*/
} /*namespace xo*/
/* end random_seed.hpp */

View file

@ -0,0 +1,169 @@
/* @file xoshiro256.hpp */
#pragma once
#include "engine_concept.hpp"
#include <iostream>
#include <array>
#include <limits>
#include <cstdint>
namespace xo {
namespace rng {
/* engine for producing 64-bit random numbers
*
* see https:/en.wikipedia.org/wiki/Xorshift#xoshiro256**
*
* - satisfies c++ UniformRandomBitGenerator
* - satisfies c++
*
* Note: zero seed --> constant output sequence {0, 0, 0, ...}
*/
class xoshiro256ss {
public:
using result_type = std::uint64_t;
using seed_type = std::array<std::uint64_t, 4>;
public:
/* null state -- generates constant stream of 0 bits */
xoshiro256ss() : xoshiro256ss(0) {}
/* copy ctor */
xoshiro256ss(xoshiro256ss const & x) = default;
xoshiro256ss(seed_type const & seed) : s_(seed) {}
/* fallback version -- deprecated */
xoshiro256ss(std::uint64_t seed)
{
this->s_[0] = 0;
this->s_[1] = seed;
this->s_[2] = 0;
this->s_[3] = 0;
generate();
}
static constexpr std::uint64_t min() { return 0; }
static constexpr std::uint64_t max() { return std::numeric_limits<std::uint64_t>::max(); }
static std::uint64_t rol64(std::uint64_t x, std::int64_t k)
{
return (x << k) | (x >> (64 - k));
}
static bool equal(xoshiro256ss const & x, xoshiro256ss const & y) {
return ((x.s_[0] == y.s_[0])
&& (x.s_[1] == y.s_[1])
&& (x.s_[2] == y.s_[2])
&& (x.s_[3] == y.s_[3]));
}
/* puts generator into null state */
void seed() { *this = xoshiro256ss(); }
void seed(std::uint64_t s) { *this = xoshiro256ss{s}; }
/* e.g. used with std::seed_seq<> */
template <typename SeedSeq>
void seed(SeedSeq & sseq) {
sseq.generate(s_.begin(), s_.end());
}
std::uint64_t generate() {
std::array<std::uint64_t, 4> & s = (this->s_);
std::uint64_t const result = rol64(s[1] * 5, 7) * 9;
std::uint64_t const t = s[1] << 17;
s[2] ^= s[0];
s[3] ^= s[1];
s[1] ^= s[2];
s[0] ^= s[3];
s[2] ^= t;
s[3] = rol64(s[3], 45);
return result;
} /*generate*/
/* advance to same state as obtained from z calls to .generate(). O(z) !
* usually better to use jump().
*
* providing .discard() to satisfy c++ named requirement _RandomNumberEngine_
*/
void discard(std::uint64_t z) {
for (std::uint64_t i=0; i<z; ++i)
this->generate();
}
/* equivalent to .discard(2^128), but uses O(1) time
*
* (may use in multithreaded program to get determinsitic non-overlapping random sequences)
*/
void jump() {
std::array<std::uint64_t, 4> const s_jump_v
= {{0x180ec6d33cfd0aba,
0xd5a61266f0c9392c,
0xa9582618e03fc9aa,
0x39abdc4529b1661c}};
std::array<std::uint64_t, 4> & s = (this->s_);
std::uint64_t s0 = 0;
std::uint64_t s1 = 0;
std::uint64_t s2 = 0;
std::uint64_t s3 = 0;
for (std::uint32_t i = 0; i < s_jump_v.size(); ++i) {
for (std::uint32_t bit = 0; bit < 64; ++bit) {
if (s_jump_v[i] & 1UL << bit) {
s0 ^= s[0];
s1 ^= s[1];
s2 ^= s[2];
s3 ^= s[3];
}
this->generate();
}
}
s[0] = s0;
s[1] = s1;
s[2] = s2;
s[3] = s3;
} /*jump*/
/* inverse of .load() */
void print(std::ostream & os) const {
os << "<xoshiro256ss " << s_[0] << " " << s_[1] << " " << s_[2] << " " << s_[3] << ">";
}
/* inverse of .print() */
void load(std::istream & is) {
std::string header, trailer;
std::array<std::uint64_t, 4> sv;
is >> header >> sv[0] >> sv[1] >> sv[2] >> sv[3] >> trailer;
if ((header != "<xoshiro256ss") || trailer != ">")
throw std::runtime_error("xoshiro256ss.load: bad input format, expecting input like <xoshiro256ss $s0 $s1 $s2 $s3>");
this->s_ = sv;
} /*load*/
std::uint64_t operator()() { return generate(); }
private:
/* state */
std::array<std::uint64_t, 4> s_;
}; /*xoshiro256ss*/
inline bool operator==(xoshiro256ss const & x, xoshiro256ss const & y) {
return xoshiro256ss::equal(x, y);
}
inline bool operator!=(xoshiro256ss const & x, xoshiro256ss const & y) {
return !xoshiro256ss::equal(x, y);
}
static_assert(engine_concept<xoshiro256ss>);
} /*namespace rng*/
} /*namespace xo*/
/* end xoshiro256.hpp */