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

// owens_t.hpp
// Owen's T function template
//
// Functions provided:
//   owensT(h, a)  — T(h, a) = (1/(2π)) ∫₀ᵃ exp(-h²(1+t²)/2)/(1+t²) dt
//
// Owen's T function is used to compute orthant probabilities of the bivariate normal distribution.
// It is also the basis for CDF computation of the noncentral t and skew-normal distributions.
//
// Supported types:
//   float, double, long double — 10-point Gauss-Legendre quadrature
//   Float                      — Cauchy-product series (arbitrary precision)
//
// References:
//   Owen, D.B. (1956) "Tables for computing bivariate normal probabilities"
//   Patefield & Tandy (2000) "Fast and accurate calculation of Owen's T function"
//   No DLMF entry (a statistics-specific function)

#ifndef SANGI_SPECIAL_OWENS_T_HPP
#define SANGI_SPECIAL_OWENS_T_HPP

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

namespace sangi {
namespace special {

// ================================================================
// Owen's T function (double)
// ================================================================
// T(h, a) = (1/(2π)) ∫₀ᵃ exp(-h²(1+t²)/2)/(1+t²) dt
//
// double: compute ∫₀ᵃ exp(-h²(1+t²)/2)/(1+t²) dt directly via Gauss-Legendre quadrature
// Float:  expand the Cauchy product Σ d_n t^{2n} of exp(-h²t²/2)/(1+t²) as a series
//         d_n = a_n − d_{n−1},  a_n = (−h²/2)^n/n!
// Transformation for a > 1 (Owen's relation):
//   T(h, a) = (Φ(h)+Φ(ah))/2 − Φ(h)·Φ(ah) − T(ah, 1/a)
//
// Implementation:
//   double — 10-point Gauss-Legendre quadrature (high accuracy for all a)
//   Float  — Cauchy-product series + the a > 1 transformation

namespace detail {

/// Owen's T — 10-point Gauss-Legendre quadrature (for native float)
/// Compute ∫₀ᵃ exp(-h²(1+t²)/2)/(1+t²) dt by mapping to [-1,1]
template<IsNativeFloat T>
[[nodiscard]] T owensT_gl(T h, T a) {
    // 10-point GL nodes/weights on [-1, 1]
    static constexpr double gl_x[] = {
        -0.97390652851717172, -0.86506336668898451,
        -0.67940956829902441, -0.43339539412924719,
        -0.14887433898163121,
         0.14887433898163121,  0.43339539412924719,
         0.67940956829902441,  0.86506336668898451,
         0.97390652851717172
    };
    static constexpr double gl_w[] = {
        0.06667134430868814, 0.14945134915058059,
        0.21908636251598204, 0.26926671930999636,
        0.29552422471475287,
        0.29552422471475287, 0.26926671930999636,
        0.21908636251598204, 0.14945134915058059,
        0.06667134430868814
    };
    // Change of variable: t = a(1+x)/2, dt = a/2 dx, x ∈ [-1,1]
    T half_a = a / T(2);
    T h2 = h * h;
    T sum = T(0);
    for (int i = 0; i < 10; i++) {
        T t = half_a * (T(1) + T(gl_x[i]));
        T one_plus_t2 = T(1) + t * t;
        sum += T(gl_w[i]) * std::exp(-h2 * one_plus_t2 / T(2)) / one_plus_t2;
    }
    return half_a * sum / (T(2) * std::numbers::pi_v<T>);
}

/// Owen's T series expansion — for Float (the Cauchy product correctly accounts for 1/(1+t²))
/// f(t) = exp(-h²t²/2)/(1+t²) = Σ d_n t^{2n}, d_n = a_n − d_{n−1}
///   (a_n = (−h²/2)^n / n! is the Taylor coefficient of exp(−h²t²/2))
/// T(h,a) = exp(−h²/2)/(2π) · Σ d_n a^{2n+1}/(2n+1)
template<typename R>
[[nodiscard]] R owensT_series(R h, R a, R eps, int max_iter, int wp) {
    R h_sq_half = h * h / R(2);
    R exp_val = sangi::exp(-h_sq_half, wp);

    R a_sq = a * a;
    R a_pow = a;               // a^{2n+1}
    R neg_h_sq_half = -h_sq_half;
    R a_n = R(1);              // (-h²/2)^n / n!
    R d_n = R(1);              // Cauchy-product coefficient, d_0 = a_0 = 1
    R sum = d_n * a;           // n=0 term

    for (int n = 1; n < max_iter; n++) {
        a_n = a_n * neg_h_sq_half / R(n);   // a_n = (-h²/2)^n / n!
        d_n = a_n - d_n;                     // d_n = a_n − d_{n−1}
        a_pow = a_pow * a_sq;
        R term = d_n * a_pow / R(2 * n + 1);
        sum = sum + term;

        if (n >= 5) {
            R abs_term = sangi::abs(term);
            R abs_sum = sangi::abs(sum);
            if (abs_term < eps * abs_sum) break;
        }
    }

    R pi_val = Float::pi(wp);          // ★ π at wp precision (was: capped at ~56 digits via defaultPrecision)
    return exp_val * sum / (R(2) * pi_val);
}

} // namespace detail

/// Owen's T function T(h, a)
template<IsNativeFloat T>
[[nodiscard]] T owensT(T h, T a) {
    if (std::isnan(h) || std::isnan(a)) return std::numeric_limits<T>::quiet_NaN();

    // Special cases
    if (a == T(0)) return T(0);
    if (h == T(0)) return std::atan(a) / (T(2) * std::numbers::pi_v<T>);

    // T(h, a) is even in h: T(-h, a) = T(h, a)
    T abs_h = std::abs(h);

    // T(h, -a) = -T(h, a)
    T sign_a = (a < T(0)) ? T(-1) : T(1);
    T abs_a = std::abs(a);

    T result;

    if (abs_a <= T(1)) {
        // 10-point GL quadrature
        result = detail::owensT_gl(abs_h, abs_a);
    } else {
        // Owen's transformation (a > 1):
        // T(h, a) = 1/2 · [Φ(h) + Φ(ah)] - Φ(h)·Φ(ah) - T(ah, 1/a)
        T inv_a = T(1) / abs_a;
        T ah = abs_h * abs_a;

        T sqrt2 = std::numbers::sqrt2_v<T>;
        T phi_h = std::erfc(-abs_h / sqrt2) / T(2);
        T phi_ah = std::erfc(-ah / sqrt2) / T(2);

        T t_ah_inv_a = detail::owensT_gl(ah, inv_a);

        result = (phi_h + phi_ah) / T(2) - phi_h * phi_ah - t_ah_inv_a;
    }

    return sign_a * result;
}

/// Owen's T function — Float (arbitrary precision)
[[nodiscard]] inline Float owensT(const Float& h, const Float& a, int precision) {
    int wp = precision + 20;
    // PrecisionGuard removed: setResultPrecision(wp) on hw/aw propagates req=wp
    // (FLOAT_PRECISION_PLAN Step 12 follow-up).

    Float hw = h; hw.setResultPrecision(wp);
    Float aw = a; aw.setResultPrecision(wp);

    // Special cases
    if (aw.isZero()) { Float z(0); z.setPrecision(precision); return z; }
    if (hw.isZero()) {
        Float pi_val = Float::pi(wp);
        Float result = sangi::atan(aw) / mulScalarF(pi_val, uint64_t(2));
        result.setPrecision(precision);
        return result;
    }

    Float abs_h = sangi::abs(hw);
    Float sign_a = aw.isNegative() ? Float(-1) : Float(1);
    Float abs_a = sangi::abs(aw);

    Float result;

    if (abs_a <= Float(1)) {
        result = detail::owensT_series<Float>(abs_h, abs_a,
            Float::epsilon(wp), 10 * wp, wp);
    } else {
        Float inv_a = Float(1) / abs_a;
        Float ah = abs_h * abs_a;

        // Φ(x) = (1 + erf(x/√2)) / 2
        Float sqrt2 = sangi::sqrt(Float(2), wp);  // ★ wp precision (was: capped at ~56 digits via 1-arg sqrt)
        Float erf_h = sangi::erf(abs_h / sqrt2, wp);
        Float erf_ah = sangi::erf(ah / sqrt2, wp);
        Float phi_h = ldexp(Float(1) + erf_h, -1);
        Float phi_ah = ldexp(Float(1) + erf_ah, -1);

        Float t_ah_inv_a = detail::owensT_series<Float>(ah, inv_a,
            Float::epsilon(wp), 10 * wp, wp);

        result = ldexp(phi_h + phi_ah, -1) - phi_h * phi_ah - t_ah_inv_a;
    }

    result = sign_a * result;
    result.setPrecision(precision);
    return result;
}

[[nodiscard]] inline Float owensT(const Float& h, const Float& a) {
    return owensT(h, a, Float::defaultPrecision());
}

} // namespace special
} // namespace sangi

#endif // SANGI_SPECIAL_OWENS_T_HPP