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

// hypergeometric_pq.hpp
// General pFq hypergeometric function and Meijer G function (partial, stocking phase initial drop)
//
// The existing hypergeometric.hpp provides ₀F₁, ₁F₁, ₂F₁ individually. This header:
//
//   pFq(a, b, z)            — general ₚFq Σ_{k} [Π(a_i)_k / Π(b_j)_k] z^k / k!
//   meijerG(m, n, a, b, z)  — Meijer G function G^{m,n}_{p,q}(z | a; b) (residue representation)
//
// Convergence conditions for pFq:
//   p ≤ q          : converges for any z (entire function)
//   p = q + 1      : converges for |z| < 1, conditional at z=1
//   p > q + 1      : formal series (asymptotic) — not handled by this implementation
//
// Current state of Meijer G (stocking phase initial):
//   - n = 0 (a parameters do not appear in the numerator) and the m ≤ q constraint
//   - series representation via the residue theorem only when all b_j are simple poles (distinct)
//   - resonance (integer differences between b_j) is not handled in this drop (returns NaN)
//   - the general (m, n, p, q) case is continued development in a future release
//
// API design:
//   - p, q are passed as std::span<const T> (C++20)
//   - polynomial termination (a_i a non-positive integer) is detected automatically
//   - non-converging cases (p > q+1, |z| ≥ 1 with p = q+1) return NaN
//
// References:
//   - DLMF §16.2 (pFq), §16.17 (Meijer G)
//   - Bühring 1992 (pFq analytic continuation) — this drop is direct Taylor summation only
//   - mpmath.hyper, mpmath.meijerg

#ifndef SANGI_SPECIAL_HYPERGEOMETRIC_PQ_HPP
#define SANGI_SPECIAL_HYPERGEOMETRIC_PQ_HPP

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

namespace sangi {
namespace special {

namespace detail_pq {

    // Detect polynomial termination: whether any a_i is a non-positive integer
    template<IsNativeFloat T>
    [[nodiscard]] int polynomial_terms(std::span<const T> a) {
        int min_terms = -1;
        for (T ai : a) {
            if (ai <= T(0) && ai == std::floor(ai)) {
                int terms = static_cast<int>(-ai) + 1;
                if (min_terms < 0 || terms < min_terms) min_terms = terms;
            }
        }
        return min_terms;  // -1 = not a polynomial
    }

