// Copyright (C) 2026 Kiyotsugu Arai
// SPDX-License-Identifier: LGPL-3.0-or-later

// distributions.hpp
// Additional probability distributions (5-3f)
//
// Continuous distributions:
//   ExponentialDistribution, GammaDistribution, BetaDistribution,
//   ChiSquaredDistribution, FDistribution, StudentTDistribution,
//   CauchyDistribution, RayleighDistribution, WeibullDistribution,
//   LognormalDistribution, LaplaceDistribution, ParetoDistribution,
//   LogisticDistribution, GumbelDistribution
//
// Discrete distributions:
//   BernoulliDistribution, PoissonDistribution, BinomialDistribution,
//   GeometricDistribution, NegativeBinomialDistribution

#ifndef SANGI_DISTRIBUTIONS_HPP
#define SANGI_DISTRIBUTIONS_HPP

#include <math/special/gamma.hpp>
#include <math/random/random_sampling.hpp>   // FloatingPointType + Float sampler foundation
#include <cmath>
#include <random>
#include <stdexcept>
#include <vector>
#include <numbers>
#include <limits>
#include <numeric>

namespace sangi {

// ================================================================
// Common utilities
// ================================================================

namespace detail {

template<typename T>
inline constexpr T pi_v = std::numbers::pi_v<T>;

// Factorial (for small integers)
inline double factorial(int n) {
    double r = 1.0;
    for (int i = 2; i <= n; ++i) r *= i;
    return r;
}

// Binomial coefficient
inline double binomCoeff(int n, int k) {
    if (k < 0 || k > n) return 0.0;
    if (k == 0 || k == n) return 1.0;
    if (k > n - k) k = n - k;
    double r = 1.0;
    for (int i = 0; i < k; ++i)
        r = r * (n - i) / (i + 1);
    return r;
}

} // namespace detail

// ================================================================
// ExponentialDistribution — exponential distribution Exp(λ)
// ================================================================

template<FloatingPointType T>
class ExponentialDistribution {
public:
    explicit ExponentialDistribution(T lambda = T(1)) : lambda_(lambda) {
        if (lambda <= T(0))
            throw std::invalid_argument("ExponentialDistribution: lambda must be > 0");
    }

    [[nodiscard]] T pdf(T x) const {
        if (x < T(0)) return T(0);
        return lambda_ * std::exp(-lambda_ * x);
    }

    [[nodiscard]] T cdf(T x) const {
        if (x < T(0)) return T(0);
        return T(1) - std::exp(-lambda_ * x);
    }

    [[nodiscard]] T quantile(T p) const {
        if (p < T(0) || p > T(1))
            throw std::invalid_argument("quantile: p must be in [0,1]");
        if (p == T(1)) return std::numeric_limits<T>::infinity();
        return -std::log(T(1) - p) / lambda_;
    }

    [[nodiscard]] T mean() const { return T(1) / lambda_; }
    [[nodiscard]] T variance() const { return T(1) / (lambda_ * lambda_); }

    [[nodiscard]] T sample(std::mt19937& gen) const {
        if constexpr (std::is_floating_point_v<T>) {
            std::exponential_distribution<T> d(lambda_);
            return d(gen);
        } else {
            return detail::sampleExponentialStd<T>(gen) / lambda_;
        }
    }

    [[nodiscard]] std::vector<T> sample(std::mt19937& gen, size_t n) const {
        std::vector<T> r(n);
        for (auto& v : r) v = sample(gen);
        return r;
    }

private:
    T lambda_;
};

// ================================================================
// GammaDistribution — gamma distribution Γ(α, β)
// ================================================================

template<FloatingPointType T>
class GammaDistribution {
public:
    // α = shape, β = rate (1/scale)
    GammaDistribution(T alpha = T(1), T beta = T(1)) : alpha_(alpha), beta_(beta) {
        if (alpha <= T(0) || beta <= T(0))
            throw std::invalid_argument("GammaDistribution: alpha, beta must be > 0");
    }

    [[nodiscard]] T pdf(T x) const {
        if (x <= T(0)) return T(0);
        return std::pow(beta_, alpha_) / detail::generic_tgamma(alpha_)
             * std::pow(x, alpha_ - T(1)) * std::exp(-beta_ * x);
    }

    [[nodiscard]] T cdf(T x) const {
        if (x <= T(0)) return T(0);
        return special::gammaP<T>(alpha_, beta_ * x);
    }

    [[nodiscard]] T mean() const { return alpha_ / beta_; }
    [[nodiscard]] T variance() const { return alpha_ / (beta_ * beta_); }

