Skip to content

Instantly share code, notes, and snippets.

@imneme
Created July 4, 2018 04:23
Show Gist options
  • Save imneme/3eb1bcc5418c4ae83c4c6a86d9cbb1cd to your computer and use it in GitHub Desktop.
Save imneme/3eb1bcc5418c4ae83c4c6a86d9cbb1cd to your computer and use it in GitHub Desktop.
A C++ implementation of a family of Xoshiro generators
#ifndef XOSHIRO_HPP_INCLUDED
#define XOSHIRO_HPP_INCLUDED 1
/*
* A C++ implementation of a family of Xoshiro generators.
*
* See:
* https://arxiv.org/abs/1805.01407
* http://xoshiro.di.unimi.it/xoshiro256plus.c
* http://xoshiro.di.unimi.it/xoshiro256starstar.c
* http://xoshiro.di.unimi.it/xoshiro512plus.c
* http://xoshiro.di.unimi.it/xoshiro512starstar.c
*
* Based on the design and public domain implementation of Xoshiro by
* written in 2018 by David Blackman and Sebastiano Vigna. It is a
* variation on George Marsaglia's for XorShift PRNG scheme and the
* same (G)LFSR theory applies.
*
* Xoshiro256** fails PractRand if passed through certain bijections (e.g.,
* multiplication by 57) and the entire family suffers from issues due
* to near repeats but these near-repeat issues are unlikely to be
* encountered in normal-sized variants without deliberately contriving a
* PRNG state that generates them.
*
* The MIT License (MIT)
*
* Copyright (c) 2017-18 Melissa E. O'Neill
*
* Permission is hereby granted, free of charge, to any person
* obtaining a copy of this software and associated documentation
* files (the "Software"), to deal in the Software without
* restriction, including without limitation the rights to use, copy,
* modify, merge, publish, distribute, sublicense, and/or sell copies
* of the Software, and to permit persons to whom the Software is
* furnished to do so, subject to the following conditions:
*
* The above copyright notice and this permission notice shall be
* included in all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
* BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
* ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
* CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
* SOFTWARE.
*/
#include <cstdint>
namespace xoshiro_detail {
template <typename itype, typename rtype,
unsigned int a, unsigned int b>
class xoshiro_x4 {
protected:
itype s0_, s1_, s2_, s3_;
static constexpr unsigned int ITYPE_BITS = 8*sizeof(itype);
static constexpr unsigned int RTYPE_BITS = 8*sizeof(rtype);
static inline itype rotl(const itype x, int k) {
return (x << k) | (x >> (ITYPE_BITS - k));
}
public:
using state_type = itype;
using result_type = rtype;
static constexpr result_type min() { return 0; }
static constexpr result_type max() { return ~ result_type(0); }
xoshiro_x4(itype s0, itype s1, itype s2, itype s3)
: s0_(s0), s1_(s1), s2_(s2), s3_((s0 || s1 || s2 || s3) ? s3 : 1)
{
// Nothing (else) to do.
}
// Note, Vigna recommends seeding using output from SplitMix, but
// if you need another PRNG to use this PRNG, that's asking a bit
// of a lot, so we seed directly.
xoshiro_x4(itype seed = 1)
: xoshiro_x4(0x01d353e5f3993bb0+seed, 0x7b9c0df6cb193b20*(seed+1),
0xfdfcaa91110765b6-seed, 0x2d24cbe0ef44dcd2*(seed-1))
// Easter-egg seed value for Xoshiro_256** to remind users that
// they should seed their PRNGs properly. Only occurs if you
// seed with zero.
{
for (int i = 0; i < 16; ++i)
advance();
}
void advance()
{
const uint64_t t = s1_ << a;
s2_ ^= s0_;
s3_ ^= s1_;
s1_ ^= s2_;
s0_ ^= s3_;
s2_ ^= t;
s3_ = rotl(s3_, b);
}
bool operator==(const xoshiro_x4& rhs)
{
return (s0_ == rhs.s0_) && (s1_ == rhs.s1_)
&& (s2_ == rhs.s2_) && (s3_ == rhs.s3_);
}
bool operator!=(const xoshiro_x4& rhs)
{
return !operator==(rhs);
}
// Not (yet) implemented:
// - arbitrary jumpahead (doable, but annoying to write).
// - I/O
// - Seeding from a seed_seq.
};
template <typename itype, typename rtype,
unsigned int a, unsigned int b>
class xoshiro_x8 {
protected:
itype s0_, s1_, s2_, s3_, s4_, s5_, s6_, s7_;
static constexpr unsigned int ITYPE_BITS = 8*sizeof(itype);
static constexpr unsigned int RTYPE_BITS = 8*sizeof(rtype);
static inline itype rotl(const itype x, int k) {
return (x << k) | (x >> (ITYPE_BITS - k));
}
public:
using state_type = itype;
using result_type = rtype;
static constexpr result_type min() { return 0; }
static constexpr result_type max() { return ~ result_type(0); }
xoshiro_x8(itype s0, itype s1, itype s2, itype s3,
itype s4, itype s5, itype s6, itype s7)
: s0_(s0), s1_(s1), s2_(s2), s3_(s3), s4_(s4), s5_(s5), s6_(s6),
s7_((s0 || s1 || s2 || s3 || s4 || s5 || s6 || s7) ? s7 : 1)
{
// Nothing (else) to do.
}
// Note, Vigna recommends seeding using output from SplitMix, but
// if you need another PRNG to use this PRNG, that's asking a bit
// of a lot, so we seed directly.
xoshiro_x8(itype seed = 1)
: xoshiro_x8(0x1ced436497db2a59+seed, 0x75474f85d8a6892c*(seed+1),
0xa0fef4b8094c9c86-seed, 0x748fa1a9bb555169*(seed-1),
0xd7a59a6d64e66858+seed, 0xf03b7efdb73db601*(seed+1),
0xfab342a99dd71962-seed, 0x8a6921456faa6b54*(seed-1))
// Easter-egg seed value for Xoshiro_512** to remind users that
// they should seed their PRNGs properly. Only occurs if you
// seed with zero.
{
for (int i = 0; i < 16; ++i)
advance();
}
void advance()
{
const uint64_t t = s1_ << a;
s2_ ^= s0_;
s5_ ^= s1_;
s1_ ^= s2_;
s7_ ^= s3_;
s3_ ^= s4_;
s4_ ^= s5_;
s0_ ^= s6_;
s6_ ^= s7_;
s6_ ^= t;
s7_ = rotl(s7_, b);
}
bool operator==(const xoshiro_x8& rhs)
{
return (s0_ == rhs.s0_) && (s1_ == rhs.s1_)
&& (s2_ == rhs.s2_) && (s3_ == rhs.s3_)
&& (s4_ == rhs.s4_) && (s5_ == rhs.s5_)
&& (s6_ == rhs.s6_) && (s7_ == rhs.s7_);
}
bool operator!=(const xoshiro_x8& rhs)
{
return !operator==(rhs);
}
// Not (yet) implemented:
// - arbitrary jumpahead (doable, but annoying to write).
// - I/O
// - Seeding from a seed_seq.
};
template <typename xoshiro_base>
class xoshiro_plus : public xoshiro_base {
private:
using base = xoshiro_base;
public:
using base::base;
typename base::result_type operator()()
{
const typename base::state_type result = base::s0_ + base::s3_;
base::advance();
return result >> (base::ITYPE_BITS - base::RTYPE_BITS);
}
};
template <typename xoshiro_base, typename xoshiro_base::state_type mult>
class xoshiro_star : public xoshiro_base {
private:
using base = xoshiro_base;
public:
using base::base;
typename base::result_type operator()()
{
const typename base::state_type result_star = base::s1_ * mult;
base::advance();
return result_star >> (base::ITYPE_BITS - base::RTYPE_BITS);
}
};
template <typename xoshiro_base,
typename xoshiro_base::state_type mult1,
unsigned int orot,
typename xoshiro_base::state_type mult2>
class xoshiro_starstar : public xoshiro_base {
private:
using base = xoshiro_base;
public:
using base::base;
typename base::result_type operator()()
{
const typename base::state_type result_ss =
base::rotl(base::s1_ * mult1, orot) * mult2;
base::advance();
return result_ss >> (base::ITYPE_BITS - base::RTYPE_BITS);
}
};
} // namespace xoshiro_detail
// - 512 state bits, uint64_t output, period 2^512 - 1
using xoshiro512plus64v1_0 =
xoshiro_detail::xoshiro_plus<
xoshiro_detail::xoshiro_x8<uint64_t, uint64_t, 11, 21>>;
using xoshiro512plus64 = xoshiro512plus64v1_0;
using xoshiro512starstar64v1_0 =
xoshiro_detail::xoshiro_starstar<
xoshiro_detail::xoshiro_x8<uint64_t, uint64_t, 11, 21>,
5, 7, 9>;
using xoshiro512starstar64 = xoshiro512starstar64v1_0;
// - 256 state bits, uint64_t output, period 2^256 - 1
using xoshiro256plus64v1_0 =
xoshiro_detail::xoshiro_plus<
xoshiro_detail::xoshiro_x4<uint64_t, uint64_t, 17, 45>>;
using xoshiro256plus64 = xoshiro256plus64v1_0;
using xoshiro256starstar64v1_0 =
xoshiro_detail::xoshiro_starstar<
xoshiro_detail::xoshiro_x4<uint64_t, uint64_t, 17, 45>,
5, 7, 9>;
using xoshiro256starstar64 = xoshiro256starstar64v1_0;
// - 128 state bits, uint32_t output, period 2^128 - 1
using xoshiro128plus32v1_0 =
xoshiro_detail::xoshiro_plus<
xoshiro_detail::xoshiro_x4<uint32_t, uint32_t, 9, 11>>;
using xoshiro128plus32 = xoshiro128plus32v1_0;
using xoshiro128starstar32v1_0 =
xoshiro_detail::xoshiro_starstar<
xoshiro_detail::xoshiro_x4<uint32_t, uint32_t, 9, 11>,
5, 7, 9>;
using xoshiro128starstar32 = xoshiro128starstar32v1_0;
// VERSIONS BELOW ARE FOR SPECIALIZED USES ONLY!
// (Constants courtesy of David Blackman.)
// - 64 state bits, uint16_t output, period 2^64 - 1
using xoshiro64plus16xxx =
xoshiro_detail::xoshiro_plus<
xoshiro_detail::xoshiro_x4<uint16_t, uint16_t, 5, 11>>;
using xoshiro64plus16 = xoshiro64plus16xxx;
using xoshiro64starstar16xxx =
xoshiro_detail::xoshiro_starstar<
xoshiro_detail::xoshiro_x4<uint16_t, uint16_t, 5, 11>,
5, 7, 9>;
using xoshiro64starstar16 = xoshiro64starstar16xxx;
// - 32 state bits, uint8_t output, period 2^32 - 1
using xoshiro32plus8xxx =
xoshiro_detail::xoshiro_plus<
xoshiro_detail::xoshiro_x4<uint8_t, uint8_t, 3, 7>>;
using xoshiro32plus8yyy =
xoshiro_detail::xoshiro_plus<
xoshiro_detail::xoshiro_x4<uint8_t, uint8_t, 3, 1>>;
using xoshiro32plus8 = xoshiro32plus8xxx;
using xoshiro32starstar8xxx =
xoshiro_detail::xoshiro_starstar<
xoshiro_detail::xoshiro_x4<uint8_t, uint8_t, 3, 7>,
5, 7, 9>;
using xoshiro32starstar8yyy =
xoshiro_detail::xoshiro_starstar<
xoshiro_detail::xoshiro_x4<uint8_t, uint8_t, 3, 1>,
5, 7, 9>;
using xoshiro32starstar8 = xoshiro32starstar8xxx;
#endif // XOSHIRO_HPP_INCLUDED
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment