libcruft-util/rand/xoshiro.cpp

/*
 * 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;
}