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

// fermi_dirac.hpp
// Fermi-Dirac integrals
//
// Provided functions:
//   fermiDirac(s, x)     — complete Fermi-Dirac integral F_s(x) = (1/Γ(s+1)) ∫₀^∞ t^s/(e^{t-x}+1) dt
//   fermiDiracHalf(x)    — F_{1/2}(x)
//   fermiDiracMHalf(x)   — F_{-1/2}(x)
//   fermiDirac3Half(x)   — F_{3/2}(x)
//
// Supported types: float, double, long double

#ifndef SANGI_SPECIAL_FERMI_DIRAC_HPP
#define SANGI_SPECIAL_FERMI_DIRAC_HPP

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

namespace sangi {
namespace special {

// ================================================================
// Fermi-Dirac integral F_s(x) — Gauss-Laguerre numerical integration
// F_s(x) = (1/Γ(s+1)) ∫₀^∞ t^s / (exp(t-x) + 1) dt
// ================================================================

namespace detail {

// Nodes and weights of 16-point Gauss-Laguerre quadrature
template<typename T>
struct GaussLaguerre16 {
    static constexpr int N = 16;
    static const T nodes[16];
    static const T weights[16];
};

template<typename T>
const T GaussLaguerre16<T>::nodes[16] = {
    T(0.08764941047892727), T(0.46269632891508116), T(1.14105777483011440),
    T(2.12928364509838070), T(3.43708663389320100), T(5.07801861543098320),
    T(7.07033853504823640), T(9.43831433639194020), T(12.21422336852287500),
    T(15.44152736879037100), T(19.18015685675218500), T(23.51590569239297000),
    T(28.57872974028905200), T(34.58339944946498500), T(41.94045264765709200),
    T(51.70116033954263500)
};

template<typename T>
const T GaussLaguerre16<T>::weights[16] = {
    T(0.20615171495780100), T(0.33105785495088000), T(0.26579577764422000),
    T(0.13629693429637800), T(0.04736616898298300), T(0.01129925104038000),
    T(0.00184907094352631), T(0.00020483687597962), T(0.00001487641200779),
    T(6.82831933087120e-07), T(1.88102484107967e-08), T(2.86235024297389e-10),
    T(2.12707903925667e-12), T(6.29796700251656e-15), T(5.05047370003552e-18),
    T(4.16146237037285e-22)
};

template<IsNativeFloat T>
T fermiDirac_quadrature(T s, T x) {
    // Substitution u = t - x (for x > 0) improves accuracy
    // F_s(x) = (1/Γ(s+1)) ∫₀^∞ t^s / (exp(t-x)+1) dt

    // For large x: Sommerfeld expansion
    // F_s(x) ≈ x^{s+1}/((s+1)Γ(s+1)) [1 + s(s+1)π²/(6x²) + ...]
    if (x > T(20)) {
        T pi = std::numbers::pi_v<T>;
        T xs1 = std::pow(x, s + T(1)) / (s + T(1));
        T correction = T(1) + s * (s + T(1)) * pi * pi / (T(6) * x * x);
        return xs1 * correction / std::tgamma(s + T(1));
    }

    // Gauss-Laguerre: ∫₀^∞ f(t) e^{-t} dt ≈ Σ w_i f(t_i)
    // F_s(x) · Γ(s+1) = ∫₀^∞ t^s / (exp(t-x)+1) dt
    //                  = ∫₀^∞ [t^s · exp(-t+x)] / (1 + exp(-t+x)) · exp(t) · exp(-t) dt
    //
    // Directly: ∫₀^∞ t^s / (exp(t-x)+1) dt = ∫₀^∞ [t^s · e^{-t}] · e^t / (e^{t-x}+1) dt
    //      = ∫₀^∞ [t^s · e^{-t}] · 1 / (1 + e^{-x} · e^{-(... no
    //
    // Direct numerical integration (trapezoidal rule + substitution)
    T gamma_s1 = std::tgamma(s + T(1));

    // Gauss-Laguerre with integrand t^s / (exp(t-x)+1) * exp(t)
    // = t^s * exp(t) / (exp(t-x)+1) = t^s * exp(x) / (1 + exp(x-t))
    // ... but this diverges. Use direct approach:
    // ∫₀^∞ t^s e^{-t} · [e^t / (e^{t-x}+1)] dt = ∫₀^∞ t^s e^{-t} · [1/(1+e^{-x})] ... no

    // Simplest: use Gauss-Laguerre with g(t) = t^s / (1 + exp(t - x))
    // ∫₀^∞ g(t) dt = ∫₀^∞ g(t)·e^t · e^{-t} dt ≈ Σ w_i · g(t_i) · e^{t_i}
    using GL = GaussLaguerre16<T>;
    T sum = T(0);
    for (int i = 0; i < GL::N; ++i) {
        T t = GL::nodes[i];
        T integrand = std::pow(t, s) / (T(1) + std::exp(t - x));
        sum += GL::weights[i] * integrand * std::exp(t);
    }

    return sum / gamma_s1;
}

} // namespace detail

/// General Fermi-Dirac integral F_s(x)
template<IsNativeFloat T>
[[nodiscard]] T fermiDirac(T s, T x) {
    if (std::isnan(s) || std::isnan(x))
        return std::numeric_limits<T>::quiet_NaN();
    return detail::fermiDirac_quadrature(s, x);
}

/// F_{1/2}(x)
template<IsNativeFloat T>
[[nodiscard]] T fermiDiracHalf(T x) {
    return fermiDirac(T(0.5), x);
}

/// F_{-1/2}(x)
template<IsNativeFloat T>
[[nodiscard]] T fermiDiracMHalf(T x) {
    return fermiDirac(T(-0.5), x);
}

/// F_{3/2}(x)
template<IsNativeFloat T>
[[nodiscard]] T fermiDirac3Half(T x) {
    return fermiDirac(T(1.5), x);
}

} // namespace special
} // namespace sangi

#endif // SANGI_SPECIAL_FERMI_DIRAC_HPP