random_rlcg

/*
Copyright (C) 2018-2024 Geoffrey Daniels. https://gpdaniels.com/

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, version 3 of the License only.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program.  If not, see <https://www.gnu.org/licenses/>.
*/

#pragma once
#ifndef GTL_RANDOM_RANDOM_RLCG_HPP
#define GTL_RANDOM_RANDOM_RLCG_HPP

// Summary: A reversible pseudo-random number generator.

#ifndef NDEBUG
#   if defined(_MSC_VER)
#       define __builtin_trap() __debugbreak()
#   endif
/// @brief A simple assert macro to break the program if the random_rlcg is misused.
#   define GTL_RANDOM_RLCG_ASSERT(ASSERTION, MESSAGE) static_cast<void>((ASSERTION) || (__builtin_trap(), 0))
#else
/// @brief At release time the assert macro is implemented as a nop.
#   define GTL_RANDOM_RLCG_ASSERT(ASSERTION, MESSAGE) static_cast<void>(0)
#endif

namespace gtl {
    /// @brief  The random_rlcg class implements a reversible linear congruential generator.
    class random_rlcg final {
    private:
        /// @brief  The modulus used to compute the inverse multiplicand.
        constexpr static const unsigned long long int modulus = 1ull << 63ull;
        static_assert((modulus & (modulus - 1)) == 0, "Modulus must be a power of two.");

        /// @brief  The factor used to progress the state by multiplication.
        constexpr static const unsigned long long int multiplicand_forward = 0x5851F42D4C957F2Dull;

        /// @brief  The inverse of the multiplicand template parameter, used for inverse steps.
        constexpr static const unsigned long long int multiplicand_inverse = [](unsigned long long int a, unsigned long long int b) {
            // This is an implementation of the extended euclidean algorithm to determine the first bezout coefficient (x).
            // of the equation defined as "ax + by = gcd(a, b)" where "gcd(a, b)" is the greatest common divisor of a and b.
            unsigned long long int bezout_coefficients[2] = { 1, 0 };
            while (b) {
                const unsigned long long int quotient = a / b;
                const unsigned long long int remainder = bezout_coefficients[0] - quotient * bezout_coefficients[1];
                bezout_coefficients[0] = bezout_coefficients[1];
                bezout_coefficients[1] = remainder;
                const unsigned long long int m = a % b;
                a = b;
                b = m;
            }
            return bezout_coefficients[0];
        }(multiplicand_forward, modulus);

    private:
        /// @brief  Current state of the random number generator.
        unsigned long long int state;

        /// @brief  Increment added to the state during pseudo-random-number generation.
        unsigned long long int increment;

    public:
        /// @brief  Defaulted destructor.
        ~random_rlcg() = default;

        /// @brief  Empty constructor.
        random_rlcg()
            : state(0x853C49E6748FEA9Bull)
            , increment(0xDA3E39CB94B95BDBull) {
        }

        /// @brief  Defaulted copy constructor.
        random_rlcg(const random_rlcg&) = default;

        /// @brief  Defaulted move constructor.
        random_rlcg(random_rlcg&&) = default;

        /// @brief  Defaulted copy assignment.
        random_rlcg& operator=(const random_rlcg&) = default;

        /// @brief  Defaulted move assignment.
        random_rlcg& operator=(random_rlcg&&) = default;

        /// @brief  Constructor with 64 bit seed.
        /// @param  seed_value Seed value used to initialise the state.
        random_rlcg(unsigned long long int seed_value) {
            this->seed(seed_value);
        }

        /// @brief  Initialise the state from a 64 bit seed.
        /// @param  seed_value Seed value used to initialise the state.
        void seed(unsigned long long int seed_value) {
            this->state = 0;
            this->increment = (seed_value << 1) | 1;
            this->get_random_raw_forward();
            this->state += seed_value;
            this->get_random_raw_forward();
        }

        /// @brief  Get the random number from the generator and advance the state forwards.
        /// @return A pseudo-random number.
        unsigned int get_random_raw_forward() {
            // Save current state for output calculation.
            unsigned long long int state_previous = this->state;
            // Advance internal state forward.
            this->state = ((multiplicand_forward * state_previous) + this->increment) & (modulus - 1);
            // Calculate output function.
            unsigned int state_shift_xor_shift = static_cast<unsigned int>(((state_previous >> 18u) ^ state_previous) >> 27u);
            int rotation = state_previous >> 59u;
            return (state_shift_xor_shift >> rotation) | (state_shift_xor_shift << ((-rotation) & 31));
        }

