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

// error_function.hpp
// Template wrappers for the error function family
//
// Functions provided:
//   erf(x)    — error function
//   erfc(x)   — complementary error function 1 - erf(x)
//   erfcx(x)  — scaled complementary error function exp(x²)·erfc(x)
//   erfi(x)   — imaginary error function erfi(x) = -i·erf(i·x) = (2/√π) ∫_0^x e^{t²} dt
//   erfInv(x) — inverse error function erfInv(erf(x)) = x, domain (-1, 1)
//
// Supported types:
//   float, double, long double — delegated to std:: / direct computation
//   Float                      — delegated to the implementation in FloatMath.cpp
//   Complex<double>            — series + Laplace continued fraction
//   Complex<Float>             — arbitrary-precision Laplace continued fraction

#ifndef SANGI_SPECIAL_ERROR_FUNCTION_HPP
#define SANGI_SPECIAL_ERROR_FUNCTION_HPP

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

namespace sangi {
namespace special {

// ================================================================
// error function erf(x) = (2/√π) ∫₀ˣ e^{-t²} dt
// ================================================================

/// Native floating-point types: delegated to std::erf
template<IsNativeFloat T>
[[nodiscard]] T erf(T x) {
    return std::erf(x);
}

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

// ================================================================
// complementary error function erfc(x) = 1 - erf(x)
// ================================================================

/// Native floating-point types: delegated to std::erfc
template<IsNativeFloat T>
[[nodiscard]] T erfc(T x) {
    return std::erfc(x);
}

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

// ================================================================
// scaled complementary error function erfcx(x) = exp(x²)·erfc(x)
// ================================================================

/// Native floating-point types: direct computation + asymptotic expansion
template<IsNativeFloat T>
[[nodiscard]] T erfcx(T x) {
    if (std::isnan(x)) return std::numeric_limits<T>::quiet_NaN();
    if (std::isinf(x)) {
        return x > T(0) ? T(0) : std::numeric_limits<T>::infinity();
    }
    if (x == T(0)) return T(1);

    // x > 0 and large: asymptotic expansion
    // erfcx(x) = (1/(x√π)) · Σ_{n=0}^{N} (-1)^n (2n-1)!! / (2x²)^n
    if (x > T(4)) {
        T t = T(1) / (T(2) * x * x);
        T series = T(1);
        T term = T(1);
        for (int k = 1; k <= 25; k++) {
            term *= -static_cast<T>(2 * k - 1) * t;
            if (std::abs(term) > std::abs(series)) break;  // divergence begins
            series += term;
            if (std::abs(term) < std::numeric_limits<T>::epsilon()) break;
        }
        const T sqrt_pi = std::sqrt(static_cast<T>(3.14159265358979323846L));
        return series / (x * sqrt_pi);
    }

    // x > 0: direct computation (range where exp(x²) does not overflow)
    if (x > T(0)) {
        return std::exp(x * x) * std::erfc(x);
    }

    // x < 0: erfcx(x) = 2·exp(x²) - erfcx(-x)
    T erfcx_pos = erfcx(-x);
    T exp_x2 = std::exp(x * x);
    if (std::isinf(exp_x2)) return exp_x2;
    return T(2) * exp_x2 - erfcx_pos;
}

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

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

// ----------------------------------------------------------------
// erf(z) = (2z/√π)·e^{-z²}·Σ_{n=0}^∞ (2z²)^n / (1·3·5···(2n+1))
// ----------------------------------------------------------------
// Series form that converges for all z (₁F₁ based).
// Accumulating only positive terms avoids cancellation.

template<SpecialFunctionScalar R>
[[nodiscard]] Complex<R> erf(Complex<R> z) {
    using C = Complex<R>;
    if (z == C(R(0))) return C(R(0));

    C z2 = z * z;
    C two_z2 = C(R(2)) * z2;
    C term(R(1));
    C sum(R(1));

    // BUGFIX (2026-05-31): getEpsilon<R>()/getPi<R>() (no argument) return
    // defaultPrecision (~55 digits), so the convergence test was cut off early at
    // ~55 digits and the constant 2/√π was also ~55 digits, capping the
    // Complex<Float> version at ~55 digits regardless of the argument's working
    // precision (detected by the Complex cross-check harness for erfcx). Take the
    // working precision of argument z and use the precision-aware version.
    const int max_iter = IsSangiFloat<R> ? 2000 : 200;
    R eps;
    if constexpr (IsSangiFloat<R>) eps = detail::getEpsilon<R>(z.re.precision());
    else                           eps = detail::getEpsilon<R>();

    for (int n = 0; n < max_iter; n++) {
        term = term * two_z2 / C(R(2 * n + 3));
        sum = sum + term;

        if (abs(term) < eps * abs(sum))
            break;
    }

    R two_over_sqrt_pi;
    if constexpr (IsSangiFloat<R>) {
        int wp = z.re.precision();
        two_over_sqrt_pi = R(2) / sqrt(detail::getPi<R>(wp), wp);
    } else
        two_over_sqrt_pi = R(1.1283791670955125738961589031215452L);
    return C(two_over_sqrt_pi) * z * exp(-z2) * sum;
}

// ----------------------------------------------------------------
// erfc(z) = 1 - erf(z)
// ----------------------------------------------------------------

template<SpecialFunctionScalar R>
[[nodiscard]] Complex<R> erfc(Complex<R> z) {
    return Complex<R>(R(1)) - erf(z);
}

// ----------------------------------------------------------------
// erfcx(z) = exp(z²)·erfc(z)  — Complex<double>
// ----------------------------------------------------------------
// Re(z) > 0 and |z| > 2: Laplace continued fraction (Modified Lentz)
//   g = z + a₁/(z + a₂/(z + ...)), a_n = n/2
//   erfcx(z) = 1/(√π · g)
// |z| ≤ 2: direct computation exp(z²)·erfc(z)
// Re(z) < 0 and |z| > 2: reflection erfcx(z) = 2·exp(z²) - erfcx(-z)

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