    [[nodiscard]] T sample(std::mt19937& gen) const {
        if constexpr (std::is_floating_point_v<T>) {
            std::gamma_distribution<T> d(alpha_, T(1) / beta_);
            return d(gen);
        } else {
            // scale-1 gamma / rate beta = gamma with scale 1/beta
            return detail::sampleGammaStd<T>(gen, alpha_) / beta_;
        }
    }

    [[nodiscard]] std::vector<T> sample(std::mt19937& gen, size_t n) const {
        std::vector<T> r(n);
        for (auto& v : r) v = sample(gen);
        return r;
    }

private:
    T alpha_, beta_;
};

// ================================================================
// BetaDistribution — beta distribution Beta(α, β)
// ================================================================

template<FloatingPointType T>
class BetaDistribution {
public:
    BetaDistribution(T alpha = T(1), T beta = T(1)) : alpha_(alpha), beta_(beta) {
        if (alpha <= T(0) || beta <= T(0))
            throw std::invalid_argument("BetaDistribution: alpha, beta must be > 0");
    }

    [[nodiscard]] T pdf(T x) const {
        if (x < T(0) || x > T(1)) return T(0);
        T B = detail::generic_tgamma(alpha_) * detail::generic_tgamma(beta_) / detail::generic_tgamma(alpha_ + beta_);
        return std::pow(x, alpha_ - T(1)) * std::pow(T(1) - x, beta_ - T(1)) / B;
    }

    [[nodiscard]] T cdf(T x) const {
        if (x <= T(0)) return T(0);
        if (x >= T(1)) return T(1);
        return special::betaRegularized<T>(x, alpha_, beta_);
    }

    [[nodiscard]] T mean() const { return alpha_ / (alpha_ + beta_); }
    [[nodiscard]] T variance() const {
        T ab = alpha_ + beta_;
        return (alpha_ * beta_) / (ab * ab * (ab + T(1)));
    }

    [[nodiscard]] T sample(std::mt19937& gen) const {
        if constexpr (std::is_floating_point_v<T>) {
            std::gamma_distribution<T> ga(alpha_, T(1));
            std::gamma_distribution<T> gb(beta_, T(1));
            T x = ga(gen), y = gb(gen);
            return x / (x + y);
        } else {
            T x = detail::sampleGammaStd<T>(gen, alpha_);
            T y = detail::sampleGammaStd<T>(gen, beta_);
            return x / (x + y);
        }
    }

    [[nodiscard]] std::vector<T> sample(std::mt19937& gen, size_t n) const {
        std::vector<T> r(n);
        for (auto& v : r) v = sample(gen);
        return r;
    }

private:
    T alpha_, beta_;
};

// ================================================================
// ChiSquaredDistribution — chi-squared distribution χ²(k)
// ================================================================

template<FloatingPointType T>
class ChiSquaredDistribution {
public:
    explicit ChiSquaredDistribution(T k = T(1)) : k_(k) {
        if (k <= T(0))
            throw std::invalid_argument("ChiSquaredDistribution: k must be > 0");
    }

    [[nodiscard]] T pdf(T x) const {
        if (x <= T(0)) return T(0);
        T half_k = k_ / T(2);
        return std::pow(x, half_k - T(1)) * std::exp(-x / T(2))
             / (std::pow(T(2), half_k) * detail::generic_tgamma(half_k));
    }

    [[nodiscard]] T cdf(T x) const {
        if (x <= T(0)) return T(0);
        return special::gammaP<T>(k_ / T(2), x / T(2));
    }

    [[nodiscard]] T mean() const { return k_; }
    [[nodiscard]] T variance() const { return T(2) * k_; }

    [[nodiscard]] T sample(std::mt19937& gen) const {
        if constexpr (std::is_floating_point_v<T>) {
            std::chi_squared_distribution<T> d(k_);
            return d(gen);
        } else {
            // χ²(k) = Gamma(shape k/2, scale 2)
            return detail::sampleGammaStd<T>(gen, k_ / T(2)) * T(2);
        }
    }

    [[nodiscard]] std::vector<T> sample(std::mt19937& gen, size_t n) const {
        std::vector<T> r(n);
        for (auto& v : r) v = sample(gen);
        return r;
    }

private:
    T k_;
};

// ================================================================
// FDistribution — F distribution F(d1, d2)
// ================================================================

template<FloatingPointType T>
class FDistribution {
public:
    FDistribution(T d1 = T(1), T d2 = T(1)) : d1_(d1), d2_(d2) {
        if (d1 <= T(0) || d2 <= T(0))
            throw std::invalid_argument("FDistribution: d1, d2 must be > 0");
    }

