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

// bessel.hpp
// Template wrappers for the Bessel function family
//
// Provided functions:
//   besselJ(nu, x), besselY(nu, x)       — ordinary Bessel of the 1st/2nd kind (real order)
//   besselI(nu, x), besselK(nu, x)       — modified Bessel of the 1st/2nd kind (real order)
//   sphericalBesselJ(n, x), sphericalBesselY(n, x) — spherical Bessel of the 1st/2nd kind
//   besselJPrime(nu, x), besselYPrime(nu, x)       — derivatives (recurrence)
//   besselIPrime(nu, x), besselKPrime(nu, x)       — modified Bessel derivatives
//   hankelH1(n, z), hankelH2(n, z)       — Hankel function of the 1st/2nd kind (DLMF 10.4.3)
//   hankelH1Prime(n, z), hankelH2Prime(n, z) — Hankel function derivatives
//   sphericalHankelH1(n, z), sphericalHankelH2(n, z) — spherical Hankel functions
//   besselJ/Y/I/K(n, z)                  — Complex<double> support
//   besselJ/Y/I/K(n, z, precision)       — Complex<Float> support
//
// Supported types:
//   float, double, long double — delegated to C++17 std::
//   Float                      — delegated to the sangi:: implementation
//   Complex<double>            — Taylor series / A&S formulas
//   Complex<Float>             — arbitrary-precision Taylor series / A&S formulas

#ifndef SANGI_SPECIAL_BESSEL_HPP
#define SANGI_SPECIAL_BESSEL_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 <limits>

