stable.hpp

Go to the documentation of this file.

#pragma once
 
#include <algorithm>
#include <array>
#include <functional>
#include <vector>
 
#include <randlib/math/constants.hpp>
#include <randlib/math/functions.hpp>
#include <randlib/math/quadrature.hpp>
#include <randlib/core/errors.hpp>
#include <randlib/core/continuous.hpp>
 
namespace randlib
{
 
namespace detail
{
 
inline double stableTanHalfPi(double alpha)
{
    const double reduced = std::remainder(0.5 * alpha, 1.0);
    if (std::fabs(reduced) == 0.25)
        return std::copysign(1.0, reduced);
    return std::tan(randlib::consts::pi * reduced);
}
 
inline double stableAtan(double value)
{
    if (std::fabs(value) == 1.0)
        return std::copysign(static_cast<double>(0.5L * randlib::consts::half_pi), value);
    return math::atan(value);
}
 
} // namespace detail
 
inline double stableTailQuantile(double p, bool upperTail, double alpha, double beta, double gamma, double mu)
{
    double weight = upperTail ? 1 + beta : 1 - beta;
    if (alpha >= 2 || weight <= 0)
        return NAN;
    double C = (alpha == 1)
                   ? randlib::consts::inv_pi
                   : std::tgamma(alpha) * std::sin(0.5 * randlib::consts::pi * alpha) * randlib::consts::inv_pi;
    double distance = gamma * std::pow(C * weight / p, 1.0 / alpha);
    return upperTail ? mu + distance : mu - distance;
}
 
template <typename RealType>
class StableDistribution : virtual public ContinuousDistribution<RealType>
{
    RealType variate(detail::EngineRef g) const override { return (*this)(g); }
    void sample(detail::EngineRef g, std::span<RealType> outputData) const override;
 
protected:
    enum class Parameterization
    {
        S0, 
        S1 
    };
 
    double alpha = 2; 
    double beta = 0; 
    double mu = 0; 
    double gamma = randlib::consts::sqrt2; 
    double logGamma = 0.5 * randlib::consts::ln2; 
 
private:
    double alphaInv = 0.5; 
    double zeta = 0; 
    double omega = 0; 
    double xi = 0; 
    double alpha_alpham1 = 2; 
    double logGammaPi_2 = randlib::consts::lnpi - 1.5 * randlib::consts::ln2; 
    Parameterization parameterization = Parameterization::S1;
    double locationParameter = 0;
 
    static constexpr double BIG_NUMBER = 1e9; 
    static constexpr double ALMOST_TWO = 1.99999; 
    static constexpr double S0_PROBABILITY_UNITY_TOLERANCE = 1.1e-8;
    static constexpr double S0_DIRECT_KERNEL_ALPHA_TOLERANCE = 1e-5;
    static constexpr double S0_DIRECT_KERNEL_RANGE = 100;
    static constexpr double S0_SAMPLING_UNITY_TOLERANCE = 1e-10;
 
    enum DISTRIBUTION_TYPE
    {
        NORMAL, 
        LEVY, 
        CAUCHY, 
        UNITY_EXPONENT, 
        GENERAL 
    };
 
    DISTRIBUTION_TYPE distributionType = NORMAL; 
 
protected:
    double pdfCoef =
        0.5 * (randlib::consts::ln2 + randlib::consts::lnpi); 
    double pdftailBound = INFINITY; 
    double cdftailBound = INFINITY; 
 
    StableDistribution(double exponent, double skewness, double scale = 1, double location = 0,
                       Parameterization selectedParameterization = Parameterization::S1);
 
public:
    virtual ~StableDistribution() {}
 
    [[nodiscard]] Support SupportType() const override;
    [[nodiscard]] RealType MinValue() const override;
    [[nodiscard]] RealType MaxValue() const override;
 
private:
    void parametersVerification(double exponent, double skewness, double scale);
    void setParametersForNormal();
    void setParametersForCauchy();
    void setParametersForLevy();
    void setParametersForUnityExponent();
    void setParametersForGeneralExponent();
    static double s0LocationShift(double exponent, double skewness, double scale);
 
protected:
    void SetParameters(double exponent, double skewness, double scale = 1, double location = 0);
 
public:
    void SetLocation(double location);
    void SetScale(double scale);
    void SetParameterization(Parameterization value) noexcept;
 
    [[nodiscard]] inline double GetExponent() const noexcept { return alpha; }
    [[nodiscard]] inline double GetSkewness() const noexcept { return beta; }
    [[nodiscard]] inline double GetScale() const noexcept { return gamma; }
    [[nodiscard]] inline double GetLocation() const noexcept { return locationParameter; }
    [[nodiscard]] inline Parameterization GetParameterization() const noexcept { return parameterization; }
    [[nodiscard]] inline double GetLogScale() const noexcept { return logGamma; }
 
protected:
    double pdfNormal(RealType x) const;
    double logpdfNormal(RealType x) const;
    double pdfCauchy(RealType x) const;
    double logpdfCauchy(RealType x) const;
    double pdfLevy(RealType x) const;
    double logpdfLevy(RealType x) const;
 
private:
    static double fastpdfExponentiation(double u);
 
    double pdfShortTailExpansionForUnityExponent(double x) const;
    std::complex<double> unityTailIntegral(int power, int logarithmPower, double distance) const;
    double unityTailExpansion(double distance, bool integrated) const;
 
    double limitCaseForIntegrandAuxForUnityExponent(double theta, double xAdj) const;
    double integrandAuxForUnityExponent(double theta, double xAdj) const;
    double integrandForUnityExponent(double theta, double xAdj) const;
    double pdfForUnityExponent(double x) const;
    static double rootOfMonotone(const std::function<double(double)> &u, double a, double b);
    static double integrateAroundRoot(const std::function<double(double)> &integrand,
                                      const std::function<double(double)> &u, double a, double b, double root);
    double peakForUnityExponent(double xAdj) const;
 
    DoublePair seriesZeroParams{};
 
    double pdfShortTailExpansionForGeneralExponent(double logX) const;
    double cdfShortTailExpansion(double logX) const;
    double pdfAtZero() const;
    double pdfSeriesExpansionAtZero(double logX, double xiAdj, int k) const;
    double pdfSeriesExpansionAtInf(double logX, double xiAdj) const;
    double pdfTaylorExpansionTailNearCauchy(double x) const;
    double limitCaseForIntegrandAuxForGeneralExponent(double theta, double xiAdj) const;
    double integrandAuxForGeneralExponent(double theta, double xAdj, double xiAdj) const;
    double integrandFoGeneralExponent(double theta, double xAdj, double xiAdj) const;
    double pdfForGeneralExponent(double x) const;
    [[nodiscard]] bool useS0UnityLimit() const noexcept;
    double pdfS0UnityLimit(double x) const;
    double cdfS0UnityLimit(double x, bool complement) const;
    [[nodiscard]] bool useS0NearUnityKernel(double x) const noexcept;
    double s0NearUnityPhase(double u, double standardizedX) const;
    double pdfS0NearUnity(double x) const;
    double cdfS0NearUnity(double x, bool complement) const;
 
public:
    [[nodiscard]] double f(const RealType &x) const override;
    [[nodiscard]] double logf(const RealType &x) const override;
 
protected:
    double cdfNormal(RealType x) const;
    double cdfNormalCompl(RealType x) const;
    double cdfCauchy(RealType x) const;
    double cdfCauchyCompl(RealType x) const;
    double cdfLevy(RealType x) const;
    double cdfLevyCompl(RealType x) const;
 
private:
    static double fastcdfExponentiation(double u);
    static double fastccdfExponentiation(double u);
    double cdfAtZero(double xiAdj) const;
    double cdfForUnityExponent(double x, bool complement) const;
    double cdfSeriesExpansionAtZero(double logX, double xiAdj, int k) const;
    double cdfSeriesExpansionAtInf(double logX, double xiAdj) const;
    double cdfIntegralRepresentation(double logX, double xiAdj, bool complement) const;
    double cdfForGeneralExponent(double x) const;
    double ccdfForGeneralExponent(double x) const;
 
public:
    [[nodiscard]] double F(const RealType &x) const override;
    [[nodiscard]] double S(const RealType &x) const override;
 
private:
    template <std::uniform_random_bit_generator G>
    double variateForUnityExponent(G &g) const;
    template <std::uniform_random_bit_generator G>
    double variateForS0UnityLimit(G &g) const;
    template <std::uniform_random_bit_generator G>
    double variateForGeneralExponent(G &g) const;
    template <std::uniform_random_bit_generator G>
    double variateForExponentEqualOneHalf(G &g) const;
 
public:
    template <std::uniform_random_bit_generator G>
    [[nodiscard]] RealType operator()(G &g) const;
 