    [[nodiscard]] T pdf(T x) const {
        if (x <= T(0)) return T(0);
        T a = d1_ / T(2), b = d2_ / T(2);
        T num = std::pow(d1_ * x, a) * std::pow(d2_, b);
        T den = std::pow(d1_ * x + d2_, a + b);
        T B = detail::generic_tgamma(a) * detail::generic_tgamma(b) / detail::generic_tgamma(a + b);
        return num / (den * x * B);
    }

    [[nodiscard]] T cdf(T x) const {
        if (x <= T(0)) return T(0);
        T t = d1_ * x / (d1_ * x + d2_);
        return special::betaRegularized<T>(t, d1_ / T(2), d2_ / T(2));
    }

    [[nodiscard]] T mean() const {
        if (d2_ <= T(2)) return std::numeric_limits<T>::infinity();
        return d2_ / (d2_ - T(2));
    }

    [[nodiscard]] T variance() const {
        if (d2_ <= T(4)) return std::numeric_limits<T>::infinity();
        return T(2) * d2_ * d2_ * (d1_ + d2_ - T(2))
             / (d1_ * (d2_ - T(2)) * (d2_ - T(2)) * (d2_ - T(4)));
    }

    [[nodiscard]] T sample(std::mt19937& gen) const {
        if constexpr (std::is_floating_point_v<T>) {
            std::fisher_f_distribution<T> d(d1_, d2_);
            return d(gen);
        } else {
            // F = (χ²(d1)/d1) / (χ²(d2)/d2),  χ²(d) = 2·Gamma(d/2)
            T g1 = detail::sampleGammaStd<T>(gen, d1_ / T(2));
            T g2 = detail::sampleGammaStd<T>(gen, d2_ / T(2));
            return (g1 * d2_) / (g2 * d1_);
        }
    }

    [[nodiscard]] std::vector<T> sample(std::mt19937& gen, size_t n) const {
        std::vector<T> r(n);
        for (auto& v : r) v = sample(gen);
        return r;
    }

private:
    T d1_, d2_;
};

// ================================================================
// StudentTDistribution — Student's t distribution t(ν)
// ================================================================

template<FloatingPointType T>
class StudentTDistribution {
public:
    explicit StudentTDistribution(T nu = T(1)) : nu_(nu) {
        if (nu <= T(0))
            throw std::invalid_argument("StudentTDistribution: nu must be > 0");
    }

    [[nodiscard]] T pdf(T x) const {
        T a = (nu_ + T(1)) / T(2);
        T B = detail::generic_tgamma(a) / (std::sqrt(nu_ * detail::generic_pi<T>()) * detail::generic_tgamma(nu_ / T(2)));
        return B * std::pow(T(1) + x * x / nu_, -a);
    }

    [[nodiscard]] T cdf(T x) const {
        T t = nu_ / (nu_ + x * x);
        T Ix = special::betaRegularized<T>(t, nu_ / T(2), T(0.5));
        return (x >= T(0)) ? T(1) - Ix / T(2) : Ix / T(2);
    }

    [[nodiscard]] T mean() const {
        if (nu_ <= T(1)) return std::numeric_limits<T>::quiet_NaN();
        return T(0);
    }

    [[nodiscard]] T variance() const {
        if (nu_ <= T(1)) return std::numeric_limits<T>::quiet_NaN();
        if (nu_ <= T(2)) return std::numeric_limits<T>::infinity();
        return nu_ / (nu_ - T(2));
    }

    [[nodiscard]] T sample(std::mt19937& gen) const {
        if constexpr (std::is_floating_point_v<T>) {
            std::student_t_distribution<T> d(nu_);
            return d(gen);
        } else {
            // t(ν) = Z / sqrt(χ²(ν)/ν),  Z~N(0,1), χ²(ν)=2·Gamma(ν/2)
            T z = detail::sampleNormalStd<T>(gen);
            T v = detail::sampleGammaStd<T>(gen, nu_ / T(2)) * T(2);
            return z / detail::generic_sqrt(v / nu_);
        }
    }

    [[nodiscard]] std::vector<T> sample(std::mt19937& gen, size_t n) const {
        std::vector<T> r(n);
        for (auto& v : r) v = sample(gen);
        return r;
    }

private:
    T nu_;
};

// ================================================================
// CauchyDistribution — Cauchy distribution Cauchy(x0, γ)
// ================================================================

template<FloatingPointType T>
class CauchyDistribution {
public:
    CauchyDistribution(T x0 = T(0), T gamma = T(1)) : x0_(x0), gamma_(gamma) {
        if (gamma <= T(0))
            throw std::invalid_argument("CauchyDistribution: gamma must be > 0");
    }

    [[nodiscard]] T pdf(T x) const {
        T z = (x - x0_) / gamma_;
        return T(1) / (detail::generic_pi<T>() * gamma_ * (T(1) + z * z));
    }

