PRNGDistribution.hpp

// INFINITY_API_PUBLIC
 
#pragma once
 
#include <Infinity/api.h>
#include <Infinity/Types/Math/Math.hpp>
 
#include <cstdint>
#include <limits>
#include <type_traits>
#include <cmath>
#include <vector>
#include <algorithm>
#include <initializer_list>
 
namespace Infinity::Engine {
 
    template<typename RealType = double>
    class INFINITY_API_TEMPLATE UniformRealDistribution
    {
        static_assert(std::is_floating_point<RealType>::value,
                    "RealType must be a floating point type");
 
        public:
            using result_type = RealType; 
 
            struct param_type
            {
                using distribution_type = UniformRealDistribution<RealType>;
 
                param_type() : param_type(RealType(0)) { }
 
                explicit param_type(RealType a, RealType b = RealType(1))
                    : m_a(a), m_b(b) {
                    #ifdef _DEBUG
                    assert(a <= b);
                    #endif
                }
 
                RealType a() const { return m_a; }
 
                RealType b() const { return m_b; }
 
                friend bool operator==(const param_type& lhs, const param_type& rhs)
                {
                    return lhs.m_a == rhs.m_a && lhs.m_b == rhs.m_b;
                }
 
                friend bool operator!=(const param_type& lhs, const param_type& rhs)
                {
                    return !(lhs == rhs);
                }
 
            private:
                RealType m_a;
                RealType m_b;
            };
 
            // Constructors
            
            UniformRealDistribution() : UniformRealDistribution(RealType(0)) { }
 
            explicit UniformRealDistribution(RealType a, RealType b = RealType(1))
                : m_param(a, b) { }
 
            explicit UniformRealDistribution(const param_type& param)
                : m_param(param) { }
 
            // Reset (no-op for uniform distribution)
            
            void reset() { }
 
            // Parameter access
            
            RealType a() const { return m_param.a(); }
 
            RealType b() const { return m_param.b(); }
 
            param_type param() const { return m_param; }
 
            void param(const param_type& param) { m_param = param; }
 
            result_type min() const { return a(); }
 
            result_type max() const { return b(); }
 
            // Generation
            
            template<typename Generator>
            result_type operator()(Generator& g)
            {
                return (*this)(g, m_param);
            }
 
            template<typename Generator>
            result_type operator()(Generator& g, const param_type& param)
            {
                uint32_t random_bits;
                using generator_result = typename Generator::result_type;
                
                if constexpr (sizeof(generator_result) >= sizeof(uint32_t))
                {
                    random_bits = static_cast<uint32_t>(g());
                } else {
                    constexpr int bits_per_call = sizeof(generator_result) * 8;
                    constexpr int calls_needed = (32 + bits_per_call - 1) / bits_per_call;
                    random_bits = 0;
                    for (int i = 0; i < calls_needed; ++i)
                    {
                        random_bits = (random_bits << bits_per_call) | static_cast<uint32_t>(g());
                    }
                }
                
                RealType unit = uint32_to_real(random_bits);
                return param.a() + unit * (param.b() - param.a());
            }
 
            // Comparison
            
            friend bool operator==(const UniformRealDistribution& lhs,
                                const UniformRealDistribution& rhs)
            {
                return lhs.m_param == rhs.m_param;
            }
 
            friend bool operator!=(const UniformRealDistribution& lhs,
                                const UniformRealDistribution& rhs)
            {
                return !(lhs == rhs);
            }
 
        private:
            param_type m_param;
 
            static result_type uint32_to_real(uint32_t value)
            {
                constexpr RealType scale = RealType(1.0) / RealType(4294967296.0);
                return static_cast<RealType>(value) * scale;
            }
    };
 
    template<typename IntType = int>
    class INFINITY_API_TEMPLATE UniformIntDistribution
    {
        static_assert(std::is_integral<IntType>::value,
                    "IntType must be an integral type");
 
        public:
            using result_type = IntType; 
 
            struct param_type
            {
                using distribution_type = UniformIntDistribution<IntType>;
 
                param_type() : param_type(0) { }
 
                explicit param_type(IntType a, IntType b = std::numeric_limits<IntType>::max())
                    : m_a(a), m_b(b) {
                    #ifdef _DEBUG
                    assert(a <= b);
                    #endif
                }
 
                IntType a() const { return m_a; }
 
                IntType b() const { return m_b; }
 
                friend bool operator==(const param_type& lhs, const param_type& rhs)
                {
                    return lhs.m_a == rhs.m_a && lhs.m_b == rhs.m_b;
                }
 
                friend bool operator!=(const param_type& lhs, const param_type& rhs)
                {
                    return !(lhs == rhs);
                }
 
            private:
                IntType m_a;
                IntType m_b;
            };
 
            // Constructors
            
            UniformIntDistribution() : UniformIntDistribution(0) { }
 
            explicit UniformIntDistribution(IntType a, IntType b = std::numeric_limits<IntType>::max())
                : m_param(a, b) { }
 
            explicit UniformIntDistribution(const param_type& param)
                : m_param(param) { }
 
            // Reset (no-op for uniform distribution)
            
            void reset() { }
 
            // Parameter access
            
            IntType a() const { return m_param.a(); }
 
            IntType b() const { return m_param.b(); }
 
            param_type param() const { return m_param; }
 
            void param(const param_type& param) { m_param = param; }
 
            result_type min() const { return a(); }
 
            result_type max() const { return b(); }
 
            // Generation
            
            template<typename Generator>
            result_type operator()(Generator& g)
            {
                return (*this)(g, m_param);
            }
 
            template<typename Generator>
            result_type operator()(Generator& g, const param_type& param)
            {
                using UIntType = typename std::make_unsigned<IntType>::type;
                
                const IntType a = param.a();
                const IntType b = param.b();
                
                if (a == b) return a;
                
                const UIntType range = static_cast<UIntType>(b) - static_cast<UIntType>(a);
                
                using generator_result = typename Generator::result_type;
                uint32_t x;
                if constexpr (sizeof(generator_result) >= sizeof(uint32_t))
                {
                    x = static_cast<uint32_t>(g());
                } else {
                    constexpr int bits_per_call = sizeof(generator_result) * 8;
                    constexpr int calls_needed = (32 + bits_per_call - 1) / bits_per_call;
                    x = 0;
                    for (int i = 0; i < calls_needed; ++i)
                    {
                        x = (x << bits_per_call) | static_cast<uint32_t>(g());
                    }
                }
                
                const UIntType range_plus_one = range + 1;
                // Fast path for power-of-2 ranges
                if ((range_plus_one & range) == 0)
                {
                    return static_cast<IntType>(static_cast<UIntType>(a) + (x & range));
                }
                
                // Lemire's algorithm for general ranges
                uint64_t m = static_cast<uint64_t>(x) * static_cast<uint64_t>(range_plus_one);
                uint32_t l = static_cast<uint32_t>(m);
                
                if (l < range_plus_one)
                {
                    uint32_t t = -range_plus_one % range_plus_one;
                    while (l < t)
                    {
                        x = static_cast<uint32_t>(g());
                        m = static_cast<uint64_t>(x) * static_cast<uint64_t>(range_plus_one);
                        l = static_cast<uint32_t>(m);
                    }
                }
                
                return static_cast<IntType>(static_cast<UIntType>(a) + (m >> 32));
            }
 
            // Comparison
            
            friend bool operator==(const UniformIntDistribution& lhs,
                                const UniformIntDistribution& rhs)
            {
                return lhs.m_param == rhs.m_param;
            }
 
            friend bool operator!=(const UniformIntDistribution& lhs,
                                const UniformIntDistribution& rhs)
            {
                return !(lhs == rhs);
            }
 
        private:
            param_type m_param;
    };
 
    template<typename RealType = double>
    class INFINITY_API_TEMPLATE NormalDistribution
    {
        static_assert(std::is_floating_point<RealType>::value,
                    "RealType must be a floating point type");
 
        public:
            using result_type = RealType; 
 
            struct param_type {
                using distribution_type = NormalDistribution<RealType>;
 
                param_type() : param_type(RealType(0)) { }
 
                explicit param_type(RealType mean, RealType stddev = RealType(1))
                    : m_mean(mean), m_stddev(stddev) { }
 
                RealType mean() const { return m_mean; }
                
                RealType stddev() const { return m_stddev; }
 
                friend bool operator==(const param_type& lhs, const param_type& rhs)
                {
                    return lhs.m_mean == rhs.m_mean && lhs.m_stddev == rhs.m_stddev;
                }
 
                friend bool operator!=(const param_type& lhs, const param_type& rhs)
                {
                    return !(lhs == rhs);
                }
 
            private:
                RealType m_mean;
                RealType m_stddev;
            };
 
            // Constructors
            
            NormalDistribution() : NormalDistribution(RealType(0)) {}
            
            explicit NormalDistribution(RealType mean, RealType stddev = RealType(1))
                : m_param(mean, stddev), m_saved_available(false) {}
                
            explicit NormalDistribution(const param_type& param)
                : m_param(param), m_saved_available(false) {}
 
            void reset() { m_saved_available = false; }
 
            // Parameter access
            
            RealType mean() const { return m_param.mean(); }
            
            RealType stddev() const { return m_param.stddev(); }
            
            param_type param() const { return m_param; }
            