        /// @brief  Get the random number from the generator and advance the state inverse.
        /// @return A pseudo-random number.
        unsigned int get_random_raw_inverse() {
            // Advance internal state inverse.
            this->state = (multiplicand_inverse * (this->state - increment)) & (modulus - 1);
            // Calculate output function.
            unsigned int state_shift_xor_shift = static_cast<unsigned int>(((this->state >> 18u) ^ this->state) >> 27u);
            int rotation = this->state >> 59u;
            return (state_shift_xor_shift >> rotation) | (state_shift_xor_shift << ((-rotation) & 31));
        }

    public:
        /// @brief  Get a random number in the open interval 0 < value < 1.
        /// @note  This function will step the state forward.
        /// @return A pseudo-random number.
        double get_random_exclusive_forward() {
            return (static_cast<double>(this->get_random_raw_forward()) + 0.5) * (1.0 / static_cast<double>(1ull << 32));
        }

        /// @brief  Get a random number in the open interval 0 < value < 1.
        /// @note  This function will step the state inverse.
        /// @return A pseudo-random number.
        double get_random_exclusive_inverse() {
            return (static_cast<double>(this->get_random_raw_inverse()) + 0.5) * (1.0 / static_cast<double>(1ull << 32));
        }

        /// @brief  Get a random number in the half-open interval 0 <= value < 1.
        /// @note  This function will step the state forward.
        /// @return A pseudo-random number.
        double get_random_exclusive_top_forward() {
            return static_cast<double>(this->get_random_raw_forward()) * (1.0 / static_cast<double>(1ull << 32));
        }

        /// @brief  Get a random number in the half-open interval 0 <= value < 1.
        /// @note  This function will step the state inverse.
        /// @return A pseudo-random number.
        double get_random_exclusive_top_inverse() {
            return static_cast<double>(this->get_random_raw_inverse()) * (1.0 / static_cast<double>(1ull << 32));
        }

        /// @brief  Get a random number in the closed interval 0 <= value <= 1.
        /// @note  This function will step the state forward.
        /// @return A pseudo-random number.
        double get_random_inclusive_forward() {
            return static_cast<double>(this->get_random_raw_forward()) * (1.0 / static_cast<double>((1ull << 32) - 1));
        }

        /// @brief  Get a random number in the closed interval 0 <= value <= 1.
        /// @note  This function will step the state inverse.
        /// @return A pseudo-random number.
        double get_random_inclusive_inverse() {
            return static_cast<double>(this->get_random_raw_inverse()) * (1.0 / static_cast<double>((1ull << 32) - 1));
        }

        /// @brief  Get a random number between two bounds.
        /// @note  This function will step the state forward.
        /// @param  inclusive_min Minimum number that can be returned.
        /// @param  inclusive_max Max number than can be returned.
        /// @return A pseudo-random number.
        unsigned int get_random_forward(unsigned int inclusive_min, unsigned int inclusive_max)  {
            GTL_RANDOM_RLCG_ASSERT(inclusive_min < inclusive_max, "Minimum bound must be lower than maximum bound.");
            return (this->get_random_raw_forward() % (1 + inclusive_max - inclusive_min)) + inclusive_min;
        }

        /// @brief  Get a random number between two bounds.
        /// @note  This function will step the state inverse.
        /// @param  inclusive_min Minimum number that can be returned.
        /// @param  inclusive_max Max number than can be returned.
        /// @return A pseudo-random number.
        unsigned int get_random_inverse(unsigned int inclusive_min, unsigned int inclusive_max)  {
            GTL_RANDOM_RLCG_ASSERT(inclusive_min < inclusive_max, "Minimum bound must be lower than maximum bound.");
            return (this->get_random_raw_inverse() % (1 + inclusive_max - inclusive_min)) + inclusive_min;
        }

        /// @brief  Get a random number between two bounds.
        /// @note  This function will step the state forward.
        /// @param  inclusive_min Minimum number that can be returned.
        /// @param  inclusive_max Max number than can be returned.
        /// @return A pseudo-random number.
        double get_random_forward(double inclusive_min, double inclusive_max) {
            GTL_RANDOM_RLCG_ASSERT(inclusive_min < inclusive_max, "Minimum bound must be lower than maximum bound.");
            return (this->get_random_inclusive_forward() * (inclusive_max - inclusive_min)) + inclusive_min;
        }

        /// @brief  Get a random number between two bounds.
        /// @note  This function will step the state inverse.
        /// @param  inclusive_min Minimum number that can be returned.
        /// @param  inclusive_max Max number than can be returned.
        /// @return A pseudo-random number.
        double get_random_inverse(double inclusive_min, double inclusive_max) {
            GTL_RANDOM_RLCG_ASSERT(inclusive_min < inclusive_max, "Minimum bound must be lower than maximum bound.");
            return (this->get_random_inclusive_inverse() * (inclusive_max - inclusive_min)) + inclusive_min;
        }
    };
}

#undef GTL_RANDOM_RLCG_ASSERT

#endif // GTL_RANDOM_RANDOM_RLCG_HPP