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

// elliptic.hpp
// Templates for elliptic integrals (Carlson symmetric form + Legendre form) and Jacobi elliptic functions
//
// Provided functions:
//   carlsonRF(x,y,z), carlsonRC(x,y), carlsonRD(x,y,z), carlsonRJ(x,y,z,p)
//   ellipticK(k), ellipticE(k), ellipticPi(n,k)        — complete elliptic integrals
//   ellipticF(phi,k), ellipticE(phi,k), ellipticPi(n,phi,k) — incomplete elliptic integrals
//   jacobiSn(u,k), jacobiCn(u,k), jacobiDn(u,k)
//
// Supported types:
//   float, double, long double — implemented directly in this header
//   Float                      — delegates to the sangi:: implementation
//   Complex<double>            — complex extension of the Carlson symmetric form
//   Complex<Float>             — arbitrary-precision Carlson symmetric form

#ifndef SANGI_SPECIAL_ELLIPTIC_HPP
#define SANGI_SPECIAL_ELLIPTIC_HPP

#include <math/special/special_common.hpp>
#include <math/special/gamma.hpp>
#include <math/core/mp/Float.hpp>
#include <math/core/Complex.hpp>
#include <cmath>
#include <limits>
#include <algorithm>
#include <vector>

namespace sangi {
namespace special {

// ================================================================
// Carlson R_C(x, y)
// R_C(x, y) = (1/2)∫₀^∞ dt / [(t+y)√(t+x)]
// Quadratically convergent iteration via the duplication theorem
// ================================================================

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

    const T eps = std::numeric_limits<T>::epsilon();

    for (int iter = 0; iter < 100; iter++) {
        T lam = T(2) * std::sqrt(x * y) + y;
        x = (x + lam) / T(4);
        y = (y + lam) / T(4);

        T A = (x + y + y) / T(3);
        T s = (y - A) / A;
        if (std::abs(s) < eps) break;
    }

    T A = (x + y + y) / T(3);
    T s = (y - A) / A;
    // 1 + (3/10)s² + (1/7)s³ + (3/8)s⁴ + (9/22)s⁵
    T poly = T(1) + s * s * (T(3) / T(10)
             + s * (T(1) / T(7)
             + s * (T(3) / T(8)
             + s * T(9) / T(22))));