    [[nodiscard]] T cdf(T x) const {
        return T(0.5) + std::atan((x - x0_) / gamma_) / detail::generic_pi<T>();
    }

    [[nodiscard]] T quantile(T p) const {
        if (p < T(0) || p > T(1))
            throw std::invalid_argument("quantile: p must be in [0,1]");
        return x0_ + gamma_ * std::tan(detail::generic_pi<T>() * (p - T(0.5)));
    }

    // Cauchy: mean/variance are undefined
    [[nodiscard]] T median() const { return x0_; }

    [[nodiscard]] T sample(std::mt19937& gen) const {
        if constexpr (std::is_floating_point_v<T>) {
            std::cauchy_distribution<T> d(x0_, gamma_);
            return d(gen);
        } else {
            return quantile(detail::uniform01<T>(gen));   // inverse CDF
        }
    }

    [[nodiscard]] std::vector<T> sample(std::mt19937& gen, size_t n) const {
        std::vector<T> r(n);
        for (auto& v : r) v = sample(gen);
        return r;
    }

private:
    T x0_, gamma_;
};

// ================================================================
// RayleighDistribution — Rayleigh distribution Rayleigh(σ)
// ================================================================

template<FloatingPointType T>
class RayleighDistribution {
public:
    explicit RayleighDistribution(T sigma = T(1)) : sigma_(sigma) {
        if (sigma <= T(0))
            throw std::invalid_argument("RayleighDistribution: sigma must be > 0");
    }

    [[nodiscard]] T pdf(T x) const {
        if (x < T(0)) return T(0);
        return (x / (sigma_ * sigma_)) * std::exp(-x * x / (T(2) * sigma_ * sigma_));
    }

    [[nodiscard]] T cdf(T x) const {
        if (x < T(0)) return T(0);
        return T(1) - std::exp(-x * x / (T(2) * sigma_ * sigma_));
    }

    [[nodiscard]] T mean() const { return sigma_ * std::sqrt(detail::generic_pi<T>() / T(2)); }
    [[nodiscard]] T variance() const {
        return (T(4) - detail::generic_pi<T>()) / T(2) * sigma_ * sigma_;
    }

    [[nodiscard]] T sample(std::mt19937& gen) const {
        if constexpr (std::is_floating_point_v<T>) {
            std::uniform_real_distribution<T> u(T(0), T(1));
            return sigma_ * std::sqrt(-T(2) * std::log(u(gen)));
        } else {
            return sigma_ * detail::generic_sqrt(-T(2) * detail::generic_log(detail::uniform01_pos<T>(gen)));
        }
    }

    [[nodiscard]] std::vector<T> sample(std::mt19937& gen, size_t n) const {
        std::vector<T> r(n);
        for (auto& v : r) v = sample(gen);
        return r;
    }

private:
    T sigma_;
};

// ================================================================
// WeibullDistribution — Weibull distribution Weibull(k, λ)
// ================================================================

template<FloatingPointType T>
class WeibullDistribution {
public:
    WeibullDistribution(T k = T(1), T lambda = T(1)) : k_(k), lambda_(lambda) {
        if (k <= T(0) || lambda <= T(0))
            throw std::invalid_argument("WeibullDistribution: k, lambda must be > 0");
    }

    [[nodiscard]] T pdf(T x) const {
        if (x < T(0)) return T(0);
        return (k_ / lambda_) * std::pow(x / lambda_, k_ - T(1))
             * std::exp(-std::pow(x / lambda_, k_));
    }

    [[nodiscard]] T cdf(T x) const {
        if (x < T(0)) return T(0);
        return T(1) - std::exp(-std::pow(x / lambda_, k_));
    }

    [[nodiscard]] T mean() const { return lambda_ * detail::generic_tgamma(T(1) + T(1) / k_); }
    [[nodiscard]] T variance() const {
        T g1 = detail::generic_tgamma(T(1) + T(1) / k_);
        T g2 = detail::generic_tgamma(T(1) + T(2) / k_);
        return lambda_ * lambda_ * (g2 - g1 * g1);
    }

    [[nodiscard]] T sample(std::mt19937& gen) const {
        if constexpr (std::is_floating_point_v<T>) {
            std::weibull_distribution<T> d(k_, lambda_);
            return d(gen);
        } else {
            // inverse CDF: λ·(-ln(1-u))^{1/k}
            return lambda_ * detail::generic_pow(
                -detail::generic_log(detail::uniform01_pos<T>(gen)), T(1) / k_);
        }
    }

