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

// gamma.hpp
// Template wrappers for the gamma-function family
//
// Functions provided:
//   gamma(x)      — gamma function Γ(x)
//   lnGamma(x)    — log-gamma function ln Γ(x)
//   beta(a, b)    — beta function B(a,b) = Γ(a)·Γ(b)/Γ(a+b)
//   digamma(x)    — digamma function ψ(x) = Γ'(x)/Γ(x)
//   trigamma(x)   — trigamma function ψ₁(x)
//
// Supported types:
//   float, double, long double — delegate to std::
//   Float                      — delegate to the implementation in FloatMath.cpp

#ifndef SANGI_SPECIAL_GAMMA_HPP
#define SANGI_SPECIAL_GAMMA_HPP

#include <math/special/special_common.hpp>
#include <math/core/mp/Float.hpp>
#include <cmath>
#include <numbers>
#include <vector>

namespace sangi {
namespace special {

// ================================================================
// gamma function Γ(x)
// ================================================================

/// native floating-point types: delegate to std::tgamma
template<IsNativeFloat T>
[[nodiscard]] T gamma(T x) {
    return std::tgamma(x);
}

// Float type: please use sangi::gamma(x, precision) directly.

// ================================================================
// log-gamma function ln Γ(x)
// ================================================================

/// native floating-point types: delegate to std::lgamma
template<IsNativeFloat T>
[[nodiscard]] T lnGamma(T x) {
    return std::lgamma(x);
}

// Float type: please use sangi::lnGamma(x, precision) directly.

// ================================================================
// beta function B(a,b) = Γ(a)·Γ(b)/Γ(a+b)
// ================================================================

/// native floating-point types: delegate to C++17 std::beta
template<IsNativeFloat T>
[[nodiscard]] T beta(T a, T b) {
#if defined(__cpp_lib_math_special_functions) || defined(_MSC_VER)
    return std::beta(a, b);
#else
    return std::exp(std::lgamma(a) + std::lgamma(b) - std::lgamma(a + b));
#endif
}

// Float type: please use sangi::beta(a, b, precision) directly.

// ================================================================
// digamma function ψ(x) = d/dx ln Γ(x)
// ================================================================

/// native floating-point types: implemented via asymptotic expansion
/// the C++ standard library has no digamma, so this is a custom implementation
template<IsNativeFloat T>
[[nodiscard]] T digamma(T x) {
    // negative argument: reflection formula ψ(1-x) - ψ(x) = π·cot(πx)
    if (x < T(0)) {
        if (x == std::floor(x)) {
            return std::numeric_limits<T>::quiet_NaN();  // pole
        }
        const T pi = static_cast<T>(3.14159265358979323846L);
        return digamma(T(1) - x) - pi / std::tan(pi * x);
    }

    // x == 0: pole
    if (x == T(0)) {
        return -std::numeric_limits<T>::infinity();
    }

    T result = T(0);

    // argument shift: ψ(x+1) = ψ(x) + 1/x → raise x up to >= 8
    T z = x;
    while (z < T(8)) {
        result -= T(1) / z;
        z += T(1);
    }

    // asymptotic expansion: ψ(z) ≈ ln(z) - 1/(2z) - Σ B_{2k}/(2k·z^{2k})
    // coefficients are 1/(2k) times the Bernoulli numbers B_{2k} (k=1..7)
    // reference: Abramowitz & Stegun 6.3.18
    T z2 = T(1) / (z * z);
    // 1/(2k) times the Bernoulli numbers B_2, B_4, ..., B_14
    // B_2=1/6, B_4=-1/30, B_6=1/42, B_8=-1/30, B_10=5/66, B_12=-691/2730, B_14=7/6
    constexpr T coeffs[] = {
        static_cast<T>(1.0L / 12.0L),         // B_2/2  = (1/6)/2  = 1/12
        static_cast<T>(-1.0L / 120.0L),        // B_4/4  = (-1/30)/4 = -1/120
        static_cast<T>(1.0L / 252.0L),         // B_6/6  = (1/42)/6  = 1/252
        static_cast<T>(-1.0L / 240.0L),        // B_8/8  = (-1/30)/8 = -1/240
        static_cast<T>(1.0L / 132.0L),         // B_10/10 = (5/66)/10 = 1/132
        static_cast<T>(-691.0L / 32760.0L),    // B_12/12 = (-691/2730)/12 = -691/32760
        static_cast<T>(1.0L / 12.0L),          // B_14/14 = (7/6)/14 = 1/12  (same value as B_2/2 but a different coefficient)
    };

    T series = coeffs[6];
    for (int i = 5; i >= 0; --i) {
        series = series * z2 + coeffs[i];
    }

    result += std::log(z) - T(0.5) / z - series * z2;
    return result;
}

// Float type: please use sangi::digamma(x, precision) directly.

// ================================================================
// trigamma function ψ₁(x) = d²/dx² ln Γ(x)
// ================================================================

/// native floating-point types: implemented via asymptotic expansion
template<IsNativeFloat T>
[[nodiscard]] T trigamma(T x) {
    // negative argument: reflection formula ψ₁(1-x) + ψ₁(x) = π²/sin²(πx)
    if (x < T(0)) {
        if (x == std::floor(x)) {
            return std::numeric_limits<T>::quiet_NaN();  // pole
        }
        const T pi = static_cast<T>(3.14159265358979323846L);
        T sinpx = std::sin(pi * x);
        return (pi * pi) / (sinpx * sinpx) - trigamma(T(1) - x);
    }

    // x == 0: pole
    if (x == T(0)) {
        return std::numeric_limits<T>::infinity();
    }

    T result = T(0);

    // argument shift: ψ₁(x) = ψ₁(x+1) + 1/x²
    T z = x;
    while (z < T(8)) {
        result += T(1) / (z * z);
        z += T(1);
    }

    // asymptotic expansion: ψ₁(z) ≈ 1/z + 1/(2z²) + Σ B_{2k}/z^{2k+1}
    T z_inv = T(1) / z;
    T z2 = z_inv * z_inv;
    // B_2, B_4, B_6, B_8, B_10, B_12
    constexpr T bernoulli[] = {
        static_cast<T>(1.0L / 6.0L),
        static_cast<T>(-1.0L / 30.0L),
        static_cast<T>(1.0L / 42.0L),
        static_cast<T>(-1.0L / 30.0L),
        static_cast<T>(5.0L / 66.0L),
        static_cast<T>(-691.0L / 2730.0L),
    };

    T series = bernoulli[5];
    for (int i = 4; i >= 0; --i) {
        series = series * z2 + bernoulli[i];
    }

    result += z_inv + T(0.5) * z2 + series * z2 * z_inv;
    return result;
}

// Float type: please use sangi::trigamma(x, precision) directly.

// ================================================================
// polygamma function ψ^(n)(x) = d^{n+1}/dx^{n+1} ln Γ(x)
// ================================================================

/// native floating-point types: implemented via asymptotic expansion
/// n=0 → delegate to digamma, n=1 → trigamma. n≥2 uses a generic asymptotic expansion.
template<IsNativeFloat T>
[[nodiscard]] T polygamma(int n, T x) {
    if (n < 0) return std::numeric_limits<T>::quiet_NaN();
    if (n == 0) return digamma(x);
    if (n == 1) return trigamma(x);

    // special values
    if (std::isnan(x)) return std::numeric_limits<T>::quiet_NaN();
    if (x == T(0)) return std::numeric_limits<T>::quiet_NaN();
    if (x < T(0) && x == std::floor(x)) return std::numeric_limits<T>::quiet_NaN();
    if (std::isinf(x)) {
        return x > T(0) ? T(0) : std::numeric_limits<T>::quiet_NaN();
    }

    // compute n!
    T n_fact = T(1);
    for (int i = 2; i <= n; i++) n_fact *= static_cast<T>(i);

    // argument shift: shift until z ≥ 10
    // ψ^(n)(z) = ψ^(n)(z+1) + (-1)^{n+1} · n! / z^{n+1}
    T z = x;
    T shift_sum = T(0);
    while (z < T(10)) {
        T zi_pow = T(1);
        for (int j = 0; j <= n; j++) zi_pow *= z;
        shift_sum += T(1) / zi_pow;
        z += T(1);
    }

    // asymptotic expansion: ψ^(n)(z) = (-1)^{n+1} · A
    // A = (n-1)!/z^n + n!/(2z^{n+1})
    //   + Σ_{k=1}^{K} B_{2k} · (2k+n-1)!/((2k)!) · z^{-(2k+n)}
    T n_minus_1_fact = n_fact / static_cast<T>(n);
    T z_inv = T(1) / z;
    T z_inv_n = T(1);
    for (int i = 0; i < n; i++) z_inv_n *= z_inv;

    T A = n_minus_1_fact * z_inv_n
        + n_fact * z_inv_n * z_inv / T(2);

    // B_{2k} (k=1..8): B_2, B_4, ..., B_16
    constexpr long double B[] = {
         1.0L / 6.0L,            // B_2
        -1.0L / 30.0L,           // B_4
         1.0L / 42.0L,           // B_6
        -1.0L / 30.0L,           // B_8
         5.0L / 66.0L,           // B_10
        -691.0L / 2730.0L,       // B_12
         7.0L / 6.0L,            // B_14
        -3617.0L / 510.0L,       // B_16
    };

    T z2 = z_inv * z_inv;
    T z_power = z_inv_n * z2;  // z^{-(n+2)}

    for (int k = 0; k < 8; k++) {
        int kk = k + 1;
        // R(n,kk) = (2kk+n-1)!/(2kk)! = Π_{j=1}^{n-1} (2kk+j)
        T R = T(1);
        for (int j = 1; j <= n - 1; j++) {
            R *= static_cast<T>(2 * kk + j);
        }
        T term = static_cast<T>(B[k]) * R * z_power;
        A += term;
        if (std::abs(term) < std::abs(A) * std::numeric_limits<T>::epsilon()) break;
        z_power *= z2;
    }

    // sign: (-1)^{n+1}
    T asymp = (n % 2 == 0) ? -A : A;

    // shift correction: (-1)^{n+1} · n! · shift_sum
    T shift_correction = n_fact * shift_sum;
    if (n % 2 == 0) shift_correction = -shift_correction;

    return asymp + shift_correction;
}

// Float type: please use sangi::polygamma(n, x, precision) directly.

// ================================================================
// regularized incomplete gamma functions — forward declarations (for mutual recursion)
// ================================================================

template<IsNativeFloat T> [[nodiscard]] T gammaP(T a, T x);
template<IsNativeFloat T> [[nodiscard]] T gammaQ(T a, T x);

// ================================================================
// regularized lower incomplete gamma P(a,x) = γ(a,x)/Γ(a)
// ================================================================

/// native floating-point types
/// x < a+1: Taylor series, x ≥ a+1: 1 - Q(a,x)
template<IsNativeFloat T>
[[nodiscard]] T gammaP(T a, T x) {
    if (std::isnan(a) || std::isnan(x)) return std::numeric_limits<T>::quiet_NaN();
    if (x < T(0)) return std::numeric_limits<T>::quiet_NaN();
    if (a <= T(0)) return std::numeric_limits<T>::quiet_NaN();
    if (x == T(0)) return T(0);
    if (std::isinf(x)) return T(1);

    // prefactor: exp(-x + a·ln(x) - lgamma(a))
    T front = std::exp(-x + a * std::log(x) - std::lgamma(a));

    if (x < a + T(1)) {
        // Taylor series: S = Σ x^n / (a·(a+1)···(a+n))
        T term = T(1) / a;
        T sum = term;
        for (int n = 1; n < 200; n++) {
            term *= x / (a + static_cast<T>(n));
            sum += term;
            if (std::abs(term) < std::abs(sum) * std::numeric_limits<T>::epsilon())
                break;
        }
        return front * sum;
    } else {
        return T(1) - gammaQ(a, x);
    }
}

// Float type: please use sangi::gammaP(a, x, precision) directly.

// ================================================================
// regularized upper incomplete gamma Q(a,x) = Γ(a,x)/Γ(a) = 1 - P(a,x)
// ================================================================

/// native floating-point types
/// x ≥ a+1: Legendre continued fraction, x < a+1: 1 - P(a,x)
template<IsNativeFloat T>
[[nodiscard]] T gammaQ(T a, T x) {
    if (std::isnan(a) || std::isnan(x)) return std::numeric_limits<T>::quiet_NaN();
    if (x < T(0)) return std::numeric_limits<T>::quiet_NaN();
    if (a <= T(0)) return std::numeric_limits<T>::quiet_NaN();
    if (x == T(0)) return T(1);
    if (std::isinf(x)) return T(0);

    T front = std::exp(-x + a * std::log(x) - std::lgamma(a));

    if (x >= a + T(1)) {
        // Legendre CF: Q = front / f
        // b_0 = x+1-a, a_n = -n(n-a), b_n = x+2n+1-a
        constexpr T tiny = std::numeric_limits<T>::min();

        T b0 = x + T(1) - a;
        T f = (std::abs(b0) < tiny) ? tiny : b0;
        T C = f;
        T D = T(0);

        for (int n = 1; n < 200; n++) {
            T an = static_cast<T>(-n) * (static_cast<T>(n) - a);
            T bn = x + static_cast<T>(2 * n + 1) - a;

            D = bn + an * D;
            if (std::abs(D) < tiny) D = tiny;
            D = T(1) / D;

            C = bn + an / C;
            if (std::abs(C) < tiny) C = tiny;

            T delta = C * D;
            f *= delta;

            if (std::abs(delta - T(1)) < std::numeric_limits<T>::epsilon())
                break;
        }
        return front / f;
    } else {
        return T(1) - gammaP(a, x);
    }
}

// Float type: please use sangi::gammaQ(a, x, precision) directly.

// ================================================================
// lower incomplete gamma γ(a,x) = P(a,x)·Γ(a)
// ================================================================

template<IsNativeFloat T>
[[nodiscard]] T gammaLower(T a, T x) {
    return gammaP(a, x) * std::tgamma(a);
}

// Float type: please use sangi::gammaLower(a, x, precision) directly.

// ================================================================
// upper incomplete gamma Γ(a,x) = Q(a,x)·Γ(a)
// ================================================================

template<IsNativeFloat T>
[[nodiscard]] T gammaUpper(T a, T x) {
    return gammaQ(a, x) * std::tgamma(a);
}

// Float type: please use sangi::gammaUpper(a, x, precision) directly.

// ================================================================
// regularized incomplete beta function I_x(a,b)
// ================================================================

/// native floating-point types: continued-fraction expansion + symmetry
template<IsNativeFloat T>
[[nodiscard]] T betaRegularized(T x, T a, T b) {
    if (std::isnan(x) || std::isnan(a) || std::isnan(b))
        return std::numeric_limits<T>::quiet_NaN();
    if (x < T(0) || x > T(1)) return std::numeric_limits<T>::quiet_NaN();
    if (a <= T(0) || b <= T(0)) return std::numeric_limits<T>::quiet_NaN();
    if (x == T(0)) return T(0);
    if (x == T(1)) return T(1);

    // symmetry: invert if x ≥ (a+1)/(a+b+2)
    bool flip = x >= (a + T(1)) / (a + b + T(2));
    T xx = flip ? (T(1) - x) : x;
    T aa = flip ? b : a;
    T bb = flip ? a : b;

    // front = exp(a·ln(x) + b·ln(1-x) + lgamma(a+b) - lgamma(a) - lgamma(b)) / a
    T ln_front = aa * std::log(xx) + bb * std::log(T(1) - xx)
               + std::lgamma(aa + bb) - std::lgamma(aa) - std::lgamma(bb);
    T front = std::exp(ln_front) / aa;

    // CF: 1 + d_1/(1 + d_2/(1 + ...))
    constexpr T tiny = std::numeric_limits<T>::min();

    T f = T(1);
    T C = T(1);
    T D = T(0);

    for (int n = 1; n < 200; n++) {
        T d;
        if (n % 2 == 1) {
            int m = (n - 1) / 2;
            T am = aa + static_cast<T>(m);
            T abm = aa + bb + static_cast<T>(m);
            T denom = (aa + static_cast<T>(2 * m)) * (aa + static_cast<T>(2 * m + 1));
            d = -(am * abm * xx) / denom;
        } else {
            int m = n / 2;
            T fm = static_cast<T>(m);
            T bm = bb - static_cast<T>(m);
            T denom = (aa + static_cast<T>(2 * m - 1)) * (aa + static_cast<T>(2 * m));
            d = (fm * bm * xx) / denom;
        }

        D = T(1) + d * D;
        if (std::abs(D) < tiny) D = tiny;
        D = T(1) / D;

        C = T(1) + d / C;
        if (std::abs(C) < tiny) C = tiny;

        T delta = C * D;
        f *= delta;

        if (std::abs(delta - T(1)) < std::numeric_limits<T>::epsilon())
            break;
    }

    T result = front / f;
    return flip ? (T(1) - result) : result;
}

// Float type: please use sangi::betaRegularized(x, a, b, precision) directly.

// ================================================================
// Complex<R> overloads
// ================================================================

// ----------------------------------------------------------------
// Γ(z) — Lanczos approximation (g=7, n=9)
// ----------------------------------------------------------------
// Re(z) ≥ 0.5: direct Lanczos
// Re(z) < 0.5: reflection formula Γ(z) = π / (sin(πz)·Γ(1-z))

template<IsNativeFloat R>
[[nodiscard]] Complex<R> gamma(Complex<R> z) {
    using C = Complex<R>;

    // Lanczos coefficients (g=7, n=9)
    static constexpr R lanczos_g = R(7);
    static constexpr R lanczos_c[] = {
        R( 0.99999999999980993),
        R( 676.5203681218851),
        R(-1259.1392167224028),
        R( 771.32342877765313),
        R(-176.61502916214059),
        R( 12.507343278686905),
        R(-0.13857109526572012),
        R( 9.9843695780195716e-6),
        R( 1.5056327351493116e-7),
    };

    // reflection formula: Re(z) < 0.5
    if (z.re < R(0.5)) {
        C pi(std::numbers::pi_v<R>);
        return pi / (sin(pi * z) * gamma(C(R(1)) - z));
    }

    // Lanczos: Γ(z) = √(2π)·t^{z-0.5}·e^{-t}·A(z),  t = z-1+g+0.5
    C zm1 = z - C(R(1));
    C x(lanczos_c[0]);
    for (int k = 1; k <= 8; k++) {
        x = x + C(lanczos_c[k]) / (zm1 + C(R(k)));
    }

    C t = zm1 + C(lanczos_g + R(0.5));
    constexpr R sqrt_2pi = R(2.5066282746310005024157652848110452L);
    return C(sqrt_2pi) * pow(t, zm1 + C(R(0.5))) * exp(-t) * x;
}

// ----------------------------------------------------------------
// lnΓ(z) — Lanczos + log
// ----------------------------------------------------------------

template<IsNativeFloat R>
[[nodiscard]] Complex<R> lnGamma(Complex<R> z) {
    using C = Complex<R>;

    static constexpr R lanczos_g = R(7);
    static constexpr R lanczos_c[] = {
        R( 0.99999999999980993),
        R( 676.5203681218851),
        R(-1259.1392167224028),
        R( 771.32342877765313),
        R(-176.61502916214059),
        R( 12.507343278686905),
        R(-0.13857109526572012),
        R( 9.9843695780195716e-6),
        R( 1.5056327351493116e-7),
    };

    // reflection: Re(z) < 0.5 → lnΓ(z) = ln(π) - ln(sin(πz)) - lnΓ(1-z)
    if (z.re < R(0.5)) {
        C pi(std::numbers::pi_v<R>);
        return C(std::log(std::numbers::pi_v<R>))
             - log(sin(pi * z))
             - lnGamma(C(R(1)) - z);
    }

    C zm1 = z - C(R(1));
    C x(lanczos_c[0]);
    for (int k = 1; k <= 8; k++) {
        x = x + C(lanczos_c[k]) / (zm1 + C(R(k)));
    }

    C t = zm1 + C(lanczos_g + R(0.5));
    constexpr R half_ln_2pi = R(0.91893853320467274178032973640561764L);
    return C(half_ln_2pi) + (zm1 + C(R(0.5))) * log(t) - t + log(x);
}

// ----------------------------------------------------------------
// B(a,b) = exp(lnΓ(a) + lnΓ(b) - lnΓ(a+b))
// ----------------------------------------------------------------

template<IsNativeFloat R>
[[nodiscard]] Complex<R> beta(Complex<R> a, Complex<R> b) {
    return exp(lnGamma(a) + lnGamma(b) - lnGamma(a + b));
}

// ----------------------------------------------------------------
// ψ(z) — asymptotic expansion + argument shift
// ----------------------------------------------------------------

template<IsNativeFloat R>
[[nodiscard]] Complex<R> digamma(Complex<R> z) {
    using C = Complex<R>;

    // reflection formula: Re(z) < 0.5 → ψ(z) = ψ(1-z) - π·cot(πz)
    if (z.re < R(0.5)) {
        C pi(std::numbers::pi_v<R>);
        return digamma(C(R(1)) - z) - pi * (cos(pi * z) / sin(pi * z));
    }

    C result(R(0));

    // argument shift: ψ(z+1) = ψ(z) + 1/z → until Re(z) ≥ 8
    C w = z;
    while (w.re < R(8)) {
        result = result - C(R(1)) / w;
        w = w + C(R(1));
    }

    // asymptotic expansion: ψ(z) ≈ ln(z) - 1/(2z) - Σ B_{2k}/(2k·z^{2k})
    C w2 = C(R(1)) / (w * w);
    constexpr R coeffs[] = {
        R( 1.0L / 12.0L),         // B_2 / 2
        R(-1.0L / 120.0L),        // B_4 / 4
        R( 1.0L / 252.0L),        // B_6 / 6
        R(-1.0L / 240.0L),        // B_8 / 8
        R( 1.0L / 132.0L),        // B_10 / 10
        R(-691.0L / 32760.0L),    // B_12 / 12
        R( 1.0L / 12.0L),         // B_14 / 14
    };

    C series(coeffs[6]);
    for (int i = 5; i >= 0; --i) {
        series = series * w2 + C(coeffs[i]);
    }

    result = result + log(w) - C(R(0.5)) / w - series * w2;
    return result;
}

// ----------------------------------------------------------------
// ψ^(n)(z) — polygamma (Complex<double>) n≥0
// ----------------------------------------------------------------
// asymptotic expansion + argument shift. n=0 delegates to digamma.

template<IsNativeFloat R>
[[nodiscard]] Complex<R> polygamma(int n, Complex<R> z) {
    using C = Complex<R>;
    if (n < 0) return C(std::numeric_limits<R>::quiet_NaN());
    if (n == 0) return digamma(z);

    // compute n!
    R n_fact = R(1);
    for (int i = 2; i <= n; i++) n_fact *= static_cast<R>(i);

    // argument shift: until Re(z) ≥ 10
    // ψ^(n)(z) = ψ^(n)(z+1) + (-1)^{n+1} · n! / z^{n+1}
    C w = z;
    C shift_sum(R(0));
    while (w.re < R(10)) {
        C wpow = C(R(1));
        for (int j = 0; j <= n; j++) wpow = wpow * w;
        shift_sum = shift_sum + C(R(1)) / wpow;
        w = w + C(R(1));
    }

    // asymptotic expansion: ψ^(n)(z) = (-1)^{n+1} · A
    // A = (n-1)!/z^n + n!/(2z^{n+1})
    //   + Σ_{k=1}^{K} B_{2k} · R(n,k) · z^{-(2k+n)}
    R n_minus_1_fact = n_fact / static_cast<R>(n);
    C w_inv = C(R(1)) / w;
    C w_inv_n = C(R(1));
    for (int i = 0; i < n; i++) w_inv_n = w_inv_n * w_inv;

    C A = C(n_minus_1_fact) * w_inv_n
        + C(n_fact) * w_inv_n * w_inv / C(R(2));

    constexpr long double B[] = {
         1.0L / 6.0L,
        -1.0L / 30.0L,
         1.0L / 42.0L,
        -1.0L / 30.0L,
         5.0L / 66.0L,
        -691.0L / 2730.0L,
         7.0L / 6.0L,
        -3617.0L / 510.0L,
    };

    C w2 = w_inv * w_inv;
    C w_power = w_inv_n * w2;  // z^{-(n+2)}

    for (int k = 0; k < 8; k++) {
        int kk = k + 1;
        R Rval = R(1);
        for (int j = 1; j <= n - 1; j++) {
            Rval *= static_cast<R>(2 * kk + j);
        }
        C term = C(static_cast<R>(B[k]) * Rval) * w_power;
        A = A + term;
        if (sangi::abs(term) < sangi::abs(A) * std::numeric_limits<R>::epsilon()) break;
        w_power = w_power * w2;
    }

    // sign: (-1)^{n+1}
    C asymp = (n % 2 == 0) ? C(R(0)) - A : A;

    // shift correction: (-1)^{n+1} · n! · shift_sum
    C shift_correction = C(n_fact) * shift_sum;
    if (n % 2 == 0) shift_correction = C(R(0)) - shift_correction;

    return asymp + shift_correction;
}

// ================================================================
// Complex<Float>-specific overloads (Stirling asymptotic expansion)
// ================================================================
// Lanczos (g=7, 9 coefficients, ~15-digit precision) is kept for IsNativeFloat only.
// Complex<Float> uses an arbitrary-precision Stirling series.

namespace detail {

// detail::PrecisionGuard was removed from all special-function headers when FLOAT_PRECISION_PLAN Step 12 completed (commit 6fa993f
// + Step 12b-j in bulk). From Phase 2 onward, the combination of
// the 1-arg ADL overload (via requestedPrecision(x)) + setResultPrecision(wp)
// ensures req=wp propagation, so there is no longer any need to temporarily
// rewrite Float::defaultPrecision.

/// Compute Bernoulli numbers B_{2k} in Float (Akiyama-Tanigawa)
/// result[k] = B_{2k}  (k = 0, 1, ..., max_k)
inline std::vector<Float> computeBernoulliFloat(int max_k, int precision) {
    int n = 2 * max_k + 1;
    std::vector<Float> a(n + 1, Float(0));
    std::vector<Float> result(max_k + 1, Float(0));
    result[0] = Float(1);

    // Mark the value 1 with a "precision-digit" request. This avoids the defaultPrecision
    // fallback when dividing by an exact (Float(m+1)), ensuring computation at the requested precision.
    // (this helper must work independently of PrecisionGuard)
    Float one_p(1);
    one_p.setResultPrecision(precision);

    for (int m = 0; m <= n; m++) {
        a[m] = one_p / Float(m + 1);
        a[m].setPrecision(precision);
        for (int j = m; j >= 1; j--) {
            a[j - 1] = Float(j) * (a[j - 1] - a[j]);
            a[j - 1].setPrecision(precision);
        }
        if (m >= 2 && m % 2 == 0)
            result[m / 2] = a[0];
    }
    return result;
}

} // namespace detail

// ----------------------------------------------------------------
// lnΓ(z) — Stirling asymptotic expansion (Complex<Float>)
// ----------------------------------------------------------------
// lnΓ(z) = (z-1/2)·ln(z) - z + ln(2π)/2 + Σ B_{2k}/(2k·(2k-1)·z^{2k-1})
// Re(z) < 0.5: reflection formula lnΓ(z) = ln(π) - ln(sin(πz)) - lnΓ(1-z)

[[nodiscard]] inline Complex<Float> lnGamma(Complex<Float> z, int precision) {
    using C = Complex<Float>;
    int wp = precision + 20;
    Float::PrecisionScope _ps(wp);   // compute exact÷exact inside the function at wp digits
    // PrecisionGuard removed: w.re.setResultPrecision(wp) propagates z's req=wp within the function.
    // Even in the reflection branch, Float::pi(wp) carries req=wp, and Complex arithmetic keeps >= wp.
    // computeBernoulliFloat is not called (used only by digamma/polygamma).

    // reflection formula: Re(z) < 0.5
    if (z.re.toDouble() < 0.5) {
        Float pi_val = Float::pi(wp);
        C ln_pi(log(pi_val, wp));
        C sin_pi_z = sin(C(pi_val) * z);
        return ln_pi - log(sin_pi_z) - lnGamma(C(Float(1)) - z, precision);
    }

    // argument shift: until Re(w) is large enough
    double shift_target = wp * 0.45 + 10;  // BUGFIX 2026-05-30: 0.35 is too small for asymptotic convergence
    C w = z;
    w.re.setResultPrecision(wp);  // pad the input to working precision (also sets eff)
    w.im.setResultPrecision(wp);

    C prod_log(Float(0));
    double w_re = w.re.toDouble();
    int m = 0;
    if (w_re < shift_target) {
        m = static_cast<int>(shift_target - w_re) + 1;
        for (int i = 0; i < m; i++) {
            C val = w + C(Float(i));
            prod_log = prod_log + log(val);
        }
        w = w + C(Float(m));
    }

    // number of Bernoulli terms: optimal truncation k* ≈ π·|w| (same as the real lnGamma/digamma).
    // BUGFIX (2026-05-30): num_terms = wp/(2·log2|w|) misused a bit-formula assuming a geometric series,
    // giving too few terms → capped at ~0.3·P digits (AUDIT_FLOAT_UNIT_MIXING pattern B).
    // the convergence test abs(term) < eps·abs(result) uses relative eps, so it is correct and fires before truncation.
    double w_mag = std::sqrt(w.re.toDouble() * w.re.toDouble()
                           + w.im.toDouble() * w.im.toDouble());
    constexpr double PI_VAL = 3.141592653589793238462643383279502884;
    int num_terms = static_cast<int>(PI_VAL * w_mag) + 10;
    if (num_terms < 5) num_terms = 5;
    int num_terms_cap = static_cast<int>(PI_VAL * shift_target * 1.2 + 100.0);
    if (num_terms_cap < 2000) num_terms_cap = 2000;
    if (num_terms > num_terms_cap) num_terms = num_terms_cap;

    auto bern = detail::computeBernoulliFloat(num_terms, wp);

    // main computation
    C ln_w = log(w);
    Float half = Float(1) / Float(2);

    // (w - 1/2)·ln(w) - w
    C result = (w - C(half)) * ln_w - w;

    // + ln(2π)/2
    Float ln2pi = log(Float(2) * Float::pi(wp), wp);
    result = result + C(ln2pi * half);

    // Stirling correction: Σ B_{2k} / (2k·(2k-1)·z^{2k-1})
    C w_inv = C(Float(1)) / w;
    C w_inv2 = w_inv * w_inv;
    C w_power = w_inv;

    Float eps = Float::epsilon(wp);

    // BUGFIX (2026-05-30): when num_terms is raised to the optimal point π·|w|, the asymptotic series diverges
    // past its smallest term, so divergence detection is required (without it, it explodes to ~10^116). Same as the real version.
    Float prev_abs = Float::positiveInfinity();
    for (int k = 1; k <= num_terms; k++) {
        Float coeff = bern[k] / Float(2 * k * (2 * k - 1));
        C term = C(coeff) * w_power;
        Float at = abs(term);
        if (k >= 3 && at > prev_abs) break;  // divergence detection: passed the smallest term
        result = result + term;
        if (k >= 3 && at < eps * abs(result)) break;
        prev_abs = at;
        w_power = w_power * w_inv2;
    }

    // shift correction
    result = result - prod_log;

    result.re.setPrecision(precision);
    result.im.setPrecision(precision);
    return result;
}

[[nodiscard]] inline Complex<Float> lnGamma(Complex<Float> z) {
    return lnGamma(std::move(z), Float::defaultPrecision());
}

// ----------------------------------------------------------------
// Γ(z) = exp(lnΓ(z)) (Complex<Float>)
// ----------------------------------------------------------------

[[nodiscard]] inline Complex<Float> gamma(Complex<Float> z, int precision) {
    using C = Complex<Float>;
    int wp = precision + 15;
    Float::PrecisionScope _ps(wp);   // compute exact÷exact inside the function at wp digits
    // PrecisionGuard removed: lnGamma(z, wp) gives a result with req=wp at explicit precision,
    // and exp(lng) is Complex<Float> exp → uses requestedPrecision = wp via the 1-arg ADL.

    C lng = lnGamma(z, wp);
    C result = exp(lng);

    result.re.setPrecision(precision);
    result.im.setPrecision(precision);
    return result;
}

[[nodiscard]] inline Complex<Float> gamma(Complex<Float> z) {
    return gamma(std::move(z), Float::defaultPrecision());
}

// ----------------------------------------------------------------
// ψ(z) — asymptotic expansion + argument shift (Complex<Float>)
// ----------------------------------------------------------------
// ψ(z) ≈ ln(z) - 1/(2z) - Σ B_{2k}/(2k·z^{2k})
// Re(z) < 0.5: reflection ψ(z) = ψ(1-z) - π·cot(πz)

[[nodiscard]] inline Complex<Float> digamma(Complex<Float> z, int precision) {
    using C = Complex<Float>;
    int wp = precision + 20;
    Float::PrecisionScope _ps(wp);   // compute exact÷exact inside the function at wp digits
    // PrecisionGuard removed: w.re.setResultPrecision(wp) propagates z's req=wp,
    // and computeBernoulliFloat is also independent of defaultPrecision via its precision argument.

    // reflection formula
    if (z.re.toDouble() < 0.5) {
        Float pi_val = Float::pi(wp);
        C pi_c(pi_val);
        C pi_z = pi_c * z;
        return digamma(C(Float(1)) - z, precision) - pi_c * (cos(pi_z) / sin(pi_z));
    }

    C result(Float(0));
    C w = z;
    w.re.setResultPrecision(wp);  // pad the input to working precision (also sets eff)
    w.im.setResultPrecision(wp);

    // argument shift: ψ(z+1) = ψ(z) + 1/z
    double shift_target = wp * 0.45 + 10;  // BUGFIX 2026-05-30: 0.35 is too small for asymptotic convergence
    double w_re = w.re.toDouble();
    if (w_re < shift_target) {
        int m_shift = static_cast<int>(shift_target - w_re) + 1;
        for (int i = 0; i < m_shift; i++) {
            C val = w + C(Float(i));
            result = result - C(Float(1)) / val;
        }
        w = w + C(Float(m_shift));
    }

    // number of Bernoulli terms: optimal truncation k* ≈ π·|w| (same as the real lnGamma/digamma).
    // BUGFIX (2026-05-30): num_terms = wp/(2·log2|w|) misused a bit-formula assuming a geometric series,
    // giving too few terms → capped at ~0.3·P digits (AUDIT_FLOAT_UNIT_MIXING pattern B).
    // the convergence test abs(term) < eps·abs(result) uses relative eps, so it is correct and fires before truncation.
    double w_mag = std::sqrt(w.re.toDouble() * w.re.toDouble()
                           + w.im.toDouble() * w.im.toDouble());
    constexpr double PI_VAL = 3.141592653589793238462643383279502884;
    int num_terms = static_cast<int>(PI_VAL * w_mag) + 10;
    if (num_terms < 5) num_terms = 5;
    int num_terms_cap = static_cast<int>(PI_VAL * shift_target * 1.2 + 100.0);
    if (num_terms_cap < 2000) num_terms_cap = 2000;
    if (num_terms > num_terms_cap) num_terms = num_terms_cap;

    auto bern = detail::computeBernoulliFloat(num_terms, wp);

    // asymptotic expansion: ψ(w) ≈ ln(w) - 1/(2w) - Σ B_{2k}/(2k·w^{2k})
    result = result + log(w) - C(Float(1) / Float(2)) / w;

    C w_inv2 = C(Float(1)) / (w * w);
    C w_power = w_inv2;

    Float eps = Float::epsilon(wp);

    // BUGFIX (2026-05-30): divergence detection (truncate when the smallest term is passed). Required for num_terms=π·|w|.
    Float prev_abs = Float::positiveInfinity();
    for (int k = 1; k <= num_terms; k++) {
        Float coeff = bern[k] / Float(2 * k);
        C term = C(coeff) * w_power;
        Float at = abs(term);
        if (k >= 3 && at > prev_abs) break;  // divergence detection
        result = result - term;
        if (k >= 3 && at < eps * abs(result)) break;
        prev_abs = at;
        w_power = w_power * w_inv2;
    }

    result.re.setPrecision(precision);
    result.im.setPrecision(precision);
    return result;
}

[[nodiscard]] inline Complex<Float> digamma(Complex<Float> z) {
    return digamma(std::move(z), Float::defaultPrecision());
}

// ----------------------------------------------------------------
// ψ^(n)(z) — polygamma (Complex<Float>) n≥0
// ----------------------------------------------------------------

[[nodiscard]] inline Complex<Float> polygamma(int n, Complex<Float> z, int precision) {
    using C = Complex<Float>;
    if (n < 0) return C(Float::nan());
    if (n == 0) return digamma(std::move(z), precision);

    int wp = precision + 20;
    Float::PrecisionScope _ps(wp);   // compute exact÷exact inside the function at wp digits
    // PrecisionGuard removed: w's setResultPrecision(wp) and computeBernoulliFloat's
    // precision-independence make manipulating defaultPrecision unnecessary.

    C w = z;
    w.re.setResultPrecision(wp);
    w.im.setResultPrecision(wp);

    // compute n!
    Float n_fact(1);
    for (int i = 2; i <= n; i++) n_fact = n_fact * Float(i);

    // argument shift: until Re(w) is large enough
    // ψ^(n)(z) = ψ^(n)(z+1) + (-1)^{n+1} · n! / z^{n+1}
    double shift_target = wp * 0.45 + 10;  // BUGFIX 2026-05-30: 0.35 is too small for asymptotic convergence
    C shift_sum(Float(0));
    double w_re = w.re.toDouble();
    if (w_re < shift_target) {
        int m_shift = static_cast<int>(shift_target - w_re) + 1;
        for (int i = 0; i < m_shift; i++) {
            C val = w + C(Float(i));
            C vpow = C(Float(1));
            for (int j = 0; j <= n; j++) vpow = vpow * val;
            shift_sum = shift_sum + C(Float(1)) / vpow;
        }
        w = w + C(Float(m_shift));
    }

    // number of Bernoulli terms: optimal truncation k* ≈ π·|w| (same as the real lnGamma/digamma).
    // BUGFIX (2026-05-30): num_terms = wp/(2·log2|w|) misused a bit-formula assuming a geometric series,
    // giving too few terms → capped at ~0.3·P digits (AUDIT_FLOAT_UNIT_MIXING pattern B).
    // the convergence test abs(term) < eps·abs(result) uses relative eps, so it is correct and fires before truncation.
    double w_mag = std::sqrt(w.re.toDouble() * w.re.toDouble()
                           + w.im.toDouble() * w.im.toDouble());
    constexpr double PI_VAL = 3.141592653589793238462643383279502884;
    int num_terms = static_cast<int>(PI_VAL * w_mag) + 10;
    if (num_terms < 5) num_terms = 5;
    int num_terms_cap = static_cast<int>(PI_VAL * shift_target * 1.2 + 100.0);
    if (num_terms_cap < 2000) num_terms_cap = 2000;
    if (num_terms > num_terms_cap) num_terms = num_terms_cap;

    auto bern = detail::computeBernoulliFloat(num_terms, wp);

    // asymptotic expansion: ψ^(n)(z) = (-1)^{n+1} · A
    // A = (n-1)!/z^n + n!/(2z^{n+1})
    //   + Σ_{k=1}^{K} B_{2k} · R(n,k) · z^{-(2k+n)}
    Float n_minus_1_fact = n_fact / Float(n);
    C w_inv = C(Float(1)) / w;
    C w_inv_n = C(Float(1));
    for (int i = 0; i < n; i++) w_inv_n = w_inv_n * w_inv;

    C A = C(n_minus_1_fact) * w_inv_n
        + C(n_fact) * w_inv_n * w_inv / C(Float(2));

    C w2 = w_inv * w_inv;
    C w_power = w_inv_n * w2;  // z^{-(n+2)}
    Float eps = Float::epsilon(wp);

    // BUGFIX (2026-05-30): divergence detection is required. The polygamma asymptotic series, due to polynomial growth of R(n,k),
    // diverges sharply past the smallest term; summing up to num_terms=π·|w| exploded to ~10^116.
    Float prev_abs = Float::positiveInfinity();
    for (int k = 1; k <= num_terms; k++) {
        // R(n,k) = Π_{j=1}^{n-1} (2k+j)
        Float Rval(1);
        for (int j = 1; j <= n - 1; j++) {
            Rval = Rval * Float(2 * k + j);
        }
        C term = C(bern[k] * Rval) * w_power;
        Float at = abs(term);
        if (k >= 3 && at > prev_abs) break;  // divergence detection: passed the smallest term
        A = A + term;
        if (k >= 3 && at < eps * abs(A)) break;
        prev_abs = at;
        w_power = w_power * w2;
    }

    // sign: (-1)^{n+1}
    C asymp = (n % 2 == 0) ? C(Float(0)) - A : A;

    // shift correction: (-1)^{n+1} · n! · shift_sum
    C shift_correction = C(n_fact) * shift_sum;
    if (n % 2 == 0) shift_correction = C(Float(0)) - shift_correction;

    C result = asymp + shift_correction;
    result.re.setPrecision(precision);
    result.im.setPrecision(precision);
    return result;
}

[[nodiscard]] inline Complex<Float> polygamma(int n, Complex<Float> z) {
    return polygamma(n, std::move(z), Float::defaultPrecision());
}

// ----------------------------------------------------------------
// B(a,b) = exp(lnΓ(a) + lnΓ(b) - lnΓ(a+b)) (Complex<Float>)
// ----------------------------------------------------------------

[[nodiscard]] inline Complex<Float> beta(Complex<Float> a, Complex<Float> b, int precision) {
    int wp = precision + 10;
    // PrecisionGuard removed: lnGamma(a, wp) etc. produce results with req=wp at explicit precision,
    // and the subsequent Complex<Float> exp uses requestedPrecision (= wp) via the 1-arg ADL,
    // so manipulating defaultPrecision is unnecessary.

    auto result = exp(lnGamma(a, wp) + lnGamma(b, wp) - lnGamma(a + b, wp));

    result.re.setPrecision(precision);
    result.im.setPrecision(precision);
    return result;
}

[[nodiscard]] inline Complex<Float> beta(Complex<Float> a, Complex<Float> b) {
    return beta(std::move(a), std::move(b), Float::defaultPrecision());
}

// ================================================================
// Complex<double> incomplete gamma functions
// ================================================================
// P(a,z) = γ(a,z)/Γ(a) — Taylor series (|z| < Re(a)+1) / 1-Q (otherwise)
// Q(a,z) = Γ(a,z)/Γ(a) — Legendre CF (|z| ≥ Re(a)+1) / 1-P (otherwise)

template<IsNativeFloat R>
[[nodiscard]] Complex<R> gammaP(Complex<R> a, Complex<R> z);
template<IsNativeFloat R>
[[nodiscard]] Complex<R> gammaQ(Complex<R> a, Complex<R> z);

template<IsNativeFloat R>
[[nodiscard]] Complex<R> gammaP(Complex<R> a, Complex<R> z) {
    using C = Complex<R>;
    if (sangi::abs(z) == R(0)) return C(R(0));

    // prefactor: exp(-z + a·ln(z) - lnΓ(a))
    C front = exp(-z + a * log(z) - lnGamma(a));

    if (sangi::abs(z) < a.re + R(1)) {
        // Taylor series: S = Σ z^n / (a·(a+1)···(a+n))
        C term = C(R(1)) / a;
        C sum = term;
        for (int n = 1; n < 500; n++) {
            term = term * z / (a + C(R(n)));
            sum = sum + term;
            if (sangi::abs(term) < sangi::abs(sum) * std::numeric_limits<R>::epsilon())
                break;
        }
        return front * sum;
    } else {
        return C(R(1)) - gammaQ(a, z);
    }
}

template<IsNativeFloat R>
[[nodiscard]] Complex<R> gammaQ(Complex<R> a, Complex<R> z) {
    using C = Complex<R>;
    if (sangi::abs(z) == R(0)) return C(R(1));

    C front = exp(-z + a * log(z) - lnGamma(a));

    if (sangi::abs(z) >= a.re + R(1)) {
        // Legendre CF: Q = front / f
        constexpr R tiny = std::numeric_limits<R>::min();

        C b0 = z + C(R(1)) - a;
        C f = (sangi::abs(b0) < tiny) ? C(tiny) : b0;
        C big_C = f;
        C D(R(0));

        for (int n = 1; n < 500; n++) {
            C an = C(R(-n)) * (C(R(n)) - a);
            C bn = z + C(R(2 * n + 1)) - a;

            D = bn + an * D;
            if (sangi::abs(D) < tiny) D = C(tiny);
            D = C(R(1)) / D;

            big_C = bn + an / big_C;
            if (sangi::abs(big_C) < tiny) big_C = C(tiny);

            C delta = big_C * D;
            f = f * delta;

            if (sangi::abs(delta - C(R(1))) < std::numeric_limits<R>::epsilon())
                break;
        }
        return front / f;
    } else {
        return C(R(1)) - gammaP(a, z);
    }
}

template<IsNativeFloat R>
[[nodiscard]] Complex<R> gammaLower(Complex<R> a, Complex<R> z) {
    return gammaP(a, z) * gamma(a);
}

template<IsNativeFloat R>
[[nodiscard]] Complex<R> gammaUpper(Complex<R> a, Complex<R> z) {
    return gammaQ(a, z) * gamma(a);
}

// ================================================================
// Complex<double> regularized incomplete beta I_z(a,b)
// ================================================================

template<IsNativeFloat R>
[[nodiscard]] Complex<R> betaRegularized(Complex<R> z, Complex<R> a, Complex<R> b) {
    using C = Complex<R>;
    if (sangi::abs(z) == R(0)) return C(R(0));
    if (sangi::abs(z - C(R(1))) < std::numeric_limits<R>::epsilon()) return C(R(1));

    // symmetry inversion: invert if |z| ≥ (Re(a)+1)/(Re(a)+Re(b)+2)
    bool flip = sangi::abs(z) >= (a.re + R(1)) / (a.re + b.re + R(2));
    C zz = flip ? (C(R(1)) - z) : z;
    C aa = flip ? b : a;
    C bb = flip ? a : b;

    // front = exp(a·ln(z) + b·ln(1-z) + lnΓ(a+b) - lnΓ(a) - lnΓ(b)) / a
    C ln_front = aa * log(zz) + bb * log(C(R(1)) - zz)
               + lnGamma(aa + bb) - lnGamma(aa) - lnGamma(bb);
    C front = exp(ln_front) / aa;

    // CF: Modified Lentz
    constexpr R tiny = std::numeric_limits<R>::min();
    C f(R(1)), big_C(R(1)), D(R(0));

    for (int n = 1; n < 500; n++) {
        C d;
        if (n % 2 == 1) {
            int m = (n - 1) / 2;
            C am = aa + C(R(m));
            C abm = aa + bb + C(R(m));
            C denom = (aa + C(R(2 * m))) * (aa + C(R(2 * m + 1)));
            d = -(am * abm * zz) / denom;
        } else {
            int m = n / 2;
            C c_fm{R(m)};
            C bm = bb - C(R(m));
            C denom = (aa + C(R(2 * m - 1))) * (aa + C(R(2 * m)));
            d = (c_fm * bm * zz) / denom;
        }

        C one_c{R(1)};
        D = one_c + d * D;
        if (sangi::abs(D) < tiny) D = C(tiny);
        D = one_c / D;

        big_C = one_c + d / big_C;
        if (sangi::abs(big_C) < tiny) big_C = C(tiny);

        C delta = big_C * D;
        f = f * delta;

        if (sangi::abs(delta - one_c) < std::numeric_limits<R>::epsilon())
            break;
    }

    C result = front / f;
    return flip ? (C{R(1)} - result) : result;
}

// ================================================================
// Complex<Float> incomplete gamma functions
// ================================================================

namespace detail {

template<typename R>
Complex<R> gammaPQ_taylor(const Complex<R>& a, const Complex<R>& z,
                           const Complex<R>& front, R eps, int max_iter) {
    using C = Complex<R>;
    C term = C(R(1)) / a;
    C sum = term;
    for (int n = 1; n < max_iter; n++) {
        term = term * z / (a + C(R(n)));
        sum = sum + term;
        if (n >= 3 && sangi::abs(term) < sangi::abs(sum) * eps)
            break;
    }
    return front * sum;
}

template<typename R>
Complex<R> gammaPQ_cf(const Complex<R>& a, const Complex<R>& z,
                       const Complex<R>& front, R eps, int max_iter) {
    using C = Complex<R>;
    R tiny = eps * eps;

    C b0 = z + C(R(1)) - a;
    C f = (sangi::abs(b0) < sangi::abs(C(tiny))) ? C(tiny) : b0;
    C big_C = f;
    C D(R(0));

    for (int n = 1; n < max_iter; n++) {
        C an = C(R(-n)) * (C(R(n)) - a);
        C bn = z + C(R(2 * n + 1)) - a;

        D = bn + an * D;
        if (sangi::abs(D) < tiny) D = C(tiny);
        D = C(R(1)) / D;

        big_C = bn + an / big_C;
        if (sangi::abs(big_C) < tiny) big_C = C(tiny);

        C delta = big_C * D;
        f = f * delta;

        if (n >= 3 && sangi::abs(delta - C(R(1))) < eps)
            break;
    }
    return front / f;
}

} // namespace detail

[[nodiscard]] inline Complex<Float> gammaP(Complex<Float> a, Complex<Float> z, int precision) {
    using C = Complex<Float>;
    int wp = precision + 25;
    // PrecisionGuard removed: a, z's setResultPrecision(wp) propagates req=wp.
    a.re.setResultPrecision(wp); a.im.setResultPrecision(wp);
    z.re.setResultPrecision(wp); z.im.setResultPrecision(wp);

    if (sangi::abs(z) == Float(0)) {
        Float zero(0); zero.setPrecision(precision);
        return C(zero);
    }

    // prefactor
    C front = exp(-z + a * log(z) - lnGamma(a, wp));
    Float eps = Float::epsilon(wp);
    int max_iter = 10 * (wp / 50 + 1);
    if (max_iter < 500) max_iter = 500;

    C result;
    if (sangi::abs(z).toDouble() < a.re.toDouble() + 1.0) {
        result = detail::gammaPQ_taylor<Float>(a, z, front, eps, max_iter);
    } else {
        result = C(Float(1)) - detail::gammaPQ_cf<Float>(a, z, front, eps, max_iter);
    }

    result.re.setPrecision(precision);
    result.im.setPrecision(precision);
    return result;
}

[[nodiscard]] inline Complex<Float> gammaQ(Complex<Float> a, Complex<Float> z, int precision) {
    using C = Complex<Float>;
    int wp = precision + 25;
    // PrecisionGuard removed: a, z's setResultPrecision(wp) propagates req=wp.
    a.re.setResultPrecision(wp); a.im.setResultPrecision(wp);
    z.re.setResultPrecision(wp); z.im.setResultPrecision(wp);

    if (sangi::abs(z) == Float(0)) {
        Float one(1); one.setPrecision(precision);
        return C(one);
    }

    C front = exp(-z + a * log(z) - lnGamma(a, wp));
    Float eps = Float::epsilon(wp);
    int max_iter = 10 * (wp / 50 + 1);
    if (max_iter < 500) max_iter = 500;

    C result;
    if (sangi::abs(z).toDouble() >= a.re.toDouble() + 1.0) {
        result = detail::gammaPQ_cf<Float>(a, z, front, eps, max_iter);
    } else {
        result = C(Float(1)) - detail::gammaPQ_taylor<Float>(a, z, front, eps, max_iter);
    }

    result.re.setPrecision(precision);
    result.im.setPrecision(precision);
    return result;
}

[[nodiscard]] inline Complex<Float> gammaLower(Complex<Float> a, Complex<Float> z, int precision) {
    using C = Complex<Float>;
    int wp = precision + 15;
    // PrecisionGuard removed: explicit wp to gammaP / gamma; result has req=wp.
    C p = gammaP(a, z, wp);
    C ga = gamma(a, wp);
    C result = p * ga;
    result.re.setPrecision(precision);
    result.im.setPrecision(precision);
    return result;
}

[[nodiscard]] inline Complex<Float> gammaUpper(Complex<Float> a, Complex<Float> z, int precision) {
    using C = Complex<Float>;
    int wp = precision + 15;
    // PrecisionGuard removed: explicit wp to gammaQ / gamma; result has req=wp.
    C q = gammaQ(a, z, wp);
    C ga = gamma(a, wp);
    C result = q * ga;
    result.re.setPrecision(precision);
    result.im.setPrecision(precision);
    return result;
}

// ================================================================
// Complex<Float> regularized incomplete beta I_z(a,b)
// ================================================================

[[nodiscard]] inline Complex<Float> betaRegularized(Complex<Float> z, Complex<Float> a,
                                                      Complex<Float> b, int precision) {
    using C = Complex<Float>;
    int wp = precision + 25;
    // PrecisionGuard removed: z, a, b's setResultPrecision(wp) propagates req=wp.
    z.re.setResultPrecision(wp); z.im.setResultPrecision(wp);
    a.re.setResultPrecision(wp); a.im.setResultPrecision(wp);
    b.re.setResultPrecision(wp); b.im.setResultPrecision(wp);

    if (sangi::abs(z) == Float(0)) {
        Float zero(0); zero.setPrecision(precision);
        return C(zero);
    }

    Float eps = Float::epsilon(wp);
    int max_iter = 10 * (wp / 50 + 1);
    if (max_iter < 500) max_iter = 500;

    // symmetry inversion
    Float one(1), zero_f(0);
    bool flip = sangi::abs(z).toDouble() >= (a.re.toDouble() + 1.0) / (a.re.toDouble() + b.re.toDouble() + 2.0);
    C zz = flip ? (C(one) - z) : z;
    C aa = flip ? b : a;
    C bb = flip ? a : b;

    // front = exp(a·ln(z) + b·ln(1-z) + lnΓ(a+b) - lnΓ(a) - lnΓ(b)) / a
    C ln_front = aa * log(zz) + bb * log(C(one) - zz)
               + lnGamma(aa + bb, wp) - lnGamma(aa, wp) - lnGamma(bb, wp);
    C front = exp(ln_front) / aa;

    // CF: Modified Lentz
    Float tiny = eps * eps;
    C f{one}, big_C{one}, D{zero_f};

    for (int n = 1; n < max_iter; n++) {
        C d;
        if (n % 2 == 1) {
            int m = (n - 1) / 2;
            Float fm(m), f2m(2 * m), f2m1(2 * m + 1);
            C am = aa + C(fm);
            C abm = aa + bb + C(fm);
            C denom = (aa + C(f2m)) * (aa + C(f2m1));
            d = -(am * abm * zz) / denom;
        } else {
            int m = n / 2;
            Float fm(m), fbm(m), f2m(2 * m), f2m_1(2 * m - 1);
            C c_fm{fm};
            C bm = bb - c_fm;
            C denom = (aa + C(f2m_1)) * (aa + C(f2m));
            d = (c_fm * bm * zz) / denom;
        }

        D = C(one) + d * D;
        if (sangi::abs(D) < tiny) D = C(tiny);
        D = C(one) / D;

        big_C = C(one) + d / big_C;
        if (sangi::abs(big_C) < tiny) big_C = C(tiny);

        C delta = big_C * D;
        f = f * delta;

        if (n >= 3 && sangi::abs(delta - C(one)) < eps)
            break;
    }

    C result = front / f;
    if (flip) result = C(one) - result;

    result.re.setPrecision(precision);
    result.im.setPrecision(precision);
    return result;
}

} // namespace special
} // namespace sangi

#endif // SANGI_SPECIAL_GAMMA_HPP