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

// dawson.hpp
// Dawson function
//
// Provided functions:
//   dawson(x) — Dawson function F(x) = exp(-x²) ∫₀ˣ exp(t²) dt
//
// Supported types: float, double, long double

#ifndef SANGI_SPECIAL_DAWSON_HPP
#define SANGI_SPECIAL_DAWSON_HPP

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

namespace sangi {
namespace special {

// ================================================================
// Dawson function F(x) = exp(-x²) ∫₀ˣ exp(t²) dt
// Algorithm of Rybicki (1989)
// ================================================================

template<IsNativeFloat T>
[[nodiscard]] T dawson(T x) {
    if (std::isnan(x)) return std::numeric_limits<T>::quiet_NaN();
    if (x == T(0)) return T(0);

    // Odd function: F(-x) = -F(x)
    T sign = T(1);
    T ax = x;
    if (x < T(0)) { sign = T(-1); ax = -x; }

    T result;

    if (ax < T(0.2)) {
        // Taylor series: F(x) = x - 2x³/3 + 4x⁵/15 - 8x⁷/105 + ...
        // F(x) = Σ_{n=0}^∞ (-1)^n · 2^n · x^{2n+1} / (1·3·5···(2n+1))
        T x2 = ax * ax;
        T term = ax;
        T sum = term;
        for (int n = 1; n < 50; ++n) {
            term *= -T(2) * x2 / T(2 * n + 1);
            sum += term;
            if (std::abs(term) < std::numeric_limits<T>::epsilon() * std::abs(sum))
                break;
        }
        result = sum;
    } else if (ax > T(5.0)) {
        // Asymptotic expansion: F(x) ≈ 1/(2x) + 1/(4x³) + 3/(8x⁵) + ...
        // F(x) = (1/(2x)) Σ_{n=0}^∞ (2n-1)!! / (2x²)^n
        T x2 = ax * ax;
        T inv_2x2 = T(1) / (T(2) * x2);
        T term = T(1);
        T sum = T(1);
        for (int n = 1; n < 30; ++n) {
            term *= T(2 * n - 1) * inv_2x2;
            sum += term;
            if (std::abs(term) < std::numeric_limits<T>::epsilon() * std::abs(sum))
                break;
        }
        result = sum / (T(2) * ax);
    } else {
        // Simpson's rule: F(x) = exp(-x²) ∫₀ˣ exp(t²) dt
        int N = static_cast<int>(ax * T(100)) + 100;
        if (N % 2 != 0) ++N;
        T dx = ax / T(N);
        T sum = T(0);
        T x2 = ax * ax;

        for (int i = 0; i <= N; ++i) {
            T t = dx * T(i);
            T w = (i == 0 || i == N) ? T(1) : (i % 2 == 0 ? T(2) : T(4));
            sum += w * std::exp(t * t - x2);
        }
        result = sum * dx / T(3);
    }

    return sign * result;
}

} // namespace special
} // namespace sangi

#endif // SANGI_SPECIAL_DAWSON_HPP