    // NaN when b_j is a non-positive integer (pole) and the corresponding a_i does not terminate earlier
    template<IsNativeFloat T>
    [[nodiscard]] bool b_is_pole(std::span<const T> a, std::span<const T> b) {
        // For each b_j, if b_j is a non-positive integer → pole candidate
        for (T bj : b) {
            if (bj <= T(0) && bj == std::floor(bj)) {
                int bj_int = static_cast<int>(bj);  // ≤ 0
                // if any a_i is a larger (or equal) non-positive integer → the polynomial stops first
                bool saved = false;
                for (T ai : a) {
                    if (ai <= T(0) && ai == std::floor(ai)) {
                        int ai_int = static_cast<int>(ai);
                        if (ai_int >= bj_int) { saved = true; break; }
                    }
                }
                if (!saved) return true;
            }
        }
        return false;
    }

}  // namespace detail_pq

// ================================================================
// general ₚFq hypergeometric function
// ================================================================
//
// pFq(a; b; z) = Σ_{k=0}^∞ [Π(a_i)_k / Π(b_j)_k] · z^k / k!
//
// Recurrence:
//   term_{k+1}/term_k = z · Π(a_i + k) / [(k+1) · Π(b_j + k)]
//
// Arguments:
//   a: numerator parameters a_1, ..., a_p
//   b: denominator parameters b_1, ..., b_q
//   z: argument
//
// Return value:
//   p ≤ q  : converges for any z
//   p=q+1  : converges for |z| < 1, |z|=1 is parameter-dependent
//   p>q+1  : NaN
//
template<IsNativeFloat T>
[[nodiscard]] T pFq(std::span<const T> a, std::span<const T> b, T z) {
    if (std::isnan(z)) return std::numeric_limits<T>::quiet_NaN();
    for (T ai : a) if (std::isnan(ai)) return std::numeric_limits<T>::quiet_NaN();
    for (T bj : b) if (std::isnan(bj)) return std::numeric_limits<T>::quiet_NaN();

    const int p = static_cast<int>(a.size());
    const int q = static_cast<int>(b.size());

    // divergent case: p > q+1
    if (p > q + 1) return std::numeric_limits<T>::quiet_NaN();

    // b has a pole and is not rescued by a polynomial
    if (detail_pq::b_is_pole(a, b)) return std::numeric_limits<T>::quiet_NaN();

    if (z == T(0)) return T(1);

    // detect polynomial termination
    int poly_terms = detail_pq::polynomial_terms(a);
    bool is_poly = (poly_terms > 0);

    // p = q+1 with |z| ≥ 1 and not a polynomial → divergent (conditional convergence not handled here)
    if (!is_poly && p == q + 1 && std::abs(z) >= T(1))
        return std::numeric_limits<T>::quiet_NaN();

    int max_iter = is_poly ? poly_terms : 5000;

    T sum = T(1);
    T term = T(1);
    const T eps = std::numeric_limits<T>::epsilon();

    for (int k = 0; k < max_iter; ++k) {
        // Pochhammer factors
        T num = T(1), den = T(k + 1);
        bool zero_term = false;
        for (T ai : a) {
            T factor = ai + T(k);
            if (factor == T(0)) { zero_term = true; break; }
            num *= factor;
        }
        if (zero_term) break;
        for (T bj : b) {
            T factor = bj + T(k);
            den *= factor;
        }
        term = term * z * num / den;
        sum += term;

        if (!is_poly && k >= 5 && std::abs(term) < eps * std::abs(sum))
            break;
    }
    return sum;
}

// initializer-list version (convenience)
template<IsNativeFloat T>
[[nodiscard]] T pFq(std::initializer_list<T> a, std::initializer_list<T> b, T z) {
    std::vector<T> av(a), bv(b);
    return pFq(std::span<const T>(av), std::span<const T>(bv), z);
}

// ================================================================
// Float arbitrary-precision ₚFq
// ================================================================
//
// Same strategy as native. Computes while truncating the Pochhammer factor each time.

[[nodiscard]] inline Float pFq(std::span<const Float> a, std::span<const Float> b,
                                const Float& z, int precision) {
    if (z.isNaN()) return Float::nan();
    for (const auto& ai : a) if (ai.isNaN()) return Float::nan();
    for (const auto& bj : b) if (bj.isNaN()) return Float::nan();

    const int p = static_cast<int>(a.size());
    const int q = static_cast<int>(b.size());
    if (p > q + 1) return Float::nan();

    int wp = precision + 30;

    if (z.isZero()) return Float::one(precision);

    // detect polynomial termination (a_i a non-positive integer)
    int poly_terms = -1;
    for (const auto& ai : a) {
        if (ai.isInteger() && (ai.isZero() || ai.isNegative())) {
            // -ai is a non-negative integer
            int terms = static_cast<int>(-ai.toDouble()) + 1;
            if (poly_terms < 0 || terms < poly_terms) poly_terms = terms;
        }
    }
    bool is_poly = (poly_terms > 0);

    Float zw = z; zw.truncateToApprox(wp);

    // p = q+1 with |z| ≥ 1 and not a polynomial → divergent
    if (!is_poly && p == q + 1) {
        Float az = abs(zw);
        if (az >= Float::one(wp)) return Float::nan();
    }

    int max_iter = is_poly ? poly_terms : (wp * 6 + 500);

    Float sum = Float::one(wp);
    Float term = Float::one(wp);
    int small_count = 0;

    for (int k = 0; k < max_iter; ++k) {
        Float num = Float::one(wp);
        Float den = Float(k + 1);
        bool zero_term = false;
        for (const auto& ai : a) {
            Float factor = ai + Float(k);
            if (factor.isZero()) { zero_term = true; break; }
            num = num * factor;
            num.truncateToApprox(wp);
        }
        if (zero_term) break;
        for (const auto& bj : b) {
            Float factor = bj + Float(k);
            den = den * factor;
            den.truncateToApprox(wp);
        }
        term = term * zw * num / den;
        term.truncateToApprox(wp);
        sum = sum + term;
        sum.truncateToApprox(wp);

        if (!is_poly && k >= 5) {
            Float at = abs(term);
            Float as = abs(sum);
            Float thresh = as * Float::epsilon(wp);
            if (at < thresh) {
                if (++small_count >= 3) break;
            } else {
                small_count = 0;
            }
        }
    }
    sum.setPrecision(precision);
    return sum;
}

[[nodiscard]] inline Float pFq(std::initializer_list<Float> a,
                                std::initializer_list<Float> b,
                                const Float& z, int precision) {
    std::vector<Float> av(a), bv(b);
    return pFq(std::span<const Float>(av), std::span<const Float>(bv), z, precision);
}

// ================================================================
// Meijer G function (residue representation, partial)
// ================================================================
//
// G^{m,n}_{p,q}(z | a_1...a_p; b_1...b_q)
//   = (1/2πi) ∫_L Π_{j=1}^m Γ(b_j-s) · Π_{j=1}^n Γ(1-a_j+s)
//                / [Π_{j=m+1}^q Γ(1-b_j+s) · Π_{j=n+1}^p Γ(a_j-s)]
//                · z^s ds
//
// Constraints of the partial implementation provided here:
//   - n = 0 : the a parameters are not split into upper/lower halves but simply collected
//             on the "denominator side" (= best tested when a is an empty vector)
//   - m ≤ q : the poles of Γ(b_j-s) lie in the left half-plane
//   - all b_j distinct (also distinct mod ℤ) : no resonance
//
// In this setting G is a sum of residues at the simple poles at b_h:
//   G(z) = Σ_{h=1}^m Π'_{j=1}^m Γ(b_j-b_h) · Π_{j=1}^n Γ(1-a_j+b_h)
//          / [Π_{j=m+1}^q Γ(1-b_j+b_h) · Π_{j=n+1}^p Γ(a_j-b_h)]
//          · z^{b_h} · ₚFq-1(...; z·(-1)^{p-m-n})
//
// For simplicity, the n=0 formula (based on DLMF 16.17.2):
//   G^{m,0}_{p,q}(z | a; b)
//     = Σ_{h=1}^m  [Π_{j≠h, j≤m} Γ(b_j-b_h)] / [Π_{j>m} Γ(1-b_j+b_h) · Π_j Γ(a_j-b_h)]
//                · z^{b_h} · pFq^*
//
// The argument of pFq^* is determined by the shift of the b parameters, but this drop
// only computes the subfamily with empty a (n=0, p=0) (= covers the Helmholtz-Green
// function family).
//
// ⚠️ More general cases are staged extensions in a future release:
//      - n > 0
//      - resonance handling (logarithmic terms)
//      - p > 0 with non-empty a (interaction between the b and a parameters)

template<IsNativeFloat T>
[[nodiscard]] T meijerG_simple_b(int m, std::span<const T> b, T z) {
    // G^{m,0}_{0,q}(z | -; b) — a empty, n=0
    // = Σ_{h=1}^m [Π_{j≠h, j≤m} Γ(b_j - b_h)] / [Π_{j=m+1}^q Γ(1 - b_j + b_h)]
    //          · z^{b_h} · ₀F_{q-1}(; b_*; z·(-1)^{q-m})
    // where b_* = the set of shifted parameters obtained from b by removing b_h
    if (std::isnan(z)) return std::numeric_limits<T>::quiet_NaN();
    if (z <= T(0)) return std::numeric_limits<T>::quiet_NaN();  // z^{b_h} could become complex

    const int q = static_cast<int>(b.size());
    if (m < 1 || m > q) return std::numeric_limits<T>::quiet_NaN();

    // resonance check: NaN in this implementation if b_j - b_h is an integer (j ≠ h, both ≤ m)
    for (int i = 0; i < m; ++i) {
        for (int j = i + 1; j < m; ++j) {
            T diff = b[i] - b[j];
            if (diff == std::floor(diff))
                return std::numeric_limits<T>::quiet_NaN();
        }
    }

    T result = T(0);
    for (int h = 0; h < m; ++h) {
        // coefficient: prefactor
        T prefac = T(1);
        for (int j = 0; j < m; ++j) {
            if (j == h) continue;
            prefac *= std::tgamma(b[j] - b[h]);
        }
        for (int j = m; j < q; ++j) {
            prefac /= std::tgamma(T(1) - b[j] + b[h]);
        }
        // pFq^* part (₀F_{q-1}): arguments b_*= 1 + b_h - b_j (j ≠ h)
        std::vector<T> b_star;
        b_star.reserve(q - 1);
        for (int j = 0; j < q; ++j) {
            if (j == h) continue;
            b_star.push_back(T(1) + b[h] - b[j]);
        }
        // general sign factor: (-1)^{p-m-n}. The partial implementation in this header has p=n=0, so (-1)^m.
        T sign = (m % 2 == 0) ? T(1) : T(-1);
        // ₀F_{q-1}(; b_*; sign · z) has p=0
        std::vector<T> a_empty;
        T hyp = pFq<T>(std::span<const T>(a_empty),
                       std::span<const T>(b_star),
                       sign * z);

        // z^{b_h}
        T z_pow = std::pow(z, b[h]);
        result += prefac * z_pow * hyp;
    }
    return result;
}

// more general stub (currently delegates only to n=0, empty a)
template<IsNativeFloat T>
[[nodiscard]] T meijerG(int m, int n,
                         std::span<const T> a, std::span<const T> b, T z) {
    if (n != 0 || !a.empty()) {
        // extension in a future release
        return std::numeric_limits<T>::quiet_NaN();
    }
    return meijerG_simple_b(m, b, z);
}

template<IsNativeFloat T>
[[nodiscard]] T meijerG(int m, int n,
                         std::initializer_list<T> a, std::initializer_list<T> b, T z) {
    std::vector<T> av(a), bv(b);
    return meijerG<T>(m, n, std::span<const T>(av), std::span<const T>(bv), z);
}

} // namespace special
} // namespace sangi

#endif // SANGI_SPECIAL_HYPERGEOMETRIC_PQ_HPP