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

// struve.hpp
//
// Struve function H_ν(z) and modified Struve function L_ν(z)
//
// Provided functions:
//   struveH(nu, z)   — Struve function H_ν(z)         (real order, real argument)
//   struveL(nu, z)   — modified Struve function L_ν(z) (real order, real argument)
//
// Supported types:
//   float, double, long double — direct implementation (undefined in std::)
//   Float                      — not yet supported (planned extension in a future release)
//
// Algorithm:
//   |z| <= cutoff (≈16): direct series from DLMF 11.2.1
//     H_ν(z) = (z/2)^{ν+1} Σ_{k=0}^∞ ((-1)^k / (Γ(k+3/2) Γ(ν+k+3/2))) (z/2)^{2k}
//     L_ν(z) = (z/2)^{ν+1} Σ_{k=0}^∞ (1     / (Γ(k+3/2) Γ(ν+k+3/2))) (z/2)^{2k}
//
//   |z| >  cutoff: asymptotic expansion from DLMF 11.6.1
//     H_ν(z) ~ Y_ν(z) + (1/π) Σ_{k=0}^N Γ(k+1/2) Γ(ν+1/2-k)^{-1} (z/2)^{ν-2k-1}
//     L_ν(z) ~ I_ν(z) - (1/π) Σ_{k=0}^N Γ(k+1/2) Γ(ν+1/2-k)^{-1} (z/2)^{ν-2k-1}
//   The asymptotic series diverges, so it is truncated once the absolute value of a term starts to increase.
//
// Constraints:
//   - For now z >= 0 only. Negative z is not supported (for integer ν the symmetry
//     H_ν(-z) = (-1)^{ν+1} H_ν(z) could be used, but for non-integer ν the solution is complex, so a separate API is needed).
//   - For half-integer ν, Γ(ν+1/2-k) may have a pole (excluded in the asymptotic series / early termination).

#ifndef SANGI_SPECIAL_STRUVE_HPP
#define SANGI_SPECIAL_STRUVE_HPP

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

namespace sangi {
namespace special {

namespace detail {

    // 1/π (constant at long double precision)
    template<class T>
    inline constexpr T struve_inv_pi() {
        return T(0.31830988618379067153776752674502872406891929148L);
    }

    // Direct series: switch between H_ν and L_ν via a sign flag
    // alt = true → H_ν (factor (-1)^k on each term), alt = false → L_ν (no sign alternation)
    template<IsNativeFloat T>
    [[nodiscard]] T struve_series(T nu, T z, bool alt) {
        if (z == T(0)) {
            // (z/2)^{ν+1} → 0 (ν > -1) or diverges (ν < -1)
            // By convention assume ν > -1: H_ν(0) = L_ν(0) = 0
            return T(0);
        }
        const T half_z = z / T(2);
        const T half_z_sq = half_z * half_z;
        const T common = std::pow(half_z, nu + T(1));

        // a_0 = 1 / (Γ(3/2) Γ(ν + 3/2))
        T a = T(1) / (std::tgamma(T(1.5)) * std::tgamma(nu + T(1.5)));
        T sum = a;

        constexpr int MAX_ITER = 2000;
        const T eps = std::numeric_limits<T>::epsilon();
        for (int k = 0; k < MAX_ITER; ++k) {
            // a_{k+1} = ±a_k * (z/2)^2 / ((k + 3/2)(ν + k + 3/2))
            const T denom = (T(k) + T(1.5)) * (nu + T(k) + T(1.5));
            a *= half_z_sq / denom;
            if (alt) a = -a;
            sum += a;
            if (std::abs(a) < std::abs(sum) * eps * T(0.5)) break;
        }
        return common * sum;
    }

