/*
 * Copyright Matt Borland 2022 - 2025.
 * Distributed under the Boost Software License, Version 1.0. (See
        * accompanying file LICENSE_1_0.txt or copy at
 * http://www.boost.org/LICENSE_1_0.txt)
 *
 * See http://www.boost.org for most recent version including documentation.
 *
 * $Id$
 */

#ifndef BOOST_RANDOM_DETAIL_XOSHIRO_BASE
#define BOOST_RANDOM_DETAIL_XOSHIRO_BASE

#include <boost/random/splitmix64.hpp>
#include <boost/random/detail/seed.hpp>
#include <boost/throw_exception.hpp>
#include <boost/config.hpp>
#include <array>
#include <utility>
#include <stdexcept>
#include <limits>
#include <initializer_list>
#include <cstdint>
#include <cstdlib>
#include <string>
#include <ios>
#include <istream>
#include <type_traits>
#include <iterator>

namespace boost {
namespace random {
namespace detail {

// N is the number of words (e.g. for xoshiro 256 N=4)
template <typename Derived, std::size_t N, typename OutputType = std::uint64_t, typename BlockType = std::uint64_t>
class xoshiro_base
{
protected:

    std::array<BlockType, N> state_;

private:

    using xoshiro_type = std::integral_constant<BlockType, N>;

    inline std::uint64_t concatenate(std::uint32_t word1, std::uint32_t word2) noexcept
    {
        return static_cast<std::uint64_t>(word1) << 32U | word2;
    }

    template <typename Sseq>
    inline void sseq_seed_64(Sseq& seq)
    {
        std::array<std::uint32_t, N * 2> seeds;
        seq.generate(seeds.begin(), seeds.end());

        for (std::size_t i = 0; i < state_.size(); ++i)
        {
            state_[i] = concatenate(seeds[2*i], seeds[2*i + 1]);
        }
    }

    template <typename Sseq>
    inline void sseq_seed_32(Sseq& seq)
    {
        seq.generate(state_.begin(), state_.end());
    }

    inline void jump_impl(const std::integral_constant<std::uint64_t, 4>&) noexcept
    {
        constexpr std::array<std::uint64_t, 4U> jump = {{ UINT64_C(0x180ec6d33cfd0aba), UINT64_C(0xd5a61266f0c9392c),
                                                          UINT64_C(0xa9582618e03fc9aa), UINT64_C(0x39abdc4529b1661c) }};

        std::uint64_t s0 = 0;
        std::uint64_t s1 = 0;
        std::uint64_t s2 = 0;
        std::uint64_t s3 = 0;

        for (std::uint64_t i = 0; i < jump.size(); i++)
        {
            for (std::size_t b = 0; b < 64U; b++)
            {
                if (jump[i] & UINT64_C(1) << b) {
                    s0 ^= state_[0];
                    s1 ^= state_[1];
                    s2 ^= state_[2];
                    s3 ^= state_[3];
                }

                next();
            }
        }

        state_[0] = s0;
        state_[1] = s1;
        state_[2] = s2;
        state_[3] = s3;
    }

    inline void jump_impl(const std::integral_constant<std::uint64_t, 8>&) noexcept
    {
        constexpr std::array<std::uint64_t, 8U> jump = {{ UINT64_C(0x33ed89b6e7a353f9), UINT64_C(0x760083d7955323be),
                                                          UINT64_C(0x2837f2fbb5f22fae), UINT64_C(0x4b8c5674d309511c),
                                                          UINT64_C(0xb11ac47a7ba28c25), UINT64_C(0xf1be7667092bcc1c),
                                                          UINT64_C(0x53851efdb6df0aaf), UINT64_C(0x1ebbc8b23eaf25db) }};

        std::array<std::uint64_t, 8U> t = {{ 0, 0, 0, 0, 0, 0, 0, 0 }};

        for (std::size_t i = 0; i < jump.size(); ++i)
        {
            for (std::size_t b = 0; b < 64U; ++b)
            {
                if (jump[i] & UINT64_C(1) << b)
                {
                    for (std::size_t w = 0; w < state_.size(); ++w)
                    {
                        t[w] ^= state_[w];
                    }
                }

                next();
            }
        }

        state_ = t;
    }

    inline void long_jump_impl(const std::integral_constant<std::uint64_t, 4>&) noexcept
    {
        constexpr std::array<std::uint64_t, 4> long_jump = {{ UINT64_C(0x76e15d3efefdcbbf), UINT64_C(0xc5004e441c522fb3),
                                                              UINT64_C(0x77710069854ee241), UINT64_C(0x39109bb02acbe635) }};

        std::uint64_t s0 = 0;
        std::uint64_t s1 = 0;
        std::uint64_t s2 = 0;
        std::uint64_t s3 = 0;

        for (std::size_t i = 0; i < long_jump.size(); i++)
        {
            for (std::size_t b = 0; b < 64; b++)
            {
                if (long_jump[i] & UINT64_C(1) << b)
                {
                    s0 ^= state_[0];
                    s1 ^= state_[1];
                    s2 ^= state_[2];
                    s3 ^= state_[3];
                }

                next();
            }
        }

        state_[0] = s0;
        state_[1] = s1;
        state_[2] = s2;
        state_[3] = s3;
    }

