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

// exponential_integral.hpp
// Templates for exponential integrals, trigonometric integrals, and the dilogarithm
//
// Provided functions:
//   expint(x)          — exponential integral Ei(x)
//   li(x)              — logarithmic integral li(x) = Ei(ln(x))
//   dilog(x)           — dilogarithm Li₂(x)
//   expintN(n, x)      — generalized exponential integral E_n(x)
//   sinIntegral(x)     — sine integral Si(x)
//   cosIntegral(x)     — cosine integral Ci(x)
//   sinhIntegral(x)    — hyperbolic sine integral Shi(x)
//   coshIntegral(x)    — hyperbolic cosine integral Chi(x)
//
// Supported types:
//   float, double, long double — series / continued fractions
//   Float                      — delegated to the sangi:: implementation

#ifndef SANGI_SPECIAL_EXPONENTIAL_INTEGRAL_HPP
#define SANGI_SPECIAL_EXPONENTIAL_INTEGRAL_HPP

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

namespace sangi {
namespace special {

// Euler-Mascheroni constant
namespace detail {
    template<typename T>
    constexpr T euler_gamma = T(0.5772156649015328606065120900824024L);
}

// ================================================================
// Exponential integral Ei(x) = γ + ln|x| + Σ_{n=1}^∞ x^n/(n·n!)
// ================================================================

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

    T sum = detail::euler_gamma<T> + std::log(std::abs(x));
    T term = x;
    sum += term;  // n=1

    for (int n = 2; n < 200; n++) {
        term *= x / T(n);
        T contribution = term / T(n);
        sum += contribution;

        if (n >= 5 && std::abs(contribution) <
            std::numeric_limits<T>::epsilon() * std::abs(sum))
            break;
    }

    return sum;
}

// ================================================================
// Logarithmic integral li(x) = Ei(ln(x))
// ================================================================

template<IsNativeFloat T>
[[nodiscard]] T li(T x) {
    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(1)) return -std::numeric_limits<T>::infinity();
    if (std::isinf(x)) return std::numeric_limits<T>::infinity();
    return expint(std::log(x));
}

// ================================================================
// Dilogarithm Li₂(x) = -∫₀ˣ ln(1-t)/t dt = Σ_{n=1}^∞ x^n/n²
// ================================================================

template<IsNativeFloat T>
[[nodiscard]] T dilog(T x) {
    if (std::isnan(x)) return std::numeric_limits<T>::quiet_NaN();
    if (x == T(0)) return T(0);
    if (x == T(1)) return std::numbers::pi_v<T> * std::numbers::pi_v<T> / T(6);
    if (x == T(-1)) return -std::numbers::pi_v<T> * std::numbers::pi_v<T> / T(12);

    // |x| > 0.5: reflection formula Li₂(x) = -Li₂(1-x) + π²/6 - ln(x)·ln(1-x)
    if (std::abs(x) > T(0.5) && std::abs(x) <= T(1)) {
        T one_minus_x = T(1) - x;
        return -dilog(one_minus_x) +
               std::numbers::pi_v<T> * std::numbers::pi_v<T> / T(6) -
               std::log(x) * std::log(one_minus_x);
    }

    // |x| <= 0.5: Taylor series
    T sum = T(0);
    T x_power = x;

    for (int n = 1; n < 200; n++) {
        T term = x_power / (T(n) * T(n));
        sum += term;

        if (n >= 5 && std::abs(term) <
            std::numeric_limits<T>::epsilon() * std::abs(sum))
            break;

        x_power *= x;
    }

    return sum;
}

// ================================================================
// Generalized exponential integral E_n(x) = ∫₁^∞ e^{-xt}/t^n dt
// ================================================================
// n ≥ 0, x > 0 (n=1 and x=0 → +∞)
// E_1(x) = -Ei(-x) (x > 0)
// recurrence: n·E_{n+1}(x) = e^{-x} - x·E_n(x)
// continued fraction (large x): E_n(x) = e^{-x} · CF