            void param(const param_type& param) { m_param = param; m_saved_available = false; }
            
            result_type min() const { return std::numeric_limits<result_type>::lowest(); }
            
            result_type max() const { return std::numeric_limits<result_type>::max(); }
 
            // Generation
            
            template<typename Generator>
            result_type operator()(Generator& g) { return (*this)(g, m_param); }
 
            template<typename Generator>
            result_type operator()(Generator& g, const param_type& param)
            {
                // Use cached value if available
                if (m_saved_available)
                {
                    m_saved_available = false;
                    return param.mean() + m_saved * param.stddev();
                }
                
                UniformRealDistribution<RealType> uniform(RealType(0), RealType(1));
                
                // Generate two uniform random values
                RealType u1, u2;
                do {
                    u1 = uniform(g);
                } while (u1 == RealType(0));  // Avoid log(0)
                
                u2 = uniform(g);
                
                // Box-Muller transform
                const RealType magnitude = std::sqrt(-RealType(2) * std::log(u1));
                const RealType z0 = magnitude * std::cos(RealType(2) * RealType(Infinity::Types::Math::PI) * u2);
                const RealType z1 = magnitude * std::sin(RealType(2) * RealType(Infinity::Types::Math::PI) * u2);
                
                // Cache second value for next call
                m_saved = z1;
                m_saved_available = true;
                
                // Return first value scaled to desired mean and stddev
                return param.mean() + z0 * param.stddev();
            }
 
            // Comparison
            
            friend bool operator==(const NormalDistribution& lhs, const NormalDistribution& rhs)
            {
                return lhs.m_param == rhs.m_param && lhs.m_saved_available == rhs.m_saved_available &&
                    (!lhs.m_saved_available || lhs.m_saved == rhs.m_saved);
            }
 
            friend bool operator!=(const NormalDistribution& lhs, const NormalDistribution& rhs)
            {
                return !(lhs == rhs);
            }
 
        private:
            param_type m_param;
            result_type m_saved = RealType(0);      
            bool m_saved_available = false;          
    };
 
    template<typename RealType = double>
    class INFINITY_API_TEMPLATE LognormalDistribution {
        static_assert(std::is_floating_point<RealType>::value,
                    "RealType must be a floating point type");
 
        public:
            using result_type = RealType; 
 
            struct param_type
            {
                using distribution_type = LognormalDistribution<RealType>;
 
                param_type() : param_type(RealType(0)) { }
                
                explicit param_type(RealType m, RealType s = RealType(1)) : m_m(m), m_s(s) { }
 
                RealType m() const { return m_m; }
                
                RealType s() const { return m_s; }
 
                friend bool operator==(const param_type& lhs, const param_type& rhs)
                {
                    return lhs.m_m == rhs.m_m && lhs.m_s == rhs.m_s;
                }
                friend bool operator!=(const param_type& lhs, const param_type& rhs)
                {
                    return !(lhs == rhs);
                }
 
            private:
                RealType m_m, m_s;
            };
 
            // Constructors
            
            LognormalDistribution() : LognormalDistribution(RealType(0)) { }
            
            explicit LognormalDistribution(RealType m, RealType s = RealType(1))
                : m_param(m, s), m_nd(m, s) { }
                
            explicit LognormalDistribution(const param_type& param)
                : m_param(param), m_nd(param.m(), param.s()) { }
 
            void reset() { m_nd.reset(); }
 
            // Parameter access
            
            RealType m() const { return m_param.m(); }
            
            RealType s() const { return m_param.s(); }
            
            param_type param() const { return m_param; }
            
            void param(const param_type& param) { m_param = param; }
            
            result_type min() const { return RealType(0); }
            
            result_type max() const { return std::numeric_limits<result_type>::max(); }
 
            // Generation
            
            template<typename Generator>
            result_type operator()(Generator& g) { return std::exp(m_nd(g)); }
 
            template<typename Generator>
            result_type operator()(Generator& g, const param_type& p)
            {
                NormalDistribution<RealType> nd(p.m(), p.s());
                return std::exp(nd(g));
            }
 
            // Comparison
            
            friend bool operator==(const LognormalDistribution& lhs, const LognormalDistribution& rhs)
            {
                return lhs.m_param == rhs.m_param && lhs.m_nd == rhs.m_nd;
            }
            friend bool operator!=(const LognormalDistribution& lhs, const LognormalDistribution& rhs)
            {
                return !(lhs == rhs);
            }
 
        private:
            param_type m_param;
            NormalDistribution<result_type> m_nd;  
    };
 
    template<typename RealType = double>
    class INFINITY_API_TEMPLATE CauchyDistribution
    {
        static_assert(std::is_floating_point<RealType>::value,
                    "RealType must be a floating point type");
 
        public:
            using result_type = RealType; 
 
            struct param_type {
                using distribution_type = CauchyDistribution<RealType>;
 
                param_type() : param_type(RealType(0)) {}
                
                explicit param_type(RealType a, RealType b = RealType(1)) : m_a(a), m_b(b) { }
 
                RealType a() const { return m_a; }
                
                RealType b() const { return m_b; }
 
                friend bool operator==(const param_type& lhs, const param_type& rhs)
                {
                    return lhs.m_a == rhs.m_a && lhs.m_b == rhs.m_b;
                }
                friend bool operator!=(const param_type& lhs, const param_type& rhs)
                {
                    return !(lhs == rhs);
                }
 
            private:
                RealType m_a, m_b;
            };
 
            // Constructors
            
            CauchyDistribution() : CauchyDistribution(RealType(0)) { }
            
            explicit CauchyDistribution(RealType a, RealType b = RealType(1)) : m_param(a, b) { }
            
            explicit CauchyDistribution(const param_type& param) : m_param(param) { }
 
            void reset() { }
 
            // Parameter access
            
            RealType a() const { return m_param.a(); }
            
            RealType b() const { return m_param.b(); }
            
            param_type param() const { return m_param; }
            
            void param(const param_type& param) { m_param = param; }
            
            result_type min() const { return std::numeric_limits<result_type>::lowest(); }
            
            result_type max() const { return std::numeric_limits<result_type>::max(); }
 
            // Generation
            
            template<typename Generator>
            result_type operator()(Generator& g) { return (*this)(g, m_param); }
 
            template<typename Generator>
            result_type operator()(Generator& g, const param_type& param)
            {
                UniformRealDistribution<RealType> uniform(RealType(0), RealType(1));
                const RealType u = uniform(g);
                return param.a() + param.b() * std::tan(RealType(Infinity::Types::Math::PI) * (u - RealType(0.5)));
            }
 
            // Comparison
            
            friend bool operator==(const CauchyDistribution& lhs, const CauchyDistribution& rhs)
            {
                return lhs.m_param == rhs.m_param;
            }
            friend bool operator!=(const CauchyDistribution& lhs, const CauchyDistribution& rhs)
            {
                return !(lhs == rhs);
            }
 
        private:
            param_type m_param;
    };
 
    template<typename RealType = double>
    class INFINITY_API_TEMPLATE GammaDistribution
    {
        static_assert(std::is_floating_point<RealType>::value,
                    "RealType must be a floating point type");
 
        public:
            using result_type = RealType; 
 
            struct param_type
            {
                using distribution_type = GammaDistribution<RealType>;
 
                param_type() : param_type(RealType(1)) {}
                
                explicit param_type(RealType alpha, RealType beta = RealType(1))
                    : m_alpha(alpha), m_beta(beta) {}
 
                RealType alpha() const { return m_alpha; }
                
                RealType beta() const { return m_beta; }
 
                friend bool operator==(const param_type& lhs, const param_type& rhs)
                {
                    return lhs.m_alpha == rhs.m_alpha && lhs.m_beta == rhs.m_beta;
                }
                friend bool operator!=(const param_type& lhs, const param_type& rhs)
                {
                    return !(lhs == rhs);
                }
 
            private:
                RealType m_alpha, m_beta;
            };
 
            // Constructors
            
            GammaDistribution() : GammaDistribution(RealType(1)) { }
            
            explicit GammaDistribution(RealType alpha, RealType beta = RealType(1))
                : m_param(alpha, beta), m_nd() { }
                
            explicit GammaDistribution(const param_type& param)
                : m_param(param), m_nd() { }
 
            void reset() { m_nd.reset(); }
 
            // Parameter access
            
            RealType alpha() const { return m_param.alpha(); }
            
            RealType beta() const { return m_param.beta(); }
            
            param_type param() const { return m_param; }
            
            void param(const param_type& param) { m_param = param; }
            
            result_type min() const { return RealType(0); }
            
            result_type max() const { return std::numeric_limits<result_type>::max(); }
 
            // Generation
            
            template<typename Generator>
            result_type operator()(Generator& g) { return (*this)(g, m_param); }
 
