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

// clausen.hpp
// Clausen functions
//
// Provided functions:
//   clausen(x)     — Clausen function Cl₂(θ) = -∫₀^θ ln|2sin(t/2)| dt
//   clausenN(n, x) — generalized Clausen function Cl_n(θ) / Sl_n(θ)
//
// Supported types: float, double, long double

#ifndef SANGI_SPECIAL_CLAUSEN_HPP
#define SANGI_SPECIAL_CLAUSEN_HPP

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

namespace sangi {
namespace special {

// ================================================================
// Clausen function Cl₂(θ) = Σ_{k=1}^∞ sin(kθ)/k²
// = -∫₀^θ ln|2sin(t/2)| dt
// ================================================================

template<IsNativeFloat T>
[[nodiscard]] T clausen(T theta) {
    if (std::isnan(theta)) return std::numeric_limits<T>::quiet_NaN();

    // Periodicity: Cl₂(θ + 2π) = Cl₂(θ)
    T pi2 = T(2) * std::numbers::pi_v<T>;
    T x = std::fmod(theta, pi2);
    if (x < T(0)) x += pi2;
    if (x > std::numbers::pi_v<T>) x = pi2 - x;  // antisymmetry: Cl₂(2π-θ) = -Cl₂(θ)
    T sign = (std::fmod(theta, pi2) > std::numbers::pi_v<T>) ? T(-1) : T(1);
    if (std::fmod(theta, pi2) < T(0))
        sign = (std::fmod(theta, pi2) + pi2 > std::numbers::pi_v<T>) ? T(-1) : T(1);

    // x ∈ [0, π]
    if (x < T(1e-15)) return T(0);
    if (std::abs(x - std::numbers::pi_v<T>) < T(1e-15)) return T(0);

    // Fourier series: Cl₂(θ) = Σ sin(kθ)/k²
    // Convergence is slow, so use acceleration via Bernoulli numbers
    // Small θ: Cl₂(θ) = θ(1 - ln|θ|) + θ³/36 - θ⁵/3600 + ...

    if (x < T(0.5)) {
        // Small-θ expansion: Cl₂(θ) = θ - θ·ln(θ) + θ³/36 - θ⁵/3600 + θ⁷/211680 - ...
        T ln_x = std::log(x);
        T x2 = x * x;
        T result = x * (T(1) - ln_x)
                   + x * x2 / T(36)
                   - x * x2 * x2 / T(3600)
                   + x * x2 * x2 * x2 / T(211680)
                   - x * x2 * x2 * x2 * x2 / T(10886400);
        return sign * result;
    }

    // General case: direct Fourier series (Kahan summation)
    T sum = T(0);
    T c = T(0);
    for (int k = 1; k <= 1000; ++k) {
        T term = std::sin(T(k) * x) / (T(k) * T(k));
        T y = term - c;
        T t = sum + y;
        c = (t - sum) - y;
        sum = t;
        if (k >= 10 && std::abs(term) <
            std::numeric_limits<T>::epsilon() * std::abs(sum))
            break;
    }
    return sign * sum;
}

// ================================================================
// Generalized Clausen: Cl_n(θ) = Σ cos(kθ)/k^n (n even)
//                  Sl_n(θ) = Σ sin(kθ)/k^n (n odd)
// (Cl₂ corresponds to Sl₂)
// ================================================================

template<IsNativeFloat T>
[[nodiscard]] T clausenCl(int n, T theta) {
    if (std::isnan(theta) || n < 1)
        return std::numeric_limits<T>::quiet_NaN();

    T sum = T(0);
    for (int k = 1; k <= 2000; ++k) {
        T term;
        if (n % 2 == 0) {
            term = std::cos(T(k) * theta) / std::pow(T(k), T(n));
        } else {
            term = std::sin(T(k) * theta) / std::pow(T(k), T(n));
        }
        sum += term;
        if (k >= 10 && std::abs(term) <
            std::numeric_limits<T>::epsilon() * std::abs(sum) * T(10))
            break;
    }
    return sum;
}

} // namespace special
} // namespace sangi

#endif // SANGI_SPECIAL_CLAUSEN_HPP