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

// hypergeometric.hpp
// Hypergeometric function templates
//
// Provided functions:
//   confHyperg(a, b, z)      — confluent hypergeometric function ₁F₁(a; b; z) (Kummer M)
//   hyperg(a, b, c, z)       — Gauss hypergeometric function ₂F₁(a, b; c; z)
//   hyperg0F1(b, z)          — confluent limit function ₀F₁(; b; z)
//
// Supported types:
//   float, double, long double — Taylor series
//   Float                      — delegated to the sangi:: implementation

#ifndef SANGI_SPECIAL_HYPERGEOMETRIC_HPP
#define SANGI_SPECIAL_HYPERGEOMETRIC_HPP

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

namespace sangi {
namespace special {

// ================================================================
// Confluent hypergeometric function ₁F₁(a; b; z) (Kummer M)
// ================================================================
// M(a, b, z) = Σ_{k=0}^∞ (a)_k z^k / ((b)_k k!)
// term_{k+1}/term_k = (a+k)·z / ((b+k)·(k+1))
// b cannot be a non-positive integer (pole)

template<IsNativeFloat T>
[[nodiscard]] T confHyperg(T a, T b, T z) {
    if (std::isnan(a) || std::isnan(b) || std::isnan(z))
        return std::numeric_limits<T>::quiet_NaN();
    // b is a non-positive integer → pole
    if (b <= T(0) && b == std::floor(b))
        return std::numeric_limits<T>::quiet_NaN();
    if (a == T(0)) return T(1);
    if (z == T(0)) return T(1);
    if (a == b) return std::exp(z);

    // Kummer transformation: for z < 0, M(a,b,z) = e^z · M(b-a, b, -z)
    if (z < T(0)) {
        return std::exp(z) * confHyperg(b - a, b, -z);
    }

    // a is a non-positive integer → polynomial (finite sum); otherwise cap the convergence truncation at 1000 terms
    // 1000 terms: for double, this is fast enough to converge when |z| ≤ 100 or so. Beware that it diverges for |z| >> 100
    bool is_poly = (a < T(0) && a == std::floor(a));
    int max_terms = is_poly ? static_cast<int>(-a) + 1 : 1000;

    T sum = T(1);
    T term = T(1);

    for (int k = 0; k < max_terms; k++) {
        term *= (a + T(k)) * z / ((b + T(k)) * T(k + 1));
        sum += term;

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

    return sum;
}

// Float type: use sangi::confHyperg(a, b, z, precision) directly.

// ================================================================
// Gauss hypergeometric function ₂F₁(a, b; c; z)
// ================================================================
// ₂F₁(a, b; c; z) = Σ_{k=0}^∞ (a)_k (b)_k z^k / ((c)_k k!)
// Converges for |z| < 1. If a or b is a non-positive integer it is a polynomial (finite sum for all z).
// If c is a non-positive integer it is a pole (unless a or b is a larger negative integer).

template<IsNativeFloat T>
[[nodiscard]] T hyperg(T a, T b, T c, T z) {
    if (std::isnan(a) || std::isnan(b) || std::isnan(c) || std::isnan(z))
        return std::numeric_limits<T>::quiet_NaN();
    if (z == T(0)) return T(1);
    if (a == T(0) || b == T(0)) return T(1);

    // c is a non-positive integer → a pole unless a or b terminates first
    if (c <= T(0) && c == std::floor(c)) {
        int ci = static_cast<int>(c);
        bool a_terminates = (a <= T(0) && a == std::floor(a) && static_cast<int>(a) >= ci);
        bool b_terminates = (b <= T(0) && b == std::floor(b) && static_cast<int>(b) >= ci);
        if (!a_terminates && !b_terminates)
            return std::numeric_limits<T>::quiet_NaN();
    }

    // Detection of the polynomial case
    bool is_poly = false;
    int poly_terms = 0;
    if (a <= T(0) && a == std::floor(a)) {
        is_poly = true;
        poly_terms = static_cast<int>(-a) + 1;
    } else if (b <= T(0) && b == std::floor(b)) {
        is_poly = true;
        poly_terms = static_cast<int>(-b) + 1;
    }

    T z_abs = std::abs(z);

    // |z| >= 1 and non-polynomial → diverges
    if (!is_poly && z_abs >= T(1))
        return std::numeric_limits<T>::quiet_NaN();

    // Pfaff transformation: z < -0.5 → (1-z)^{-a} · ₂F₁(a, c-b; c; z/(z-1))
    if (!is_poly && z < T(-0.5)) {
        T one_minus_z = T(1) - z;
        T z_mapped = z / (z - T(1));
        return std::pow(one_minus_z, -a) * hyperg(a, c - b, c, z_mapped);
    }

    // Taylor series
    T sum = T(1);
    T term = T(1);
    int max_iter = is_poly ? poly_terms : 1000;

    for (int k = 0; k < max_iter; k++) {
        // Pochhammer (a)_k or (b)_k is zero → all subsequent terms are 0
        if (a + T(k) == T(0) || b + T(k) == T(0)) break;

        term *= (a + T(k)) * (b + T(k)) * z / ((c + T(k)) * T(k + 1));
        sum += term;

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

    return sum;
}

// Float type: use sangi::hyperg(a, b, c, z, precision) directly.

// ================================================================
// Confluent limit function ₀F₁(; b; z)
// ================================================================
// ₀F₁(; b; z) = Σ_{k=0}^∞ z^k / ((b)_k k!)
// Converges for all z. b cannot be a non-positive integer.
// Related: ₀F₁(;1;-z²/4) = J₀(z), ₀F₁(;1/2;-z²/4) = cos(z)

template<IsNativeFloat T>
[[nodiscard]] T hyperg0F1(T b, T z) {
    if (std::isnan(b) || std::isnan(z))
        return std::numeric_limits<T>::quiet_NaN();
    if (b <= T(0) && b == std::floor(b))
        return std::numeric_limits<T>::quiet_NaN();

    T sum = T(1);
    T term = T(1);

    for (int k = 0; k < 500; k++) {
        term *= z / ((b + T(k)) * T(k + 1));
        sum += term;

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

    return sum;
}

// Float type: use sangi::hyperg0F1(b, z, precision) directly.

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

// ----------------------------------------------------------------
// ₁F₁(a; b; z) — Taylor series
// ----------------------------------------------------------------

template<SpecialFunctionScalar R>
[[nodiscard]] Complex<R> confHyperg(Complex<R> a, Complex<R> b, Complex<R> z) {
    using C = Complex<R>;
    if (z == C(R(0))) return C(R(1));

    C sum(R(1));
    C term(R(1));

    const int max_iter = IsSangiFloat<R> ? 5000 : 1000;
    R eps = detail::getEpsilon<R>();
    if constexpr (IsSangiFloat<R>) {
        eps = detail::getEpsilon<R>(z.re.precision());  // converge using the argument's working precision (avoid the default ~55 digits)
    }

    for (int k = 0; k < max_iter; k++) {
        term = term * (a + C(R(k))) * z / ((b + C(R(k))) * C(R(k + 1)));
        sum = sum + term;

        if (k >= 3 && abs(term) < eps * abs(sum))
            break;
    }

    return sum;
}

// ----------------------------------------------------------------
// ₂F₁(a, b; c; z) — Taylor series (|z| < 1)
// ----------------------------------------------------------------

template<SpecialFunctionScalar R>
[[nodiscard]] Complex<R> hyperg(Complex<R> a, Complex<R> b, Complex<R> c, Complex<R> z) {
    using C = Complex<R>;
    if (z == C(R(0))) return C(R(1));

    // |z| >= 1 → diverges (the polynomial case is not handled)
    if (abs(z) >= R(1))
        return C(detail::getNaN<R>());

    C sum(R(1));
    C term(R(1));

    const int max_iter = IsSangiFloat<R> ? 5000 : 1000;
    R eps = detail::getEpsilon<R>();
    if constexpr (IsSangiFloat<R>) {
        eps = detail::getEpsilon<R>(z.re.precision());  // converge using the argument's working precision (avoid the default ~55 digits)
    }

    for (int k = 0; k < max_iter; k++) {
        term = term * (a + C(R(k))) * (b + C(R(k))) * z
             / ((c + C(R(k))) * C(R(k + 1)));
        sum = sum + term;

        if (k >= 3 && abs(term) < eps * abs(sum))
            break;
    }

    return sum;
}

// ----------------------------------------------------------------
// ₀F₁(; b; z) — Taylor series
// ----------------------------------------------------------------

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

    C sum(R(1));
    C term(R(1));

    const int max_iter = IsSangiFloat<R> ? 2000 : 500;
    R eps = detail::getEpsilon<R>();
    if constexpr (IsSangiFloat<R>) {
        eps = detail::getEpsilon<R>(z.re.precision());  // converge using the argument's working precision (avoid the default ~55 digits)
    }

    for (int k = 0; k < max_iter; k++) {
        term = term * z / ((b + C(R(k))) * C(R(k + 1)));
        sum = sum + term;

        if (k >= 3 && abs(term) < eps * abs(sum))
            break;
    }

    return sum;
}

} // namespace special
} // namespace sangi

#endif // SANGI_SPECIAL_HYPERGEOMETRIC_HPP