            template<typename Generator>
            result_type operator()(Generator& g, const param_type& param)
            {
                const RealType alpha = param.alpha();
                const RealType beta = param.beta();
                
                if (alpha < RealType(1))
                {
                    // Transformation method for alpha < 1
                    // Generate Gamma(alpha+1) and scale by U^(1/alpha)
                    RealType result = (*this)(g, param_type(alpha + RealType(1), RealType(1)));
                    UniformRealDistribution<RealType> uniform(RealType(0), RealType(1));
                    return result * std::pow(uniform(g), RealType(1) / alpha) * beta;
                }
                
                // Marsaglia & Tsang (2000) method for alpha >= 1
                const RealType d = alpha - RealType(1) / RealType(3);
                const RealType c = RealType(1) / std::sqrt(RealType(9) * d);
                
                UniformRealDistribution<RealType> uniform(RealType(0), RealType(1));
                
                while (true)
                {
                    RealType x, v;
                    // Generate candidate from transformed normal
                    do {
                        x = m_nd(g);
                        v = RealType(1) + c * x;
                    } while (v <= RealType(0));  // Ensure v > 0
                    
                    v = v * v * v;
                    const RealType u = uniform(g);
                    const RealType x_sq = x * x;
                    
                    // Squeeze test (fast acceptance)
                    if (u < RealType(1) - RealType(0.0331) * x_sq * x_sq)
                    {
                        return d * v * beta;
                    }
                    
                    // Log test (slower but exact)
                    if (std::log(u) < RealType(0.5) * x_sq + d * (RealType(1) - v + std::log(v)))
                    {
                        return d * v * beta;
                    }
                    
                    // Rejection - loop continues
                }
            }
 
            // Comparison
            
            friend bool operator==(const GammaDistribution& lhs, const GammaDistribution& rhs)
            {
                return lhs.m_param == rhs.m_param && lhs.m_nd == rhs.m_nd;
            }
            friend bool operator!=(const GammaDistribution& lhs, const GammaDistribution& rhs)
            {
                return !(lhs == rhs);
            }
 
        private:
            param_type m_param;
            NormalDistribution<result_type> m_nd;  
    };
 
    template<typename RealType = double>
    class INFINITY_API_TEMPLATE ChiSquaredDistribution
    {
        static_assert(std::is_floating_point<RealType>::value,
                    "RealType must be a floating point type");
 
        public:
            using result_type = RealType; 
 
            struct param_type {
                using distribution_type = ChiSquaredDistribution<RealType>;
 
                param_type() : param_type(RealType(1)) { }
                
                explicit param_type(RealType n) : m_n(n) { }
 
                RealType n() const { return m_n; }
 
                friend bool operator==(const param_type& lhs, const param_type& rhs)
                {
                    return lhs.m_n == rhs.m_n;
                }
                friend bool operator!=(const param_type& lhs, const param_type& rhs)
                {
                    return !(lhs == rhs);
                }
 
            private:
                RealType m_n;
            };
 
            // Constructors
            
            ChiSquaredDistribution() : ChiSquaredDistribution(RealType(1)) { }
            
            explicit ChiSquaredDistribution(RealType n)
                : m_param(n)
            {
                m_gd = new GammaDistribution<result_type>(n / RealType(2), RealType(2));
            }
            
            explicit ChiSquaredDistribution(const param_type& param)
                : m_param(param)
            {
                m_gd = new GammaDistribution<result_type>(param.n() / RealType(2), RealType(2));
            }
 
            void reset()
            {
                if (m_gd) m_gd->reset();
            }
 
            // Parameter access
            
            RealType n() const { return m_param.n(); }
            
            param_type param() const { return m_param; }
            
            void param(const param_type& param)
            {
                m_param = param;
                delete m_gd;
                m_gd = new GammaDistribution<result_type>(param.n() / RealType(2), RealType(2));
            }
            
            result_type min() const { return RealType(0); }
            
            result_type max() const { return std::numeric_limits<result_type>::max(); }
 
            // Generation
            
            template<typename Generator>
            result_type operator()(Generator& g)
            {
                return RealType(2) * (*m_gd)(g);
            }
 
            template<typename Generator>
            result_type operator()(Generator& g, const param_type& param)
            {
                GammaDistribution<result_type> gd(param.n() / RealType(2), RealType(2));
                return RealType(2) * gd(g);
            }
 
            // Comparison
            
            friend bool operator==(const ChiSquaredDistribution& lhs, const ChiSquaredDistribution& rhs)
            {
                if (lhs.m_param != rhs.m_param) return false;
                if ((lhs.m_gd == nullptr) != (rhs.m_gd == nullptr)) return false;
                if (lhs.m_gd && rhs.m_gd && *lhs.m_gd != *rhs.m_gd) return false;
                return true;
            }
            friend bool operator!=(const ChiSquaredDistribution& lhs, const ChiSquaredDistribution& rhs)
            {
                return !(lhs == rhs);
            }
 
        private:
            param_type m_param;
            class GammaDistribution<result_type>* m_gd;  
    };
 
    template<typename RealType = double>
    class INFINITY_API_TEMPLATE FisherFDistribution
    {
        static_assert(std::is_floating_point<RealType>::value,
                    "RealType must be a floating point type");
 
        public:
            using result_type = RealType; 
 
            struct param_type
            {
                using distribution_type = FisherFDistribution<RealType>;
 
                param_type() : param_type(RealType(1)) { }
                
                explicit param_type(RealType m, RealType n = RealType(1)) : m_m(m), m_n(n) { }
 
                RealType m() const { return m_m; }
                
                RealType n() const { return m_n; }
 
                friend bool operator==(const param_type& lhs, const param_type& rhs)
                {
                    return lhs.m_m == rhs.m_m && lhs.m_n == rhs.m_n;
                }
                friend bool operator!=(const param_type& lhs, const param_type& rhs)
                {
                    return !(lhs == rhs);
                }
 
            private:
                RealType m_m, m_n;
            };
 
            // Constructors
            
            FisherFDistribution() : FisherFDistribution(RealType(1)) { }
            
            explicit FisherFDistribution(RealType m, RealType n = RealType(1))
                : m_param(m, n)
            {
                m_gd_x = new GammaDistribution<result_type>(m / RealType(2), RealType(2));
                m_gd_y = new GammaDistribution<result_type>(n / RealType(2), RealType(2));
            }
            
            explicit FisherFDistribution(const param_type& param)
                : m_param(param)
            {
                m_gd_x = new GammaDistribution<result_type>(param.m() / RealType(2), RealType(2));
                m_gd_y = new GammaDistribution<result_type>(param.n() / RealType(2), RealType(2));
            }
 
            void reset()
            {
                if (m_gd_x) m_gd_x->reset();
                if (m_gd_y) m_gd_y->reset();
            }
 
            // Parameter access
            
            RealType m() const { return m_param.m(); }
            
            RealType n() const { return m_param.n(); }
            
            param_type param() const { return m_param; }
            
            void param(const param_type& param)
            {
                m_param = param;
                delete m_gd_x;
                delete m_gd_y;
                m_gd_x = new GammaDistribution<result_type>(param.m() / RealType(2), RealType(2));
                m_gd_y = new GammaDistribution<result_type>(param.n() / RealType(2), RealType(2));
            }
            
            result_type min() const { return RealType(0); }
            
            result_type max() const { return std::numeric_limits<result_type>::max(); }
 
            // Generation
            
            template<typename Generator>
            result_type operator()(Generator& g)
            {
                const RealType x = (*m_gd_x)(g);
                const RealType y = (*m_gd_y)(g);
                return (x * m_param.n()) / (y * m_param.m());
            }
 
            template<typename Generator>
            result_type operator()(Generator& g, const param_type& param)
            {
                GammaDistribution<result_type> gd_x(param.m() / RealType(2), RealType(2));
                GammaDistribution<result_type> gd_y(param.n() / RealType(2), RealType(2));
                const RealType x = gd_x(g);
                const RealType y = gd_y(g);
                return (x * param.n()) / (y * param.m());
            }
 
            // Comparison
            
            friend bool operator==(const FisherFDistribution& lhs, const FisherFDistribution& rhs)
            {
                if (lhs.m_param != rhs.m_param) return false;
                if ((lhs.m_gd_x == nullptr) != (rhs.m_gd_x == nullptr)) return false;
                if ((lhs.m_gd_y == nullptr) != (rhs.m_gd_y == nullptr)) return false;
                if (lhs.m_gd_x && rhs.m_gd_x && *lhs.m_gd_x != *rhs.m_gd_x) return false;
                if (lhs.m_gd_y && rhs.m_gd_y && *lhs.m_gd_y != *rhs.m_gd_y) return false;
                return true;
            }
            friend bool operator!=(const FisherFDistribution& lhs, const FisherFDistribution& rhs)
            {
                return !(lhs == rhs);
            }
 
        private:
            param_type m_param;
            class GammaDistribution<result_type>* m_gd_x;  
            class GammaDistribution<result_type>* m_gd_y;  
    };
 
    template<typename RealType = double>
    class INFINITY_API_TEMPLATE StudentTDistribution
    {
        static_assert(std::is_floating_point<RealType>::value,
                    "RealType must be a floating point type");
 
        public:
            using result_type = RealType; 
 
            struct param_type
            {
                using distribution_type = StudentTDistribution<RealType>;
 
                param_type() : param_type(RealType(1)) { }
                
                explicit param_type(RealType n) : m_n(n) { }
 
                RealType n() const { return m_n; }
 
                friend bool operator==(const param_type& lhs, const param_type& rhs)
                {
                    return lhs.m_n == rhs.m_n;
                }
                friend bool operator!=(const param_type& lhs, const param_type& rhs)
                {
                    return !(lhs == rhs);
                }
 
            private:
                RealType m_n;
            };
 
            // Constructors
            