    [[nodiscard]] std::vector<T> sample(std::mt19937& gen, size_t n) const {
        std::vector<T> r(n);
        for (auto& v : r) v = sample(gen);
        return r;
    }

private:
    T k_, lambda_;
};

// ================================================================
// LognormalDistribution — log-normal distribution LogNormal(μ, σ)
// ================================================================

template<FloatingPointType T>
class LognormalDistribution {
public:
    LognormalDistribution(T mu = T(0), T sigma = T(1)) : mu_(mu), sigma_(sigma) {
        if (sigma <= T(0))
            throw std::invalid_argument("LognormalDistribution: sigma must be > 0");
    }

    [[nodiscard]] T pdf(T x) const {
        if (x <= T(0)) return T(0);
        T z = (std::log(x) - mu_) / sigma_;
        return std::exp(-T(0.5) * z * z) / (x * sigma_ * std::sqrt(T(2) * detail::generic_pi<T>()));
    }

    [[nodiscard]] T cdf(T x) const {
        if (x <= T(0)) return T(0);
        return T(0.5) * detail::generic_erfc(-(std::log(x) - mu_) / (sigma_ * std::sqrt(T(2))));
    }

    [[nodiscard]] T mean() const { return std::exp(mu_ + sigma_ * sigma_ / T(2)); }
    [[nodiscard]] T variance() const {
        T s2 = sigma_ * sigma_;
        return (std::exp(s2) - T(1)) * std::exp(T(2) * mu_ + s2);
    }

    [[nodiscard]] T sample(std::mt19937& gen) const {
        if constexpr (std::is_floating_point_v<T>) {
            std::lognormal_distribution<T> d(mu_, sigma_);
            return d(gen);
        } else {
            // exp(μ + σ·N(0,1))
            return detail::generic_exp(mu_ + sigma_ * detail::sampleNormalStd<T>(gen));
        }
    }

    [[nodiscard]] std::vector<T> sample(std::mt19937& gen, size_t n) const {
        std::vector<T> r(n);
        for (auto& v : r) v = sample(gen);
        return r;
    }

private:
    T mu_, sigma_;
};

// ================================================================
// LaplaceDistribution — Laplace distribution Laplace(μ, b)
// ================================================================

template<FloatingPointType T>
class LaplaceDistribution {
public:
    LaplaceDistribution(T mu = T(0), T b = T(1)) : mu_(mu), b_(b) {
        if (b <= T(0))
            throw std::invalid_argument("LaplaceDistribution: b must be > 0");
    }

    [[nodiscard]] T pdf(T x) const {
        return std::exp(-std::abs(x - mu_) / b_) / (T(2) * b_);
    }

    [[nodiscard]] T cdf(T x) const {
        if (x < mu_)
            return T(0.5) * std::exp((x - mu_) / b_);
        else
            return T(1) - T(0.5) * std::exp(-(x - mu_) / b_);
    }

    [[nodiscard]] T quantile(T p) const {
        if (p < T(0) || p > T(1))
            throw std::invalid_argument("quantile: p must be in [0,1]");
        if (p < T(0.5))
            return mu_ + b_ * std::log(T(2) * p);
        else
            return mu_ - b_ * std::log(T(2) * (T(1) - p));
    }

    [[nodiscard]] T mean() const { return mu_; }
    [[nodiscard]] T variance() const { return T(2) * b_ * b_; }

    [[nodiscard]] T sample(std::mt19937& gen) const {
        if constexpr (std::is_floating_point_v<T>) {
            std::uniform_real_distribution<T> u(T(0), T(1));
            T p = u(gen) - T(0.5);
            return mu_ - b_ * std::copysign(std::log(T(1) - T(2) * std::abs(p)), p);
        } else {
            return quantile(detail::uniform01<T>(gen));   // inverse CDF
        }
    }

    [[nodiscard]] std::vector<T> sample(std::mt19937& gen, size_t n) const {
        std::vector<T> r(n);
        for (auto& v : r) v = sample(gen);
        return r;
    }

private:
    T mu_, b_;
};

// ================================================================
// ParetoDistribution — Pareto distribution Pareto(α, xm)
// ================================================================

template<FloatingPointType T>
class ParetoDistribution {
public:
    ParetoDistribution(T alpha = T(1), T xm = T(1)) : alpha_(alpha), xm_(xm) {
        if (alpha <= T(0) || xm <= T(0))
            throw std::invalid_argument("ParetoDistribution: alpha, xm must be > 0");
    }

    [[nodiscard]] T pdf(T x) const {
        if (x < xm_) return T(0);
        return alpha_ * std::pow(xm_, alpha_) / std::pow(x, alpha_ + T(1));
    }

    [[nodiscard]] T cdf(T x) const {
        if (x < xm_) return T(0);
        return T(1) - std::pow(xm_ / x, alpha_);
    }