    // Asymptotic series: the (z/2)^{ν-2k-1} series for large |z|
    // The return value is Σ_{k=0}^N Γ(k+1/2) Γ(ν+1/2-k)^{-1} (z/2)^{ν-2k-1}
    // Recurrence: term_{k+1} = term_k * (k + 1/2)(ν - 1/2 - k) / (z/2)^2
    template<IsNativeFloat T>
    [[nodiscard]] T struve_asymp_series(T nu, T z) {
        const T half_z = z / T(2);
        const T inv_half_z_sq = T(1) / (half_z * half_z);

        // term_0 = Γ(1/2)/Γ(ν+1/2) * (z/2)^{ν-1}
        // = √π / Γ(ν+1/2) * (z/2)^{ν-1}
        T term = std::pow(half_z, nu - T(1)) / std::tgamma(nu + T(0.5))
                 * std::tgamma(T(0.5));
        T sum = term;

        constexpr int MAX_TERMS = 50;
        T prev_abs = std::abs(term);
        for (int k = 0; k < MAX_TERMS; ++k) {
            // next term = current * (k+1/2)(ν-1/2-k) / (z/2)^2
            T factor = (T(k) + T(0.5)) * (nu - T(0.5) - T(k));
            term *= factor * inv_half_z_sq;
            T abs_term = std::abs(term);
            // Divergence detection: since this is an asymptotic series, stop once the term's absolute value starts to increase
            if (abs_term > prev_abs) break;
            sum += term;
            prev_abs = abs_term;
            if (abs_term < std::abs(sum) * std::numeric_limits<T>::epsilon()) break;
        }
        return sum;
    }

    // |z| switchover threshold: direct series vs asymptotic series
    // Switch to the asymptotic series just before the loss of significance from intermediate-term growth becomes noticeable.
    template<IsNativeFloat T>
    inline constexpr T struve_series_cutoff() {
        return T(16);
    }

} // namespace detail

// ================================================================
// Struve function H_ν(z)
// ================================================================

/// Native floating-point types: series + asymptotic expansion
template<IsNativeFloat T>
[[nodiscard]] T struveH(T nu, T z) {
    if (std::isnan(nu) || std::isnan(z)) {
        return std::numeric_limits<T>::quiet_NaN();
    }
    if (z < T(0)) {
        // For integer ν this could be handled via H_ν(-z) = (-1)^{ν+1} H_ν(z),
        // but for general ν the solution is complex, so it is not supported for now.
        return std::numeric_limits<T>::quiet_NaN();
    }
    if (z == T(0)) {
        // ν > -1: H_ν(0) = 0, ν <= -1: diverges
        if (nu > T(-1)) return T(0);
        return std::numeric_limits<T>::infinity();
    }

    if (z <= detail::struve_series_cutoff<T>()) {
        return detail::struve_series(nu, z, /*alt=*/true);
    }
    // Asymptotic: H_ν(z) ~ Y_ν(z) + (1/π) Σ
    const T y_part = std::cyl_neumann(nu, z);
    const T r_part = detail::struve_asymp_series(nu, z);
    return y_part + detail::struve_inv_pi<T>() * r_part;
}

/// Convenience overload for integer order
template<IsNativeFloat T>
[[nodiscard]] T struveH(int n, T z) {
    return struveH(static_cast<T>(n), z);
}

// ================================================================
// Modified Struve function L_ν(z)
// ================================================================

/// Native floating-point types: series + asymptotic expansion
/// Relation: L_ν(z) = -i e^{-iνπ/2} H_ν(iz), but for real arguments a different asymptotic form is used.
template<IsNativeFloat T>
[[nodiscard]] T struveL(T nu, T z) {
    if (std::isnan(nu) || std::isnan(z)) {
        return std::numeric_limits<T>::quiet_NaN();
    }
    if (z < T(0)) {
        return std::numeric_limits<T>::quiet_NaN();
    }
    if (z == T(0)) {
        if (nu > T(-1)) return T(0);
        return std::numeric_limits<T>::infinity();
    }

    if (z <= detail::struve_series_cutoff<T>()) {
        return detail::struve_series(nu, z, /*alt=*/false);
    }
    // Asymptotic: L_ν(z) ~ I_ν(z) - (1/π) Σ
    const T i_part = std::cyl_bessel_i(nu, z);
    const T r_part = detail::struve_asymp_series(nu, z);
    return i_part - detail::struve_inv_pi<T>() * r_part;
}

/// Convenience overload for integer order
template<IsNativeFloat T>
[[nodiscard]] T struveL(int n, T z) {
    return struveL(static_cast<T>(n), z);
}

} // namespace special
} // namespace sangi

#endif // SANGI_SPECIAL_STRUVE_HPP