            StudentTDistribution() : StudentTDistribution(RealType(1)) { }
            
            explicit StudentTDistribution(RealType n)
                : m_param(n), m_nd()
            {
                m_gd = new GammaDistribution<result_type>(n / RealType(2), RealType(2));
            }
            
            explicit StudentTDistribution(const param_type& param)
                : m_param(param), m_nd()
            {
                m_gd = new GammaDistribution<result_type>(param.n() / RealType(2), RealType(2));
            }
 
            void reset()
            {
                m_nd.reset();
                if (m_gd) m_gd->reset();
            }
 
            // Parameter access
            
            RealType n() const { return m_param.n(); }
            
            param_type param() const { return m_param; }
            
            void param(const param_type& param)
            {
                m_param = param;
                delete m_gd;
                m_gd = new GammaDistribution<result_type>(param.n() / RealType(2), RealType(2));
            }
            
            result_type min() const { return std::numeric_limits<result_type>::lowest(); }
            
            result_type max() const { return std::numeric_limits<result_type>::max(); }
 
            // Generation
            
            template<typename Generator>
            result_type operator()(Generator& g)
            {
                const RealType z = m_nd(g);
                const RealType v = (*m_gd)(g);
                return z * std::sqrt(m_param.n() / v);
            }
 
            template<typename Generator>
            result_type operator()(Generator& g, const param_type& param)
            {
                NormalDistribution<result_type> nd;
                GammaDistribution<result_type> gd(param.n() / RealType(2), RealType(2));
                const RealType z = nd(g);
                const RealType v = gd(g);
                return z * std::sqrt(param.n() / v);
            }
 
            // Comparison
            
            friend bool operator==(const StudentTDistribution& lhs, const StudentTDistribution& rhs)
            {
                if (lhs.m_param != rhs.m_param) return false;
                if (lhs.m_nd != rhs.m_nd) return false;
                if ((lhs.m_gd == nullptr) != (rhs.m_gd == nullptr)) return false;
                if (lhs.m_gd && rhs.m_gd && *lhs.m_gd != *rhs.m_gd) return false;
                return true;
            }
            friend bool operator!=(const StudentTDistribution& lhs, const StudentTDistribution& rhs)
            {
                return !(lhs == rhs);
            }
 
        private:
            param_type m_param;
            NormalDistribution<result_type> m_nd;              
            class GammaDistribution<result_type>* m_gd;       
    };
 
    template<typename RealType = double>
    class INFINITY_API_TEMPLATE ExponentialDistribution
    {
        static_assert(std::is_floating_point<RealType>::value,
                    "RealType must be a floating point type");
 
        public:
            using result_type = RealType; 
 
            struct param_type
            {
                using distribution_type = ExponentialDistribution<RealType>;
 
                param_type() : param_type(RealType(1)) { }
                
                explicit param_type(RealType lambda) : m_lambda(lambda) { }
 
                RealType lambda() const { return m_lambda; }
 
                friend bool operator==(const param_type& lhs, const param_type& rhs)
                {
                    return lhs.m_lambda == rhs.m_lambda;
                }
                friend bool operator!=(const param_type& lhs, const param_type& rhs)
                {
                    return !(lhs == rhs);
                }
 
            private:
                RealType m_lambda;
            };
 
            // Constructors
            
            ExponentialDistribution() : ExponentialDistribution(RealType(1)) { }
            
            explicit ExponentialDistribution(RealType lambda) : m_param(lambda) { }
            
            explicit ExponentialDistribution(const param_type& param) : m_param(param) { }
 
            void reset() { }
 
            // Parameter access
            
            RealType lambda() const { return m_param.lambda(); }
            
            param_type param() const { return m_param; }
            
            void param(const param_type& param) { m_param = param; }
            
            result_type min() const { return RealType(0); }
            
            result_type max() const { return std::numeric_limits<result_type>::max(); }
 
            // Generation
            
            template<typename Generator>
            result_type operator()(Generator& g) { return (*this)(g, m_param); }
 
            template<typename Generator>
            result_type operator()(Generator& g, const param_type& param)
            {
                UniformRealDistribution<RealType> uniform(RealType(0), RealType(1));
                RealType u;
                do {
                    u = uniform(g);
                } while (u == RealType(0) || u == RealType(1));
                
                return -std::log(RealType(1) - u) / param.lambda();
            }
 
            // Comparison
            
            friend bool operator==(const ExponentialDistribution& lhs, const ExponentialDistribution& rhs)
            {
                return lhs.m_param == rhs.m_param;
            }
            friend bool operator!=(const ExponentialDistribution& lhs, const ExponentialDistribution& rhs)
            {
                return !(lhs == rhs);
            }
 
        private:
            param_type m_param;
    };
 
    template<typename RealType = double>
    class INFINITY_API_TEMPLATE WeibullDistribution
    {
        static_assert(std::is_floating_point<RealType>::value,
                    "RealType must be a floating point type");
 
        public:
            using result_type = RealType; 
 
            struct param_type
            {
                using distribution_type = WeibullDistribution<RealType>;
 
                param_type() : param_type(RealType(1)) {}
                
                explicit param_type(RealType a, RealType b = RealType(1)) : m_a(a), m_b(b) { }
 
                RealType a() const { return m_a; }
                
                RealType b() const { return m_b; }
 
                friend bool operator==(const param_type& lhs, const param_type& rhs)
                {
                    return lhs.m_a == rhs.m_a && lhs.m_b == rhs.m_b;
                }
                friend bool operator!=(const param_type& lhs, const param_type& rhs)
                {
                    return !(lhs == rhs);
                }
 
            private:
                RealType m_a, m_b;
            };
 
            // Constructors
            
            WeibullDistribution() : WeibullDistribution(RealType(1)) { }
            
            explicit WeibullDistribution(RealType a, RealType b = RealType(1)) : m_param(a, b) { }
            
            explicit WeibullDistribution(const param_type& param) : m_param(param) { }
 
            void reset() { }
 
            // Parameter access
            
            RealType a() const { return m_param.a(); }
            
            RealType b() const { return m_param.b(); }
            
            param_type param() const { return m_param; }
            
            void param(const param_type& param) { m_param = param; }
            
            result_type min() const { return RealType(0); }
            
            result_type max() const { return std::numeric_limits<result_type>::max(); }
 
            // Generation
            
            template<typename Generator>
            result_type operator()(Generator& g) { return (*this)(g, m_param); }
 
            template<typename Generator>
            result_type operator()(Generator& g, const param_type& param)
            {
                ExponentialDistribution<RealType> exp(RealType(1));
                return param.b() * std::pow(exp(g), RealType(1) / param.a());
            }
 
            // Comparison
            
            friend bool operator==(const WeibullDistribution& lhs, const WeibullDistribution& rhs)
            {
                return lhs.m_param == rhs.m_param;
            }
            friend bool operator!=(const WeibullDistribution& lhs, const WeibullDistribution& rhs)
            {
                return !(lhs == rhs);
            }
 
        private:
            param_type m_param;
    };
 
    template<typename RealType = double>
    class INFINITY_API_TEMPLATE ExtremeValueDistribution
    {
        static_assert(std::is_floating_point<RealType>::value,
                    "RealType must be a floating point type");
 
        public:
            using result_type = RealType; 
 
            struct param_type
            {
                using distribution_type = ExtremeValueDistribution<RealType>;
 
                param_type() : param_type(RealType(0)) { }
                
                explicit param_type(RealType a, RealType b = RealType(1)) : m_a(a), m_b(b) { }
 
                RealType a() const { return m_a; }
                
                RealType b() const { return m_b; }
 
                friend bool operator==(const param_type& lhs, const param_type& rhs)
                {
                    return lhs.m_a == rhs.m_a && lhs.m_b == rhs.m_b;
                }
                friend bool operator!=(const param_type& lhs, const param_type& rhs)
                {
                    return !(lhs == rhs);
                }
 
            private:
                RealType m_a, m_b;
            };
 
            // Constructors
            
            ExtremeValueDistribution() : ExtremeValueDistribution(RealType(0)) { }
            
            explicit ExtremeValueDistribution(RealType a, RealType b = RealType(1)) : m_param(a, b) { }
            
            explicit ExtremeValueDistribution(const param_type& param) : m_param(param) { }
 
            void reset() { }
 
            // Parameter access
            
            RealType a() const { return m_param.a(); }
            
            RealType b() const { return m_param.b(); }
            
            param_type param() const { return m_param; }
            
            void param(const param_type& param) { m_param = param; }
            
            result_type min() const { return std::numeric_limits<result_type>::lowest(); }
            
            result_type max() const { return std::numeric_limits<result_type>::max(); }
 
            // Generation
            
            template<typename Generator>
            result_type operator()(Generator& g) { return (*this)(g, m_param); }
 
            template<typename Generator>
            result_type operator()(Generator& g, const param_type& param)
            {
                UniformRealDistribution<RealType> uniform(RealType(0), RealType(1));
                RealType u;
                do {
                    u = uniform(g);
                } while (u == RealType(0) || u == RealType(1));
                
                return param.a() - param.b() * std::log(-std::log(u));
            }
 
            // Comparison
            
            friend bool operator==(const ExtremeValueDistribution& lhs, const ExtremeValueDistribution& rhs)
            {
                return lhs.m_param == rhs.m_param;
            }
            friend bool operator!=(const ExtremeValueDistribution& lhs, const ExtremeValueDistribution& rhs)
            {
                return !(lhs == rhs);
            }
 