    // small |z|: direct computation (exp(z²) does not overflow)
    R z_abs = sangi::abs(z);
    if (z_abs <= R(2)) {
        C z2 = z * z;
        return exp(z2) * erfc(z);
    }

    // Re(z) ≥ 0: Laplace continued fraction
    if (z.re >= R(0)) {
        R sqrt_pi = std::sqrt(std::numbers::pi_v<R>);
        R eps = std::numeric_limits<R>::epsilon();
        R tiny = R(1e-300);

        C f = z;
        C big_C = z;
        C D(R(0));

        for (int n = 1; n <= 300; n++) {
            C a_n{R(n) / R(2)};

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

            big_C = z + a_n / 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))) < eps)
                break;
        }
        return C(R(1)) / (C(sqrt_pi) * f);
    }

    // Re(z) < 0: reflection erfcx(z) = 2·exp(z²) - erfcx(-z)
    C z2 = z * z;
    return C(R(2)) * exp(z2) - erfcx(C(R(0)) - z);
}

// ----------------------------------------------------------------
// erfcx(z) = exp(z²)·erfc(z)  — Complex<Float>
// ----------------------------------------------------------------

[[nodiscard]] inline Complex<Float> erfcx(Complex<Float> z, int precision) {
    using C = Complex<Float>;
    int wp = precision + 20;
    // PrecisionGuard removed: setResultPrecision(wp) on w.re/im propagates req=wp
    // (FLOAT_PRECISION_PLAN Step 12 follow-up, same pattern as gamma.hpp 6fa993f).

    C w = z;
    w.re.setResultPrecision(wp);
    w.im.setResultPrecision(wp);
    // BUGFIX (2026-05-31): for exact input (e.g. 1+1i, re/im being finite binary
    // fractions), the internal exact/exact division drops to default_precision
    // (~55 digits) (detected by the Complex cross-check harness). Add working
    // precision to make it non-exact.
    w.re.setEffectiveBits(Float::precisionToBits(wp));
    w.im.setEffectiveBits(Float::precisionToBits(wp));

    double z_abs = std::sqrt(w.re.toDouble() * w.re.toDouble()
                           + w.im.toDouble() * w.im.toDouble());

    // small |z|: direct computation erfcx = exp(z²)·erfc(z).
    // BUGFIX (2026-05-31): raised the threshold from 2 to 4. The Laplace continued
    // fraction converges extremely slowly near |z|≈2 and capped out at about 18
    // digits (detected by harness erfcx(0.3+2i)). Actually apply to the erfc call
    // the wp2 guard that compensates for erfc's cancellation (z²/ln2 bits)
    // (previously wp2 was computed but left unused).
    if (z_abs <= 4.0) {
        int extra = static_cast<int>(std::ceil(z_abs * z_abs / std::log(2.0))) + 10;
        int wp2 = wp + extra;
        C ww = z;
        ww.re.setResultPrecision(wp2);
        ww.im.setResultPrecision(wp2);
        ww.re.setEffectiveBits(Float::precisionToBits(wp2));
        ww.im.setEffectiveBits(Float::precisionToBits(wp2));
        C erfc_val = erfc(ww);
        C z2 = ww * ww;
        C result = exp(z2) * erfc_val;
        result.re.setPrecision(precision);
        result.im.setPrecision(precision);
        return result;
    }

    // Re(z) < 0: reflection erfcx(z) = 2·exp(z²) - erfcx(-z)
    if (w.re.toDouble() < 0.0) {
        C neg_w = -w;
        C erfcx_pos = erfcx(neg_w, precision + 10);
        erfcx_pos.re.setResultPrecision(wp);
        erfcx_pos.im.setResultPrecision(wp);
        C z2 = w * w;
        C result = C(Float(2)) * exp(z2) - erfcx_pos;
        result.re.setPrecision(precision);
        result.im.setPrecision(precision);
        return result;
    }

    // Re(z) ≥ 0, |z| > 2: Laplace continued fraction (Modified Lentz)
    Float two(2);
    Float eps = Float::epsilon(wp);

    C f = w;
    C big_C = w;
    C D;

    for (int n = 1; n < 10 * wp; n++) {
        C a_n{Float(n) / two};

        D = w + a_n * D;
        if (abs(D) < eps) D = C(eps);
        D = C(Float(1)) / D;

        big_C = w + a_n / big_C;
        if (abs(big_C) < eps) big_C = C(eps);

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

        // convergence test
        C diff = delta - C(Float(1));
        Float diff_abs = abs(diff);
        if (diff_abs.isZero() || diff_abs < eps) break;
    }

    Float sqrtpi = sangi::sqrt(Float::pi(wp), wp);
    C result = C(Float(1)) / (C(sqrtpi) * f);
    result.re.setPrecision(precision);
    result.im.setPrecision(precision);
    return result;
}

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