template<IsNativeFloat T>
[[nodiscard]] T expintN(int n, T x) {
    if (std::isnan(x)) return std::numeric_limits<T>::quiet_NaN();
    if (n < 0) return std::numeric_limits<T>::quiet_NaN();
    if (x < T(0)) return std::numeric_limits<T>::quiet_NaN();

    // x = 0: E_n(0) = 1/(n-1) for n > 1, E_1(0) = +∞, E_0(0) = +∞
    if (x == T(0)) {
        if (n <= 1) return std::numeric_limits<T>::infinity();
        return T(1) / T(n - 1);
    }

    // E_0(x) = e^{-x}/x
    if (n == 0) return std::exp(-x) / x;

    // Continued fraction (Modified Lentz) — converges fast when x is large enough
    // CF: E_n(x) = e^{-x} / (x + n / (1 + 1 / (x + (n+1) / (1 + 2 / (x + ...)))))
    // Lentz form: b_0 = x+n, a_1 = 1·n, b_1 = 1, a_2 = 1·(n+1), b_2 = x+n+2, ...
    // (the Thompson-Barnett CF is more stable)
    if (x > T(1) || n > 10) {
        // Lentz: E_n(x) = e^{-x} · 1/(x + n/(1 + 1/(x + (n+1)/(1 + ...))))
        // implementation: standard Gauss CF
        constexpr T tiny = T(1e-30);
        T b = x + T(n);
        T c = T(1) / tiny;
        T d = T(1) / b;
        T h = d;

        for (int i = 1; i <= 200; i++) {
            T a = -T(i) * T(n - 1 + i);
            b += T(2);
            d = T(1) / (a * d + b);
            c = b + a / c;
            T delta = c * d;
            h *= delta;

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

        return std::exp(-x) * h;
    }

    // Series: E_n(x) = [(-x)^{n-1}/(n-1)!][ψ(n) - ln(x)] - Σ_{m=0,m≠n-1}^∞ (-x)^m/((m-n+1)·m!)
    // For E_1: E_1(x) = -γ - ln(x) - Σ_{k=1}^∞ (-x)^k/(k·k!)
    // For general n use the recurrence: first compute E_1, then the upward recurrence
    T e1;
    {
        // E_1(x) = -Ei(-x) = -γ - ln(x) + Σ_{k=1}^∞ (-1)^{k+1} x^k/(k·k!)
        T sum = -detail::euler_gamma<T> - std::log(x);
        T term = T(1);
        for (int k = 1; k < 200; k++) {
            term *= -x / T(k);
            T contribution = -term / T(k);
            sum += contribution;

            if (k >= 5 && std::abs(contribution) <
                std::numeric_limits<T>::epsilon() * std::abs(sum))
                break;
        }
        e1 = sum;
    }

    if (n == 1) return e1;

    // Upward recurrence: n·E_{n+1}(x) = e^{-x} - x·E_n(x)
    T emx = std::exp(-x);
    T en = e1;
    for (int k = 1; k < n; k++) {
        en = (emx - x * en) / T(k);
    }

    return en;
}

// ================================================================
// Sine integral Si(x) = ∫₀ˣ sin(t)/t dt
// ================================================================
// Si(x) = Σ_{k=0}^∞ (-1)^k · x^{2k+1} / ((2k+1)·(2k+1)!)
// Si(-x) = -Si(x) (odd function)
// Si(∞) = π/2

template<IsNativeFloat T>
[[nodiscard]] T sinIntegral(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::numbers::pi_v<T> / T(2)
                         : -std::numbers::pi_v<T> / T(2);
    }

    // odd function
    T sign = T(1);
    T ax = x;
    if (x < T(0)) {
        sign = T(-1);
        ax = -x;
    }

    // Taylor: term_0 = x, term_{k+1} = term_k · (-x²) / ((2k+2)(2k+3))
    // sum_k = term_k / (2k+1)
    T x2 = ax * ax;
    T term = ax;     // the x^1 / 1! part of x^{2k+1} / (2k+1)!
    T sum = ax;      // k=0: x / 1

    for (int k = 0; k < 200; k++) {
        term *= -x2 / (T(2 * k + 2) * T(2 * k + 3));
        T contribution = term / T(2 * k + 3);
        sum += contribution;

        if (k >= 3 && std::abs(contribution) <
            std::numeric_limits<T>::epsilon() * std::abs(sum))
            break;
    }

    return sign * sum;
}