        private:
            param_type m_param;
    };
 
    class INFINITY_API_TEMPLATE BernoulliDistribution
    {
        public:
            using result_type = bool; 
 
            struct param_type {
                using distribution_type = BernoulliDistribution;
 
                param_type() : param_type(0.5) { }
                
                explicit param_type(double p) : m_p(p) { }
 
                double p() const { return m_p; }
 
                friend bool operator==(const param_type& lhs, const param_type& rhs)
                {
                    return lhs.m_p == rhs.m_p;
                }
                friend bool operator!=(const param_type& lhs, const param_type& rhs)
                {
                    return !(lhs == rhs);
                }
 
            private:
                double m_p;
            };
 
            // Constructors
            
            BernoulliDistribution() : BernoulliDistribution(0.5) { }
            
            explicit BernoulliDistribution(double p) : m_param(p) { }
            
            explicit BernoulliDistribution(const param_type& param) : m_param(param) { }
 
            void reset() { }
 
            // Parameter access
            
            double p() const { return m_param.p(); }
            
            param_type param() const { return m_param; }
            
            void param(const param_type& param) { m_param = param; }
            
            result_type min() const { return false; }
            
            result_type max() const { return true; }
 
            // Generation
            
            template<typename Generator>
            result_type operator()(Generator& g) { return (*this)(g, m_param); }
 
            template<typename Generator>
            result_type operator()(Generator& g, const param_type& param)
            {
                UniformRealDistribution<double> uniform(0.0, 1.0);
                return uniform(g) < param.p();
            }
 
            // Comparison
            
            friend bool operator==(const BernoulliDistribution& lhs, const BernoulliDistribution& rhs)
            {
                return lhs.m_param == rhs.m_param;
            }
            friend bool operator!=(const BernoulliDistribution& lhs, const BernoulliDistribution& rhs)
            {
                return !(lhs == rhs);
            }
 
        private:
            param_type m_param;
    };
 
    template<typename IntType = int>
    class INFINITY_API_TEMPLATE BinomialDistribution
    {
        static_assert(std::is_integral<IntType>::value,
                    "IntType must be an integral type");
 
        public:
            using result_type = IntType; 
 
            struct param_type {
                using distribution_type = BinomialDistribution<IntType>;
 
                param_type() : param_type(1) { }
                
                explicit param_type(IntType t, double p = 0.5) : m_t(t), m_p(p) { }
 
                IntType t() const { return m_t; }
                
                double p() const { return m_p; }
 
                friend bool operator==(const param_type& lhs, const param_type& rhs)
                {
                    return lhs.m_t == rhs.m_t && lhs.m_p == rhs.m_p;
                }
                friend bool operator!=(const param_type& lhs, const param_type& rhs)
                {
                    return !(lhs == rhs);
                }
 
            private:
                IntType m_t;
                double m_p;
            };
 
            // Constructors
            
            BinomialDistribution() : BinomialDistribution(1) { }
            
            explicit BinomialDistribution(IntType t, double p = 0.5)
                : m_param(t, p), m_nd(nullptr)
            {
                if (t * p >= 5.0 && t * (1.0 - p) >= 5.0)
                {
                    double mean = t * p;
                    double stddev = std::sqrt(t * p * (1.0 - p));
                    m_nd = new NormalDistribution<double>(mean, stddev);
                }
            }
            
            explicit BinomialDistribution(const param_type& param)
                : m_param(param), m_nd(nullptr)
            {
                IntType t = param.t();
                double p = param.p();
                if (t * p >= 5.0 && t * (1.0 - p) >= 5.0)
                {
                    double mean = t * p;
                    double stddev = std::sqrt(t * p * (1.0 - p));
                    m_nd = new NormalDistribution<double>(mean, stddev);
                }
            }
 
            void reset()
            {
                if (m_nd) m_nd->reset();
            }
 
            // Parameter access
            
            IntType t() const { return m_param.t(); }
            
            double p() const { return m_param.p(); }
            
            param_type param() const { return m_param; }
            
            void param(const param_type& param)
            {
                m_param = param;
                delete m_nd;
                m_nd = nullptr;
                IntType t = param.t();
                double p = param.p();
                if (t * p >= 5.0 && t * (1.0 - p) >= 5.0)
                {
                    double mean = t * p;
                    double stddev = std::sqrt(t * p * (1.0 - p));
                    m_nd = new NormalDistribution<double>(mean, stddev);
                }
            }
            
            result_type min() const { return 0; }
            
            result_type max() const { return m_param.t(); }
 
            // Generation
            
            template<typename Generator>
            result_type operator()(Generator& g)
            {
                return (*this)(g, m_param);
            }
 
            template<typename Generator>
            result_type operator()(Generator& g, const param_type& param)
            {
                const IntType t = param.t();
                const double p = param.p();
                
                if (t * p >= 5.0 && t * (1.0 - p) >= 5.0)
                {
                    // Normal approximation with continuity correction
                    double mean = t * p;
                    double stddev = std::sqrt(t * p * (1.0 - p));
                    NormalDistribution<double> nd(mean, stddev);
                    IntType result;
                    do {
                        result = static_cast<IntType>(std::floor(nd(g) + 0.5));
                    } while (result < 0 || result > t);
                    return result;
                } else {
                    // Direct method: simulate each trial
                    BernoulliDistribution bern(p);
                    IntType sum = 0;
                    for (IntType i = 0; i < t; ++i)
                    {
                        if (bern(g)) ++sum;
                    }
                    return sum;
                }
            }
 
            // Comparison
            
            friend bool operator==(const BinomialDistribution& lhs, const BinomialDistribution& rhs)
            {
                if (lhs.m_param != rhs.m_param) return false;
                if ((lhs.m_nd == nullptr) != (rhs.m_nd == nullptr)) return false;
                if (lhs.m_nd && rhs.m_nd && *lhs.m_nd != *rhs.m_nd) return false;
                return true;
            }
            friend bool operator!=(const BinomialDistribution& lhs, const BinomialDistribution& rhs)
            {
                return !(lhs == rhs);
            }
 
        private:
            param_type m_param;
            class NormalDistribution<double>* m_nd;  
    };
 
    template<typename IntType = int>
    class INFINITY_API_TEMPLATE GeometricDistribution
    {
        static_assert(std::is_integral<IntType>::value,
                    "IntType must be an integral type");
 
        public:
            using result_type = IntType; 
 
            struct param_type {
                using distribution_type = GeometricDistribution<IntType>;
 
                param_type() : param_type(0.5) { }
                
                explicit param_type(double p) : m_p(p) { }
 
                double p() const { return m_p; }
 
                friend bool operator==(const param_type& lhs, const param_type& rhs)
                {
                    return lhs.m_p == rhs.m_p;
                }
                friend bool operator!=(const param_type& lhs, const param_type& rhs)
                {
                    return !(lhs == rhs);
                }
 
            private:
                double m_p;
            };
 
            // Constructors
            
            GeometricDistribution() : GeometricDistribution(0.5) { }
            
            explicit GeometricDistribution(double p) : m_param(p) { }
            
            explicit GeometricDistribution(const param_type& param) : m_param(param) { }
 
            void reset() { }
 
            // Parameter access
            
            double p() const { return m_param.p(); }
            
            param_type param() const { return m_param; }
            
            void param(const param_type& param) { m_param = param; }
            
            result_type min() const { return 0; }
            
            result_type max() const { return std::numeric_limits<result_type>::max(); }
 
            // Generation
            
            template<typename Generator>
            result_type operator()(Generator& g) { return (*this)(g, m_param); }
 
            template<typename Generator>
            result_type operator()(Generator& g, const param_type& param)
            {
                UniformRealDistribution<double> uniform(0.0, 1.0);
                double u;
                do {
                    u = uniform(g);
                } while (u == 0.0 || u == 1.0);
                
                const double log_1_p = std::log(1.0 - param.p());
                return static_cast<IntType>(std::ceil(std::log(u) / log_1_p) - 1.0);
            }
 
            // Comparison
            
            friend bool operator==(const GeometricDistribution& lhs, const GeometricDistribution& rhs)
            {
                return lhs.m_param == rhs.m_param;
            }
            friend bool operator!=(const GeometricDistribution& lhs, const GeometricDistribution& rhs)
            {
                return !(lhs == rhs);
            }
 
        private:
            param_type m_param;
    };
 
    template<typename IntType = int>
    class INFINITY_API_TEMPLATE PoissonDistribution
    {
        static_assert(std::is_integral<IntType>::value,
                    "IntType must be an integral type");
 
        public:
            using result_type = IntType; 
 
            struct param_type {
                using distribution_type = PoissonDistribution<IntType>;
 
                param_type() : param_type(1.0) { }
                
                explicit param_type(double mean) : m_mean(mean) { }
 
                double mean() const { return m_mean; }
 
                friend bool operator==(const param_type& lhs, const param_type& rhs)
                {
                    return lhs.m_mean == rhs.m_mean;
                }
                friend bool operator!=(const param_type& lhs, const param_type& rhs)
                {
                    return !(lhs == rhs);
                }
 
            private:
                double m_mean;
            };
 
            // Constructors
            