    inline void long_jump_impl(const std::integral_constant<std::uint64_t, 8>&) noexcept
    {
        constexpr std::array<std::uint64_t, 8U> long_jump = {{ UINT64_C(0x11467fef8f921d28), UINT64_C(0xa2a819f2e79c8ea8),
                                                               UINT64_C(0xa8299fc284b3959a), UINT64_C(0xb4d347340ca63ee1),
                                                               UINT64_C(0x1cb0940bedbff6ce), UINT64_C(0xd956c5c4fa1f8e17),
                                                               UINT64_C(0x915e38fd4eda93bc), UINT64_C(0x5b3ccdfa5d7daca5) }};

        std::array<std::uint64_t, 8U> t = {{ 0, 0, 0, 0, 0, 0, 0, 0 }};

        for (std::size_t i = 0; i < long_jump.size(); ++i)
        {
            for (std::size_t b = 0; b < 64U; ++b)
            {
                if (long_jump[i] & UINT64_C(1) << b)
                {
                    for (std::size_t w = 0; w < state_.size(); ++w)
                    {
                        t[w] ^= state_[w];
                    }
                }

                next();
            }
        }

        state_ = t;
    }

    inline void jump_impl(const std::integral_constant<std::uint32_t, 4>&) noexcept
    {
        constexpr std::array<std::uint32_t, 4> jump = {{ UINT32_C(0x8764000b), UINT32_C(0xf542d2d3), 
                                                         UINT32_C(0x6fa035c3), UINT32_C(0x77f2db5b) }};

        std::uint32_t s0 = 0;
        std::uint32_t s1 = 0;
        std::uint32_t s2 = 0;
        std::uint32_t s3 = 0;

        for (std::size_t i = 0; i < jump.size(); i++)
        {
            for (std::size_t b = 0; b < 32U; b++)
            {
                if (jump[i] & UINT32_C(1) << b)
                {
                    s0 ^= state_[0];
                    s1 ^= state_[1];
                    s2 ^= state_[2];
                    s3 ^= state_[3];
                }

                next();
            }
        }

        state_[0] = s0;
        state_[1] = s1;
        state_[2] = s2;
        state_[3] = s3;
    }

    inline void long_jump_impl(const std::integral_constant<std::uint32_t, 4>&) noexcept
    {
        constexpr std::array<std::uint32_t, 4> jump = {{ UINT32_C(0xb523952e), UINT32_C(0x0b6f099f),
                                                         UINT32_C(0xccf5a0ef), UINT32_C(0x1c580662) }};

        std::uint32_t s0 = 0;
        std::uint32_t s1 = 0;
        std::uint32_t s2 = 0;
        std::uint32_t s3 = 0;

        for (std::size_t i = 0; i < jump.size(); i++)
        {
            for (std::size_t b = 0; b < 32; b++)
            {
                if (jump[i] & UINT32_C(1) << b)
                {
                    s0 ^= state_[0];
                    s1 ^= state_[1];
                    s2 ^= state_[2];
                    s3 ^= state_[3];
                }

                next();
            }
        }

        state_[0] = s0;
        state_[1] = s1;
        state_[2] = s2;
        state_[3] = s3;
    }

public:

    using result_type = OutputType;
    using seed_type = BlockType;

    static constexpr bool has_fixed_range {false};

    /** Seeds the generator using the default seed of boost::random::splitmix64 */
    void seed()
    {
        splitmix64 gen;
        for (auto& i : state_)
        {
            i = static_cast<seed_type>(gen());
        }
    }

    /** Seeds the generator with a user provided seed. */
    void seed(const seed_type value)
    {
        splitmix64 gen(value);
        for (auto& i : state_)
        {
            i = static_cast<seed_type>(gen());
        }
    }

    /**
     * Seeds the generator with 32-bit values produced by @c seq.generate().
     */
    template <typename Sseq, typename std::enable_if<!std::is_convertible<Sseq, seed_type>::value, bool>::type = true>
    void seed(Sseq& seq)
    {
        BOOST_IF_CONSTEXPR (std::is_same<BlockType, std::uint64_t>::value)
        {
            sseq_seed_64(seq);
        }
        else
        {
            sseq_seed_32(seq);
        }
    }