// ================================================================
// Cosine integral Ci(x) = γ + ln|x| + ∫₀ˣ (cos(t)-1)/t dt
// ================================================================
// Ci(x) = γ + ln|x| + Σ_{k=1}^∞ (-1)^k · x^{2k} / (2k·(2k)!)
// Defined only for x > 0 (real range). Ci(0) = -∞.

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

    // Ci(x) = γ + ln(x) + Σ_{k=1}^∞ (-1)^k · x^{2k} / (2k·(2k)!)
    // term tracks (-1)^k · x^{2k} / (2k)!
    // term_{k} = term_{k-1} · (-x²) / ((2k-1)·2k)
    T sum = detail::euler_gamma<T> + std::log(x);
    T x2 = x * x;
    T term = T(1);  // k=0 initial value

    for (int k = 1; k < 200; k++) {
        term *= -x2 / (T(2 * k - 1) * T(2 * k));
        T contribution = term / T(2 * k);
        sum += contribution;

        if (k >= 3 && std::abs(contribution) <
            std::numeric_limits<T>::epsilon() * std::abs(sum))
            break;
    }

    return sum;
}

// Float type: sangi::expint, sangi::li, sangi::dilog already exist.
// Use sangi::expintN, sangi::sinIntegral, sangi::cosIntegral directly.

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

// ----------------------------------------------------------------
// Ei(z) = γ + ln(z) + Σ_{n=1}^∞ z^n / (n·n!)
// ----------------------------------------------------------------

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

    if (z == C(R(0)))
        return C(-detail::getInfinity<R>());

    // The no-precision-argument getEulerGamma<R>()/getEpsilon<R>() return defaultPrecision (~55 digits),
    // so take the working precision of the argument z and use the precision-aware versions (same fix as complex erf).
    R gam = detail::getEulerGamma<R>();
    R eps = detail::getEpsilon<R>();
    if constexpr (IsSangiFloat<R>) {
        const int wp = z.re.precision();
        gam = R::euler(wp);                 // detail::getEulerGamma has no precision-argument version, so use Float::euler
        eps = detail::getEpsilon<R>(wp);
    }

    C sum = C(gam) + log(z);
    C term = z;
    sum = sum + term;  // n=1

    const int max_iter = IsSangiFloat<R> ? 2000 : 200;

    for (int n = 2; n < max_iter; n++) {
        // term = term*z/n, contribution = term/n. The generic Complex/Complex division
        // loses precision on all-dyadic (exact) inputs (exact/exact poison), so divide
        // the real/imaginary parts by Float/Float(n) (same approach as dilog; avoids the ~55-digit cap on e.g. 0.5+0.5i).
        term = term * z;
        term = C(term.re / R(n), term.im / R(n));
        C contribution(term.re / R(n), term.im / R(n));
        sum = sum + contribution;

        if (n >= 5 && abs(contribution) < eps * abs(sum))
            break;
    }

    return sum;
}

// ----------------------------------------------------------------
// Li₂(z) — Taylor + analytic continuation
// ----------------------------------------------------------------
// |z| ≤ 0.5: Σ z^n/n²
// 0.5 < |z| ≤ 1: reflection Li₂(z) = -Li₂(1-z) + π²/6 - ln(z)·ln(1-z)
// |z| > 1: inversion Li₂(z) = -Li₂(1/z) - π²/6 - (ln(-z))²/2

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

    if (z == C(R(0))) return C(R(0));

    R pi_val = detail::getPi<R>();
    if constexpr (IsSangiFloat<R>) {
        pi_val = detail::getPi<R>(z.re.precision());
    }
    const R pi2_6 = pi_val * pi_val / R(6);

    // z = 1: Li₂(1) = π²/6
    if (z.re == R(1) && z.im == R(0)) {
        return C(pi2_6);
    }
    // z = -1: Li₂(-1) = -π²/12
    if (z.re == R(-1) && z.im == R(0)) {
        return C(-pi_val * pi_val / R(12));
    }

    R az = abs(z);

    // |z| > 1: inversion formula (1/|z| < 1 always holds, so it is safe)
    if (az > R(1)) {
        C lnmz = log(-z);
        return -dilog(C(R(1)) / z) - C(pi2_6) - lnmz * lnmz / C(R(2));
    }

    // |z| > 0.5: reflection formula — but only when |1-z| < |z|
    // (e.g. z=0.3+0.4i has |z|=0.5, |1-z|=0.806, so reflecting would recurse infinitely)
    if (az > R(0.5)) {
        C one_minus_z = C(R(1)) - z;
        R abs_omz = abs(one_minus_z);
        if (abs_omz < az) {
            return -dilog(one_minus_z) + C(pi2_6) - log(z) * log(one_minus_z);
        }
        // when |1-z| ≥ |z|, fall back to Taylor (slow convergence but safe)
    }

    // Taylor Σ z^n/n² — pass n² directly as Float(n*n) (exact within INT_MAX),
    // avoiding the trap where C(R(n) * R(n)) does exact*exact → falls to defaultPrecision
    C sum(R(0));
    C z_power = z;

    const int max_iter = IsSangiFloat<R> ? 2000 : 200;
    R eps = detail::getEpsilon<R>();
    if constexpr (IsSangiFloat<R>) {
        eps = detail::getEpsilon<R>(z.re.precision());
    }

    for (int n = 1; n < max_iter; n++) {
        // term = z_power / (n²) — Complex / scalar (re/n², im/n² separately)
        // R(int64_t(n) * n) is exact, but on the Complex/scalar path it mixes with the
        // high-precision z_power and min/max propagation preserves the precision of z
        const int64_t n_sq = static_cast<int64_t>(n) * static_cast<int64_t>(n);
        C term(z_power.re / R(n_sq), z_power.im / R(n_sq));
        sum = sum + term;

        if (n >= 5 && abs(term) < eps * abs(sum))
            break;

        z_power = z_power * z;
    }

    return sum;
}