            PoissonDistribution() : PoissonDistribution(1.0) { }
            
            explicit PoissonDistribution(double mean)
                : m_param(mean), m_nd(nullptr)
            {
                if (mean >= 10.0)
                {
                    m_nd = new NormalDistribution<double>(mean, std::sqrt(mean));
                }
            }
            
            explicit PoissonDistribution(const param_type& param)
                : m_param(param), m_nd(nullptr)
            {
                if (param.mean() >= 10.0)
                {
                    m_nd = new NormalDistribution<double>(param.mean(), std::sqrt(param.mean()));
                }
            }
 
            void reset()
            {
                if (m_nd) m_nd->reset();
            }
 
            // Parameter access
            
            double mean() const { return m_param.mean(); }
            
            param_type param() const { return m_param; }
            
            void param(const param_type& param)
            {
                m_param = param;
                delete m_nd;
                m_nd = nullptr;
                if (param.mean() >= 10.0)
                {
                    m_nd = new NormalDistribution<double>(param.mean(), std::sqrt(param.mean()));
                }
            }
            
            result_type min() const { return 0; }
            
            result_type max() const { return std::numeric_limits<result_type>::max(); }
 
            // Generation
            
            template<typename Generator>
            result_type operator()(Generator& g)
            {
                return (*this)(g, m_param);
            }
 
            template<typename Generator>
            result_type operator()(Generator& g, const param_type& param)
            {
                const double mean = param.mean();
                
                if (mean >= 10.0)
                {
                    // Normal approximation for large mean
                    NormalDistribution<double> nd(mean, std::sqrt(mean));
                    IntType result;
                    do {
                        result = static_cast<IntType>(std::floor(nd(g) + 0.5));
                    } while (result < 0);
                    return result;
                } else {
                    // Knuth method for small mean
                    UniformRealDistribution<double> uniform(0.0, 1.0);
                    const double l = std::exp(-mean);
                    IntType k = 0;
                    double p = 1.0;
                    
                    do {
                        ++k;
                        p *= uniform(g);
                    } while (p > l);
                    
                    return k - 1;
                }
            }
 
            // Comparison
            
            friend bool operator==(const PoissonDistribution& lhs, const PoissonDistribution& rhs)
            {
                if (lhs.m_param != rhs.m_param) return false;
                if ((lhs.m_nd == nullptr) != (rhs.m_nd == nullptr)) return false;
                if (lhs.m_nd && rhs.m_nd && *lhs.m_nd != *rhs.m_nd) return false;
                return true;
            }
            friend bool operator!=(const PoissonDistribution& lhs, const PoissonDistribution& rhs)
            {
                return !(lhs == rhs);
            }
 
        private:
            param_type m_param;
            class NormalDistribution<double>* m_nd;  
    };
 
    template<typename IntType = int>
    class INFINITY_API_TEMPLATE NegativeBinomialDistribution
    {
        static_assert(std::is_integral<IntType>::value,
                    "IntType must be an integral type");
 
        public:
            using result_type = IntType; 
 
            struct param_type {
                using distribution_type = NegativeBinomialDistribution<IntType>;
 
                param_type() : param_type(1) { }
                
                explicit param_type(IntType k, double p = 0.5) : m_k(k), m_p(p) {}
 
                IntType k() const { return m_k; }
                
                double p() const { return m_p; }
 
                friend bool operator==(const param_type& lhs, const param_type& rhs)
                {
                    return lhs.m_k == rhs.m_k && lhs.m_p == rhs.m_p;
                }
                friend bool operator!=(const param_type& lhs, const param_type& rhs)
                {
                    return !(lhs == rhs);
                }
 
            private:
                IntType m_k;
                double m_p;
            };
 
            // Constructors
            
            NegativeBinomialDistribution() : NegativeBinomialDistribution(1) { }
            
            explicit NegativeBinomialDistribution(IntType k, double p = 0.5)
                : m_param(k, p)
            {
                m_gd = new GammaDistribution<double>(k, (1.0 - p) / p);
            }
            
            explicit NegativeBinomialDistribution(const param_type& param)
                : m_param(param)
            {
                m_gd = new GammaDistribution<double>(param.k(), (1.0 - param.p()) / param.p());
            }
 
            void reset()
            {
                if (m_gd) m_gd->reset();
            }
 
            // Parameter access
            
            IntType k() const { return m_param.k(); }
            
            double p() const { return m_param.p(); }
            
            param_type param() const { return m_param; }
            
            void param(const param_type& param)
            {
                m_param = param;
                delete m_gd;
                m_gd = new GammaDistribution<double>(param.k(), (1.0 - param.p()) / param.p());
            }
            
            result_type min() const { return 0; }
            
            result_type max() const { return std::numeric_limits<result_type>::max(); }
 
            // Generation
            
            template<typename Generator>
            result_type operator()(Generator& g)
            {
                return (*this)(g, m_param);
            }
 
            template<typename Generator>
            result_type operator()(Generator& g, const param_type& param)
            {
                GammaDistribution<double> gd(param.k(), (1.0 - param.p()) / param.p());
                double gamma_variate = gd(g);
                PoissonDistribution<IntType> pd(gamma_variate);
                return pd(g);
            }
 
            // Comparison
            
            friend bool operator==(const NegativeBinomialDistribution& lhs, const NegativeBinomialDistribution& rhs)
            {
                if (lhs.m_param != rhs.m_param) return false;
                if ((lhs.m_gd == nullptr) != (rhs.m_gd == nullptr)) return false;
                if (lhs.m_gd && rhs.m_gd && *lhs.m_gd != *rhs.m_gd) return false;
                return true;
            }
            friend bool operator!=(const NegativeBinomialDistribution& lhs, const NegativeBinomialDistribution& rhs)
            {
                return !(lhs == rhs);
            }
 
        private:
            param_type m_param;
            class GammaDistribution<double>* m_gd;  
    };
 
    template<typename IntType = int>
    class INFINITY_API_TEMPLATE DiscreteDistribution
    {
        static_assert(std::is_integral<IntType>::value,
                    "IntType must be an integral type");
 
        public:
            using result_type = IntType; 
 
            struct param_type
            {
                using distribution_type = DiscreteDistribution<IntType>;
 
                param_type() : m_prob(1, 1.0), m_cp(1, 1.0) {}
 
                template<typename InputIterator>
                param_type(InputIterator first, InputIterator last)
                    : m_prob(first, last), m_cp()
                {
                    normalize();
                }
 
                param_type(std::initializer_list<double> wl)
                    : m_prob(wl.begin(), wl.end()), m_cp()
                {
                    normalize();
                }
 
                template<typename UnaryOperation>
                param_type(size_t count, double xmin, double xmax, UnaryOperation fw)
                    : m_prob(count), m_cp()
                {
                    const double delta = (xmax - xmin) / count;
                    for (size_t i = 0; i < count; ++i)
                    {
                        m_prob[i] = fw(xmin + (i + 0.5) * delta);
                    }
                    normalize();
                }
 
                std::vector<double> probabilities() const
                {
                    return m_prob.empty() ? std::vector<double>(1, 1.0) : m_prob;
                }
 
                friend bool operator==(const param_type& lhs, const param_type& rhs)
                {
                    return lhs.m_prob == rhs.m_prob;
                }
                friend bool operator!=(const param_type& lhs, const param_type& rhs)
                {
                    return !(lhs == rhs);
                }
 
                std::vector<double> m_prob;  
                std::vector<double> m_cp;    
 
            private:
                void normalize()
                {
                    if (m_prob.empty())
                    {
                        m_prob.push_back(1.0);
                        m_cp.push_back(1.0);
                        return;
                    }
 
                    // Normalize probabilities
                    double sum = 0.0;
                    for (double p : m_prob)
                    {
                        sum += p;
                    }
 
                    if (sum == 0.0)
                    {
                        // All weights are zero, use uniform distribution
                        const double uniform = 1.0 / m_prob.size();
                        for (auto& p : m_prob)
                        {
                            p = uniform;
                        }
                    } else {
                        for (auto& p : m_prob)
                        {
                            p /= sum;
                        }
                    }
 
                    // Build cumulative probability table
                    m_cp.resize(m_prob.size());
                    m_cp[0] = m_prob[0];
                    for (size_t i = 1; i < m_prob.size(); ++i)
                    {
                        m_cp[i] = m_cp[i - 1] + m_prob[i];
                    }
                    m_cp.back() = 1.0; // Ensure last element is exactly 1.0
                }
                friend class DiscreteDistribution<IntType>;
            };
 
            // Constructors
            
            DiscreteDistribution() : m_param() { }
 
            template<typename InputIterator>
            DiscreteDistribution(InputIterator first, InputIterator last)
                : m_param(first, last) { }
 
            DiscreteDistribution(std::initializer_list<double> wl)
                : m_param(wl) { }
 
            template<typename UnaryOperation>
            DiscreteDistribution(size_t count, double xmin, double xmax, UnaryOperation fw)
                : m_param(count, xmin, xmax, fw) { }
 
            explicit DiscreteDistribution(const param_type& param) : m_param(param) { }
 
            void reset() { }
 
            // Parameter access
            
            std::vector<double> probabilities() const { return m_param.probabilities(); }
            
            param_type param() const { return m_param; }
            
            void param(const param_type& param) { m_param = param; }
            
            result_type min() const { return 0; }
            