    [[nodiscard]] long double Mean() const override;
    [[nodiscard]] long double Variance() const override;
    [[nodiscard]] RealType Median() const override;
    [[nodiscard]] RealType Mode() const override;
    [[nodiscard]] long double Skewness() const override;
    [[nodiscard]] long double ExcessKurtosis() const override;
 
protected:
    RealType quantileNormal(double p) const;
    RealType quantileNormal1m(double p) const;
    RealType quantileCauchy(double p) const;
    RealType quantileCauchy1m(double p) const;
    RealType quantileLevy(double p) const;
    RealType quantileLevy1m(double p) const;
 
private:
    RealType searchTailQuantile(double p, bool upperTail) const;
    RealType quantileImpl(double p) const override;
    RealType quantileImpl1m(double p) const override;
 
protected:
    std::complex<double> cfNormal(double t) const;
    std::complex<double> cfCauchy(double t) const;
    std::complex<double> cfLevy(double t) const;
 
private:
    std::complex<double> CFImpl(double t) const override;
};
 
template <typename RealType>
class StableRand final : public StableDistribution<RealType>
{
public:
    using Parameterization = typename StableDistribution<RealType>::Parameterization;
 
    StableRand(double exponent = 2, double skewness = 0, double scale = 1, double location = 0,
               Parameterization selectedParameterization = Parameterization::S1)
        : StableDistribution<RealType>(exponent, skewness, scale, location, selectedParameterization)
    {
    }
    [[nodiscard]] std::string Name() const override;
 
    void SetExponent(double exponent);
    void SetSkewness(double skewness);
};
 
template <typename RealType>
class HoltsmarkRand final : public StableDistribution<RealType>
{
public:
    HoltsmarkRand(double scale = 1, double location = 0) : StableDistribution<RealType>(1.5, 0.0, scale, location) {}
    [[nodiscard]] std::string Name() const override;
};
 
template <typename RealType>
class LandauRand final : public StableDistribution<RealType>
{
public:
    LandauRand(double scale = 1, double location = 0) : StableDistribution<RealType>(1.0, 1.0, scale, location) {}
    [[nodiscard]] std::string Name() const override;
};
 
} // namespace randlib
 
#include <randlib/continuous/levy.hpp>
#include <randlib/continuous/cauchy.hpp>
#include <randlib/continuous/normal.hpp>
#include <randlib/continuous/uniform.hpp>
#include <randlib/continuous/exponential.hpp>
 