// ================================================================
// imaginary error function erfi(x) = (2/√π) ∫_0^x e^{t²} dt
// ================================================================
// Properties: odd function, erfi(0) = 0, erfi(∞) = +∞ (diverges on the order of e^{x²})
// Relation: erfi(x) = -i · erf(i·x) = imag(erf(0 + i·x))
//
// Implementation: reuses the existing sangi::special::erf(Complex<double>) (Taylor +
//       Laplace continued fraction, ~1e-15 accuracy across the whole double range).
//       Since erf(i·x) is purely imaginary, just returning its imag part is OK.
//
// Overflow: in double precision, e^{x²} overflows at x ≈ 27 → ±Inf.

template<IsNativeFloat T>
[[nodiscard]] T erfi(T x) {
    if (std::isnan(x)) return std::numeric_limits<T>::quiet_NaN();
    if (x == T(0)) return T(0);
    if (std::isinf(x)) {
        return x > T(0) ? std::numeric_limits<T>::infinity()
                        : -std::numeric_limits<T>::infinity();
    }

    double xd = static_cast<double>(x);
    if (std::abs(xd) > 27.0) {
        return xd > 0.0 ? std::numeric_limits<T>::infinity()
                        : -std::numeric_limits<T>::infinity();
    }

    // erfi(x) = imag(erf(i x))
    Complex<double> r = erf(Complex<double>(0.0, xd));
    return static_cast<T>(r.im);
}

// ================================================================
// inverse error function erfInv(p), p ∈ (-1, 1)
// ================================================================
// Solution: y = erfInv(p) is the real number satisfying erf(y) = p
//
// Implementation: Newton's method (machine precision in 3-4 iterations)
//   y_{n+1} = y_n - (erf(y_n) - p) / (2/√π · e^{-y_n²})
// Initial guess:
//   |p| ≤ 0.7  : y_0 = p · (a + b·p²) (Acklam-style minimax)
//   |p| > 0.7  : y_0 = sign(p) · √(-ln((1-|p|)(1+|p|)))
//
// Degenerate: p = ±1 → ±∞, |p| > 1 → NaN

template<IsNativeFloat T>
[[nodiscard]] T erfInv(T p) {
    if (std::isnan(p)) return std::numeric_limits<T>::quiet_NaN();
    if (p == T(0)) return T(0);
    if (p == T(1)) return std::numeric_limits<T>::infinity();
    if (p == T(-1)) return -std::numeric_limits<T>::infinity();
    if (std::abs(p) > T(1)) return std::numeric_limits<T>::quiet_NaN();

    constexpr T sqrt_pi  = T(1.7724538509055160272981674833);
    constexpr T half_sqrt_pi = sqrt_pi / T(2);

    T ap = std::abs(p);
    T y;
    if (ap <= T(0.7)) {
        // first-order initial guess: erfInv(p) ≈ (√π/2) · p
        y = half_sqrt_pi * p;
    } else {
        // tail approximation: erfInv(p) ~ √(-ln((1-|p|)(1+|p|))) · sign(p)
        T t = std::sqrt(-std::log((T(1) - ap) * (T(1) + ap)));
        y = (p > T(0) ? T(1) : T(-1)) * t;
    }

    // Newton iteration (up to 8 times until convergence)
    for (int iter = 0; iter < 8; ++iter) {
        T erf_y = std::erf(y);
        T err = erf_y - p;
        if (std::abs(err) <
            std::numeric_limits<T>::epsilon() * (T(1) + std::abs(y)))
            break;
        T dy = err / (T(2) / sqrt_pi * std::exp(-y * y));
        y -= dy;
        if (std::abs(dy) <
            std::numeric_limits<T>::epsilon() * (T(1) + std::abs(y)))
            break;
    }
    return y;
}

} // namespace special
} // namespace sangi

#endif // SANGI_SPECIAL_ERROR_FUNCTION_HPP