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

// orthogonal_classical.hpp
// Classical orthogonal polynomials (Hermite / Laguerre / Chebyshev)
//
// Functions provided:
//   hermiteH(n, x)            — Hermite polynomial H_n(x) (physicists' version, weight e^{-x²})
//   hermiteHe(n, x)           — Hermite polynomial He_n(x) (probabilists' version, weight e^{-x²/2})
//   laguerreL(n, x)           — Laguerre polynomial L_n(x) (weight e^{-x}, [0,∞))
//   assocLaguerreL(n, α, x)   — generalized Laguerre polynomial L_n^α(x) (weight x^α e^{-x})
//   chebyshevT(n, x)          — Chebyshev of the first kind T_n(x) (weight 1/√(1-x²))
//   chebyshevU(n, x)          — Chebyshev of the second kind U_n(x) (weight √(1-x²))
//
// All computed stably via three-term recurrences (Bonnet/Favard).
//
// Relation to the existing jacobiP:
//   T_n(x) ∝ P_n^{(-1/2,-1/2)}(x), U_n(x) ∝ P_n^{(1/2,1/2)}(x)
// The standalone implementation is faster with lower overhead (about one order of magnitude in the double-precision range).
//
// Supported types:
//   float, double, long double — direct recurrence
//   Float<P>                   — TODO: arbitrary-precision version (native only for now)
// Audit: TODO/PLAN_CAS_PHASE2C_ALGEBRAIC_20260507.md is a separate matter; this implementation is
//       priority item 1 of the audit (project_sangi_special_functions_audit_20260507).

#ifndef SANGI_SPECIAL_ORTHOGONAL_CLASSICAL_HPP
#define SANGI_SPECIAL_ORTHOGONAL_CLASSICAL_HPP

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