            result_type max() const {
                return m_param.m_prob.empty() ? 0 : 
                    static_cast<result_type>(m_param.m_prob.size() - 1);
            }
 
            // Generation
            
            template<typename Generator>
            result_type operator()(Generator& g) { return (*this)(g, m_param); }
 
            template<typename Generator>
            result_type operator()(Generator& g, const param_type& param)
            {
                UniformRealDistribution<double> uniform(0.0, 1.0);
                const double u = uniform(g);
 
                // Binary search in cumulative probability table
                auto it = std::lower_bound(param.m_cp.begin(), param.m_cp.end(), u);
                return static_cast<IntType>(it - param.m_cp.begin());
            }
 
            // Comparison
            
            friend bool operator==(const DiscreteDistribution& lhs, const DiscreteDistribution& rhs)
            {
                return lhs.m_param == rhs.m_param;
            }
            friend bool operator!=(const DiscreteDistribution& lhs, const DiscreteDistribution& rhs)
            {
                return !(lhs == rhs);
            }
 
        private:
            param_type m_param;
    };
 
    template<typename RealType = double>
    class INFINITY_API_TEMPLATE PiecewiseConstantDistribution
    {
        static_assert(std::is_floating_point<RealType>::value,
                    "RealType must be a floating point type");
 
        public:
            using result_type = RealType; 
 
            struct param_type
            {
                using distribution_type = PiecewiseConstantDistribution<RealType>;
 
                param_type()
                    : m_int{RealType(0), RealType(1)}, m_den{1.0}, m_cp{1.0}
                { }
 
                template<typename InputIteratorB, typename InputIteratorW>
                param_type(InputIteratorB first_b, InputIteratorB last_b,
                        InputIteratorW first_w)
                    : m_int(first_b, last_b), m_den(), m_cp()
                {
                    
                    if (m_int.size() < 2)
                    {
                        m_int = {RealType(0), RealType(1)};
                        m_den = {1.0};
                        m_cp = {1.0};
                        return;
                    }
 
                    m_den.resize(m_int.size() - 1);
                    for (size_t i = 0; i < m_den.size(); ++i, ++first_w)
                    {
                        m_den[i] = *first_w;
                    }
                    
                    initialize();
                }
 
                template<typename UnaryOperation>
                param_type(std::initializer_list<RealType> bl, UnaryOperation fw)
                    : m_int(bl.begin(), bl.end()), m_den(), m_cp()
                {
                    
                    if (m_int.size() < 2)
                    {
                        m_int = {RealType(0), RealType(1)};
                        m_den = {1.0};
                        m_cp = {1.0};
                        return;
                    }
 
                    m_den.resize(m_int.size() - 1);
                    for (size_t i = 0; i < m_den.size(); ++i)
                    {
                        const RealType x = (m_int[i] + m_int[i + 1]) / RealType(2);
                        m_den[i] = fw(x);
                    }
                    
                    initialize();
                }
 
                template<typename UnaryOperation>
                param_type(size_t nw, RealType xmin, RealType xmax, UnaryOperation fw)
                    : m_int(nw + 1), m_den(nw), m_cp()
                {
                    
                    const RealType delta = (xmax - xmin) / nw;
                    for (size_t i = 0; i <= nw; ++i)
                    {
                        m_int[i] = xmin + i * delta;
                    }
 
                    for (size_t i = 0; i < nw; ++i)
                    {
                        const RealType x = (m_int[i] + m_int[i + 1]) / RealType(2);
                        m_den[i] = fw(x);
                    }
                    
                    initialize();
                }
 
                std::vector<RealType> intervals() const
                {
                    if (m_int.empty()) {
                        std::vector<RealType> result(2);
                        result[0] = RealType(0);
                        result[1] = RealType(1);
                        return result;
                    }
                    return m_int;
                }
                
                std::vector<double> densities() const
                {
                    return m_den.empty() ? std::vector<double>(1, 1.0) : m_den;
                }
 
                friend bool operator==(const param_type& lhs, const param_type& rhs)
                {
                    return lhs.m_int == rhs.m_int && lhs.m_den == rhs.m_den;
                }
                friend bool operator!=(const param_type& lhs, const param_type& rhs)
                {
                    return !(lhs == rhs);
                }
 
                std::vector<RealType> m_int;  
                std::vector<double> m_den;    
                std::vector<double> m_cp;     
 
            private:
                void initialize()
                {
                    // Calculate area under each piece
                    std::vector<double> areas(m_den.size());
                    double total_area = 0.0;
 
                    for (size_t i = 0; i < m_den.size(); ++i)
                    {
                        const double width = static_cast<double>(m_int[i + 1] - m_int[i]);
                        areas[i] = m_den[i] * width;
                        total_area += areas[i];
                    }
 
                    // Normalize densities
                    if (total_area == 0.0)
                    {
                        for (auto& d : m_den)
                        {
                            d = 1.0;
                        }
                        total_area = static_cast<double>(m_int.back() - m_int.front());
                    }
 
                    for (auto& d : m_den)
                    {
                        d /= total_area;
                    }
 
                    // Build cumulative probability table
                    m_cp.resize(m_den.size());
                    m_cp[0] = m_den[0] * static_cast<double>(m_int[1] - m_int[0]);
                    for (size_t i = 1; i < m_den.size(); ++i)
                    {
                        const double width = static_cast<double>(m_int[i + 1] - m_int[i]);
                        m_cp[i] = m_cp[i - 1] + m_den[i] * width;
                    }
                    m_cp.back() = 1.0;
                }
                friend class PiecewiseConstantDistribution<RealType>;
            };
 
            // Constructors
            
            PiecewiseConstantDistribution() : m_param() { }
 
            template<typename InputIteratorB, typename InputIteratorW>
            PiecewiseConstantDistribution(InputIteratorB first_b, InputIteratorB last_b,
                                        InputIteratorW first_w)
                : m_param(first_b, last_b, first_w) { }
 
            template<typename UnaryOperation>
            PiecewiseConstantDistribution(std::initializer_list<RealType> bl, UnaryOperation fw)
                : m_param(bl, fw) { }
 
            template<typename UnaryOperation>
            PiecewiseConstantDistribution(size_t nw, RealType xmin, RealType xmax, UnaryOperation fw)
                : m_param(nw, xmin, xmax, fw) { }
 
            explicit PiecewiseConstantDistribution(const param_type& param) : m_param(param) { }
 
            void reset() { }
 
            // Parameter access
            
            std::vector<RealType> intervals() const { return m_param.intervals(); }
            
            std::vector<double> densities() const { return m_param.densities(); }
            
            param_type param() const { return m_param; }
            
            void param(const param_type& param) { m_param = param; }
            
            result_type min() const
            {
                return m_param.m_int.empty() ? RealType(0) : m_param.m_int.front();
            }
            
            result_type max() const
            {
                return m_param.m_int.empty() ? RealType(1) : m_param.m_int.back();
            }
 
            // Generation
            
            template<typename Generator>
            result_type operator()(Generator& g) { return (*this)(g, m_param); }
 
            template<typename Generator>
            result_type operator()(Generator& g, const param_type& param)
            {
                UniformRealDistribution<double> uniform(0.0, 1.0);
                const double u = uniform(g);
 
                // Find which interval using binary search
                auto it = std::lower_bound(param.m_cp.begin(), param.m_cp.end(), u);
                const size_t idx = it - param.m_cp.begin();
 
                // Generate uniform within that interval
                UniformRealDistribution<RealType> interval_uniform(param.m_int[idx], param.m_int[idx + 1]);
                return interval_uniform(g);
            }
 
            // Comparison
            
            friend bool operator==(const PiecewiseConstantDistribution& lhs,
                                const PiecewiseConstantDistribution& rhs)
            {
                return lhs.m_param == rhs.m_param;
            }
            friend bool operator!=(const PiecewiseConstantDistribution& lhs,
                                const PiecewiseConstantDistribution& rhs)
            {
                return !(lhs == rhs);
            }
 
        private:
            param_type m_param;
    };
 
    template<typename RealType = double>
    class INFINITY_API_TEMPLATE PiecewiseLinearDistribution
    {
        static_assert(std::is_floating_point<RealType>::value,
                    "RealType must be a floating point type");
 
        public:
            using result_type = RealType; 
 
            struct param_type
            {
                using distribution_type = PiecewiseLinearDistribution<RealType>;
 
                param_type()
                    : m_int{RealType(0), RealType(1)}, m_den{1.0, 1.0}, m_cp{1.0}, m_m{0.0}
                { }
 
                template<typename InputIteratorB, typename InputIteratorW>
                param_type(InputIteratorB first_b, InputIteratorB last_b,
                        InputIteratorW first_w)
                    : m_int(first_b, last_b), m_den(), m_cp(), m_m() 
                {
                    
                    if (m_int.size() < 2)
                    {
                        m_int = {RealType(0), RealType(1)};
                        m_den = {1.0, 1.0};
                        m_cp = {1.0};
                        m_m = {0.0};
                        return;
                    }
 
                    m_den.resize(m_int.size());
                    for (size_t i = 0; i < m_den.size(); ++i, ++first_w)
                    {
                        m_den[i] = *first_w;
                    }
                    
                    initialize();
                }
 