    [[nodiscard]] T mean() const {
        if (alpha_ <= T(1)) return std::numeric_limits<T>::infinity();
        return alpha_ * xm_ / (alpha_ - T(1));
    }

    [[nodiscard]] T variance() const {
        if (alpha_ <= T(2)) return std::numeric_limits<T>::infinity();
        return xm_ * xm_ * alpha_ / ((alpha_ - T(1)) * (alpha_ - T(1)) * (alpha_ - T(2)));
    }

    [[nodiscard]] T sample(std::mt19937& gen) const {
        if constexpr (std::is_floating_point_v<T>) {
            std::uniform_real_distribution<T> u(T(0), T(1));
            return xm_ / std::pow(u(gen), T(1) / alpha_);
        } else {
            return xm_ / detail::generic_pow(detail::uniform01_pos<T>(gen), T(1) / alpha_);
        }
    }

    [[nodiscard]] std::vector<T> sample(std::mt19937& gen, size_t n) const {
        std::vector<T> r(n);
        for (auto& v : r) v = sample(gen);
        return r;
    }

private:
    T alpha_, xm_;
};

// ================================================================
// LogisticDistribution — logistic distribution Logistic(μ, s)
// ================================================================

template<FloatingPointType T>
class LogisticDistribution {
public:
    LogisticDistribution(T mu = T(0), T s = T(1)) : mu_(mu), s_(s) {
        if (s <= T(0))
            throw std::invalid_argument("LogisticDistribution: s must be > 0");
    }

    [[nodiscard]] T pdf(T x) const {
        T z = std::exp(-(x - mu_) / s_);
        return z / (s_ * (T(1) + z) * (T(1) + z));
    }

    [[nodiscard]] T cdf(T x) const {
        return T(1) / (T(1) + std::exp(-(x - mu_) / s_));
    }

    [[nodiscard]] T quantile(T p) const {
        if (p <= T(0) || p >= T(1))
            throw std::invalid_argument("quantile: p must be in (0,1)");
        return mu_ + s_ * std::log(p / (T(1) - p));
    }

    [[nodiscard]] T mean() const { return mu_; }
    [[nodiscard]] T variance() const {
        return s_ * s_ * detail::generic_pi<T>() * detail::generic_pi<T>() / T(3);
    }

    [[nodiscard]] T sample(std::mt19937& gen) const {
        if constexpr (std::is_floating_point_v<T>) {
            std::uniform_real_distribution<T> u(T(0), T(1));
            T p = u(gen);
            return mu_ + s_ * std::log(p / (T(1) - p));
        } else {
            T p = detail::uniform01_pos<T>(gen);   // (0,1]
            return mu_ + s_ * detail::generic_log(p / (T(1) - p));
        }
    }

    [[nodiscard]] std::vector<T> sample(std::mt19937& gen, size_t n) const {
        std::vector<T> r(n);
        for (auto& v : r) v = sample(gen);
        return r;
    }

private:
    T mu_, s_;
};

// ================================================================
// GumbelDistribution — Gumbel distribution (Type-1 extreme-value distribution)
// ================================================================

template<FloatingPointType T>
class GumbelDistribution {
public:
    GumbelDistribution(T mu = T(0), T beta = T(1)) : mu_(mu), beta_(beta) {
        if (beta <= T(0))
            throw std::invalid_argument("GumbelDistribution: beta must be > 0");
    }

    [[nodiscard]] T pdf(T x) const {
        T z = (x - mu_) / beta_;
        return std::exp(-(z + std::exp(-z))) / beta_;
    }

    [[nodiscard]] T cdf(T x) const {
        T z = (x - mu_) / beta_;
        return std::exp(-std::exp(-z));
    }

    [[nodiscard]] T quantile(T p) const {
        if (p <= T(0) || p >= T(1))
            throw std::invalid_argument("quantile: p must be in (0,1)");
        return mu_ - beta_ * std::log(-std::log(p));
    }

    // Euler-Mascheroni constant γ ≈ 0.5772
    [[nodiscard]] T mean() const { return mu_ + beta_ * T(0.5772156649015329); }
    [[nodiscard]] T variance() const {
        return detail::generic_pi<T>() * detail::generic_pi<T>() * beta_ * beta_ / T(6);
    }

    [[nodiscard]] T sample(std::mt19937& gen) const {
        if constexpr (std::is_floating_point_v<T>) {
            std::extreme_value_distribution<T> d(mu_, beta_);
            return d(gen);
        } else {
            // inverse CDF: μ - β·ln(-ln(u)),  u∈(0,1]
            T u = detail::uniform01_pos<T>(gen);
            return mu_ - beta_ * detail::generic_log(-detail::generic_log(u));
        }
    }