namespace sangi {
namespace special {

// ================================================================
// Hermite polynomial (physicists' version) H_n(x)
// Recurrence: H_{n+1}(x) = 2x·H_n(x) - 2n·H_{n-1}(x)
// H_0 = 1, H_1 = 2x
// Orthogonality: ∫_{-∞}^{∞} H_m(x) H_n(x) e^{-x²} dx = √π · 2^n · n! · δ_{mn}
// Eigenfunctions of the quantum harmonic oscillator
// ================================================================

template<IsNativeFloat T>
[[nodiscard]] T hermiteH(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 (n == 0)          return T(1);
    T two_x = T(2) * x;
    if (n == 1)          return two_x;

    T h_prev = T(1);    // H_0
    T h_curr = two_x;   // H_1
    for (int k = 1; k < n; ++k) {
        T h_next = two_x * h_curr - T(2 * k) * h_prev;
        h_prev = h_curr;
        h_curr = h_next;
    }
    return h_curr;
}

// ================================================================
// Hermite polynomial (probabilists' version) He_n(x)
// Recurrence: He_{n+1}(x) = x·He_n(x) - n·He_{n-1}(x)
// He_0 = 1, He_1 = x
// Orthogonality: ∫ He_m He_n e^{-x²/2} dx = √(2π) · n! · δ_{mn}
// Relation: H_n(x) = 2^{n/2} · He_n(x·√2)
// ================================================================

template<IsNativeFloat T>
[[nodiscard]] T hermiteHe(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 (n == 0)          return T(1);
    if (n == 1)          return x;

    T h_prev = T(1);    // He_0
    T h_curr = x;       // He_1
    for (int k = 1; k < n; ++k) {
        T h_next = x * h_curr - T(k) * h_prev;
        h_prev = h_curr;
        h_curr = h_next;
    }
    return h_curr;
}

// ================================================================
// Laguerre polynomial L_n(x)
// Recurrence: (k+1) L_{k+1}(x) = (2k+1-x) L_k(x) - k L_{k-1}(x)
// L_0 = 1, L_1 = 1 - x
// Orthogonality: ∫_0^∞ L_m L_n e^{-x} dx = δ_{mn}
// ================================================================

template<IsNativeFloat T>
[[nodiscard]] T laguerreL(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 (n == 0)          return T(1);
    T l1 = T(1) - x;
    if (n == 1)          return l1;

    T l_prev = T(1);
    T l_curr = l1;
    for (int k = 1; k < n; ++k) {
        T kf = T(k);
        T l_next = ((T(2) * kf + T(1) - x) * l_curr - kf * l_prev) / (kf + T(1));
        l_prev = l_curr;
        l_curr = l_next;
    }
    return l_curr;
}

// ================================================================
// Generalized Laguerre polynomial L_n^α(x)
// Recurrence: (k+1) L_{k+1}^α(x) = (2k+1+α-x) L_k^α(x) - (k+α) L_{k-1}^α(x)
// L_0^α = 1, L_1^α = 1 + α - x
// Orthogonality: ∫_0^∞ L_m^α L_n^α x^α e^{-x} dx = Γ(n+α+1)/n! · δ_{mn}
// Used in the radial wavefunction of the hydrogen atom (α depends on the azimuthal quantum number)
// ================================================================

template<IsNativeFloat T>
[[nodiscard]] T assocLaguerreL(int n, T alpha, T x) {
    if (std::isnan(x) || std::isnan(alpha))
        return std::numeric_limits<T>::quiet_NaN();
    if (n < 0) return std::numeric_limits<T>::quiet_NaN();
    if (n == 0) return T(1);
    T l1 = T(1) + alpha - x;
    if (n == 1) return l1;

    T l_prev = T(1);
    T l_curr = l1;
    for (int k = 1; k < n; ++k) {
        T kf = T(k);
        T l_next = ((T(2) * kf + T(1) + alpha - x) * l_curr
                    - (kf + alpha) * l_prev) / (kf + T(1));
        l_prev = l_curr;
        l_curr = l_next;
    }
    return l_curr;
}

// ================================================================
// Chebyshev of the first kind T_n(x)
// Recurrence: T_{n+1}(x) = 2x·T_n(x) - T_{n-1}(x)
// T_0 = 1, T_1 = x
// Orthogonality: ∫_{-1}^{1} T_m T_n / √(1-x²) dx = δ_{mn} · π/(2 if n>0, 1 if n=0)
// Relation: T_n(cos θ) = cos(nθ)
// Stable for |x| ≤ 1; the recurrence still works for |x| > 1 (can also be computed via cosh)
// ================================================================

template<IsNativeFloat T>
[[nodiscard]] T chebyshevT(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 (n == 0)          return T(1);
    if (n == 1)          return x;

    T t_prev = T(1);    // T_0
    T t_curr = x;       // T_1
    T two_x = T(2) * x;
    for (int k = 1; k < n; ++k) {
        T t_next = two_x * t_curr - t_prev;
        t_prev = t_curr;
        t_curr = t_next;
    }
    return t_curr;
}

// ================================================================
// Chebyshev of the second kind U_n(x)
// Recurrence: U_{n+1}(x) = 2x·U_n(x) - U_{n-1}(x)
// U_0 = 1, U_1 = 2x
// Orthogonality: ∫_{-1}^{1} U_m U_n √(1-x²) dx = (π/2)·δ_{mn}
// Relation: U_n(cos θ) = sin((n+1)θ) / sin θ
// ================================================================

template<IsNativeFloat T>
[[nodiscard]] T chebyshevU(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 (n == 0)          return T(1);
    T two_x = T(2) * x;
    if (n == 1)          return two_x;

    T u_prev = T(1);    // U_0
    T u_curr = two_x;   // U_1
    for (int k = 1; k < n; ++k) {
        T u_next = two_x * u_curr - u_prev;
        u_prev = u_curr;
        u_curr = u_next;
    }
    return u_curr;
}

} // namespace special
} // namespace sangi

#endif // SANGI_SPECIAL_ORTHOGONAL_CLASSICAL_HPP