    return poly / std::sqrt(A);
}

// Float type: use sangi::carlsonRC(x, y, precision) directly.

// ================================================================
// Carlson R_F(x, y, z)
// R_F(x, y, z) = (1/2)∫₀^∞ dt / √[(t+x)(t+y)(t+z)]
// x, y, z ≥ 0, at most one is 0
// ================================================================

template<IsNativeFloat T>
[[nodiscard]] T carlsonRF(T x, T y, T z) {
    if (std::isnan(x) || std::isnan(y) || std::isnan(z))
        return std::numeric_limits<T>::quiet_NaN();
    if (x < T(0) || y < T(0) || z < T(0))
        return std::numeric_limits<T>::quiet_NaN();

    const T eps = std::numeric_limits<T>::epsilon();

    for (int iter = 0; iter < 100; iter++) {
        T sqx = std::sqrt(x), sqy = std::sqrt(y), sqz = std::sqrt(z);
        T lam = sqx * sqy + sqy * sqz + sqz * sqx;
        x = (x + lam) / T(4);
        y = (y + lam) / T(4);
        z = (z + lam) / T(4);

        T A = (x + y + z) / T(3);
        T dx = (A - x) / A;
        T dy = (A - y) / A;
        T dz = (A - z) / A;
        if (std::max({std::abs(dx), std::abs(dy), std::abs(dz)}) < eps) break;
    }

    T A = (x + y + z) / T(3);
    T X = (A - x) / A;
    T Y = (A - y) / A;
    T Z = -(X + Y);

    T E2 = X * Y - Z * Z;
    T E3 = X * Y * Z;
    T poly = T(1) - E2 / T(10) + E3 / T(14) + E2 * E2 / T(24)
             - T(3) * E2 * E3 / T(44);

    return poly / std::sqrt(A);
}

// Float type: use sangi::carlsonRF(x, y, z, precision) directly.

// ================================================================
// Carlson R_D(x, y, z)
// R_D(x, y, z) = (3/2)∫₀^∞ dt / [(t+z)^{3/2}√((t+x)(t+y))]
// x, y ≥ 0 (at most one is 0), z > 0
// ================================================================

template<IsNativeFloat T>
[[nodiscard]] T carlsonRD(T x, T y, T z) {
    if (std::isnan(x) || std::isnan(y) || std::isnan(z))
        return std::numeric_limits<T>::quiet_NaN();
    if (x < T(0) || y < T(0) || z <= T(0))
        return std::numeric_limits<T>::quiet_NaN();

    const T eps = std::numeric_limits<T>::epsilon();
    T sigma = T(0);
    T fac = T(1);

    for (int iter = 0; iter < 100; iter++) {
        T sqx = std::sqrt(x), sqy = std::sqrt(y), sqz = std::sqrt(z);
        T lam = sqx * sqy + sqy * sqz + sqz * sqx;
        sigma += fac / (sqz * (z + lam));
        fac /= T(4);
        x = (x + lam) / T(4);
        y = (y + lam) / T(4);
        z = (z + lam) / T(4);

        T A = (x + y + T(3) * z) / T(5);
        T dx = (A - x) / A;
        T dy = (A - y) / A;
        T dz = (A - z) / A;
        if (std::max({std::abs(dx), std::abs(dy), std::abs(dz)}) < eps) break;
    }

    T A = (x + y + T(3) * z) / T(5);
    T dx = (A - x) / A;
    T dy = (A - y) / A;
    T dz = (A - z) / A;
    T ea = dx * dy;
    T eb = dz * dz;
    T ec = ea - eb;
    T ed = ea - T(6) * eb;
    T ee = ed + ec + ec;
    T poly = T(1)
        + ed * (T(-3) / T(14) + T(9) / T(88) * ed - T(9) / T(52) * dz * ee)
        + dz * (T(1) / T(6) * ee + dz * (T(-9) / T(22) * ec + dz * T(3) / T(26) * ea));

    return T(3) * sigma + fac * poly / (A * std::sqrt(A));
}

// Float type: use sangi::carlsonRD(x, y, z, precision) directly.

// ================================================================
// Carlson R_J(x, y, z, p)
// R_J(x, y, z, p) = (3/2)∫₀^∞ dt / [(t+p)√((t+x)(t+y)(t+z))]
// x, y, z ≥ 0 (at most one is 0), p > 0
// ================================================================

template<IsNativeFloat T>
[[nodiscard]] T carlsonRJ(T x, T y, T z, T p) {
    if (std::isnan(x) || std::isnan(y) || std::isnan(z) || std::isnan(p))
        return std::numeric_limits<T>::quiet_NaN();
    if (x < T(0) || y < T(0) || z < T(0) || p <= T(0))
        return std::numeric_limits<T>::quiet_NaN();

    const T eps = std::numeric_limits<T>::epsilon();
    T sigma = T(0);
    T fac = T(1);

    for (int iter = 0; iter < 100; iter++) {
        T sqx = std::sqrt(x), sqy = std::sqrt(y), sqz = std::sqrt(z), sqp = std::sqrt(p);
        T lam = sqx * sqy + sqy * sqz + sqz * sqx;
        T alpha = p * (sqx + sqy + sqz) + sqx * sqy * sqz;
        alpha *= alpha;
        T beta = p * (p + lam) * (p + lam);
        sigma += fac * carlsonRC(alpha, beta);
        fac /= T(4);
        x = (x + lam) / T(4);
        y = (y + lam) / T(4);
        z = (z + lam) / T(4);
        p = (p + lam) / T(4);

        T A = (x + y + z + p + p) / T(5);
        T dx = (A - x) / A;
        T dy = (A - y) / A;
        T dz = (A - z) / A;
        T dp = (A - p) / A;
        if (std::max({std::abs(dx), std::abs(dy), std::abs(dz), std::abs(dp)}) < eps) break;
    }

    T A = (x + y + z + p + p) / T(5);
    T dx = (A - x) / A;
    T dy = (A - y) / A;
    T dz = (A - z) / A;
    T dp = (A - p) / A;
    T ea = dx * (dy + dz) + dy * dz;
    T eb = dx * dy * dz;
    T ec = dp * dp;
    T ed = ea - T(3) * ec;
    T ee = eb + T(2) * dp * (ea - ec);
    T poly = T(1)
        + ed * (T(-3) / T(14) + T(9) / T(88) * ed - T(9) / T(52) * ee)
        + eb * (T(1) / T(6) + dp * (T(-6) / T(22) + dp * T(3) / T(26)))
        + dp * ea * (T(1) / T(3) - dp * T(3) / T(22))
        - T(1) / T(3) * dp * ec;

    return T(3) * sigma + fac * poly / (A * std::sqrt(A));
}

// Float type: use sangi::carlsonRJ(x, y, z, p, precision) directly.

// ================================================================
// Complete elliptic integral of the first kind K(k)
// K(k) = R_F(0, 1-k², 1)
// ================================================================

template<IsNativeFloat T>
[[nodiscard]] T ellipticK(T k) {
    if (std::isnan(k)) return std::numeric_limits<T>::quiet_NaN();
    T k2 = k * k;
    if (k2 >= T(1)) return std::numeric_limits<T>::infinity(); // K(±1) = ∞
    return carlsonRF(T(0), T(1) - k2, T(1));
}

// Float type: use sangi::ellipticK(k, precision) directly.

// ================================================================
// Complete elliptic integral of the second kind E(k)
// E(k) = R_F(0, 1-k², 1) - (k²/3)·R_D(0, 1-k², 1)
// ================================================================

template<IsNativeFloat T>
[[nodiscard]] T ellipticE(T k) {
    if (std::isnan(k)) return std::numeric_limits<T>::quiet_NaN();
    T k2 = k * k;
    if (k2 >= T(1)) return T(1);  // E(±1) = 1
    T kp2 = T(1) - k2;
    T rf = carlsonRF(T(0), kp2, T(1));
    T rd = carlsonRD(T(0), kp2, T(1));
    return rf - k2 * rd / T(3);
}

// Float type: use sangi::ellipticE(k, precision) directly.

// ================================================================
// Complete elliptic integral of the third kind Π(n, k)
// Π(n, k) = R_F(0, 1-k², 1) + (n/3)·R_J(0, 1-k², 1, 1-n)
// ================================================================

template<IsNativeFloat T>
[[nodiscard]] T ellipticPi(T n, T k) {
    if (std::isnan(n) || std::isnan(k)) return std::numeric_limits<T>::quiet_NaN();
    T k2 = k * k;
    T kp2 = T(1) - k2;
    if (kp2 < T(0)) return std::numeric_limits<T>::quiet_NaN();
    T rf = carlsonRF(T(0), kp2, T(1));
    T rj = carlsonRJ(T(0), kp2, T(1), T(1) - n);
    return rf + n * rj / T(3);
}

// Float type: use sangi::ellipticPi(n, k, precision) directly.

// ================================================================
// Incomplete elliptic integral of the first kind F(φ, k)
// F(φ, k) = sin(φ)·R_F(cos²φ, 1-k²sin²φ, 1)
// ================================================================

template<IsNativeFloat T>
[[nodiscard]] T ellipticF(T phi, T k) {
    if (std::isnan(phi) || std::isnan(k)) return std::numeric_limits<T>::quiet_NaN();
    T s = std::sin(phi);
    T c = std::cos(phi);
    T k2 = k * k;
    return s * carlsonRF(c * c, T(1) - k2 * s * s, T(1));
}

// Float type: use sangi::ellipticF(phi, k, precision) directly.

// ================================================================
// Incomplete elliptic integral of the second kind E(φ, k)
// E(φ, k) = sin(φ)·R_F(...) - (k²/3)sin³φ·R_D(...)
// ================================================================

template<IsNativeFloat T>
[[nodiscard]] T ellipticE(T phi, T k) {
    if (std::isnan(phi) || std::isnan(k)) return std::numeric_limits<T>::quiet_NaN();
    T s = std::sin(phi);
    T c = std::cos(phi);
    T k2 = k * k;
    T s2 = s * s;
    T c2 = c * c;
    T w = T(1) - k2 * s2;
    T rf = carlsonRF(c2, w, T(1));
    T rd = carlsonRD(c2, w, T(1));
    return s * rf - k2 * s * s2 * rd / T(3);
}

// Float type: use sangi::ellipticE(phi, k, precision) directly.

// ================================================================
// Incomplete elliptic integral of the third kind Π(n, φ, k)
// Π(n, φ, k) = sin(φ)·R_F(...) + (n/3)sin³φ·R_J(...)
// ================================================================

template<IsNativeFloat T>
[[nodiscard]] T ellipticPi(T n, T phi, T k) {
    if (std::isnan(n) || std::isnan(phi) || std::isnan(k))
        return std::numeric_limits<T>::quiet_NaN();
    T s = std::sin(phi);
    T c = std::cos(phi);
    T k2 = k * k;
    T s2 = s * s;
    T c2 = c * c;
    T w = T(1) - k2 * s2;
    T rf = carlsonRF(c2, w, T(1));
    T rj = carlsonRJ(c2, w, T(1), T(1) - n * s2);
    return s * rf + n * s * s2 * rj / T(3);
}

// Float type: use sangi::ellipticPi(n, phi, k, precision) directly.

// ================================================================
// Jacobi elliptic functions — AGM descending Landen transform
// ================================================================

namespace detail {

/// Compute the amplitude am(u, k) via AGM iteration + arcsin inverse transform
template<IsNativeFloat T>
[[nodiscard]] T jacobiAmplitude(T u, T k) {
    // k = 0: am(u, 0) = u
    if (k == T(0)) return u;
    // k = 1: am(u, 1) = gd(u) = 2·arctan(tanh(u/2))
    if (k == T(1)) return T(2) * std::atan(std::tanh(u / T(2)));

    T kp = std::sqrt(T(1) - k * k);

    // AGM forward iteration: a₀=1, b₀=k', c₀=k
    std::vector<T> a_seq, c_seq;
    T a = T(1), b = kp, c = k;
    a_seq.push_back(a);
    c_seq.push_back(c);

    const T eps = std::numeric_limits<T>::epsilon();
    for (int n = 0; n < 100; n++) {
        T a_new = (a + b) / T(2);
        T c_new = (a - b) / T(2);
        T b_new = std::sqrt(a * b);
        a = a_new; b = b_new; c = c_new;
        a_seq.push_back(a);
        c_seq.push_back(c);
        if (std::abs(c) < eps) break;
    }

    int N = static_cast<int>(a_seq.size()) - 1;

    // φ_N = 2^N · a_N · u
    T phi = u * a_seq[N];
    for (int i = 0; i < N; i++) {
        phi *= T(2);
    }

    // Inverse transform: φ_{n-1} = (φ_n + arcsin(c_n/a_n · sin(φ_n))) / 2
    for (int n = N; n >= 1; n--) {
        T sinphi = std::sin(phi);
        T arg = c_seq[n] * sinphi / a_seq[n];
        // Clamp arg (in case numerical error makes |arg| > 1)
        arg = std::max(T(-1), std::min(T(1), arg));
        phi = (phi + std::asin(arg)) / T(2);
    }

    return phi;
}

} // namespace detail

/// sn(u, k) = sin(am(u, k))
template<IsNativeFloat T>
[[nodiscard]] T jacobiSn(T u, T k) {
    if (std::isnan(u) || std::isnan(k)) return std::numeric_limits<T>::quiet_NaN();
    if (u == T(0)) return T(0);
    return std::sin(detail::jacobiAmplitude(u, k));
}

/// cn(u, k) = cos(am(u, k))
template<IsNativeFloat T>
[[nodiscard]] T jacobiCn(T u, T k) {
    if (std::isnan(u) || std::isnan(k)) return std::numeric_limits<T>::quiet_NaN();
    if (u == T(0)) return T(1);
    return std::cos(detail::jacobiAmplitude(u, k));
}

/// dn(u, k) = √(1 - k²·sn²(u, k))
template<IsNativeFloat T>
[[nodiscard]] T jacobiDn(T u, T k) {
    if (std::isnan(u) || std::isnan(k)) return std::numeric_limits<T>::quiet_NaN();
    if (u == T(0)) return T(1);
    T sn = jacobiSn(u, k);
    return std::sqrt(T(1) - k * k * sn * sn);
}

// Float type: use sangi::jacobiSn/Cn/Dn(u, k, precision) directly.

// ================================================================
// Complex Carlson symmetric integrals — internal templates
// ================================================================
// Same iterative algorithm as the real version, but with sqrt, abs replaced by complex ones.
// The convergence test is |d| < eps (scalar comparison).

namespace detail {

/// Helper that produces an integer constant at the appropriate precision
template<typename R>
[[nodiscard]] R ellipConst(int val) {
    if constexpr (IsSangiFloat<R>) {
        Float f(val);
        f.setResultPrecision(Float::defaultPrecision());
        return f;
    } else {
        return R(val);
    }
}

/// Helper that produces a fractional constant at the appropriate precision
template<typename R>
[[nodiscard]] R ellipFrac(int num, int den) {
    if constexpr (IsSangiFloat<R>) {
        Float n(num); n.setResultPrecision(Float::defaultPrecision());
        Float d(den); d.setResultPrecision(Float::defaultPrecision());
        return n / d;
    } else {
        return R(num) / R(den);
    }
}

template<typename R>
[[nodiscard]] Complex<R> carlsonRC_complex(Complex<R> x, Complex<R> y, R eps, int max_iter) {
    using C = Complex<R>;
    R _2 = ellipConst<R>(2), _3 = ellipConst<R>(3), _4 = ellipConst<R>(4);
    for (int iter = 0; iter < max_iter; iter++) {
        C lam = _2 * sangi::sqrt(x * y) + y;
        x = (x + lam) / _4;
        y = (y + lam) / _4;

        C A = (x + y + y) / _3;
        C s = (y - A) / A;
        if (sangi::abs(s) < eps) break;
    }
    C A = (x + y + y) / _3;
    C s = (y - A) / A;
    C s2 = s * s;
    C poly = C(ellipConst<R>(1)) + s2 * (C(ellipFrac<R>(3, 10))
             + s * (C(ellipFrac<R>(1, 7))
             + s * (C(ellipFrac<R>(3, 8))
             + s * C(ellipFrac<R>(9, 22)))));
    return poly / sangi::sqrt(A);
}

template<typename R>
[[nodiscard]] Complex<R> carlsonRF_complex(Complex<R> x, Complex<R> y, Complex<R> z,
                                            R eps, int max_iter) {
    using C = Complex<R>;
    R _3 = ellipConst<R>(3), _4 = ellipConst<R>(4);
    for (int iter = 0; iter < max_iter; iter++) {
        C sqx = sangi::sqrt(x), sqy = sangi::sqrt(y), sqz = sangi::sqrt(z);
        C lam = sqx * sqy + sqy * sqz + sqz * sqx;
        x = (x + lam) / _4;
        y = (y + lam) / _4;
        z = (z + lam) / _4;

        C A = (x + y + z) / _3;
        R mx = std::max({sangi::abs((A - x) / A), sangi::abs((A - y) / A), sangi::abs((A - z) / A)});
        if (mx < eps) break;
    }
    C A = (x + y + z) / _3;
    C X = (A - x) / A;
    C Y = (A - y) / A;
    C Z = -(X + Y);
    C E2 = X * Y - Z * Z;
    C E3 = X * Y * Z;
    R _1 = ellipConst<R>(1);
    C poly = C(_1) - E2 / ellipConst<R>(10) + E3 / ellipConst<R>(14)
             + E2 * E2 / ellipConst<R>(24)
             - C(_3) * E2 * E3 / ellipConst<R>(44);
    return poly / sangi::sqrt(A);
}

template<typename R>
[[nodiscard]] Complex<R> carlsonRD_complex(Complex<R> x, Complex<R> y, Complex<R> z,
                                            R eps, int max_iter) {
    using C = Complex<R>;
    R _3 = ellipConst<R>(3), _4 = ellipConst<R>(4), _5 = ellipConst<R>(5);
    C sigma(ellipConst<R>(0));
    R fac = ellipConst<R>(1);

    for (int iter = 0; iter < max_iter; iter++) {
        C sqx = sangi::sqrt(x), sqy = sangi::sqrt(y), sqz = sangi::sqrt(z);
        C lam = sqx * sqy + sqy * sqz + sqz * sqx;
        sigma = sigma + C(fac) / (sqz * (z + lam));
        fac = fac / _4;
        x = (x + lam) / _4;
        y = (y + lam) / _4;
        z = (z + lam) / _4;

        C A = (x + y + C(_3) * z) / _5;
        R mx = std::max({sangi::abs((A - x) / A), sangi::abs((A - y) / A), sangi::abs((A - z) / A)});
        if (mx < eps) break;
    }

    C A = (x + y + C(_3) * z) / _5;
    C dx = (A - x) / A;
    C dy = (A - y) / A;
    C dz = (A - z) / A;
    C ea = dx * dy;
    C eb = dz * dz;
    C ec = ea - eb;
    C ed = ea - C(ellipConst<R>(6)) * eb;
    C ee = ed + ec + ec;
    C poly = C(ellipConst<R>(1))
        + ed * (C(ellipFrac<R>(-3, 14)) + C(ellipFrac<R>(9, 88)) * ed
                - C(ellipFrac<R>(9, 52)) * dz * ee)
        + dz * (C(ellipFrac<R>(1, 6)) * ee
                + dz * (C(ellipFrac<R>(-9, 22)) * ec
                        + dz * C(ellipFrac<R>(3, 26)) * ea));

    return C(_3) * sigma + C(fac) * poly / (A * sangi::sqrt(A));
}

template<typename R>
[[nodiscard]] Complex<R> carlsonRJ_complex(Complex<R> x, Complex<R> y, Complex<R> z, Complex<R> p,
                                            R eps, int max_iter) {
    using C = Complex<R>;
    R _3 = ellipConst<R>(3), _4 = ellipConst<R>(4), _5 = ellipConst<R>(5);
    C sigma(ellipConst<R>(0));
    R fac = ellipConst<R>(1);

    for (int iter = 0; iter < max_iter; iter++) {
        C sqx = sangi::sqrt(x), sqy = sangi::sqrt(y), sqz = sangi::sqrt(z), sqp = sangi::sqrt(p);
        C lam = sqx * sqy + sqy * sqz + sqz * sqx;
        C alpha = p * (sqx + sqy + sqz) + sqx * sqy * sqz;
        alpha = alpha * alpha;
        C beta = p * (p + lam) * (p + lam);
        sigma = sigma + C(fac) * carlsonRC_complex<R>(alpha, beta, eps, max_iter);
        fac = fac / _4;
        x = (x + lam) / _4;
        y = (y + lam) / _4;
        z = (z + lam) / _4;
        p = (p + lam) / _4;

        C A = (x + y + z + p + p) / _5;
        R mx = std::max({sangi::abs((A - x) / A), sangi::abs((A - y) / A),
                         sangi::abs((A - z) / A), sangi::abs((A - p) / A)});
        if (mx < eps) break;
    }

    C A = (x + y + z + p + p) / _5;
    C dx = (A - x) / A;
    C dy = (A - y) / A;
    C dz = (A - z) / A;
    C dp = (A - p) / A;
    C ea = dx * (dy + dz) + dy * dz;
    C eb = dx * dy * dz;
    C ec = dp * dp;
    C ed = ea - C(_3) * ec;
    C ee = eb + C(ellipConst<R>(2)) * dp * (ea - ec);
    C poly = C(ellipConst<R>(1))
        + ed * (C(ellipFrac<R>(-3, 14)) + C(ellipFrac<R>(9, 88)) * ed
                - C(ellipFrac<R>(9, 52)) * ee)
        + eb * (C(ellipFrac<R>(1, 6)) + dp * (C(ellipFrac<R>(-6, 22))
                + dp * C(ellipFrac<R>(3, 26))))
        + dp * ea * (C(ellipFrac<R>(1, 3)) - dp * C(ellipFrac<R>(3, 22)))
        - C(ellipFrac<R>(1, 3)) * dp * ec;

    return C(_3) * sigma + C(fac) * poly / (A * sangi::sqrt(A));
}

} // namespace detail

// ================================================================
// Complex<double> elliptic integral overloads
// ================================================================

// Carlson symmetric integrals
[[nodiscard]] inline Complex<double> carlsonRC(const Complex<double>& x, const Complex<double>& y) {
    return detail::carlsonRC_complex<double>(x, y,
        std::numeric_limits<double>::epsilon(), 100);
}

[[nodiscard]] inline Complex<double> carlsonRF(const Complex<double>& x,
                                                const Complex<double>& y,
                                                const Complex<double>& z) {
    return detail::carlsonRF_complex<double>(x, y, z,
        std::numeric_limits<double>::epsilon(), 100);
}

[[nodiscard]] inline Complex<double> carlsonRD(const Complex<double>& x,
                                                const Complex<double>& y,
                                                const Complex<double>& z) {
    return detail::carlsonRD_complex<double>(x, y, z,
        std::numeric_limits<double>::epsilon(), 100);
}

[[nodiscard]] inline Complex<double> carlsonRJ(const Complex<double>& x,
                                                const Complex<double>& y,
                                                const Complex<double>& z,
                                                const Complex<double>& p) {
    return detail::carlsonRJ_complex<double>(x, y, z, p,
        std::numeric_limits<double>::epsilon(), 100);
}

// Complete elliptic integrals
[[nodiscard]] inline Complex<double> ellipticK(const Complex<double>& k) {
    using C = Complex<double>;
    C k2 = k * k;
    return carlsonRF(C(0.0), C(1.0) - k2, C(1.0));
}

[[nodiscard]] inline Complex<double> ellipticE(const Complex<double>& k) {
    using C = Complex<double>;
    C k2 = k * k;
    C kp2 = C(1.0) - k2;
    return carlsonRF(C(0.0), kp2, C(1.0)) - k2 * carlsonRD(C(0.0), kp2, C(1.0)) / C(3.0);
}

[[nodiscard]] inline Complex<double> ellipticPi(const Complex<double>& n,
                                                 const Complex<double>& k) {
    using C = Complex<double>;
    C k2 = k * k;
    C kp2 = C(1.0) - k2;
    return carlsonRF(C(0.0), kp2, C(1.0))
         + n * carlsonRJ(C(0.0), kp2, C(1.0), C(1.0) - n) / C(3.0);
}

// ================================================================
// Complex<Float> elliptic integral overloads
// ================================================================

[[nodiscard]] inline Complex<Float> ellipticK(const Complex<Float>& k, int precision) {
    using C = Complex<Float>;
    int wp = precision + 20;
    Float::PrecisionScope _ps(wp);   // compute the exact÷exact (/4,/6,…) inside Carlson RF to wp digits
    // PrecisionGuard removed: setResultPrecision(wp) on kw propagates req=wp.
    C kw = k;
    kw.re.setResultPrecision(wp);
    kw.im.setResultPrecision(wp);

    Float one(1);  // setResultPrecision not needed: PrecisionScope promotes the constant exact÷exact to wp
    C k2 = kw * kw;
    C kp2 = C(one) - k2;
    C result = detail::carlsonRF_complex<Float>(C(), kp2, C(one),
                                                 Float::epsilon(wp), 10 * wp);
    result.re.setPrecision(precision);
    result.im.setPrecision(precision);
    return result;
}

[[nodiscard]] inline Complex<Float> ellipticE(const Complex<Float>& k, int precision) {
    using C = Complex<Float>;
    int wp = precision + 20;
    Float::PrecisionScope _ps(wp);   // compute the exact÷exact inside Carlson RF/RD to wp digits
    // PrecisionGuard removed: setResultPrecision(wp) on kw propagates req=wp.
    C kw = k;
    kw.re.setResultPrecision(wp);
    kw.im.setResultPrecision(wp);

    Float one(1);  // setResultPrecision not needed: PrecisionScope promotes the constant exact÷exact to wp
    Float three(3);
    C k2 = kw * kw;
    C kp2 = C(one) - k2;
    C zero;
    C rf = detail::carlsonRF_complex<Float>(zero, kp2, C(one), Float::epsilon(wp), 10 * wp);
    C rd = detail::carlsonRD_complex<Float>(zero, kp2, C(one), Float::epsilon(wp), 10 * wp);
    C result = rf - k2 * rd / three;
    result.re.setPrecision(precision);
    result.im.setPrecision(precision);
    return result;
}

[[nodiscard]] inline Complex<Float> ellipticPi(const Complex<Float>& n,
                                                const Complex<Float>& k, int precision) {
    using C = Complex<Float>;
    int wp = precision + 20;
    Float::PrecisionScope _ps(wp);   // compute the exact÷exact inside Carlson RF/RJ to wp digits
    // PrecisionGuard removed: setResultPrecision(wp) on nw/kw propagates req=wp.
    C nw = n; nw.re.setResultPrecision(wp); nw.im.setResultPrecision(wp);
    C kw = k; kw.re.setResultPrecision(wp); kw.im.setResultPrecision(wp);

    Float one(1);  // setResultPrecision not needed: PrecisionScope promotes the constant exact÷exact to wp
    Float three(3);
    C k2 = kw * kw;
    C kp2 = C(one) - k2;
    C zero;
    C rf = detail::carlsonRF_complex<Float>(zero, kp2, C(one), Float::epsilon(wp), 10 * wp);
    C rj = detail::carlsonRJ_complex<Float>(zero, kp2, C(one), C(one) - nw,
                                             Float::epsilon(wp), 10 * wp);
    C result = rf + nw * rj / three;
    result.re.setPrecision(precision);
    result.im.setPrecision(precision);
    return result;
}

} // namespace special
} // namespace sangi

#endif // SANGI_SPECIAL_ELLIPTIC_HPP