namespace sangi {
namespace special {

// ================================================================
// Ordinary Bessel function of the 1st kind J_ν(x)
// ================================================================

/// Native floating-point types: delegate to std::cyl_bessel_j
template<IsNativeFloat T>
[[nodiscard]] T besselJ(T nu, T x) {
    if (std::isnan(nu) || std::isnan(x)) return std::numeric_limits<T>::quiet_NaN();
    return std::cyl_bessel_j(nu, x);
}

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

// Float type: use sangi::besselJ(n, x, precision) directly.

// ================================================================
// Ordinary Bessel function of the 2nd kind Y_ν(x) (Neumann function)
// ================================================================

/// Native floating-point types: delegate to std::cyl_neumann
template<IsNativeFloat T>
[[nodiscard]] T besselY(T nu, T x) {
    if (std::isnan(nu) || std::isnan(x)) return std::numeric_limits<T>::quiet_NaN();
    return std::cyl_neumann(nu, x);
}

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

// Float type: use sangi::besselY(n, x, precision) directly.

// ================================================================
// Modified Bessel function of the 1st kind I_ν(x)
// ================================================================

/// Native floating-point types: delegate to std::cyl_bessel_i
template<IsNativeFloat T>
[[nodiscard]] T besselI(T nu, T x) {
    if (std::isnan(nu) || std::isnan(x)) return std::numeric_limits<T>::quiet_NaN();
    return std::cyl_bessel_i(nu, x);
}

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

// Float type: use sangi::besselI(n, x, precision) directly.

// ================================================================
// Modified Bessel function of the 2nd kind K_ν(x)
// ================================================================

/// Native floating-point types: delegate to std::cyl_bessel_k
template<IsNativeFloat T>
[[nodiscard]] T besselK(T nu, T x) {
    if (std::isnan(nu) || std::isnan(x)) return std::numeric_limits<T>::quiet_NaN();
    if (x <= T(0)) return std::numeric_limits<T>::quiet_NaN();
    return std::cyl_bessel_k(nu, x);
}

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

// Float type: use sangi::besselK(n, x, precision) directly.

// ================================================================
// Spherical Bessel function of the 1st kind j_n(x)
// j_n(x) = √(π/(2x)) · J_{n+1/2}(x)
// ================================================================

/// Native floating-point types: delegate to std::sph_bessel
template<IsNativeFloat T>
[[nodiscard]] T sphericalBesselJ(unsigned int n, T x) {
    if (std::isnan(x)) return std::numeric_limits<T>::quiet_NaN();
    return std::sph_bessel(n, x);
}

// Float type: use sangi::sphericalBesselJ(n, x, precision) directly.

// ================================================================
// Spherical Bessel function of the 2nd kind y_n(x) (spherical Neumann function)
// y_n(x) = √(π/(2x)) · Y_{n+1/2}(x)
// ================================================================

/// Native floating-point types: delegate to std::sph_neumann
template<IsNativeFloat T>
[[nodiscard]] T sphericalBesselY(unsigned int n, T x) {
    if (std::isnan(x)) return std::numeric_limits<T>::quiet_NaN();
    return std::sph_neumann(n, x);
}

// Float type: use sangi::sphericalBesselY(n, x, precision) directly.

// ================================================================
// Derivatives (via recurrence relations)
// DLMF 10.6.1: J'_ν(x) = (J_{ν-1}(x) - J_{ν+1}(x)) / 2
// DLMF 10.6.3: Y'_ν(x) = (Y_{ν-1}(x) - Y_{ν+1}(x)) / 2
// DLMF 10.29.1: I'_ν(x) = (I_{ν-1}(x) + I_{ν+1}(x)) / 2
// DLMF 10.29.3: K'_ν(x) = -(K_{ν-1}(x) + K_{ν+1}(x)) / 2
// ================================================================

template<IsNativeFloat T>
[[nodiscard]] T besselJPrime(T nu, T x) {
    return (besselJ(nu - T(1), x) - besselJ(nu + T(1), x)) / T(2);
}

template<IsNativeFloat T>
[[nodiscard]] T besselYPrime(T nu, T x) {
    return (besselY(nu - T(1), x) - besselY(nu + T(1), x)) / T(2);
}

template<IsNativeFloat T>
[[nodiscard]] T besselIPrime(T nu, T x) {
    return (besselI(nu - T(1), x) + besselI(nu + T(1), x)) / T(2);
}

template<IsNativeFloat T>
[[nodiscard]] T besselKPrime(T nu, T x) {
    return -(besselK(nu - T(1), x) + besselK(nu + T(1), x)) / T(2);
}

// ================================================================
// Complex Bessel functions (Complex<double> and Complex<Float>)
// ================================================================
// J_n(z): Taylor series  (z/2)^n Σ (-1)^m (z/2)^{2m} / (m!(m+n)!)
// I_n(z): Taylor series  (z/2)^n Σ (z/2)^{2m} / (m!(m+n)!)
// Y_n(z): A&S 9.1.11  (2/π)J_n·ln(z/2) - (1/π)Σ... - (1/π)Σ...
// K_n(z): A&S 9.6.11  (-1)^{n+1}·I_n·ln(z/2) + ...
// ================================================================

namespace detail {

/// J_n(z) — Taylor series (integer order, complex argument)
template<typename R>
[[nodiscard]] Complex<R> besselJ_complex(int n, const Complex<R>& z,
                                          R eps, int max_iter) {
    using C = Complex<R>;
    int abs_n = (n >= 0) ? n : -n;

    if (sangi::abs(z) == R(0))
        return (abs_n == 0) ? C(R(1)) : C(R(0));

    C z_half = z / R(2);
    C z_half_sq = z_half * z_half;
    C neg_z_half_sq = -z_half_sq;

    // (z/2)^n
    C z_half_n(R(1));
    for (int i = 0; i < abs_n; i++) z_half_n = z_half_n * z_half;

    // term_0 = 1/n!
    R fact_inv(1);
    for (int i = 1; i <= abs_n; i++) fact_inv = fact_inv / R(i);
    C term(fact_inv);
    C sum = term;

    for (int m = 1; m < max_iter; m++) {
        // term_m = term_{m-1} · (-(z/2)²) / (m·(m+n))
        term = term * neg_z_half_sq / R(static_cast<int64_t>(m) * (m + abs_n));
        sum = sum + term;
        if (m >= 3 && sangi::abs(term) < eps * sangi::abs(sum)) break;
    }

    C result = z_half_n * sum;
    if (n < 0 && abs_n % 2 != 0) result = -result;
    return result;
}

/// I_n(z) — Taylor series (integer order, complex argument)
template<typename R>
[[nodiscard]] Complex<R> besselI_complex(int n, const Complex<R>& z,
                                          R eps, int max_iter) {
    using C = Complex<R>;
    int abs_n = (n >= 0) ? n : -n;

    if (sangi::abs(z) == R(0))
        return (abs_n == 0) ? C(R(1)) : C(R(0));

    C z_half = z / R(2);
    C z_half_sq = z_half * z_half;

    // (z/2)^n
    C z_half_n(R(1));
    for (int i = 0; i < abs_n; i++) z_half_n = z_half_n * z_half;

    // term_0 = 1/n!
    R fact_inv(1);
    for (int i = 1; i <= abs_n; i++) fact_inv = fact_inv / R(i);
    C term(fact_inv);
    C sum = term;

    for (int m = 1; m < max_iter; m++) {
        // term_m = term_{m-1} · (z/2)² / (m·(m+n))  — no sign flip
        term = term * z_half_sq / R(static_cast<int64_t>(m) * (m + abs_n));
        sum = sum + term;
        if (m >= 3 && sangi::abs(term) < eps * sangi::abs(sum)) break;
    }

    C result = z_half_n * sum;
    // I_{-n}(z) = I_n(z) for integer n — unchanged
    return result;
}

/// Y_n(z) — A&S 9.1.11 (integer order, complex argument)
/// Y_n(z) = (2/π)·J_n(z)·ln(z/2)
///        - (1/π)·(z/2)^{-n}·Σ_{k=0}^{n-1} (n-k-1)!/k! · (z²/4)^k
///        - (1/π)·(z/2)^n · Σ_{k=0}^{∞} (-1)^k·(H_k+H_{k+n})/(k!·(k+n)!) · (z²/4)^k
template<typename R>
[[nodiscard]] Complex<R> besselY_complex(int n, const Complex<R>& z,
                                          R eps, int max_iter,
                                          R pi_val, R euler_gamma) {
    using C = Complex<R>;
    int abs_n = (n >= 0) ? n : -n;

    C z_half = z / R(2);
    C z_half_sq = z_half * z_half;
    C ln_z_half = sangi::log(z_half);
    R two_over_pi = R(2) / pi_val;
    R one_over_pi = R(1) / pi_val;

    // Compute J_n(z)
    C jn = besselJ_complex(abs_n, z, eps, max_iter);

    // Part A: (2/π) · J_n(z) · (γ + ln(z/2))  [DLMF 10.8.1]
    C result = C(two_over_pi) * jn * (C(euler_gamma) + ln_z_half);

    // Part B: -(1/π) · (z/2)^{-n} · Σ_{k=0}^{n-1} (n-k-1)!/k! · (z²/4)^k
    if (abs_n > 0) {
        C z_half_neg_n(R(1));
        for (int i = 0; i < abs_n; i++) z_half_neg_n = z_half_neg_n / z_half;

        // k=0: (n-1)!
        R bk_coeff(1);
        for (int i = 1; i < abs_n; i++) bk_coeff = bk_coeff * R(i);
        C sumB(bk_coeff);

        for (int k = 1; k < abs_n; k++) {
            // recurrence: coeff_k = coeff_{k-1} · (z²/4) / (k·(n-k))
            bk_coeff = bk_coeff / R(static_cast<int64_t>(k) * (abs_n - k));
            sumB = sumB + C(bk_coeff) * sangi::pow(z_half_sq, k);
        }

        result = result - C(one_over_pi) * z_half_neg_n * sumB;
    }

    // Part C: -(1/π) · (z/2)^n · Σ_{k=0}^{∞} (-1)^k (H_k+H_{k+n})/(k!·(k+n)!) · (z²/4)^k
    C z_half_n(R(1));
    for (int i = 0; i < abs_n; i++) z_half_n = z_half_n * z_half;

    // Avoid exact/exact poison: for R=Float, divide harmonic numbers / factorial
    // inverses by a working-precision unit. R(1)/R(j) with both operands exact
    // falls back to default precision and tops out at ~57 digits
    // (AUDIT_FLOAT_UNIT_MIXING root cause #1). For R=double, unit=1 is the identity.
    R unit(1);
    if constexpr (std::is_same_v<R, sangi::Float>) {
        unit.setEffectiveBits(sangi::Float::precisionToBits(pi_val.precision()));
    }
    R H_k(0);       // H_0 = 0
    R H_kn(0);      // H_n = Σ_{j=1}^{n} 1/j
    for (int j = 1; j <= abs_n; j++) H_kn = H_kn + unit / R(j);

    // k=0: 1/n!
    R ck_coeff = unit;
    for (int i = 1; i <= abs_n; i++) ck_coeff = ck_coeff / R(i);

    C neg_z_half_sq = -z_half_sq;
    C ck_term(ck_coeff);  // track the z^0 term
    C sumC = ck_term * C(H_k + H_kn);

    for (int k = 1; k < max_iter; k++) {
        H_k = H_k + unit / R(k);
        H_kn = H_kn + unit / R(k + abs_n);

        // recurrence: ck_term_{k} = ck_term_{k-1} · (-(z²/4)) / (k·(k+n))
        ck_term = ck_term * neg_z_half_sq / R(static_cast<int64_t>(k) * (k + abs_n));

        C contrib = ck_term * C(H_k + H_kn);
        sumC = sumC + contrib;

        if (k >= 5 && sangi::abs(contrib) < eps * sangi::abs(sumC)) break;
    }

    result = result - C(one_over_pi) * z_half_n * sumC;

    // Y_{-n}(z) = (-1)^n · Y_n(z)
    if (n < 0 && abs_n % 2 != 0) result = -result;
    return result;
}

/// K_n(z) — A&S 9.6.11 (integer order, complex argument)
/// K_n(z) = (-1)^{n+1}·I_n(z)·ln(z/2)
///        + (1/2)·(z/2)^{-n}·Σ_{k=0}^{n-1} (-1)^k·(n-k-1)!/k! · (z²/4)^k
///        + (-1)^n·(1/2)·(z/2)^n·Σ_{k=0}^{∞} (H_k+H_{k+n}-2γ)/(k!·(k+n)!) · (z²/4)^k
template<typename R>
[[nodiscard]] Complex<R> besselK_complex(int n, const Complex<R>& z,
                                          R eps, int max_iter,
                                          R pi_val, R euler_gamma) {
    using C = Complex<R>;
    int abs_n = (n >= 0) ? n : -n;

    C z_half = z / R(2);
    C z_half_sq = z_half * z_half;
    C ln_z_half = sangi::log(z_half);
    R two_gamma = R(2) * euler_gamma;

    // Compute I_n(z)
    C in = besselI_complex(abs_n, z, eps, max_iter);

    // Part A: (-1)^{n+1} · I_n(z) · ln(z/2)
    C partA = in * ln_z_half;
    if (abs_n % 2 == 0) partA = -partA;  // (-1)^{n+1}

    // Part B: (1/2)·(z/2)^{-n} · Σ_{k=0}^{n-1} (-1)^k·(n-k-1)!/k! · (z²/4)^k
    C partB(R(0));
    if (abs_n > 0) {
        C z_half_neg_n(R(1));
        for (int i = 0; i < abs_n; i++) z_half_neg_n = z_half_neg_n / z_half;

        // k=0: (n-1)!
        R bk_coeff(1);
        for (int i = 1; i < abs_n; i++) bk_coeff = bk_coeff * R(i);
        C neg_z_half_sq = -z_half_sq;
        C bk_term(bk_coeff);
        C sumB = bk_term;

        for (int k = 1; k < abs_n; k++) {
            // recurrence: bk_term_{k} = bk_term_{k-1} · (-(z²/4)) / (k·(n-k))
            bk_term = bk_term * neg_z_half_sq / R(static_cast<int64_t>(k) * (abs_n - k));
            sumB = sumB + bk_term;
        }
        partB = C(R(1) / R(2)) * z_half_neg_n * sumB;
    }

    // Part C: (-1)^n · (1/2) · (z/2)^n · Σ_{k=0}^{∞} (H_k+H_{k+n}-2γ)/(k!·(k+n)!) · (z²/4)^k
    C z_half_n(R(1));
    for (int i = 0; i < abs_n; i++) z_half_n = z_half_n * z_half;

    // Avoid exact/exact poison (same as besselY): divide harmonic numbers / factorial inverses by the working-precision unit.
    R unit(1);
    if constexpr (std::is_same_v<R, sangi::Float>) {
        unit.setEffectiveBits(sangi::Float::precisionToBits(pi_val.precision()));
    }
    R H_k(0);
    R H_kn(0);
    for (int j = 1; j <= abs_n; j++) H_kn = H_kn + unit / R(j);

    // k=0: 1/n!
    R ck_coeff = unit;
    for (int i = 1; i <= abs_n; i++) ck_coeff = ck_coeff / R(i);

    C ck_term(ck_coeff);
    C sumC = ck_term * C(H_k + H_kn - two_gamma);

    for (int k = 1; k < max_iter; k++) {
        H_k = H_k + unit / R(k);
        H_kn = H_kn + unit / R(k + abs_n);

        // recurrence: (z²/4) / (k·(k+n))
        ck_term = ck_term * z_half_sq / R(static_cast<int64_t>(k) * (k + abs_n));
        C contrib = ck_term * C(H_k + H_kn - two_gamma);
        sumC = sumC + contrib;

        if (k >= 5 && sangi::abs(contrib) < eps * sangi::abs(sumC)) break;
    }

    R sign_n = (abs_n % 2 == 0) ? R(1) : R(-1);
    C partC = C(sign_n * R(1) / R(2)) * z_half_n * sumC;

    C result = partA + partB + partC;
    // K_{-n}(z) = K_n(z) for integer n
    return result;
}

} // namespace detail

// ================================================================
// Complex<double> overloads
// ================================================================

[[nodiscard]] inline Complex<double> besselJ(int n, const Complex<double>& z) {
    return detail::besselJ_complex<double>(n, z,
        std::numeric_limits<double>::epsilon(), 300);
}

[[nodiscard]] inline Complex<double> besselI(int n, const Complex<double>& z) {
    return detail::besselI_complex<double>(n, z,
        std::numeric_limits<double>::epsilon(), 300);
}

[[nodiscard]] inline Complex<double> besselY(int n, const Complex<double>& z) {
    return detail::besselY_complex<double>(n, z,
        std::numeric_limits<double>::epsilon(), 300,
        std::numbers::pi, 0.5772156649015328606065120900824024310421593359);
}

[[nodiscard]] inline Complex<double> besselK(int n, const Complex<double>& z) {
    return detail::besselK_complex<double>(n, z,
        std::numeric_limits<double>::epsilon(), 300,
        std::numbers::pi, 0.5772156649015328606065120900824024310421593359);
}

// ================================================================
// Complex<Float> overloads
// ================================================================

[[nodiscard]] inline Complex<Float> besselJ(int n, const Complex<Float>& z, int precision) {
    using C = Complex<Float>;
    int wp = precision + 20;
    Float::PrecisionScope _ps(wp);   // compute in-function exact÷exact (harmonic numbers H_k etc.) at wp digits
    // PrecisionGuard removed: setResultPrecision(wp) on zw propagates req=wp.
    C zw = z;
    zw.re.setResultPrecision(wp);
    zw.im.setResultPrecision(wp);
    C result = detail::besselJ_complex<Float>(n, zw,
        Float::epsilon(wp), 10 * wp);
    result.re.setPrecision(precision);
    result.im.setPrecision(precision);
    return result;
}

[[nodiscard]] inline Complex<Float> besselJ(int n, const Complex<Float>& z) {
    return besselJ(n, z, Float::defaultPrecision());
}

[[nodiscard]] inline Complex<Float> besselI(int n, const Complex<Float>& z, int precision) {
    using C = Complex<Float>;
    int wp = precision + 20;
    Float::PrecisionScope _ps(wp);   // compute in-function exact÷exact (harmonic numbers H_k etc.) at wp digits
    // PrecisionGuard removed: setResultPrecision(wp) on zw propagates req=wp.
    C zw = z;
    zw.re.setResultPrecision(wp);
    zw.im.setResultPrecision(wp);
    C result = detail::besselI_complex<Float>(n, zw,
        Float::epsilon(wp), 10 * wp);
    result.re.setPrecision(precision);
    result.im.setPrecision(precision);
    return result;
}

[[nodiscard]] inline Complex<Float> besselI(int n, const Complex<Float>& z) {
    return besselI(n, z, Float::defaultPrecision());
}

[[nodiscard]] inline Complex<Float> besselY(int n, const Complex<Float>& z, int precision) {
    using C = Complex<Float>;
    int wp = precision + 30;
    Float::PrecisionScope _ps(wp);   // compute in-function exact÷exact (harmonic numbers H_k etc.) at wp digits
    // PrecisionGuard removed: setResultPrecision(wp) on zw propagates req=wp.
    C zw = z;
    zw.re.setResultPrecision(wp);
    zw.im.setResultPrecision(wp);
    Float pi_val = Float::pi(wp);
    Float euler_gamma = Float::euler(wp);
    C result = detail::besselY_complex<Float>(n, zw,
        Float::epsilon(wp), 10 * wp, pi_val, euler_gamma);
    result.re.setPrecision(precision);
    result.im.setPrecision(precision);
    return result;
}

[[nodiscard]] inline Complex<Float> besselY(int n, const Complex<Float>& z) {
    return besselY(n, z, Float::defaultPrecision());
}

[[nodiscard]] inline Complex<Float> besselK(int n, const Complex<Float>& z, int precision) {
    using C = Complex<Float>;
    int wp = precision + 30;
    Float::PrecisionScope _ps(wp);   // compute in-function exact÷exact (harmonic numbers H_k etc.) at wp digits
    // PrecisionGuard removed: setResultPrecision(wp) on zw propagates req=wp.
    C zw = z;
    zw.re.setResultPrecision(wp);
    zw.im.setResultPrecision(wp);
    Float pi_val = Float::pi(wp);
    Float euler_gamma = Float::euler(wp);
    C result = detail::besselK_complex<Float>(n, zw,
        Float::epsilon(wp), 10 * wp, pi_val, euler_gamma);
    result.re.setPrecision(precision);
    result.im.setPrecision(precision);
    return result;
}

[[nodiscard]] inline Complex<Float> besselK(int n, const Complex<Float>& z) {
    return besselK(n, z, Float::defaultPrecision());
}

// ================================================================
// Hankel functions (1st/2nd kind)
// DLMF 10.4.3:
//   H_n^{(1)}(z) = J_n(z) + i·Y_n(z)
//   H_n^{(2)}(z) = J_n(z) - i·Y_n(z)
// ================================================================

/// Hankel function of the 1st kind H_n^{(1)}(z) — Complex<double>
[[nodiscard]] inline Complex<double> hankelH1(int n, const Complex<double>& z) {
    using C = Complex<double>;
    C jn = besselJ(n, z);
    C yn = besselY(n, z);
    // H^(1) = J + i·Y → (Jre - Yim, Jim + Yre)
    return C(jn.re - yn.im, jn.im + yn.re);
}

/// Hankel function of the 2nd kind H_n^{(2)}(z) — Complex<double>
[[nodiscard]] inline Complex<double> hankelH2(int n, const Complex<double>& z) {
    using C = Complex<double>;
    C jn = besselJ(n, z);
    C yn = besselY(n, z);
    // H^(2) = J - i·Y → (Jre + Yim, Jim - Yre)
    return C(jn.re + yn.im, jn.im - yn.re);
}

/// Hankel function of the 1st kind H_n^{(1)}(z) — Complex<Float>
[[nodiscard]] inline Complex<Float> hankelH1(int n, const Complex<Float>& z, int precision) {
    using C = Complex<Float>;
    int wp = precision + 30;
    Float::PrecisionScope _ps(wp);   // compute in-function exact÷exact at wp digits
    // PrecisionGuard removed: setResultPrecision(wp) on zw propagates req=wp.
    C zw = z;
    zw.re.setResultPrecision(wp);
    zw.im.setResultPrecision(wp);
    Float pi_val = Float::pi(wp);
    Float euler_gamma = Float::euler(wp);
    Float eps = Float::epsilon(wp);
    int max_iter = 10 * wp;

    C jn = detail::besselJ_complex<Float>(n, zw, eps, max_iter);
    C yn = detail::besselY_complex<Float>(n, zw, eps, max_iter, pi_val, euler_gamma);
    // H^(1) = J + i·Y → (Jre - Yim, Jim + Yre)
    C result(jn.re - yn.im, jn.im + yn.re);
    result.re.setPrecision(precision);
    result.im.setPrecision(precision);
    return result;
}

[[nodiscard]] inline Complex<Float> hankelH1(int n, const Complex<Float>& z) {
    return hankelH1(n, z, Float::defaultPrecision());
}

/// Hankel function of the 2nd kind H_n^{(2)}(z) — Complex<Float>
[[nodiscard]] inline Complex<Float> hankelH2(int n, const Complex<Float>& z, int precision) {
    using C = Complex<Float>;
    int wp = precision + 30;
    Float::PrecisionScope _ps(wp);   // compute in-function exact÷exact at wp digits
    // PrecisionGuard removed: setResultPrecision(wp) on zw propagates req=wp.
    C zw = z;
    zw.re.setResultPrecision(wp);
    zw.im.setResultPrecision(wp);
    Float pi_val = Float::pi(wp);
    Float euler_gamma = Float::euler(wp);
    Float eps = Float::epsilon(wp);
    int max_iter = 10 * wp;

    C jn = detail::besselJ_complex<Float>(n, zw, eps, max_iter);
    C yn = detail::besselY_complex<Float>(n, zw, eps, max_iter, pi_val, euler_gamma);
    // H^(2) = J - i·Y → (Jre + Yim, Jim - Yre)
    C result(jn.re + yn.im, jn.im - yn.re);
    result.re.setPrecision(precision);
    result.im.setPrecision(precision);
    return result;
}

[[nodiscard]] inline Complex<Float> hankelH2(int n, const Complex<Float>& z) {
    return hankelH2(n, z, Float::defaultPrecision());
}

// ================================================================
// Hankel function derivatives (recurrence)
// DLMF 10.6.5:
//   H'^{(k)}_n(z) = (H^{(k)}_{n-1}(z) - H^{(k)}_{n+1}(z)) / 2
// ================================================================

/// H'^{(1)}_n(z) — Complex<double>
[[nodiscard]] inline Complex<double> hankelH1Prime(int n, const Complex<double>& z) {
    return (hankelH1(n - 1, z) - hankelH1(n + 1, z)) / Complex<double>(2.0);
}

/// H'^{(2)}_n(z) — Complex<double>
[[nodiscard]] inline Complex<double> hankelH2Prime(int n, const Complex<double>& z) {
    return (hankelH2(n - 1, z) - hankelH2(n + 1, z)) / Complex<double>(2.0);
}

/// H'^{(1)}_n(z) — Complex<Float>
[[nodiscard]] inline Complex<Float> hankelH1Prime(int n, const Complex<Float>& z, int precision) {
    return (hankelH1(n - 1, z, precision) - hankelH1(n + 1, z, precision)) / Complex<Float>(Float(2));
}

[[nodiscard]] inline Complex<Float> hankelH1Prime(int n, const Complex<Float>& z) {
    return hankelH1Prime(n, z, Float::defaultPrecision());
}

/// H'^{(2)}_n(z) — Complex<Float>
[[nodiscard]] inline Complex<Float> hankelH2Prime(int n, const Complex<Float>& z, int precision) {
    return (hankelH2(n - 1, z, precision) - hankelH2(n + 1, z, precision)) / Complex<Float>(Float(2));
}

[[nodiscard]] inline Complex<Float> hankelH2Prime(int n, const Complex<Float>& z) {
    return hankelH2Prime(n, z, Float::defaultPrecision());
}

// ================================================================
// Spherical Hankel functions
// DLMF 10.47.5:
//   h_n^{(1)}(z) = √(π/(2z)) · H_{n+1/2}^{(1)}(z)
//                = j_n(z) + i·y_n(z)
//   h_n^{(2)}(z) = √(π/(2z)) · H_{n+1/2}^{(2)}(z)
//                = j_n(z) - i·y_n(z)
// ================================================================

/// Complex version of the spherical Bessel functions j_n(z), y_n(z) (recurrence)
/// DLMF 10.51.1: j_{n+1}(z) = (2n+1)/z · j_n(z) - j_{n-1}(z)
namespace detail {

template<typename R>
[[nodiscard]] Complex<R> sphericalBesselJ_complex(int n, const Complex<R>& z) {
    using C = Complex<R>;
    if (sangi::abs(z) == R(0))
        return (n == 0) ? C(R(1)) : C(R(0));

    // j_0(z) = sin(z)/z, j_1(z) = sin(z)/z² - cos(z)/z
    C sinz = sangi::sin(z);
    C cosz = sangi::cos(z);
    C inv_z = C(R(1)) / z;

    C j0 = sinz * inv_z;
    if (n == 0) return j0;

    C j1 = sinz * inv_z * inv_z - cosz * inv_z;
    if (n == 1) return j1;

    C j_prev = j0;
    C j_curr = j1;
    for (int k = 1; k < n; k++) {
        C j_next = C(R(2 * k + 1)) * inv_z * j_curr - j_prev;
        j_prev = j_curr;
        j_curr = j_next;
    }
    return j_curr;
}

template<typename R>
[[nodiscard]] Complex<R> sphericalBesselY_complex(int n, const Complex<R>& z) {
    using C = Complex<R>;
    // y_0(z) = -cos(z)/z, y_1(z) = -cos(z)/z² - sin(z)/z
    C sinz = sangi::sin(z);
    C cosz = sangi::cos(z);
    C inv_z = C(R(1)) / z;

    C y0 = -cosz * inv_z;
    if (n == 0) return y0;

    C y1 = -cosz * inv_z * inv_z - sinz * inv_z;
    if (n == 1) return y1;

    C y_prev = y0;
    C y_curr = y1;
    for (int k = 1; k < n; k++) {
        C y_next = C(R(2 * k + 1)) * inv_z * y_curr - y_prev;
        y_prev = y_curr;
        y_curr = y_next;
    }
    return y_curr;
}

} // namespace detail

/// Spherical Hankel function of the 1st kind h_n^{(1)}(z) = j_n(z) + i·y_n(z)
[[nodiscard]] inline Complex<double> sphericalHankelH1(int n, const Complex<double>& z) {
    using C = Complex<double>;
    C jn = detail::sphericalBesselJ_complex<double>(n, z);
    C yn = detail::sphericalBesselY_complex<double>(n, z);
    // j + i·y → (j.re - y.im, j.im + y.re)
    return C(jn.re - yn.im, jn.im + yn.re);
}

/// Spherical Hankel function of the 2nd kind h_n^{(2)}(z) = j_n(z) - i·y_n(z)
[[nodiscard]] inline Complex<double> sphericalHankelH2(int n, const Complex<double>& z) {
    using C = Complex<double>;
    C jn = detail::sphericalBesselJ_complex<double>(n, z);
    C yn = detail::sphericalBesselY_complex<double>(n, z);
    // j - i·y → (j.re + y.im, j.im - y.re)
    return C(jn.re + yn.im, jn.im - yn.re);
}

} // namespace special
} // namespace sangi

#endif // SANGI_SPECIAL_BESSEL_HPP