    /** Sets the state of the generator using values from an iterator range. */
    template <typename FIter>
    void seed(FIter first, FIter last)
    {
        static_assert(std::is_integral<typename std::iterator_traits<FIter>::value_type>::value,
                      "Value type must be a built-in integer type" );

        std::size_t offset = 0;
        while (first != last && offset < state_.size())
        {
            state_[offset++] = static_cast<seed_type>(*first++);
        }

        if (offset != state_.size())
        {
            boost::throw_exception(std::invalid_argument("Not enough elements in call to seed."));
        }
    }

    /**
     * Constructs a @c xoshiro and calls @c seed().
     */
    xoshiro_base() { seed(); }

    /** Seeds the generator with a user provided seed. */
    explicit xoshiro_base(const seed_type value)
    {
        seed(value);
    }

    template <typename FIter>
    xoshiro_base(FIter& first, FIter last) { seed(first, last); }

    /**
     * Seeds the generator with 64-bit values produced by @c seq.generate().
     *
     * @xmlnote
     * The copy constructor will always be preferred over
     * the templated constructor.
     * @endxmlnote
     */
    template <typename Sseq, typename std::enable_if<!std::is_convertible<Sseq, xoshiro_base>::value, bool>::type = true>
    explicit xoshiro_base(Sseq& seq)
    {
        seed(seq);
    }

    // Hit all of our rule of 5 explicitly to ensure old platforms work correctly
    ~xoshiro_base() = default;
    xoshiro_base(const xoshiro_base& other) noexcept { state_ = other.state(); }
    xoshiro_base& operator=(const xoshiro_base& other) noexcept { state_ = other.state(); return *this; }
    xoshiro_base(xoshiro_base&& other) noexcept { state_ = other.state(); }
    xoshiro_base& operator=(xoshiro_base&& other) noexcept { state_ = other.state(); return *this; }

    inline result_type next() noexcept
    {
        return static_cast<Derived*>(this)->next();
    }

    /** This is the jump function for the generator. It is equivalent
     *  to 2^128 calls to next(); it can be used to generate 2^128
     *  non-overlapping subsequences for parallel computations. */
    inline void jump() noexcept
    {
        jump_impl(xoshiro_type());
    }

    /** This is the long-jump function for the generator. It is equivalent to
     *  2^192 calls to next(); it can be used to generate 2^64 starting points,
     *  from each of which jump() will generate 2^64 non-overlapping
     *  subsequences for parallel distributed computations. */
    inline void long_jump() noexcept
    {
        long_jump_impl(xoshiro_type());
    }

    /**  Returns the next value of the generator. */
    inline result_type operator()() noexcept
    {
        return next();
    }

    /** Advances the state of the generator by @c z. */
    inline void discard(const std::uint64_t z) noexcept
    {
        for (std::uint64_t i {}; i < z; ++i)
        {
            next();
        }
    }

    /**
     * Returns true if the two generators will produce identical
     * sequences of values.
     */
    inline friend bool operator==(const xoshiro_base& lhs, const xoshiro_base& rhs) noexcept
    {
        return lhs.state_ == rhs.state_;
    }

    /**
     * Returns true if the two generators will produce different
     * sequences of values.
     */
    inline friend bool operator!=(const xoshiro_base& lhs, const xoshiro_base& rhs) noexcept
    {
        return lhs.state_ != rhs.state_;
    }

    /**  Writes a @c xorshiro to a @c std::ostream. */
    template <typename CharT, typename Traits>
    inline friend std::basic_ostream<CharT, Traits>& operator<<(std::basic_ostream<CharT, Traits>& ost,
                                                                const xoshiro_base& e)
    {
        for (std::size_t i {}; i < e.state_.size(); ++i)
        {
            ost << e.state_[i] << ' ';
        }

        return ost;
    }

    /**  Reads a @c xorshiro from a @c std::istream. */
    template <typename CharT, typename Traits>
    inline friend std::basic_istream<CharT, Traits>& operator>>(std::basic_istream<CharT, Traits>& ist,
                                                                xoshiro_base& e)
    {
        for (std::size_t i {}; i < e.state_.size(); ++i)
        {
            ist >> e.state_[i] >> std::ws;
        }

        return ist;
    }

    /** Fills a range with random values */
    template <typename FIter>
    inline void generate(FIter first, FIter last) noexcept
    {
        using iter_type = typename std::iterator_traits<FIter>::value_type;

        while (first != last)
        {
            *first++ = static_cast<iter_type>(next());
        }
    }

    /**
     * Returns the largest value that the @c xorshiro
     * can produce.
     */
    static constexpr result_type (max)() noexcept
    {
        return (std::numeric_limits<result_type>::max)();
    }

    /**
     * Returns the smallest value that the @c xorshiro
     * can produce.
     */
    static constexpr result_type (min)() noexcept
    {
        return (std::numeric_limits<result_type>::min)();
    }

    inline std::array<BlockType, N> state() const noexcept
    {
        return state_;
    }
};

} // namespace detail
} // namespace random
} // namespace boost

#endif