namespace randlib
{
 
template <typename RealType>
StableDistribution<RealType>::StableDistribution(double exponent, double skewness, double scale, double location,
                                                 Parameterization selectedParameterization)
    : parameterization(selectedParameterization)
{
    SetParameters(exponent, skewness, scale, location);
}
 
template <typename RealType>
Support StableDistribution<RealType>::SupportType() const
{
    if (alpha < 1) {
        if (beta == 1)
            return Support::RightSemiFinite;
        if (beta == -1)
            return Support::LeftSemiFinite;
    }
    return Support::Infinite;
}
 
template <typename RealType>
RealType StableDistribution<RealType>::MinValue() const
{
    return (alpha < 1 && beta == 1) ? mu : -INFINITY;
}
 
template <typename RealType>
RealType StableDistribution<RealType>::MaxValue() const
{
    return (alpha < 1 && beta == -1) ? mu : INFINITY;
}
 
template <typename RealType>
void StableDistribution<RealType>::parametersVerification(double exponent, double skewness, double scale)
{
    if (!(exponent > 0.0 && exponent <= 2.0))
        throw std::invalid_argument(
            "Stable distribution: exponent should be in the interval (0, 2], but it's equal to " +
            std::to_string(exponent));
    if (!(std::fabs(skewness) <= 1.0))
        throw std::invalid_argument(
            "Stable distribution: skewness should be in the interval [-1, 1], but it's equal to " +
            std::to_string(skewness));
    if (!(scale > 0.0))
        throw std::invalid_argument("Stable distribution: scale should be positive, but it's equal to " +
                                    std::to_string(scale));
}
 
template <typename RealType>
void StableDistribution<RealType>::setParametersForNormal()
{
    distributionType = NORMAL;
    pdfCoef = randlib::consts::ln2 + logGamma + 0.5 * randlib::consts::lnpi;
}
 
template <typename RealType>
void StableDistribution<RealType>::setParametersForCauchy()
{
    distributionType = CAUCHY;
    pdfCoef = -logGamma - randlib::consts::lnpi;
}
 
template <typename RealType>
void StableDistribution<RealType>::setParametersForLevy()
{
    distributionType = LEVY;
    pdfCoef = logGamma - randlib::consts::ln2 - randlib::consts::lnpi;
}
 
template <typename RealType>
void StableDistribution<RealType>::setParametersForUnityExponent()
{
    distributionType = UNITY_EXPONENT;
    pdfCoef = 0.5 / (gamma * std::fabs(beta));
    pdftailBound = 0; // not in the use for now
    logGammaPi_2 = randlib::consts::lnpi - randlib::consts::ln2 - logGamma;
}
 
template <typename RealType>
void StableDistribution<RealType>::setParametersForGeneralExponent()
{
    distributionType = GENERAL;
    if (beta != 0.0) {
        zeta = -beta * detail::stableTanHalfPi(alpha);
        const double zetaSquared = zeta * zeta;
        omega = (zetaSquared == 1.0) ? 0.5 * alphaInv * static_cast<double>(randlib::consts::ln2)
                                     : 0.5 * alphaInv * std::log1pl(zetaSquared);
        xi = alphaInv * detail::stableAtan(-zeta);
    } else {
        zeta = omega = xi = 0.0;
    }
    pdfCoef = randlib::consts::inv_pi * std::fabs(alpha_alpham1) / gamma;
    pdftailBound = 3.0 / (1.0 + alpha) * randlib::consts::ln10;
    cdftailBound = 3.0 / alpha * randlib::consts::ln10;
    if (alpha <= ALMOST_TWO) {
        seriesZeroParams.first = std::round(std::min(alpha * alpha * 40 + 1, 10.0));
        seriesZeroParams.second =
            -(alphaInv * 1.5 + 0.5) * randlib::consts::ln10; 
    } else {
        seriesZeroParams.first = 85;
        seriesZeroParams.second = randlib::consts::ln2 + randlib::consts::ln3; 
    }
}
 
template <typename RealType>
double StableDistribution<RealType>::s0LocationShift(double exponent, double skewness, double scale)
{
    if (skewness == 0.0 || exponent == 2.0)
        return 0.0;
    if (exponent == 1.0)
        return randlib::consts::two_inv_pi * skewness * scale * std::log(scale);
    return skewness * scale * detail::stableTanHalfPi(exponent);
}
 
template <typename RealType>
void StableDistribution<RealType>::SetParameters(double exponent, double skewness, double scale, double location)
{
    if (!std::isfinite(location))
        throw std::invalid_argument("Stable distribution: location should be finite");
    parametersVerification(exponent, skewness, scale);
    alpha = exponent;
    alphaInv = 1.0 / alpha;
    beta = skewness;
    locationParameter = location;
    mu = locationParameter -
         ((parameterization == Parameterization::S0) ? s0LocationShift(exponent, skewness, scale) : 0.0);
    gamma = scale;
    logGamma = std::log(gamma);
    alpha_alpham1 = alpha / (alpha - 1.0);
 
    if (alpha == 2.0)
        return setParametersForNormal();
    if (alpha == 1.0)
        return (beta == 0.0) ? setParametersForCauchy() : setParametersForUnityExponent();
    if (alpha == 0.5 && std::fabs(beta) == 1.0)
        return setParametersForLevy();
    return setParametersForGeneralExponent();
}
 
template <typename RealType>
void StableDistribution<RealType>::SetLocation(double location)
{
    SetParameters(alpha, beta, gamma, location);
}
 
template <typename RealType>
void StableDistribution<RealType>::SetScale(double scale)
{
    SetParameters(alpha, beta, scale, GetLocation());
}
 
template <typename RealType>
void StableDistribution<RealType>::SetParameterization(Parameterization value) noexcept
{
    parameterization = value;
    locationParameter = mu + ((parameterization == Parameterization::S0) ? s0LocationShift(alpha, beta, gamma) : 0.0);
}
 
template <typename RealType>
double StableDistribution<RealType>::pdfNormal(RealType x) const
{
    return std::exp(logpdfNormal(x));
}
 
template <typename RealType>
double StableDistribution<RealType>::logpdfNormal(RealType x) const
{
    double y = x - mu;
    y *= 0.5 / gamma;
    y *= y;
    y += pdfCoef;
    return -y;
}
 
template <typename RealType>
double StableDistribution<RealType>::pdfCauchy(RealType x) const
{
    double y = x - mu;
    y *= y;
    y /= gamma;
    y += gamma;
    return randlib::consts::inv_pi / y;
}
 
template <typename RealType>
double StableDistribution<RealType>::logpdfCauchy(RealType x) const
{
    double x0 = x - mu;
    x0 /= gamma;
    double xSq = x0 * x0;
    return pdfCoef - std::log1pl(xSq);
}
 
template <typename RealType>
double StableDistribution<RealType>::pdfLevy(RealType x) const
{
    return (x <= mu) ? 0.0 : std::exp(logpdfLevy(x));
}
 
template <typename RealType>
double StableDistribution<RealType>::logpdfLevy(RealType x) const
{
    double x0 = x - mu;
    if (x0 <= 0.0)
        return -INFINITY;
    double y = gamma / x0;
    y += 3 * std::log(x0);
    y -= pdfCoef;
    return -0.5 * y;
}
 
template <typename RealType>
double StableDistribution<RealType>::fastpdfExponentiation(double u)
{
    return (u < -700) ? std::exp(u) : std::exp(u - std::exp(u));
}
 
template <typename RealType>
double StableDistribution<RealType>::pdfShortTailExpansionForUnityExponent(double x) const
{
    if (x > 10)
        return 0.0;
    double xm1 = x - 1.0;
    double y = 0.5 * xm1 - std::exp(xm1);
    y -= 0.5 * (randlib::consts::ln2 + randlib::consts::lnpi);
    return std::exp(y);
}
 
template <typename RealType>
std::complex<double> StableDistribution<RealType>::unityTailIntegral(int power, int logarithmPower,
                                                                     double distance) const
{
    const double absoluteDistance = std::fabs(distance);
    const double direction = std::copysign(1.0, distance);
    const double halfPi = static_cast<double>(randlib::consts::half_pi);
    const std::complex<double> imaginary(0, 1);
    const std::complex<double> integral = std::tgamma(power + 1.0) * std::pow(absoluteDistance, -power - 1.0) *
                                          std::exp(imaginary * direction * halfPi * (power + 1.0));
    if (logarithmPower == 0)
        return integral;
 
    std::array<std::complex<double>, 9> derivatives{};
    derivatives[1] = math::digamma(power + 1.0) - std::log(absoluteDistance) + imaginary * direction * halfPi;
    for (int order = 2; order <= logarithmPower; ++order) {
        double harmonic = 0;
        for (int term = 1; term <= power; ++term)
            harmonic += std::pow(static_cast<double>(term), -order);
        const double factorial = std::tgamma(static_cast<double>(order));
        const double sign = (order % 2 == 0) ? 1.0 : -1.0;
        derivatives[order] = sign * factorial * (math::riemannZeta(static_cast<double>(order)) - harmonic);
    }
 
    std::array<std::complex<double>, 9> bell{};
    bell[0] = 1;
    for (int order = 1; order <= logarithmPower; ++order) {
        double combination = 1;
        for (int part = 1; part <= order; ++part) {
            bell[order] += combination * bell[order - part] * derivatives[part];
            combination *= static_cast<double>(order - part) / part;
        }
    }
    return integral * bell[logarithmPower];
}
 
template <typename RealType>
double StableDistribution<RealType>::unityTailExpansion(double distance, bool integrated) const
{
    static constexpr int ORDER = 8;
    const std::complex<double> imaginary(0, 1);
    const double logarithmicCoefficient = randlib::consts::two_inv_pi * beta * gamma;
    std::complex<double> sum = 0;
    for (int totalOrder = 1; totalOrder <= ORDER; ++totalOrder) {
        for (int logarithmPower = 0; logarithmPower <= totalOrder; ++logarithmPower) {
            const int exponentialPower = totalOrder - logarithmPower;
            const std::complex<double> coefficient =
                std::pow(-gamma, exponentialPower) / std::tgamma(exponentialPower + 1.0) *
                std::pow(imaginary * logarithmicCoefficient, logarithmPower) / std::tgamma(logarithmPower + 1.0);
            const int power = totalOrder - (integrated ? 1 : 0);
            sum += coefficient * unityTailIntegral(power, logarithmPower, distance);
        }
    }
    const double inversePi = static_cast<double>(randlib::consts::inv_pi);
    if (!integrated)
        return std::max(0.0, sum.real() * inversePi);
    return std::clamp(-std::copysign(1.0, distance) * sum.imag() * inversePi, 0.0, 1.0);
}
 
template <typename RealType>
double StableDistribution<RealType>::limitCaseForIntegrandAuxForUnityExponent(double theta, double xAdj) const
{
    if (theta > 0.0) {
        if (beta > 0.0)
            return BIG_NUMBER;
        return (beta == -1) ? xAdj - 1.0 : -BIG_NUMBER;
    }
    if (beta < 0.0)
        return BIG_NUMBER;
    return (beta == 1) ? xAdj - 1.0 : -BIG_NUMBER;
}
 
template <typename RealType>
double StableDistribution<RealType>::integrandAuxForUnityExponent(double theta, double xAdj) const
{
    if (std::fabs(theta) >= randlib::consts::half_pi)
        return limitCaseForIntegrandAuxForUnityExponent(theta, xAdj);
    if (theta == 0.0)
        return xAdj + randlib::consts::lnpi - randlib::consts::ln2;
    double thetaAdj = (randlib::consts::half_pi + beta * theta) / std::cos(theta);
    double u = std::log(thetaAdj);
    u += thetaAdj * std::sin(theta) / beta;
    return std::isfinite(u) ? u + xAdj : limitCaseForIntegrandAuxForUnityExponent(theta, xAdj);
}
 
template <typename RealType>
double StableDistribution<RealType>::integrandForUnityExponent(double theta, double xAdj) const
{
    if (std::fabs(theta) >= randlib::consts::half_pi)
        return 0.0;
    double u = integrandAuxForUnityExponent(theta, xAdj);
    return fastpdfExponentiation(u);
}
 
template <typename RealType>
double StableDistribution<RealType>::rootOfMonotone(const std::function<double(double)> &u, double a, double b)
{
    bool increasing = u(b) > u(a);
    for (int i = 0; i < 200; ++i) {
        double middle = 0.5 * (a + b);
        if (!(middle > a && middle < b))
            break;
        bool belowRoot = (u(middle) < 0) == increasing;
        (belowRoot ? a : b) = middle;
    }
    return 0.5 * (a + b);
}
 
template <typename RealType>
double StableDistribution<RealType>::integrateAroundRoot(const std::function<double(double)> &integrand,
                                                         const std::function<double(double)> &u, double a, double b,
                                                         double root)
{
    double h = 1e-7 * std::max(std::min(root - a, b - root), 1e-300);
    double slope = (u(root + h) - u(root - h)) / (2 * h);
    double width = (std::isfinite(slope) && slope != 0.0) ? 1.0 / std::fabs(slope) : b - a;
    std::vector<double> points = {a, root, b};
    for (double scale : {1e-3, 1e-2, 1e-1, 1.0, 10.0, 100.0, 1000.0}) {
        points.push_back(root - scale * width);
        points.push_back(root + scale * width);
    }
    std::sort(points.begin(), points.end());
    double sum = 0.0, left = a;
    for (double point : points) {
        if (point <= left || point > b)
            continue;
        sum += math::integrate(integrand, left, point, 1e-11, 262144);
        left = point;
    }
    return sum;
}
 
template <typename RealType>
double StableDistribution<RealType>::peakForUnityExponent(double xAdj) const
{
    return rootOfMonotone([this, xAdj](double theta) { return integrandAuxForUnityExponent(theta, xAdj); },
                          -randlib::consts::half_pi, randlib::consts::half_pi);
}
 
template <typename RealType>
double StableDistribution<RealType>::pdfForUnityExponent(double x) const
{
    if (std::fabs(x - mu) >= 100 * gamma)
        return unityTailExpansion(x - mu, false);
    double xSt = (x - mu) / gamma;
    double xAdj = -randlib::consts::half_pi * xSt / beta - logGammaPi_2;
    double theta0 = peakForUnityExponent(xAdj);
    auto integrand = [this, xAdj](double theta) { return integrandForUnityExponent(theta, xAdj); };
    auto u = [this, xAdj](double theta) { return integrandAuxForUnityExponent(theta, xAdj); };
    return pdfCoef * integrateAroundRoot(integrand, u, -randlib::consts::half_pi, randlib::consts::half_pi, theta0);
}
 
template <typename RealType>
double StableDistribution<RealType>::pdfShortTailExpansionForGeneralExponent(double logX) const
{
    double logAlpha = std::log(alpha);
    double log1mAlpha = (alpha < 1) ? std::log1pl(-alpha) : std::log(alpha - 1);
    double temp = logX - logAlpha;
    double y = std::exp(alpha_alpham1 * temp);
    y *= -std::fabs(1.0 - alpha);
    double z = (1.0 - 0.5 * alpha) / (alpha - 1) * temp;
    z -= 0.5 * (randlib::consts::ln2 + randlib::consts::lnpi + logAlpha + log1mAlpha);
    return std::exp(y + z);
}
 
template <typename RealType>
double StableDistribution<RealType>::cdfShortTailExpansion(double logX) const
{
    double logU = alpha / (alpha - 1.0) * (logX - std::log(alpha));
    double u = std::exp(logU);
    double y = -(1.0 - alpha) * u -
               0.5 * (randlib::consts::ln2 + randlib::consts::lnpi + std::log(alpha) + std::log1p(-alpha) + logU);
    return std::exp(y);
}
 
template <typename RealType>
double StableDistribution<RealType>::pdfAtZero() const
{
    double y0 = 0.0;
    if (beta == 0.0)
        y0 = std::tgammal(alphaInv);
    else {
        y0 = std::lgammal(alphaInv) - omega;
        y0 = std::exp(y0) * std::cos(xi);
    }
    return y0 * randlib::consts::inv_pi / alpha;
}
 
template <typename RealType>
double StableDistribution<RealType>::pdfSeriesExpansionAtZero(double logX, double xiAdj, int k) const
{
    double y0 = pdfAtZero();
    double sum = 0.0;
    if (beta == 0.0) {
        for (int n = 1; n <= k; ++n) {
            int n2 = n + n;
            double term = std::lgammal((n2 + 1) / alpha);
            term += n2 * logX;
            term -= math::lfact(n2);
            term = std::exp(term);
            sum += (n & 1) ? -term : term;
        }
    } else {
        double rhoPi_alpha = randlib::consts::half_pi + xiAdj;
        for (int n = 1; n <= k; ++n) {
            int np1 = n + 1;
            double term = std::lgammal(np1 * alphaInv);
            term += n * logX;
            term -= math::lfact(n);
            term = std::exp(term - omega);
            term *= std::sin(np1 * rhoPi_alpha);
            sum += (n & 1) ? -term : term;
        }
    }
    return y0 + sum * randlib::consts::inv_pi / alpha;
}
 
template <typename RealType>
double StableDistribution<RealType>::pdfSeriesExpansionAtInf(double logX, double xiAdj) const
{
    static constexpr int k = 10; 
    double rhoPi = randlib::consts::half_pi + xiAdj;
    rhoPi *= alpha;
    double sum = 0.0;
    for (int n = 1; n <= k; ++n) {
        double aux = n * alpha + 1.0;
        double term = std::lgammal(aux);
        term -= aux * logX;
        term -= math::lfact(n);
        term = std::exp(term - omega);
        term *= std::sin(rhoPi * n);
        sum += (n & 1) ? term : -term;
    }
    return randlib::consts::inv_pi * sum;
}
 
template <typename RealType>
double StableDistribution<RealType>::pdfTaylorExpansionTailNearCauchy(double x) const
{
    double xSq = x * x;
    double y = 1.0 + xSq;
    double ySq = y * y;
    double z = math::atan(x);
    double zSq = z * z;
    double logY = std::log1pl(xSq);
    double alpham1 = alpha - 1.0;
    double temp = 1.0 - randlib::consts::egamma - 0.5 * logY;
    double f_a = temp;
    f_a *= xSq - 1.0;
    f_a += 2 * x * z;
    f_a /= ySq;
    static constexpr long double M_PI_SQ_6 = 1.64493406684822643647l; 
    double f_aa1 = M_PI_SQ_6;
    f_aa1 += temp * temp;
    f_aa1 -= 1.0 + z * z;
    f_aa1 *= xSq * xSq - 6.0 * xSq + 1.0;
    double f_aa2 = 0.5 + temp;
    f_aa2 *= z;
    f_aa2 *= 8 * x * (xSq - 1.0);
    double f_aa3 = (1.0 - 3 * xSq) * temp;
    f_aa3 -= x * y * z;
    f_aa3 += f_aa3;
    double f_aa = f_aa1 + f_aa2 + f_aa3;
    f_aa /= std::pow(y, 3);
    static constexpr int gammaTable[] = {1, 2, 6};
    static constexpr long double gammaDerTable[] = {1.0 - randlib::consts::egamma, 3.0 - 2.0 * randlib::consts::egamma,
                                                    11.0 - 6.0 * randlib::consts::egamma};
    static constexpr long double gammaSecDerTable[] = {0.82368066085287938958l, 2.49292999190269305794l,
                                                       11.1699273161019477314l};
    static constexpr long double digammaTable[] = {1.0 - randlib::consts::egamma, 1.5 - randlib::consts::egamma,
                                                   11.0 / 6 - randlib::consts::egamma};
    static constexpr long double digammaDerTable[] = {M_PI_SQ_6 - 1.0, M_PI_SQ_6 - 1.25, M_PI_SQ_6 - 49.0 / 36};
    static constexpr long double digammaSecDerTable[] = {-0.40411380631918857080l, -0.15411380631918857080l,
                                                         -0.08003973224511449673l};
    double gTable[] = {0, 0, 0};
    for (int i = 0; i < 3; ++i) {
        double g_11 = 0.25 * gammaTable[i] * logY * logY;
        g_11 -= gammaDerTable[i] * logY;
        g_11 += gammaSecDerTable[i];
        double aux = digammaTable[i] - 0.5 * logY;
        double g_12 = aux;
        double zip2 = z * (i + 2);
        double cosZNu = std::cos(zip2), zSinZNu = z * std::sin(zip2);
        g_12 *= cosZNu;
        g_12 -= zSinZNu;
        double g_1 = g_11 * g_12;
        double g_21 = -gammaTable[i] * logY + 2 * gammaDerTable[i];
        double g_22 = -zSinZNu * aux;
        g_22 -= zSq * cosZNu;
        g_22 += cosZNu * digammaDerTable[i];
        double g_2 = g_21 * g_22;
        double g_3 = -zSq * cosZNu * aux;
        g_3 -= 2 * zSinZNu * digammaDerTable[i];
        g_3 += zSq * zSinZNu;
        g_3 += cosZNu * digammaSecDerTable[i];
        g_3 *= gammaTable[i];
        double g = g_1 + g_2 + g_3;
        g *= std::pow(y, -0.5 * i - 1);
        gTable[i] = g;
    }
    double f_aaa = -gTable[0] + 3 * gTable[1] - gTable[2];
    double tail = f_a * alpham1;
    tail += 0.5 * f_aa * alpham1 * alpham1;
    tail += std::pow(alpham1, 3) * f_aaa / 6.0;
    tail /= randlib::consts::pi;
    return tail;
}
 
template <typename RealType>
double StableDistribution<RealType>::limitCaseForIntegrandAuxForGeneralExponent(double theta, double xiAdj) const
{
    if (theta < 0.5 * (randlib::consts::half_pi - xiAdj))
        return alpha < 1 ? -BIG_NUMBER : BIG_NUMBER;
    return alpha < 1 ? BIG_NUMBER : -BIG_NUMBER;
}
 
template <typename RealType>
double StableDistribution<RealType>::integrandAuxForGeneralExponent(double theta, double xAdj, double xiAdj) const
{
    const double halfPi = static_cast<double>(randlib::consts::half_pi);
    if (theta >= halfPi && alpha > 1.0 && std::fabs(beta) == 1.0 && xiAdj > 0.0)
        return xAdj + std::log(alpha - 1.0) - alpha_alpham1 * std::log(alpha);
    if (std::fabs(theta) >= halfPi || theta <= -xiAdj)
        return limitCaseForIntegrandAuxForGeneralExponent(theta, xiAdj);
    double thetaAdj = alpha * (theta + xiAdj);
    double sinThetaAdj = std::sin(thetaAdj);
    double y = std::log(std::cos(theta));
    y -= alpha * std::log(sinThetaAdj);
    y /= alpha - 1.0;
    y += std::log(std::cos(thetaAdj - theta));
    y += xAdj;
    return std::isfinite(y) ? y : limitCaseForIntegrandAuxForGeneralExponent(theta, xiAdj);
}
 
template <typename RealType>
double StableDistribution<RealType>::integrandFoGeneralExponent(double theta, double xAdj, double xiAdj) const
{
    if (std::fabs(theta) >= randlib::consts::half_pi)
        return 0.0;
    if (theta <= -xiAdj)
        return 0.0;
    double u = integrandAuxForGeneralExponent(theta, xAdj, xiAdj);
    return fastpdfExponentiation(u);
}
 
template <typename RealType>
double StableDistribution<RealType>::pdfForGeneralExponent(double x) const
{
    double xSt = (x - mu) / gamma;
    double absXSt = xSt;
    double xiAdj = xi;
    if (xSt > 0) {
        if (alpha < 1 && beta == -1)
            return 0.0;
    } else {
        if (alpha < 1 && beta == 1)
            return 0.0;
        absXSt = -xSt;
        xiAdj = -xi;
    }
 
    if (beta == 0.0 && std::fabs(alpha - 1.0) < 0.01)
        return pdfCauchy(x) + pdfTaylorExpansionTailNearCauchy(absXSt) / gamma;
 
    if (xSt == 0.0)
        return pdfAtZero() / gamma;
 
    double logAbsX = std::log(absXSt) - omega;
 
    if (logAbsX < seriesZeroParams.second) {
        if (alpha < 1 && std::fabs(beta) == 1)
            return pdfShortTailExpansionForGeneralExponent(logAbsX);
        return pdfSeriesExpansionAtZero(logAbsX, xiAdj, seriesZeroParams.first) / gamma;
    }
 
    if (logAbsX > pdftailBound && alpha <= ALMOST_TWO) {
        double skewnessTowardsX = (xSt > 0) ? beta : -beta;
        if (alpha > 1 && skewnessTowardsX == -1)
            return pdfShortTailExpansionForGeneralExponent(logAbsX);
        return pdfSeriesExpansionAtInf(logAbsX, xiAdj) / gamma;
    }
 
    double xAdj = alpha_alpham1 * logAbsX;
 
    auto u = [this, xAdj, xiAdj](double theta) { return integrandAuxForGeneralExponent(theta, xAdj, xiAdj); };
    auto integrand = [this, xAdj, xiAdj](double theta) { return integrandFoGeneralExponent(theta, xAdj, xiAdj); };
    double theta0 = rootOfMonotone(u, -xiAdj, randlib::consts::half_pi);
    double integral = integrateAroundRoot(integrand, u, -xiAdj, randlib::consts::half_pi, theta0);
    double res = pdfCoef * integral / absXSt;
 
    if (alpha <= ALMOST_TWO)
        return res;
 
    double alphap1 = alpha + 1.0;
    double tail = std::lgammal(alphap1);
    tail -= alphap1 * logAbsX;
    tail = std::exp(tail);
    tail *= (1.0 - 0.5 * alpha) / gamma;
    return std::max(tail, res);
}
 
template <typename RealType>
bool StableDistribution<RealType>::useS0UnityLimit() const noexcept
{
    return parameterization == Parameterization::S0 && alpha != 1.0 && beta != 0.0 &&
           std::fabs(alpha - 1.0) <= S0_PROBABILITY_UNITY_TOLERANCE;
}
 
template <typename RealType>
double StableDistribution<RealType>::pdfS0UnityLimit(double x) const
{
    const double unityMu = locationParameter - randlib::consts::two_inv_pi * beta * gamma * logGamma;
    if (std::fabs(x - unityMu) >= 100 * gamma)
        return unityTailExpansion(x - unityMu, false);
    const double standardizedX = (x - unityMu) / gamma;
    const double unityLogGammaPi2 = randlib::consts::lnpi - randlib::consts::ln2 - logGamma;
    const double xAdjusted = -randlib::consts::half_pi * standardizedX / beta - unityLogGammaPi2;
    const double theta0 = peakForUnityExponent(xAdjusted);
    const auto integrand = [this, xAdjusted](double theta) { return integrandForUnityExponent(theta, xAdjusted); };
    const auto auxiliary = [this, xAdjusted](double theta) { return integrandAuxForUnityExponent(theta, xAdjusted); };
    const double coefficient = 0.5 / (gamma * std::fabs(beta));
    return coefficient *
           integrateAroundRoot(integrand, auxiliary, -randlib::consts::half_pi, randlib::consts::half_pi, theta0);
}
 
template <typename RealType>
double StableDistribution<RealType>::cdfS0UnityLimit(double x, bool complement) const
{
    const double unityMu = locationParameter - randlib::consts::two_inv_pi * beta * gamma * logGamma;
    if (std::fabs(x - unityMu) >= 100 * gamma) {
        const double tail = unityTailExpansion(x - unityMu, true);
        return complement == (x < unityMu) ? 1.0 - tail : tail;
    }
    const double standardizedX = (x - unityMu) / gamma;
    const double unityLogGammaPi2 = randlib::consts::lnpi - randlib::consts::ln2 - logGamma;
    const double xAdjusted = -randlib::consts::half_pi * standardizedX / beta - unityLogGammaPi2;
    const bool useComplement = complement == (beta > 0);
    const auto integrand = [this, xAdjusted, useComplement](double theta) {
        const double value = integrandAuxForUnityExponent(theta, xAdjusted);
        return useComplement ? fastccdfExponentiation(value) : fastcdfExponentiation(value);
    };
    const double theta0 = peakForUnityExponent(xAdjusted);
    const auto auxiliary = [this, xAdjusted](double theta) { return integrandAuxForUnityExponent(theta, xAdjusted); };
    return randlib::consts::inv_pi *
           integrateAroundRoot(integrand, auxiliary, -randlib::consts::half_pi, randlib::consts::half_pi, theta0);
}
 
template <typename RealType>
bool StableDistribution<RealType>::useS0NearUnityKernel(double x) const noexcept
{
    const double standardizedX = (x - locationParameter) / gamma;
    return parameterization == Parameterization::S0 && alpha != 1.0 && beta != 0.0 &&
           std::fabs(alpha - 1.0) < S0_DIRECT_KERNEL_ALPHA_TOLERANCE &&
           std::fabs(standardizedX) <= S0_DIRECT_KERNEL_RANGE;
}
 
template <typename RealType>
double StableDistribution<RealType>::s0NearUnityPhase(double u, double standardizedX) const
{
    const double epsilon = alpha - 1.0;
    const double difference = std::expm1(epsilon * std::log(u));
    const double skewPhase = -beta * u * difference / std::tan(randlib::consts::half_pi * epsilon);
    return skewPhase - standardizedX * u;
}
 
template <typename RealType>
double StableDistribution<RealType>::pdfS0NearUnity(double x) const
{
    const double standardizedX = (x - locationParameter) / gamma;
    const double cutoff = std::pow(36.0, 1.0 / alpha);
    const double integral = math::integrateTanhSinh(
        [this, standardizedX](double u) {
            const double attenuation = std::exp(-std::pow(u, alpha));
            return attenuation * std::cos(s0NearUnityPhase(u, standardizedX));
        },
        0.0, cutoff, 2e-13);
    return std::max(0.0, integral * static_cast<double>(randlib::consts::inv_pi) / gamma);
}
 
template <typename RealType>
double StableDistribution<RealType>::cdfS0NearUnity(double x, bool complement) const
{
    const double standardizedX = (x - locationParameter) / gamma;
    const double cutoff = std::pow(36.0, 1.0 / alpha);
    const double integral = math::integrateTanhSinh(
        [this, standardizedX](double u) {
            const double attenuation = std::exp(-std::pow(u, alpha));
            return attenuation * std::sin(s0NearUnityPhase(u, standardizedX)) / u;
        },
        0.0, cutoff, 2e-13);
    const double correction = integral * static_cast<double>(randlib::consts::inv_pi);
    return std::clamp(complement ? 0.5 + correction : 0.5 - correction, 0.0, 1.0);
}
 
template <typename RealType>
double StableDistribution<RealType>::f(const RealType &x) const
{
    if (useS0UnityLimit())
        return pdfS0UnityLimit(x);
    if (useS0NearUnityKernel(x))
        return pdfS0NearUnity(x);
    switch (distributionType) {
    case NORMAL:
        return pdfNormal(x);
    case CAUCHY:
        return pdfCauchy(x);
    case LEVY:
        return (beta > 0) ? pdfLevy(x) : pdfLevy(2 * mu - x);
    case UNITY_EXPONENT:
        return pdfForUnityExponent(x);
    case GENERAL:
        return pdfForGeneralExponent(x);
    default:
        throw std::runtime_error("Stable distribution: invalid distribution type");
    }
}
 
template <typename RealType>
double StableDistribution<RealType>::logf(const RealType &x) const
{
    if (useS0UnityLimit())
        return std::log(pdfS0UnityLimit(x));
    if (useS0NearUnityKernel(x))
        return std::log(pdfS0NearUnity(x));
    switch (distributionType) {
    case NORMAL:
        return logpdfNormal(x);
    case CAUCHY:
        return logpdfCauchy(x);
    case LEVY:
        return (beta > 0) ? logpdfLevy(x) : logpdfLevy(2 * mu - x);
    case UNITY_EXPONENT:
        return std::log(pdfForUnityExponent(x));
    case GENERAL:
        return std::log(pdfForGeneralExponent(x));
    default:
        throw std::runtime_error("Stable distribution: invalid distribution type");
    }
}
 
template <typename RealType>
double StableDistribution<RealType>::cdfNormal(RealType x) const
{
    double y = mu - x;
    y *= 0.5 / gamma;
    return 0.5 * std::erfc(y);
}
 
template <typename RealType>
double StableDistribution<RealType>::cdfNormalCompl(RealType x) const
{
    double y = x - mu;
    y *= 0.5 / gamma;
    return 0.5 * std::erfc(y);
}
 
template <typename RealType>
double StableDistribution<RealType>::cdfCauchy(RealType x) const
{
    double x0 = x - mu;
    x0 /= gamma;
    return (x0 < -1.0) ? randlib::consts::inv_pi * std::atan(-1.0 / x0) : 0.5 + randlib::consts::inv_pi * std::atan(x0);
}
 
template <typename RealType>
double StableDistribution<RealType>::cdfCauchyCompl(RealType x) const
{
    double x0 = mu - x;
    x0 /= gamma;
    return (x0 < -1.0) ? randlib::consts::inv_pi * std::atan(-1.0 / x0) : 0.5 + randlib::consts::inv_pi * std::atan(x0);
}
 
template <typename RealType>
double StableDistribution<RealType>::cdfLevy(RealType x) const
{
    if (x <= mu)
        return 0;
    double y = x - mu;
    y += y;
    y = gamma / y;
    y = std::sqrt(y);
    return std::erfc(y);
}
 
template <typename RealType>
double StableDistribution<RealType>::cdfLevyCompl(RealType x) const
{
    if (x <= mu)
        return 1.0;
    double y = x - mu;
    y += y;
    y = gamma / y;
    y = std::sqrt(y);
    return std::erf(y);
}
 
template <typename RealType>
double StableDistribution<RealType>::fastcdfExponentiation(double u)
{
    return std::exp(-std::exp(u));
}
 
template <typename RealType>
double StableDistribution<RealType>::fastccdfExponentiation(double u)
{
    return (u > 5.0) ? 1.0 : -std::expm1(-std::exp(u));
}
 
template <typename RealType>
double StableDistribution<RealType>::cdfAtZero(double xiAdj) const
{
    return 0.5 - randlib::consts::inv_pi * xiAdj;
}
 
template <typename RealType>
double StableDistribution<RealType>::cdfForUnityExponent(double x, bool complement) const
{
    if (std::fabs(x - mu) >= 100 * gamma) {
        const double tail = unityTailExpansion(x - mu, true);
        return complement == (x < mu) ? 1.0 - tail : tail;
    }
    double xSt = (x - mu) / gamma;
    double xAdj = -randlib::consts::half_pi * xSt / beta - logGammaPi_2;
    bool useComplement = complement == (beta > 0);
    auto integrand = [this, xAdj, useComplement](double theta) {
        double u = integrandAuxForUnityExponent(theta, xAdj);
        return useComplement ? fastccdfExponentiation(u) : fastcdfExponentiation(u);
    };
    double theta0 = peakForUnityExponent(xAdj);
    auto u = [this, xAdj](double theta) { return integrandAuxForUnityExponent(theta, xAdj); };
    return randlib::consts::inv_pi *
           integrateAroundRoot(integrand, u, -randlib::consts::half_pi, randlib::consts::half_pi, theta0);
}
 
template <typename RealType>
double StableDistribution<RealType>::cdfSeriesExpansionAtZero(double logX, double xiAdj, int k) const
{
    double y0 = cdfAtZero(xiAdj);
    double sum = 0.0;
    if (beta == 0.0) {
        for (int m = 0; m <= k; ++m) {
            int m2p1 = 2 * m + 1;
            double term = std::lgammal(m2p1 * alphaInv);
            term += m2p1 * logX;
            term -= math::lfact(m2p1);
            term = std::exp(term);
            sum += (m & 1) ? -term : term;
        }
    } else {
        double rhoPi_alpha = randlib::consts::half_pi + xiAdj;
        for (int n = 1; n <= k; ++n) {
            double term = std::lgammal(n * alphaInv);
            term += n * logX;
            term -= math::lfact(n);
            term = std::exp(term);
            term *= std::sin(n * rhoPi_alpha);
            sum += (n & 1) ? term : -term;
        }
    }
    return y0 + sum * randlib::consts::inv_pi * alphaInv;
}
 
template <typename RealType>
double StableDistribution<RealType>::cdfSeriesExpansionAtInf(double logX, double xiAdj) const
{
    static constexpr int k = 10; 
    double rhoPi = randlib::consts::half_pi + xiAdj;
    rhoPi *= alpha;
    double sum = 0.0;
    for (int n = 1; n <= k; ++n) {
        double aux = n * alpha;
        double term = std::lgammal(aux);
        term -= aux * logX;
        term -= math::lfact(n);
        term = std::exp(term);
        term *= std::sin(rhoPi * n);
        sum += (n & 1) ? term : -term;
    }
    return randlib::consts::inv_pi * sum;
}
 
template <typename RealType>
double StableDistribution<RealType>::cdfIntegralRepresentation(double logX, double xiAdj, bool complement) const
{
    double xAdj = alpha_alpham1 * logX;
    auto u = [this, xAdj, xiAdj](double theta) { return integrandAuxForGeneralExponent(theta, xAdj, xiAdj); };
    auto integrand = [&u, complement](double theta) {
        return complement ? fastccdfExponentiation(u(theta)) : fastcdfExponentiation(u(theta));
    };
    double theta0 = rootOfMonotone(u, -xiAdj, randlib::consts::half_pi);
    return static_cast<double>(randlib::consts::inv_pi) *
           integrateAroundRoot(integrand, u, -xiAdj, randlib::consts::half_pi, theta0);
}
 
template <typename RealType>
double StableDistribution<RealType>::cdfForGeneralExponent(double x) const
{
    double xSt = (x - mu) / gamma; 
    if (xSt == 0) {
        if (alpha < 1 && std::fabs(beta) == 1)
            return (beta > 0) ? 0.0 : 1.0;
        return cdfAtZero(xi);
    }
    if (xSt > 0.0) {
        if (alpha > 1.0)
            return 1.0 - ccdfForGeneralExponent(x);
        if (beta == -1.0)
            return 1.0;
        double logAbsX = std::log(xSt) - omega;
        if (logAbsX > cdftailBound)
            return 1.0 - cdfSeriesExpansionAtInf(logAbsX, xi);
        if (logAbsX < seriesZeroParams.second)
            return (beta == 1.0) ? cdfShortTailExpansion(logAbsX)
                                 : cdfSeriesExpansionAtZero(logAbsX, xi, seriesZeroParams.first);
        double atZero = (beta == 1.0) ? 0.0 : cdfAtZero(xi);
        return atZero + cdfIntegralRepresentation(logAbsX, xi, false);
    }
    double logAbsX = std::log(-xSt) - omega;
    if (logAbsX < seriesZeroParams.second) {
        if (alpha < 1.0 && beta == -1.0)
            return 1.0 - cdfShortTailExpansion(logAbsX);
        return 1.0 - cdfSeriesExpansionAtZero(logAbsX, -xi, seriesZeroParams.first);
    }
    if (logAbsX > cdftailBound)
        return cdfSeriesExpansionAtInf(logAbsX, -xi);
    if (alpha > 1.0)
        return cdfIntegralRepresentation(logAbsX, -xi, false);
    return (beta == 1.0) ? 0.0 : cdfIntegralRepresentation(logAbsX, -xi, true);
}
 
template <typename RealType>
double StableDistribution<RealType>::ccdfForGeneralExponent(double x) const
{
    double xSt = (x - mu) / gamma; 
    if (xSt <= 0)
        return 1.0 - cdfForGeneralExponent(x);
    double logAbsX = std::log(xSt) - omega;
    if (logAbsX < seriesZeroParams.second) {
        if (alpha < 1.0 && beta == 1.0)
            return 1.0 - cdfShortTailExpansion(logAbsX);
        if (alpha < 1.0 && beta == -1.0)
            return cdfShortTailExpansion(logAbsX);
        return 1.0 - cdfSeriesExpansionAtZero(logAbsX, xi, seriesZeroParams.first);
    }
    if (logAbsX > cdftailBound)
        return cdfSeriesExpansionAtInf(logAbsX, xi);
    if (alpha > 1.0)
        return cdfIntegralRepresentation(logAbsX, xi, false);
    return (beta == -1.0) ? 0.0 : cdfIntegralRepresentation(logAbsX, xi, true);
}
 
template <typename RealType>
double StableDistribution<RealType>::F(const RealType &x) const
{
    if (useS0UnityLimit())
        return cdfS0UnityLimit(x, false);
    if (useS0NearUnityKernel(x))
        return cdfS0NearUnity(x, false);
    switch (distributionType) {
    case NORMAL:
        return cdfNormal(x);
    case CAUCHY:
        return cdfCauchy(x);
    case LEVY:
        return (beta > 0) ? cdfLevy(x) : cdfLevyCompl(2 * mu - x);
    case UNITY_EXPONENT:
        return cdfForUnityExponent(x, false);
    case GENERAL:
        return cdfForGeneralExponent(x);
    default:
        throw std::runtime_error("Stable distribution: invalid distribution type");
    }
}
 
template <typename RealType>
double StableDistribution<RealType>::S(const RealType &x) const
{
    if (useS0UnityLimit())
        return cdfS0UnityLimit(x, true);
    if (useS0NearUnityKernel(x))
        return cdfS0NearUnity(x, true);
    switch (distributionType) {
    case NORMAL:
        return cdfNormalCompl(x);
    case CAUCHY:
        return cdfCauchyCompl(x);
    case LEVY:
        return (beta > 0) ? cdfLevyCompl(x) : cdfLevy(2 * mu - x);
    case UNITY_EXPONENT:
        return cdfForUnityExponent(x, true);
    case GENERAL:
        return ccdfForGeneralExponent(x);
    default:
        throw std::runtime_error("Stable distribution: invalid distribution type");
    }
}
 
template <typename RealType>
template <std::uniform_random_bit_generator G>
double StableDistribution<RealType>::variateForUnityExponent(G &g) const
{
    RealType U = randlib::consts::pi * detail::uniform01<RealType>(g) - randlib::consts::half_pi;
    RealType W = detail::standardExponential<RealType>(g);
    RealType pi_2pBetaU = randlib::consts::half_pi + beta * U;
    RealType Y = W * std::cos(U) / pi_2pBetaU;
    RealType X = std::log(Y);
    X += logGammaPi_2;
    X *= -beta;
    X += pi_2pBetaU * std::tan(U);
    X *= randlib::consts::two_inv_pi;
    return mu + gamma * X;
}
 
template <typename RealType>
template <std::uniform_random_bit_generator G>
double StableDistribution<RealType>::variateForS0UnityLimit(G &g) const
{
    const RealType u = randlib::consts::pi * detail::uniform01<RealType>(g) - randlib::consts::half_pi;
    const RealType w = detail::standardExponential<RealType>(g);
    const RealType angle = randlib::consts::half_pi + beta * u;
    const RealType logRatio = std::log(w * std::cos(u) / angle);
    const RealType standardized =
        randlib::consts::two_inv_pi *
        (angle * std::tan(u) - beta * (logRatio + randlib::consts::lnpi - randlib::consts::ln2));
    return locationParameter + gamma * standardized;
}
 
template <typename RealType>
template <std::uniform_random_bit_generator G>
double StableDistribution<RealType>::variateForGeneralExponent(G &g) const
{
    RealType U = randlib::consts::pi * detail::uniform01<RealType>(g) - randlib::consts::half_pi;
    RealType W = detail::standardExponential<RealType>(g);
    RealType alphaUpxi = alpha * (U + xi);
    RealType X = std::sin(alphaUpxi);
    RealType W_adj = W / std::cos(U - alphaUpxi);
    X *= W_adj;
    RealType R = omega - alphaInv * std::log(W_adj * std::cos(U));
    X *= std::exp(R);
    return mu + gamma * X;
}
 
template <typename RealType>
template <std::uniform_random_bit_generator G>
double StableDistribution<RealType>::variateForExponentEqualOneHalf(G &g) const
{
    RealType Z1 = detail::standardNormal<RealType>(g);
    RealType Z2 = detail::standardNormal<RealType>(g);
    RealType temp1 = (1.0 + beta) / Z1, temp2 = (1.0 - beta) / Z2;
    RealType var = temp1 - temp2;
    var *= temp1 + temp2;
    var *= 0.25;
    return mu + gamma * var;
}
 
template <typename RealType>
template <std::uniform_random_bit_generator G>
RealType StableDistribution<RealType>::operator()(G &g) const
{
    if (parameterization == Parameterization::S0 && alpha != 1.0 &&
        std::fabs(alpha - 1.0) < S0_SAMPLING_UNITY_TOLERANCE)
        return variateForS0UnityLimit(g);
    switch (distributionType) {
    case NORMAL:
        return mu + randlib::consts::sqrt2 * gamma * detail::standardNormal<RealType>(g);
    case CAUCHY:
        return mu + gamma * CauchyRand<RealType>::StandardVariate(g);
    case LEVY:
        return mu + math::sign(beta) * gamma * LevyRand<RealType>::StandardVariate(g);
    case UNITY_EXPONENT:
        return variateForUnityExponent(g);
    case GENERAL:
        return (alpha == 0.5) ? variateForExponentEqualOneHalf(g) : variateForGeneralExponent(g);
    default:
        throw std::runtime_error("Stable distribution: invalid distribution type");
    }
}
 
template <typename RealType>
void StableDistribution<RealType>::sample(detail::EngineRef g, std::span<RealType> outputData) const
{
    if (parameterization == Parameterization::S0 && alpha != 1.0 &&
        std::fabs(alpha - 1.0) < S0_SAMPLING_UNITY_TOLERANCE) {
        for (RealType &var : outputData)
            var = variateForS0UnityLimit(g);
        return;
    }
    switch (distributionType) {
    case NORMAL: {
        double stdev = randlib::consts::sqrt2 * gamma;
        for (RealType &var : outputData)
            var = mu + stdev * detail::standardNormal<RealType>(g);
    } break;
    case CAUCHY: {
        for (RealType &var : outputData)
            var = mu + gamma * CauchyRand<RealType>::StandardVariate(g);
    } break;
    case LEVY: {
        if (beta > 0) {
            for (RealType &var : outputData)
                var = mu + gamma * LevyRand<RealType>::StandardVariate(g);
        } else {
            for (RealType &var : outputData)
                var = mu - gamma * LevyRand<RealType>::StandardVariate(g);
        }
    } break;
    case UNITY_EXPONENT: {
        for (RealType &var : outputData)
            var = variateForUnityExponent(g);
    } break;
    case GENERAL: {
        if (alpha == 0.5) {
            for (RealType &var : outputData)
                var = variateForExponentEqualOneHalf(g);
        } else {
            for (RealType &var : outputData)
                var = variateForGeneralExponent(g);
        }
    } break;
    default:
        break;
    }
}
 
template <typename RealType>
long double StableDistribution<RealType>::Mean() const
{
    if (alpha > 1)
        return mu;
    if (beta == 1)
        return INFINITY;
    return (beta == -1) ? -INFINITY : NAN;
}
 
template <typename RealType>
long double StableDistribution<RealType>::Variance() const
{
    return (distributionType == NORMAL) ? 2 * gamma * gamma : INFINITY;
}
 
template <typename RealType>
RealType StableDistribution<RealType>::Mode() const
{
    if (beta == 0 || distributionType == NORMAL)
        return mu;
    if (distributionType == LEVY)
        return mu + beta * gamma / 3.0;
    return ContinuousDistribution<RealType>::Mode();
}
 
template <typename RealType>
RealType StableDistribution<RealType>::Median() const
{
    if (beta == 0 || distributionType == NORMAL)
        return mu;
    return ContinuousDistribution<RealType>::Median();
}
 
template <typename RealType>
long double StableDistribution<RealType>::Skewness() const
{
    return (distributionType == NORMAL) ? 0 : NAN;
}
 
template <typename RealType>
long double StableDistribution<RealType>::ExcessKurtosis() const
{
    return (distributionType == NORMAL) ? 0 : NAN;
}
 
template <typename RealType>
RealType StableDistribution<RealType>::quantileNormal(double p) const
{
    return mu - 2 * gamma * math::erfcinv(2 * p);
}
 
template <typename RealType>
RealType StableDistribution<RealType>::quantileNormal1m(double p) const
{
    return mu + 2 * gamma * math::erfcinv(2 * p);
}
 
template <typename RealType>
RealType StableDistribution<RealType>::quantileCauchy(double p) const
{
    return mu - gamma / std::tan(randlib::consts::pi * p);
}
 
template <typename RealType>
RealType StableDistribution<RealType>::quantileCauchy1m(double p) const
{
    return mu + gamma / std::tan(randlib::consts::pi * p);
}
 
template <typename RealType>
RealType StableDistribution<RealType>::quantileLevy(double p) const
{
    double y = math::erfcinv(p);
    return mu + 0.5 * gamma / (y * y);
}
 
template <typename RealType>
RealType StableDistribution<RealType>::quantileLevy1m(double p) const
{
    double y = math::erfinv(p);
    return mu + 0.5 * gamma / (y * y);
}
 
template <typename RealType>
RealType StableDistribution<RealType>::searchTailQuantile(double p, bool upperTail) const
{
    double tail = std::min(p, 1.0 - p);
    double guess =
        (tail < 0.05) ? stableTailQuantile(tail, upperTail != (p > 0.5), alpha, beta, gamma, GetLocation()) : NAN;
    return this->searchQuantile(p, upperTail, guess, 0.05);
}
 
template <typename RealType>
RealType StableDistribution<RealType>::quantileImpl(double p) const
{
    switch (distributionType) {
    case NORMAL:
        return quantileNormal(p);
    case CAUCHY:
        return quantileCauchy(p);
    case LEVY:
        return (beta > 0) ? quantileLevy(p) : 2 * mu - quantileLevy1m(p);
    default:
        return searchTailQuantile(p, false);
    }
}
 
template <typename RealType>
RealType StableDistribution<RealType>::quantileImpl1m(double p) const
{
    switch (distributionType) {
    case NORMAL:
        return quantileNormal1m(p);
    case CAUCHY:
        return quantileCauchy1m(p);
    case LEVY:
        return (beta > 0) ? quantileLevy1m(p) : 2 * mu - quantileLevy(p);
    default:
        return searchTailQuantile(p, true);
    }
}
 
template <typename RealType>
std::complex<double> StableDistribution<RealType>::cfNormal(double t) const
{
    double gammaT = gamma * t;
    std::complex<double> y(-gammaT * gammaT, mu * t);
    return std::exp(y);
}
 
template <typename RealType>
std::complex<double> StableDistribution<RealType>::cfCauchy(double t) const
{
    std::complex<double> y(-gamma * t, mu * t);
    return std::exp(y);
}
 
template <typename RealType>
std::complex<double> StableDistribution<RealType>::cfLevy(double t) const
{
    std::complex<double> y(0.0, -2 * gamma * t);
    y = -std::sqrt(y);
    y += std::complex<double>(0.0, mu * t);
    return std::exp(y);
}
 
template <typename RealType>
std::complex<double> StableDistribution<RealType>::CFImpl(double t) const
{
    if (parameterization == Parameterization::S0) {
        const double scaledT = gamma * t;
        const double attenuation = std::pow(scaledT, alpha);
        double phase = locationParameter * t;
        if (beta != 0.0) {
            if (alpha == 1.0)
                phase -= beta * randlib::consts::two_inv_pi * scaledT * std::log(scaledT);
            else if (alpha != 2.0)
                phase -= beta * scaledT * std::expm1((alpha - 1.0) * std::log(scaledT)) /
                         std::tan(randlib::consts::half_pi * (alpha - 1.0));
        }
        return std::exp(std::complex<double>(-attenuation, phase));
    }
    double x = 0;
    switch (distributionType) {
    case NORMAL:
        return cfNormal(t);
    case CAUCHY:
        return cfCauchy(t);
    case LEVY: {
        std::complex<double> phi = cfLevy(t);
        return (beta > 0) ? phi : std::conj(phi);
    }
    case UNITY_EXPONENT:
        x = beta * randlib::consts::two_inv_pi * std::log(t);
        break;
    default:
        x = zeta;
    }
    double re = std::pow(gamma * t, alpha);
    std::complex<double> psi = std::complex<double>(re, re * x - mu * t);
    return std::exp(-psi);
}
template <typename RealType>
std::string StableRand<RealType>::Name() const
{
    std::string name = "Stable(" + this->toStringWithPrecision(this->GetExponent()) + ", " +
                       this->toStringWithPrecision(this->GetSkewness()) + ", " +
                       this->toStringWithPrecision(this->GetScale()) + ", " +
                       this->toStringWithPrecision(this->GetLocation());
    if (this->GetParameterization() == Parameterization::S0)
        name += ", S0";
    return name + ")";
}
 
template <typename RealType>
void StableRand<RealType>::SetExponent(double exponent)
{
    this->SetParameters(exponent, this->GetSkewness(), this->GetScale(), this->GetLocation());
}
 
template <typename RealType>
void StableRand<RealType>::SetSkewness(double skewness)
{
    this->SetParameters(this->GetExponent(), skewness, this->GetScale(), this->GetLocation());
}
template <typename RealType>
std::string HoltsmarkRand<RealType>::Name() const
{
    return "Holtsmark(" + this->toStringWithPrecision(this->GetScale()) + ", " +
           this->toStringWithPrecision(this->GetLocation()) + ")";
}
template <typename RealType>
std::string LandauRand<RealType>::Name() const
{
    return "Landau(" + this->toStringWithPrecision(this->GetScale()) + ", " +
           this->toStringWithPrecision(this->GetLocation()) + ")";
}
 
} // namespace randlib