    [[nodiscard]] std::vector<T> sample(std::mt19937& gen, size_t n) const {
        std::vector<T> r(n);
        for (auto& v : r) v = sample(gen);
        return r;
    }

private:
    T mu_, beta_;
};

// ================================================================
// Discrete distributions
// ================================================================

// ================================================================
// BernoulliDistribution — Bernoulli distribution Bernoulli(p)
// ================================================================

template<FloatingPointType T>
class BernoulliDistribution {
public:
    explicit BernoulliDistribution(T p = T(0.5)) : p_(p) {
        if (p < T(0) || p > T(1))
            throw std::invalid_argument("BernoulliDistribution: p must be in [0,1]");
    }

    [[nodiscard]] T pmf(int k) const {
        if (k == 0) return T(1) - p_;
        if (k == 1) return p_;
        return T(0);
    }

    [[nodiscard]] T cdf(int k) const {
        if (k < 0) return T(0);
        if (k < 1) return T(1) - p_;
        return T(1);
    }

    [[nodiscard]] T mean() const { return p_; }
    [[nodiscard]] T variance() const { return p_ * (T(1) - p_); }

    [[nodiscard]] int sample(std::mt19937& gen) const {
        if constexpr (std::is_floating_point_v<T>) {
            std::bernoulli_distribution d(static_cast<double>(p_));
            return d(gen) ? 1 : 0;
        } else {
            return (detail::uniform01<T>(gen) < p_) ? 1 : 0;
        }
    }

    [[nodiscard]] std::vector<int> sample(std::mt19937& gen, size_t n) const {
        std::vector<int> r(n);
        for (auto& v : r) v = sample(gen);
        return r;
    }

private:
    T p_;
};

// ================================================================
// PoissonDistribution — Poisson distribution Poisson(λ)
// ================================================================

template<FloatingPointType T>
class PoissonDistribution {
public:
    explicit PoissonDistribution(T lambda = T(1)) : lambda_(lambda) {
        if (lambda <= T(0))
            throw std::invalid_argument("PoissonDistribution: lambda must be > 0");
    }

    [[nodiscard]] T pmf(int k) const {
        if (k < 0) return T(0);
        return std::exp(-lambda_ + k * std::log(lambda_) - detail::generic_lgamma(T(k + 1)));
    }

    [[nodiscard]] T cdf(int k) const {
        if (k < 0) return T(0);
        // P(X ≤ k) = Q(k+1, λ) = 1 - P(k+1, λ) (upper regularized incomplete gamma)
        return T(1) - special::gammaP<T>(T(k + 1), lambda_);
    }

    [[nodiscard]] T mean() const { return lambda_; }
    [[nodiscard]] T variance() const { return lambda_; }

    [[nodiscard]] int sample(std::mt19937& gen) const {
        if constexpr (std::is_floating_point_v<T>) {
            std::poisson_distribution<int> d(static_cast<double>(lambda_));
            return d(gen);
        } else {
            return detail::samplePoissonKnuth<T>(gen, lambda_);
        }
    }

    [[nodiscard]] std::vector<int> sample(std::mt19937& gen, size_t n) const {
        std::vector<int> r(n);
        for (auto& v : r) v = sample(gen);
        return r;
    }

private:
    T lambda_;
};

// ================================================================
// BinomialDistribution — binomial distribution Binomial(n, p)
// ================================================================

template<FloatingPointType T>
class BinomialDistribution {
public:
    BinomialDistribution(int n = 1, T p = T(0.5)) : n_(n), p_(p) {
        if (n < 0)
            throw std::invalid_argument("BinomialDistribution: n must be >= 0");
        if (p < T(0) || p > T(1))
            throw std::invalid_argument("BinomialDistribution: p must be in [0,1]");
    }

    [[nodiscard]] T pmf(int k) const {
        if (k < 0 || k > n_) return T(0);
        return T(detail::binomCoeff(n_, k))
             * std::pow(p_, T(k)) * std::pow(T(1) - p_, T(n_ - k));
    }

    [[nodiscard]] T cdf(int k) const {
        if (k < 0) return T(0);
        if (k >= n_) return T(1);
        // I_{1-p}(n-k, k+1)
        return special::betaRegularized<T>(T(1) - p_, T(n_ - k), T(k + 1));
    }

    [[nodiscard]] T mean() const { return T(n_) * p_; }
    [[nodiscard]] T variance() const { return T(n_) * p_ * (T(1) - p_); }