                template<typename UnaryOperation>
                param_type(std::initializer_list<RealType> bl, UnaryOperation fw)
                    : m_int(bl.begin(), bl.end()), m_den(), m_cp(), m_m()
                {
                    
                    if (m_int.size() < 2)
                    {
                        m_int = {RealType(0), RealType(1)};
                        m_den = {1.0, 1.0};
                        m_cp = {1.0};
                        m_m = {0.0};
                        return;
                    }
 
                    m_den.resize(m_int.size());
                    for (size_t i = 0; i < m_den.size(); ++i)
                    {
                        m_den[i] = fw(m_int[i]);
                    }
                    
                    initialize();
                }
 
                template<typename UnaryOperation>
                param_type(size_t nw, RealType xmin, RealType xmax, UnaryOperation fw)
                    : m_int(nw + 1), m_den(nw + 1), m_cp(), m_m()
                {
                    
                    const RealType delta = (xmax - xmin) / nw;
                    for (size_t i = 0; i <= nw; ++i)
                    {
                        m_int[i] = xmin + i * delta;
                        m_den[i] = fw(m_int[i]);
                    }
                    
                    initialize();
                }
 
                std::vector<RealType> intervals() const
                {
                    if (m_int.empty())
                    {
                        std::vector<RealType> result(2);
                        result[0] = RealType(0);
                        result[1] = RealType(1);
                        return result;
                    }
                    return m_int;
                }
 
                std::vector<double> densities() const
                {
                    return m_den.empty() ? std::vector<double>(2, 1.0) : m_den;
                }
 
                friend bool operator==(const param_type& lhs, const param_type& rhs)
                {
                    return lhs.m_int == rhs.m_int && lhs.m_den == rhs.m_den;
                }
                friend bool operator!=(const param_type& lhs, const param_type& rhs)
                {
                    return !(lhs == rhs);
                }
 
                std::vector<RealType> m_int;  
                std::vector<double> m_den;    
                std::vector<double> m_cp;     
                std::vector<double> m_m;      
 
            private:
                void initialize()
                {
                    const size_t n = m_int.size() - 1;
                    
                    // Calculate area under each trapezoid
                    std::vector<double> areas(n);
                    double total_area = 0.0;
 
                    for (size_t i = 0; i < n; ++i)
                    {
                        const double width = static_cast<double>(m_int[i + 1] - m_int[i]);
                        areas[i] = 0.5 * width * (m_den[i] + m_den[i + 1]);
                        total_area += areas[i];
                    }
 
                    // Normalize densities
                    if (total_area == 0.0)
                    {
                        for (auto& d : m_den)
                        {
                            d = 1.0;
                        }
                        for (size_t i = 0; i < n; ++i)
                        {
                            const double width = static_cast<double>(m_int[i + 1] - m_int[i]);
                            areas[i] = width;
                            total_area += areas[i];
                        }
                    }
 
                    for (auto& d : m_den)
                    {
                        d /= total_area;
                    }
                    for (auto& a : areas)
                    {
                        a /= total_area;
                    }
 
                    // Build cumulative probability table and slopes
                    m_cp.resize(n);
                    m_m.resize(n);
                    
                    m_cp[0] = areas[0];
                    for (size_t i = 1; i < n; ++i)
                    {
                        m_cp[i] = m_cp[i - 1] + areas[i];
                    }
                    m_cp.back() = 1.0;
 
                    // Calculate slopes for each piece
                    for (size_t i = 0; i < n; ++i)
                    {
                        const double width = static_cast<double>(m_int[i + 1] - m_int[i]);
                        m_m[i] = (m_den[i + 1] - m_den[i]) / width;
                    }
                }
                
                friend class PiecewiseLinearDistribution<RealType>;
            };
 
            // Constructors
            
            PiecewiseLinearDistribution() : m_param() { }
 
            template<typename InputIteratorB, typename InputIteratorW>
            PiecewiseLinearDistribution(InputIteratorB first_b, InputIteratorB last_b,
                                        InputIteratorW first_w)
                : m_param(first_b, last_b, first_w) { }
 
            template<typename UnaryOperation>
            PiecewiseLinearDistribution(std::initializer_list<RealType> bl, UnaryOperation fw)
                : m_param(bl, fw) { }
 
            template<typename UnaryOperation>
            PiecewiseLinearDistribution(size_t nw, RealType xmin, RealType xmax, UnaryOperation fw)
                : m_param(nw, xmin, xmax, fw) { }
 
            explicit PiecewiseLinearDistribution(const param_type& param) : m_param(param) { }
 
            void reset() { }
 
            // Parameter access
            
            std::vector<RealType> intervals() const { return m_param.intervals(); }
            
            std::vector<double> densities() const { return m_param.densities(); }
            
            param_type param() const { return m_param; }
            
            void param(const param_type& param) { m_param = param; }
            
            result_type min() const
            {
                return m_param.m_int.empty() ? RealType(0) : m_param.m_int.front();
            }
            
            result_type max() const
            {
                return m_param.m_int.empty() ? RealType(1) : m_param.m_int.back();
            }
 
            // Generation
            
            template<typename Generator>
            result_type operator()(Generator& g) { return (*this)(g, m_param); }
 
            template<typename Generator>
            result_type operator()(Generator& g, const param_type& param)
            {
                UniformRealDistribution<double> uniform(0.0, 1.0);
                const double u = uniform(g);
 
                // Find which interval using binary search
                auto it = std::lower_bound(param.m_cp.begin(), param.m_cp.end(), u);
                const size_t idx = it - param.m_cp.begin();
 
                // Calculate cumulative probability at start of this interval
                const double cp_start = (idx == 0) ? 0.0 : param.m_cp[idx - 1];
                const double u_local = u - cp_start;
 
                // Solve for x in the linear piece using quadratic formula
                // For linear density d(x) = d0 + m*(x-x0), the CDF is:
                // F(x) = d0*(x-x0) + 0.5*m*(x-x0)^2
                // Solving F(x) = u_local for x
                
                const double x0 = static_cast<double>(param.m_int[idx]);
                const double d0 = param.m_den[idx];
                const double m = param.m_m[idx];
 
                if (std::abs(m) < 1e-10)
                {
                    // Nearly constant density in this interval
                    return static_cast<RealType>(x0 + u_local / d0);
                } else {
                    // Quadratic formula: 0.5*m*dx^2 + d0*dx - u_local = 0
                    const double discriminant = d0 * d0 + 2.0 * m * u_local;
                    const double dx = (-d0 + std::sqrt(discriminant)) / m;
                    return static_cast<RealType>(x0 + dx);
                }
            }
 
            // Comparison
            
            friend bool operator==(const PiecewiseLinearDistribution& lhs,
                                const PiecewiseLinearDistribution& rhs)
            {
                return lhs.m_param == rhs.m_param;
            }
            friend bool operator!=(const PiecewiseLinearDistribution& lhs,
                                const PiecewiseLinearDistribution& rhs)
            {
                return !(lhs == rhs);
            }
 
        private:
            param_type m_param;
        };
 
        extern template class UniformRealDistribution<float>;
        extern template class UniformRealDistribution<double>;
 
        extern template class UniformIntDistribution<int>;
        extern template class UniformIntDistribution<unsigned int>;
        extern template class UniformIntDistribution<long>;
        extern template class UniformIntDistribution<unsigned long>;
        extern template class UniformIntDistribution<long long>;
        extern template class UniformIntDistribution<unsigned long long>;
 
        extern template class NormalDistribution<float>;
        extern template class NormalDistribution<double>;
 
        extern template class GammaDistribution<float>;
        extern template class GammaDistribution<double>;
 
        extern template class ExponentialDistribution<float>;
        extern template class ExponentialDistribution<double>;
 
        extern template class WeibullDistribution<float>;
        extern template class WeibullDistribution<double>;
 
        extern template class ExtremeValueDistribution<float>;
        extern template class ExtremeValueDistribution<double>;
 
        extern template class CauchyDistribution<float>;
        extern template class CauchyDistribution<double>;
 
        extern template class LognormalDistribution<float>;
        extern template class LognormalDistribution<double>;
 
        extern template class GeometricDistribution<int>;
        extern template class GeometricDistribution<unsigned int>;
        extern template class GeometricDistribution<long>;
        extern template class GeometricDistribution<unsigned long>;
 
        extern template class PoissonDistribution<int>;
        extern template class PoissonDistribution<unsigned int>;
        extern template class PoissonDistribution<long>;
 
        extern template class BinomialDistribution<int>;
        extern template class BinomialDistribution<unsigned int>;
 
        extern template class NegativeBinomialDistribution<int>;
        extern template class NegativeBinomialDistribution<unsigned int>;
 
        extern template class ChiSquaredDistribution<float>;
        extern template class ChiSquaredDistribution<double>;
 
        extern template class FisherFDistribution<float>;
        extern template class FisherFDistribution<double>;
 
        extern template class StudentTDistribution<float>;
        extern template class StudentTDistribution<double>;
 
        extern template class DiscreteDistribution<int>;
        extern template class DiscreteDistribution<unsigned int>;
        extern template class DiscreteDistribution<long>;
 
        extern template class PiecewiseConstantDistribution<float>;
        extern template class PiecewiseConstantDistribution<double>;
 
        extern template class PiecewiseLinearDistribution<float>;
        extern template class PiecewiseLinearDistribution<double>;
 
}