libcruft-util/cruft/util/rand/xoshiro.cpp

139 lines
3.2 KiB
C++

/*
* 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.
*/
/*
* This Source Code Form is subject to the terms of the Mozilla Public
* License, v. 2.0. If a copy of the MPL was not distributed with this
* file, You can obtain one at http://mozilla.org/MPL/2.0/.
*
* Copyright 2020, Danny Robson <danny@nerdcruft.net>
*/
#include "xoshiro.hpp"
#include "splitmix64.hpp"
#include "../bitwise.hpp"
#include <algorithm>
using cruft::rand::xoshiro256plusplus;
///////////////////////////////////////////////////////////////////////////////
xoshiro256plusplus::xoshiro256plusplus (u64 _seed)
{
splitmix64 gen (_seed);
std::generate (
std::begin (m_state),
std::end (m_state),
gen
);
}
//-----------------------------------------------------------------------------
xoshiro256plusplus::xoshiro256plusplus (std::seed_seq &seq)
{
u32 parts[8];
seq.generate (std::begin (parts), std::end (parts));
for (int i = 0; i < std::ssize (m_state); ++i)
m_state[i] = u64 (parts[i * 2]) << 32 | parts[i * 2 + 1];
}
///////////////////////////////////////////////////////////////////////////////
u64
xoshiro256plusplus::operator() (void)
{
u64 const result = rotatel (m_state[0] + m_state[3], 23) + m_state[0];
u64 const 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] = rotatel (m_state[3], 45);
return result;
}
///////////////////////////////////////////////////////////////////////////////
void
xoshiro256plusplus::jump (void)
{
static u64 constinit JUMP[] = {
0x180ec6d33cfd0aba,
0xd5a61266f0c9392c,
0xa9582618e03fc9aa,
0x39abdc4529b1661c
};
u64 s0 = 0;
u64 s1 = 0;
u64 s2 = 0;
u64 s3 = 0;
for (std::size_t i = 0; i < std::size (JUMP); ++i) {
for (int b = 0; b < 64; b++) {
if (JUMP[i] & u64 (1) << b) {
s0 ^= m_state[0];
s1 ^= m_state[1];
s2 ^= m_state[2];
s3 ^= m_state[3];
}
(*this) ();
}
}
m_state[0] = s0;
m_state[1] = s1;
m_state[2] = s2;
m_state[3] = s3;
}
//-----------------------------------------------------------------------------
void
xoshiro256plusplus::long_jump (void)
{
static constinit u64 LONG_JUMP[] = {
0x76e15d3efefdcbbf,
0xc5004e441c522fb3,
0x77710069854ee241,
0x39109bb02acbe635
};
u64 s0 = 0;
u64 s1 = 0;
u64 s2 = 0;
u64 s3 = 0;
for (std::size_t i = 0; i < std::size (LONG_JUMP); ++i) {
for (int b = 0; b < 64; b++) {
if (LONG_JUMP[i] & u64{1} << b) {
s0 ^= m_state[0];
s1 ^= m_state[1];
s2 ^= m_state[2];
s3 ^= m_state[3];
}
(*this) ();
}
}
m_state[0] = s0;
m_state[1] = s1;
m_state[2] = s2;
m_state[3] = s3;
}