// ================================================================
// Hyperbolic sine integral Shi(x) = ∫_0^x sinh(t)/t dt
// ================================================================
// Shi(x) = Σ_{k=0}^∞ x^{2k+1} / ((2k+1)·(2k+1)!)   ← Si(x) with the alternating sign removed
// odd function: Shi(-x) = -Shi(x)
// relation: Shi(x) = (Ei(x) - Ei(-x)) / 2, but it diverges as x → ∞ (order of Ei(x))
// for large x the series converges slowly, so switch to the Ei route

template<IsNativeFloat T>
[[nodiscard]] T sinhIntegral(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();
    }

    T sign = T(1);
    T ax = x;
    if (x < T(0)) { sign = T(-1); ax = -x; }

    // large ax: Ei route (Shi(x) = (Ei(x) - Ei(-x)) / 2)
    if (ax > T(20.0)) {
        // Ei(-x) is small and negligible (below 1e-9 for x>20)
        T result = (expint(ax) - expint(-ax)) / T(2);
        return sign * result;
    }

    // Taylor: term_k = x^{2k+1} / (2k+1)!, summand = term_k / (2k+1)
    // term_{k+1}/term_k = x² / ((2k+2)(2k+3))
    T x2 = ax * ax;
    T term = ax;        // x^1 / 1!
    T sum  = ax;        // k=0
    for (int k = 0; k < 200; ++k) {
        term *= x2 / (T(2 * k + 2) * T(2 * k + 3));
        T contribution = term / T(2 * k + 3);
        sum += contribution;
        if (k >= 3 && std::abs(contribution) <
            std::numeric_limits<T>::epsilon() * std::abs(sum))
            break;
    }
    return sign * sum;
}

// ================================================================
// Hyperbolic cosine integral Chi(x) = γ + ln|x| + ∫_0^x (cosh(t)-1)/t dt   (x > 0)
// ================================================================
// Chi(x) = γ + ln(x) + Σ_{k=1}^∞ x^{2k} / (2k·(2k)!)   ← Ci with the alternating sign removed
// defined for x > 0, Chi(0) = -∞, Chi(∞) = +∞
// for x < 0 the solution is complex (the real version returns NaN)
// relation: Chi(x) = (Ei(x) + Ei(-x)) / 2

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

    // large x: Ei route (Chi(x) = (Ei(x) + Ei(-x)) / 2)
    if (x > T(20.0)) {
        return (expint(x) + expint(-x)) / T(2);
    }

    // Chi(x) = γ + ln(x) + Σ_{k=1}^∞ x^{2k} / (2k·(2k)!)
    T sum = detail::euler_gamma<T> + std::log(x);
    T x2 = x * x;
    T term = T(1);    // k=0 initial (unused, kept as a pivot)
    for (int k = 1; k < 200; ++k) {
        term *= x2 / (T(2 * k - 1) * T(2 * k));
        T contribution = term / T(2 * k);
        sum += contribution;
        if (k >= 3 && std::abs(contribution) <
            std::numeric_limits<T>::epsilon() * std::abs(sum))
            break;
    }
    return sum;
}

} // namespace special
} // namespace sangi

#endif // SANGI_SPECIAL_EXPONENTIAL_INTEGRAL_HPP