    [[nodiscard]] int sample(std::mt19937& gen) const {
        if constexpr (std::is_floating_point_v<T>) {
            std::binomial_distribution<int> d(n_, static_cast<double>(p_));
            return d(gen);
        } else {
            // sum of n Bernoulli trials (O(n))
            int count = 0;
            for (int i = 0; i < n_; ++i)
                if (detail::uniform01<T>(gen) < p_) ++count;
            return count;
        }
    }

    [[nodiscard]] std::vector<int> sample(std::mt19937& gen, size_t n) const {
        std::vector<int> r(n);
        for (auto& v : r) v = sample(gen);
        return r;
    }

private:
    int n_;
    T p_;
};

// ================================================================
// GeometricDistribution — geometric distribution Geometric(p)
// ================================================================

// P(X=k) = (1-p)^k * p, k = 0, 1, 2, ... (number of failures)
template<FloatingPointType T>
class GeometricDistribution {
public:
    explicit GeometricDistribution(T p = T(0.5)) : p_(p) {
        if (p <= T(0) || p > T(1))
            throw std::invalid_argument("GeometricDistribution: p must be in (0,1]");
    }

    [[nodiscard]] T pmf(int k) const {
        if (k < 0) return T(0);
        return std::pow(T(1) - p_, T(k)) * p_;
    }

    [[nodiscard]] T cdf(int k) const {
        if (k < 0) return T(0);
        return T(1) - std::pow(T(1) - p_, T(k + 1));
    }

    [[nodiscard]] T mean() const { return (T(1) - p_) / p_; }
    [[nodiscard]] T variance() const { return (T(1) - p_) / (p_ * p_); }

    [[nodiscard]] int sample(std::mt19937& gen) const {
        if constexpr (std::is_floating_point_v<T>) {
            std::geometric_distribution<int> d(static_cast<double>(p_));
            return d(gen);
        } else {
            if (p_ >= T(1)) return 0;
            // inverse CDF: k = floor( ln(u) / ln(1-p) ),  u∈(0,1]
            T k = detail::generic_floor(
                detail::generic_log(detail::uniform01_pos<T>(gen)) /
                detail::generic_log(T(1) - p_));
            return static_cast<int>(detail::to_int64(k));
        }
    }

    [[nodiscard]] std::vector<int> sample(std::mt19937& gen, size_t n) const {
        std::vector<int> r(n);
        for (auto& v : r) v = sample(gen);
        return r;
    }

private:
    T p_;
};

// ================================================================
// NegativeBinomialDistribution — negative binomial distribution NB(r, p)
// ================================================================

// P(X=k) = C(k+r-1, k) * p^r * (1-p)^k, k = 0, 1, 2, ...
template<FloatingPointType T>
class NegativeBinomialDistribution {
public:
    NegativeBinomialDistribution(T r = T(1), T p = T(0.5)) : r_(r), p_(p) {
        if (r <= T(0))
            throw std::invalid_argument("NegativeBinomialDistribution: r must be > 0");
        if (p <= T(0) || p > T(1))
            throw std::invalid_argument("NegativeBinomialDistribution: p must be in (0,1]");
    }

    [[nodiscard]] T pmf(int k) const {
        if (k < 0) return T(0);
        return std::exp(detail::generic_lgamma(T(k) + r_) - detail::generic_lgamma(T(k + 1)) - detail::generic_lgamma(r_)
                      + r_ * std::log(p_) + T(k) * std::log(T(1) - p_));
    }

    [[nodiscard]] T cdf(int k) const {
        if (k < 0) return T(0);
        return special::betaRegularized<T>(p_, r_, T(k + 1));
    }

    [[nodiscard]] T mean() const { return r_ * (T(1) - p_) / p_; }
    [[nodiscard]] T variance() const { return r_ * (T(1) - p_) / (p_ * p_); }

    [[nodiscard]] int sample(std::mt19937& gen) const {
        if constexpr (std::is_floating_point_v<T>) {
            std::negative_binomial_distribution<int> d(
                static_cast<int>(std::round(static_cast<double>(r_))),
                static_cast<double>(p_));
            return d(gen);
        } else {
            // gamma-Poisson mixture: λ ~ Gamma(r, (1-p)/p), then Poisson(λ).
            // (correct even for real r, and more exact than STL's integer-r rounding)
            T lambda = detail::sampleGammaStd<T>(gen, r_) * (T(1) - p_) / p_;
            return detail::samplePoissonKnuth<T>(gen, lambda);
        }
    }

    [[nodiscard]] std::vector<int> sample(std::mt19937& gen, size_t n) const {
        std::vector<int> r(n);
        for (auto& v : r) v = sample(gen);
        return r;
    }

private:
    T r_, p_;
};

} // namespace sangi

#endif // SANGI_DISTRIBUTIONS_HPP