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

// fresnel.hpp
// Fresnel integrals S(x), C(x)
//
// Definitions (DLMF §7.2.7-8, same normalization as SciPy/Boost):
//   C(x) = ∫_0^x cos(π t² / 2) dt
//   S(x) = ∫_0^x sin(π t² / 2) dt
//
// Implementation:
//   |x| ≤ 4.0  — Taylor series
//   |x| > 4.0  — asymptotic expansion (DLMF 7.12) evaluated with optimal truncation (at the smallest term)
//                C(x) = 1/2 + f(x) sin(πx²/2) - g(x) cos(πx²/2)
//                S(x) = 1/2 - f(x) cos(πx²/2) - g(x) sin(πx²/2)
//                f(x) = (1/(πx)) · Σ t_k, t_{k+1}/t_k = -(4k+1)(4k+3)/(πx²)²
//                g(x) = (1/(π²x³)) · Σ v_k, v_{k+1}/v_k = -(4k+3)(4k+5)/(πx²)²
//
// Supported types: float, double, long double

#ifndef SANGI_SPECIAL_FRESNEL_HPP
#define SANGI_SPECIAL_FRESNEL_HPP

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

namespace sangi {
namespace special {

namespace detail_fresnel {

// Compute (C, S) simultaneously via the Taylor series for |x| ≤ 4
//   C(x) = Σ_{k=0}^∞ (-1)^k (π/2)^{2k}   x^{4k+1} / ((2k)!   · (4k+1))
//   S(x) = Σ_{k=0}^∞ (-1)^k (π/2)^{2k+1} x^{4k+3} / ((2k+1)! · (4k+3))
template<IsNativeFloat T>
inline std::pair<T, T> fresnelTaylor(T x) {
    constexpr T half_pi = std::numbers::pi_v<T> / T(2);

    T x2 = x * x;
    T pix2_2 = half_pi * x2;                  // (π/2) x²
    T q = -pix2_2 * pix2_2;                   // -(π/2)² x⁴

    T t_C = x;
    T sum_C = x;
    for (int k = 0; k < 200; ++k) {
        T denom = T((2 * k + 1)) * T((2 * k + 2));
        t_C *= q / denom;
        T contribution = t_C / T(4 * k + 5);
        sum_C += contribution;
        if (std::abs(contribution) <
            std::numeric_limits<T>::epsilon() * std::abs(sum_C))
            break;
    }

    T t_S = half_pi * x * x2;                 // (π/2) x³
    T sum_S = t_S / T(3);
    for (int k = 0; k < 200; ++k) {
        T denom = T((2 * k + 2)) * T((2 * k + 3));
        t_S *= q / denom;
        T contribution = t_S / T(4 * k + 7);
        sum_S += contribution;
        if (std::abs(contribution) <
            std::numeric_limits<T>::epsilon() * std::abs(sum_S))
            break;
    }
    return { sum_C, sum_S };
}

// Asymptotic expansion for |x| > 4 — truncated at the smallest term (optimal truncation)
template<IsNativeFloat T>
inline std::pair<T, T> fresnelAsymptotic(T x) {
    constexpr T pi = std::numbers::pi_v<T>;
    constexpr T half = T(1) / T(2);

    T ax = std::abs(x);
    T x2 = ax * ax;
    T pix2 = pi * x2;
    T u = T(1) / (pix2 * pix2);    // 1 / (π x²)²

    // f-series: t_0 = 1, t_{k+1} = t_k · -(4k+1)(4k+3) · u
    T tf = T(1);
    T fsum = T(1);
    T tf_prev_abs = T(1);
    for (int k = 0; k < 50; ++k) {
        T tf_new = tf * (-T((4 * k + 1)) * T((4 * k + 3)) * u);
        T abs_new = std::abs(tf_new);
        if (abs_new > tf_prev_abs) break;   // optimal truncation
        tf = tf_new;
        fsum += tf;
        tf_prev_abs = abs_new;
        if (abs_new < std::numeric_limits<T>::epsilon() * std::abs(fsum)) break;
    }
    T f = fsum / (pi * ax);

    // g-series: v_0 = 1, v_{k+1} = v_k · -(4k+3)(4k+5) · u
    T tg = T(1);
    T gsum = T(1);
    T tg_prev_abs = T(1);
    for (int k = 0; k < 50; ++k) {
        T tg_new = tg * (-T((4 * k + 3)) * T((4 * k + 5)) * u);
        T abs_new = std::abs(tg_new);
        if (abs_new > tg_prev_abs) break;
        tg = tg_new;
        gsum += tg;
        tg_prev_abs = abs_new;
        if (abs_new < std::numeric_limits<T>::epsilon() * std::abs(gsum)) break;
    }
    T g = gsum / (pi * pi * ax * ax * ax);

    T phase = pi * x2 * half;          // π x² / 2
    T s = std::sin(phase);
    T c = std::cos(phase);

    T C = half + f * s - g * c;
    T S = half - f * c - g * s;

    if (x < T(0)) { C = -C; S = -S; }
    return { C, S };
}

}  // namespace detail_fresnel

template<IsNativeFloat T>
[[nodiscard]] T fresnelC(T x) {
    if (std::isnan(x)) return std::numeric_limits<T>::quiet_NaN();
    if (std::isinf(x)) return x > T(0) ? T(0.5) : T(-0.5);

    constexpr T threshold = T(4.0);
    if (std::abs(x) <= threshold)
        return detail_fresnel::fresnelTaylor(x).first;
    return detail_fresnel::fresnelAsymptotic(x).first;
}

template<IsNativeFloat T>
[[nodiscard]] T fresnelS(T x) {
    if (std::isnan(x)) return std::numeric_limits<T>::quiet_NaN();
    if (std::isinf(x)) return x > T(0) ? T(0.5) : T(-0.5);

    constexpr T threshold = T(4.0);
    if (std::abs(x) <= threshold)
        return detail_fresnel::fresnelTaylor(x).second;
    return detail_fresnel::fresnelAsymptotic(x).second;
}

template<IsNativeFloat T>
[[nodiscard]] std::pair<T, T> fresnelCS(T x) {
    if (std::isnan(x))
        return { std::numeric_limits<T>::quiet_NaN(),
                 std::numeric_limits<T>::quiet_NaN() };
    if (std::isinf(x))
        return { x > T(0) ? T(0.5) : T(-0.5),
                 x > T(0) ? T(0.5) : T(-0.5) };

    constexpr T threshold = T(4.0);
    if (std::abs(x) <= threshold) return detail_fresnel::fresnelTaylor(x);
    return detail_fresnel::fresnelAsymptotic(x);
}

} // namespace special
} // namespace sangi

#endif // SANGI_SPECIAL_FRESNEL_HPP