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

// theta.hpp
// Jacobi theta function templates
//
// Functions provided:
//   jacobiTheta1(z, q)  — θ₁(z, q) = 2 Σ_{n=0}^{∞} (-1)^n q^{(n+1/2)²} sin((2n+1)z)
//   jacobiTheta2(z, q)  — θ₂(z, q) = 2 Σ_{n=0}^{∞} q^{(n+1/2)²} cos((2n+1)z)
//   jacobiTheta3(z, q)  — θ₃(z, q) = 1 + 2 Σ_{n=1}^{∞} q^{n²} cos(2nz)
//   jacobiTheta4(z, q)  — θ₄(z, q) = 1 + 2 Σ_{n=1}^{∞} (-1)^n q^{n²} cos(2nz)
//   qFromTau(tau)       — q = exp(iπτ) (nome from half-period ratio)
//
// Here q = nome (|q| < 1). Relation to τ (half-period ratio): q = exp(iπτ)
//
// Supported types:
//   double                 — direct evaluation of the q-series
//   Complex<double>        — complex z, real q
//   Complex<Float>         — arbitrary precision (real q or Complex<Float> q)
//
// SPECIAL-X4 (TODO_BEYOND_MPFR): the Complex<Float> q overload and
// qFromTau(Complex<Float>) were added on 2026-05-08, fulfilling the TODO API
// `jacobiThetaN(Complex<Float> z, Complex<Float> q, int precision)`.
//
// Future work (stocking phase): AGM acceleration (real z, |q|<1),
// modular reduction of τ (convergence acceleration for small Im(τ)), a direct
// Float version (fast path for real arguments that bypasses Complex).
//
// References:
//   DLMF §20.2, A&S §16.27
//   NIST Digital Library of Mathematical Functions: https://dlmf.nist.gov/20

#ifndef SANGI_SPECIAL_THETA_HPP
#define SANGI_SPECIAL_THETA_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 <numbers>

namespace sangi {
namespace special {

// ================================================================
// Jacobi theta functions (double, real arguments)
// ================================================================
// q = nome, requires |q| < 1. The series diverges for q ≥ 1.

/// θ₁(z, q) = 2 Σ_{n=0}^{∞} (-1)^n q^{(n+1/2)²} sin((2n+1)z)
template<IsNativeFloat T>
[[nodiscard]] T jacobiTheta1(T z, T q) {
    if (std::isnan(z) || std::isnan(q)) return std::numeric_limits<T>::quiet_NaN();
    if (std::abs(q) >= T(1)) return std::numeric_limits<T>::quiet_NaN();
    if (q == T(0)) return T(0);

    T sum = T(0);
    T q_pow = std::sqrt(q);  // q^{1/4} · q^{0} → q^{(0+1/2)²} = q^{1/4}
    T q_sq = q * q;
    T sign = T(1);

    for (int n = 0; n < 300; n++) {
        // q^{(n+1/2)²}
        T qn;
        if (n == 0) {
            qn = std::sqrt(std::sqrt(q));  // q^{1/4}
        } else {
            // q^{(n+1/2)²} = q^{n²+n+1/4} = q^{n²} · q^{n} · q^{1/4}
            // Recurrence: q^{((n+1)+1/2)²} / q^{(n+1/2)²} = q^{2n+2}
            qn = q_pow * q * q;  // q^{(n+1/2)²} = previous term · q^{2n}
            // Direct computation is more stable
            qn = std::pow(q, (n + T(0.5)) * (n + T(0.5)));
        }

        T term = sign * qn * std::sin((2 * n + 1) * z);
        sum += term;

        if (n >= 2 && std::abs(term) < std::numeric_limits<T>::epsilon() * std::abs(sum))
            break;

        sign = -sign;
    }

    return T(2) * sum;
}

/// θ₂(z, q) = 2 Σ_{n=0}^{∞} q^{(n+1/2)²} cos((2n+1)z)
template<IsNativeFloat T>
[[nodiscard]] T jacobiTheta2(T z, T q) {
    if (std::isnan(z) || std::isnan(q)) return std::numeric_limits<T>::quiet_NaN();
    if (std::abs(q) >= T(1)) return std::numeric_limits<T>::quiet_NaN();
    if (q == T(0)) return T(0);

    T sum = T(0);

    for (int n = 0; n < 300; n++) {
        T qn = std::pow(q, (n + T(0.5)) * (n + T(0.5)));
        T term = qn * std::cos((2 * n + 1) * z);
        sum += term;

        if (n >= 2 && std::abs(term) < std::numeric_limits<T>::epsilon() * std::abs(sum))
            break;
    }

    return T(2) * sum;
}

/// θ₃(z, q) = 1 + 2 Σ_{n=1}^{∞} q^{n²} cos(2nz)
template<IsNativeFloat T>
[[nodiscard]] T jacobiTheta3(T z, T q) {
    if (std::isnan(z) || std::isnan(q)) return std::numeric_limits<T>::quiet_NaN();
    if (std::abs(q) >= T(1)) return std::numeric_limits<T>::quiet_NaN();
    if (q == T(0)) return T(1);

    T sum = T(0);
    T q_pow = q;  // q^{n²}: n=1 → q

    for (int n = 1; n < 300; n++) {
        T qn = std::pow(q, T(n) * T(n));
        T term = qn * std::cos(T(2 * n) * z);
        sum += term;

        if (n >= 2 && std::abs(term) < std::numeric_limits<T>::epsilon() * (T(1) + std::abs(sum)))
            break;
    }

    return T(1) + T(2) * sum;
}

/// θ₄(z, q) = 1 + 2 Σ_{n=1}^{∞} (-1)^n q^{n²} cos(2nz)
template<IsNativeFloat T>
[[nodiscard]] T jacobiTheta4(T z, T q) {
    if (std::isnan(z) || std::isnan(q)) return std::numeric_limits<T>::quiet_NaN();
    if (std::abs(q) >= T(1)) return std::numeric_limits<T>::quiet_NaN();
    if (q == T(0)) return T(1);

    T sum = T(0);
    T sign = T(-1);

    for (int n = 1; n < 300; n++) {
        T qn = std::pow(q, T(n) * T(n));
        T term = sign * qn * std::cos(T(2 * n) * z);
        sum += term;

        if (n >= 2 && std::abs(term) < std::numeric_limits<T>::epsilon() * (T(1) + std::abs(sum)))
            break;

        sign = -sign;
    }

    return T(1) + T(2) * sum;
}

// ================================================================
// Jacobi theta functions (Complex<double>)
// ================================================================

namespace detail {

/// θ₁(z, q) — complex version
template<typename R>
[[nodiscard]] Complex<R> jacobiTheta1_complex(const Complex<R>& z, R q,
                                               R eps, int max_iter) {
    using C = Complex<R>;
    if (q == R(0)) return C(R(0));

    C sum(R(0));
    R sign = R(1);

    for (int n = 0; n < max_iter; n++) {
        R nh = R(n) + R(1) / R(2);
        R qn;
        if constexpr (IsSangiFloat<R>)
            qn = sangi::pow(q, nh * nh);
        else
            qn = std::pow(q, nh * nh);
        C arg = C(R(2 * n + 1)) * z;
        C term = C(sign * qn) * sangi::sin(arg);
        sum = sum + term;

        if (n >= 2 && sangi::abs(term) < eps * sangi::abs(sum))
            break;

        sign = -sign;
    }

    return C(R(2)) * sum;
}

/// θ₂(z, q) — complex version
template<typename R>
[[nodiscard]] Complex<R> jacobiTheta2_complex(const Complex<R>& z, R q,
                                               R eps, int max_iter) {
    using C = Complex<R>;
    if (q == R(0)) return C(R(0));

    C sum(R(0));

    for (int n = 0; n < max_iter; n++) {
        R nh = R(n) + R(1) / R(2);
        R qn;
        if constexpr (IsSangiFloat<R>)
            qn = sangi::pow(q, nh * nh);
        else
            qn = std::pow(q, nh * nh);
        C arg = C(R(2 * n + 1)) * z;
        C term = C(qn) * sangi::cos(arg);
        sum = sum + term;

        if (n >= 2 && sangi::abs(term) < eps * sangi::abs(sum))
            break;
    }

    return C(R(2)) * sum;
}

/// θ₃(z, q) — complex version
template<typename R>
[[nodiscard]] Complex<R> jacobiTheta3_complex(const Complex<R>& z, R q,
                                               R eps, int max_iter) {
    using C = Complex<R>;
    if (q == R(0)) return C(R(1));

    C sum(R(0));

    for (int n = 1; n < max_iter; n++) {
        R nn = R(n);
        R qn;
        if constexpr (IsSangiFloat<R>)
            qn = sangi::pow(q, nn * nn);
        else
            qn = std::pow(q, nn * nn);
        C arg = C(R(2 * n)) * z;
        C term = C(qn) * sangi::cos(arg);
        sum = sum + term;

        if (n >= 2 && sangi::abs(term) < eps * (R(1) + sangi::abs(sum)))
            break;
    }

    return C(R(1)) + C(R(2)) * sum;
}

/// θ₄(z, q) — complex version
template<typename R>
[[nodiscard]] Complex<R> jacobiTheta4_complex(const Complex<R>& z, R q,
                                               R eps, int max_iter) {
    using C = Complex<R>;
    if (q == R(0)) return C(R(1));

    C sum(R(0));
    R sign = R(-1);

    for (int n = 1; n < max_iter; n++) {
        R nn = R(n);
        R qn;
        if constexpr (IsSangiFloat<R>)
            qn = sangi::pow(q, nn * nn);
        else
            qn = std::pow(q, nn * nn);
        C arg = C(R(2 * n)) * z;
        C term = C(sign * qn) * sangi::cos(arg);
        sum = sum + term;

        if (n >= 2 && sangi::abs(term) < eps * (R(1) + sangi::abs(sum)))
            break;

        sign = -sign;
    }

    return C(R(1)) + C(R(2)) * sum;
}

// ----------------------------------------------------------------
// Complex q (= arbitrary upper half-plane τ) version
// ----------------------------------------------------------------
// When q is complex, q^{(n+1/2)²} and q^{n²} are evaluated as
//   q^k = exp(k · log q)  (k a non-integer real)
//   q^k = pow(q, int k)   (k an integer)
// Integer powers are stable (pow(Complex<R>, int) is binary iteration),
// while half-integer powers cost one exp(log).

/// θ₁(z, q) — complex z + complex q. Decomposed into q^{1/4} and q^{n(n+1)} for iteration.
template<typename R>
[[nodiscard]] Complex<R> jacobiTheta1_cqcz(const Complex<R>& z, const Complex<R>& q,
                                           R eps, int max_iter) {
    using C = Complex<R>;
    if (q.re == R(0) && q.im == R(0)) return C(R(0));

    // q^{1/4} = exp((1/4) log q)
    C log_q = sangi::log(q);
    C q_quarter = sangi::exp(C(R(1) / R(4)) * log_q);

    C sum(R(0));
    R sign = R(1);

    for (int n = 0; n < max_iter; n++) {
        // q^{(n+1/2)²} = q^{1/4} · q^{n(n+1)} (integer power)
        C qn = q_quarter * sangi::pow(q, n * (n + 1));
        C arg = C(R(2 * n + 1)) * z;
        C term = C(sign) * qn * sangi::sin(arg);
        sum = sum + term;

        if (n >= 2 && sangi::abs(term) < eps * sangi::abs(sum))
            break;

        sign = -sign;
    }

    return C(R(2)) * sum;
}

/// θ₂(z, q) — complex z + complex q
template<typename R>
[[nodiscard]] Complex<R> jacobiTheta2_cqcz(const Complex<R>& z, const Complex<R>& q,
                                           R eps, int max_iter) {
    using C = Complex<R>;
    if (q.re == R(0) && q.im == R(0)) return C(R(0));

    C log_q = sangi::log(q);
    C q_quarter = sangi::exp(C(R(1) / R(4)) * log_q);

    C sum(R(0));

    for (int n = 0; n < max_iter; n++) {
        C qn = q_quarter * sangi::pow(q, n * (n + 1));
        C arg = C(R(2 * n + 1)) * z;
        C term = qn * sangi::cos(arg);
        sum = sum + term;

        if (n >= 2 && sangi::abs(term) < eps * sangi::abs(sum))
            break;
    }

    return C(R(2)) * sum;
}

/// θ₃(z, q) — complex z + complex q. q^{n²} uses only integer powers.
template<typename R>
[[nodiscard]] Complex<R> jacobiTheta3_cqcz(const Complex<R>& z, const Complex<R>& q,
                                           R eps, int max_iter) {
    using C = Complex<R>;
    if (q.re == R(0) && q.im == R(0)) return C(R(1));

    C sum(R(0));

    for (int n = 1; n < max_iter; n++) {
        C qn = sangi::pow(q, n * n);
        C arg = C(R(2 * n)) * z;
        C term = qn * sangi::cos(arg);
        sum = sum + term;

        if (n >= 2 && sangi::abs(term) < eps * (R(1) + sangi::abs(sum)))
            break;
    }

    return C(R(1)) + C(R(2)) * sum;
}

/// θ₄(z, q) — complex z + complex q
template<typename R>
[[nodiscard]] Complex<R> jacobiTheta4_cqcz(const Complex<R>& z, const Complex<R>& q,
                                           R eps, int max_iter) {
    using C = Complex<R>;
    if (q.re == R(0) && q.im == R(0)) return C(R(1));

    C sum(R(0));
    R sign = R(-1);

    for (int n = 1; n < max_iter; n++) {
        C qn = sangi::pow(q, n * n);
        C arg = C(R(2 * n)) * z;
        C term = C(sign) * qn * sangi::cos(arg);
        sum = sum + term;

        if (n >= 2 && sangi::abs(term) < eps * (R(1) + sangi::abs(sum)))
            break;

        sign = -sign;
    }

    return C(R(1)) + C(R(2)) * sum;
}

} // namespace detail

// ----------------------------------------------------------------
// Direct Float version — real z + real q → real Float (fast-path, SciPy-compatible API)
// ----------------------------------------------------------------
//
// For real-only arguments, compute directly in real Float arithmetic, bypassing Complex<Float>.
// About 2x faster since it does not keep both real and imaginary parts in intermediate computation.

namespace detail {

/// θ_1(z, q) — Float real fast-path
inline Float jacobiTheta1_real_(const Float& z, const Float& q, int wp) {
    if (q.isNaN() || z.isNaN()) return Float::nan();
    if (q == Float(0)) return Float(0);

    // q^{1/4} via exp((1/4) log q). For q < 0 the log is NaN, but since we require |q|<1,
    // allowing q < 0 would need separate handling. In practice only q ∈ (0,1).
    Float zw = z; zw.setResultPrecision(wp);
    Float qw = q; qw.setResultPrecision(wp);

    Float log_q = sangi::log(qw);
    Float quarter = Float(1) / Float(4); quarter.setResultPrecision(wp);
    Float q_quarter = sangi::exp(quarter * log_q);

    Float sum(Float(0)); sum.setResultPrecision(wp);
    int sign = 1;
    Float eps = Float::epsilon(wp);

    for (int n = 0; n < 10 * wp; n++) {
        // q^{(n+1/2)²} = q^{1/4} · q^{n(n+1)}
        Float qn;
        if (n == 0) {
            qn = q_quarter;
        } else {
            qn = q_quarter * sangi::pow(qw, n * (n + 1));
        }
        Float arg = Float(2 * n + 1) * zw;
        Float term = Float(sign) * qn * sangi::sin(arg);
        sum = sum + term;

        if (n >= 2 && sangi::abs(term) < eps * sangi::abs(sum)) break;
        sign = -sign;
    }

    return Float(2) * sum;
}

inline Float jacobiTheta2_real_(const Float& z, const Float& q, int wp) {
    if (q.isNaN() || z.isNaN()) return Float::nan();
    if (q == Float(0)) return Float(0);

    Float zw = z; zw.setResultPrecision(wp);
    Float qw = q; qw.setResultPrecision(wp);

    Float log_q = sangi::log(qw);
    Float quarter = Float(1) / Float(4); quarter.setResultPrecision(wp);
    Float q_quarter = sangi::exp(quarter * log_q);

    Float sum(Float(0)); sum.setResultPrecision(wp);
    Float eps = Float::epsilon(wp);

    for (int n = 0; n < 10 * wp; n++) {
        Float qn;
        if (n == 0) {
            qn = q_quarter;
        } else {
            qn = q_quarter * sangi::pow(qw, n * (n + 1));
        }
        Float arg = Float(2 * n + 1) * zw;
        Float term = qn * sangi::cos(arg);
        sum = sum + term;

        if (n >= 2 && sangi::abs(term) < eps * sangi::abs(sum)) break;
    }

    return Float(2) * sum;
}

inline Float jacobiTheta3_real_(const Float& z, const Float& q, int wp) {
    if (q.isNaN() || z.isNaN()) return Float::nan();
    if (q == Float(0)) return Float(1);

    Float zw = z; zw.setResultPrecision(wp);
    Float qw = q; qw.setResultPrecision(wp);

    Float sum(Float(0)); sum.setResultPrecision(wp);
    Float eps = Float::epsilon(wp);

    for (int n = 1; n < 10 * wp; n++) {
        Float qn = sangi::pow(qw, n * n);
        Float arg = Float(2 * n) * zw;
        Float term = qn * sangi::cos(arg);
        sum = sum + term;

        if (n >= 2 && sangi::abs(term) < eps * (Float(1) + sangi::abs(sum))) break;
    }

    return Float(1) + Float(2) * sum;
}

inline Float jacobiTheta4_real_(const Float& z, const Float& q, int wp) {
    if (q.isNaN() || z.isNaN()) return Float::nan();
    if (q == Float(0)) return Float(1);

    Float zw = z; zw.setResultPrecision(wp);
    Float qw = q; qw.setResultPrecision(wp);

    Float sum(Float(0)); sum.setResultPrecision(wp);
    int sign = -1;
    Float eps = Float::epsilon(wp);

    for (int n = 1; n < 10 * wp; n++) {
        Float qn = sangi::pow(qw, n * n);
        Float arg = Float(2 * n) * zw;
        Float term = Float(sign) * qn * sangi::cos(arg);
        sum = sum + term;

        if (n >= 2 && sangi::abs(term) < eps * (Float(1) + sangi::abs(sum))) break;
        sign = -sign;
    }

    return Float(1) + Float(2) * sum;
}

} // namespace detail

/// θ_n(z, q) — Float real-argument fast-path.
/// Assumes real q ∈ (0,1) with |q| < 1 (negative q is future work).
///
/// wp = precision + 20: standard guard bits.
/// The old version used wp = 4·precision to avoid sparse-mantissa precision loss, but
/// after the BUG_FLOAT_SPARSE_MANTISSA_PRECISION_20260509 H4 fix (commit 9e79cd5) resolved
/// the root cause on the Float library side (misjudged MSB distance in addUnsigned/subtractUnsigned),
/// the workaround is no longer needed.
namespace detail {
inline int floatRealWp_(int precision) {
    return precision + 20;
}
} // namespace detail

[[nodiscard]] inline Float jacobiTheta1(const Float& z, const Float& q, int precision) {
    Float::PrecisionScope _ps(precision + 64);   // compute exact÷exact inside the _real_ helper at working precision
    if (sangi::abs(q) >= Float(1)) return Float::nan();
    int wp = detail::floatRealWp_(precision);
    Float result = detail::jacobiTheta1_real_(z, q, wp);
    result.setPrecision(precision);
    return result;
}
[[nodiscard]] inline Float jacobiTheta1(const Float& z, const Float& q) {
    return jacobiTheta1(z, q, Float::defaultPrecision());
}

[[nodiscard]] inline Float jacobiTheta2(const Float& z, const Float& q, int precision) {
    Float::PrecisionScope _ps(precision + 64);   // compute exact÷exact inside the _real_ helper at working precision
    if (sangi::abs(q) >= Float(1)) return Float::nan();
    int wp = detail::floatRealWp_(precision);
    Float result = detail::jacobiTheta2_real_(z, q, wp);
    result.setPrecision(precision);
    return result;
}
[[nodiscard]] inline Float jacobiTheta2(const Float& z, const Float& q) {
    return jacobiTheta2(z, q, Float::defaultPrecision());
}

[[nodiscard]] inline Float jacobiTheta3(const Float& z, const Float& q, int precision) {
    if (sangi::abs(q) >= Float(1)) return Float::nan();
    int wp = detail::floatRealWp_(precision);
    Float result = detail::jacobiTheta3_real_(z, q, wp);
    result.setPrecision(precision);
    return result;
}
[[nodiscard]] inline Float jacobiTheta3(const Float& z, const Float& q) {
    return jacobiTheta3(z, q, Float::defaultPrecision());
}

[[nodiscard]] inline Float jacobiTheta4(const Float& z, const Float& q, int precision) {
    if (sangi::abs(q) >= Float(1)) return Float::nan();
    int wp = detail::floatRealWp_(precision);
    Float result = detail::jacobiTheta4_real_(z, q, wp);
    result.setPrecision(precision);
    return result;
}
[[nodiscard]] inline Float jacobiTheta4(const Float& z, const Float& q) {
    return jacobiTheta4(z, q, Float::defaultPrecision());
}

// ---- Complex<double> overloads ----

[[nodiscard]] inline Complex<double> jacobiTheta1(const Complex<double>& z, double q) {
    if (std::isnan(q) || std::abs(q) >= 1.0) return Complex<double>(std::numeric_limits<double>::quiet_NaN());
    return detail::jacobiTheta1_complex<double>(z, q,
        std::numeric_limits<double>::epsilon(), 300);
}

[[nodiscard]] inline Complex<double> jacobiTheta2(const Complex<double>& z, double q) {
    if (std::isnan(q) || std::abs(q) >= 1.0) return Complex<double>(std::numeric_limits<double>::quiet_NaN());
    return detail::jacobiTheta2_complex<double>(z, q,
        std::numeric_limits<double>::epsilon(), 300);
}

[[nodiscard]] inline Complex<double> jacobiTheta3(const Complex<double>& z, double q) {
    if (std::isnan(q) || std::abs(q) >= 1.0) return Complex<double>(std::numeric_limits<double>::quiet_NaN());
    return detail::jacobiTheta3_complex<double>(z, q,
        std::numeric_limits<double>::epsilon(), 300);
}

[[nodiscard]] inline Complex<double> jacobiTheta4(const Complex<double>& z, double q) {
    if (std::isnan(q) || std::abs(q) >= 1.0) return Complex<double>(std::numeric_limits<double>::quiet_NaN());
    return detail::jacobiTheta4_complex<double>(z, q,
        std::numeric_limits<double>::epsilon(), 300);
}

// ---- Complex<double> + Complex<double> q (complex nome) ----
// Requires |q| < 1 (q is the image of upper half-plane τ: q = exp(iπτ)).

namespace detail {
inline bool complex_q_valid_double(const Complex<double>& q) {
    if (std::isnan(q.re) || std::isnan(q.im)) return false;
    double mod_sq = q.re * q.re + q.im * q.im;
    return mod_sq < 1.0;
}
} // namespace detail

[[nodiscard]] inline Complex<double> jacobiTheta1(const Complex<double>& z, const Complex<double>& q) {
    if (!detail::complex_q_valid_double(q))
        return Complex<double>(std::numeric_limits<double>::quiet_NaN());
    return detail::jacobiTheta1_cqcz<double>(z, q,
        std::numeric_limits<double>::epsilon(), 300);
}

[[nodiscard]] inline Complex<double> jacobiTheta2(const Complex<double>& z, const Complex<double>& q) {
    if (!detail::complex_q_valid_double(q))
        return Complex<double>(std::numeric_limits<double>::quiet_NaN());
    return detail::jacobiTheta2_cqcz<double>(z, q,
        std::numeric_limits<double>::epsilon(), 300);
}

[[nodiscard]] inline Complex<double> jacobiTheta3(const Complex<double>& z, const Complex<double>& q) {
    if (!detail::complex_q_valid_double(q))
        return Complex<double>(std::numeric_limits<double>::quiet_NaN());
    return detail::jacobiTheta3_cqcz<double>(z, q,
        std::numeric_limits<double>::epsilon(), 300);
}

[[nodiscard]] inline Complex<double> jacobiTheta4(const Complex<double>& z, const Complex<double>& q) {
    if (!detail::complex_q_valid_double(q))
        return Complex<double>(std::numeric_limits<double>::quiet_NaN());
    return detail::jacobiTheta4_cqcz<double>(z, q,
        std::numeric_limits<double>::epsilon(), 300);
}

// ---- Complex<Float> overloads ----

[[nodiscard]] inline Complex<Float> jacobiTheta1(const Complex<Float>& z, const Float& q, int precision) {
    using C = Complex<Float>;
    int wp = precision + 20;
    Float::PrecisionScope _ps(wp);   // compute exact÷exact inside the q-series at wp digits
    // PrecisionGuard removed: setResultPrecision(wp) on zw / qw propagates req=wp.
    C zw = z;
    zw.re.setResultPrecision(wp);
    zw.im.setResultPrecision(wp);
    Float qw = q;
    qw.setResultPrecision(wp);
    C result = detail::jacobiTheta1_complex<Float>(zw, qw,
        Float::epsilon(wp), 10 * wp);
    result.re.setPrecision(precision);
    result.im.setPrecision(precision);
    return result;
}

[[nodiscard]] inline Complex<Float> jacobiTheta1(const Complex<Float>& z, const Float& q) {
    return jacobiTheta1(z, q, Float::defaultPrecision());
}

[[nodiscard]] inline Complex<Float> jacobiTheta2(const Complex<Float>& z, const Float& q, int precision) {
    using C = Complex<Float>;
    int wp = precision + 20;
    Float::PrecisionScope _ps(wp);   // compute exact÷exact inside the q-series at wp digits
    // PrecisionGuard removed: setResultPrecision(wp) on zw / qw propagates req=wp.
    C zw = z;
    zw.re.setResultPrecision(wp);
    zw.im.setResultPrecision(wp);
    Float qw = q;
    qw.setResultPrecision(wp);
    C result = detail::jacobiTheta2_complex<Float>(zw, qw,
        Float::epsilon(wp), 10 * wp);
    result.re.setPrecision(precision);
    result.im.setPrecision(precision);
    return result;
}

[[nodiscard]] inline Complex<Float> jacobiTheta2(const Complex<Float>& z, const Float& q) {
    return jacobiTheta2(z, q, Float::defaultPrecision());
}

[[nodiscard]] inline Complex<Float> jacobiTheta3(const Complex<Float>& z, const Float& q, int precision) {
    using C = Complex<Float>;
    int wp = precision + 20;
    Float::PrecisionScope _ps(wp);   // compute exact÷exact inside the q-series at wp digits
    // PrecisionGuard removed: setResultPrecision(wp) on zw / qw propagates req=wp.
    C zw = z;
    zw.re.setResultPrecision(wp);
    zw.im.setResultPrecision(wp);
    Float qw = q;
    qw.setResultPrecision(wp);
    C result = detail::jacobiTheta3_complex<Float>(zw, qw,
        Float::epsilon(wp), 10 * wp);
    result.re.setPrecision(precision);
    result.im.setPrecision(precision);
    return result;
}

[[nodiscard]] inline Complex<Float> jacobiTheta3(const Complex<Float>& z, const Float& q) {
    return jacobiTheta3(z, q, Float::defaultPrecision());
}

[[nodiscard]] inline Complex<Float> jacobiTheta4(const Complex<Float>& z, const Float& q, int precision) {
    using C = Complex<Float>;
    int wp = precision + 20;
    Float::PrecisionScope _ps(wp);   // compute exact÷exact inside the q-series at wp digits
    // PrecisionGuard removed: setResultPrecision(wp) on zw / qw propagates req=wp.
    C zw = z;
    zw.re.setResultPrecision(wp);
    zw.im.setResultPrecision(wp);
    Float qw = q;
    qw.setResultPrecision(wp);
    C result = detail::jacobiTheta4_complex<Float>(zw, qw,
        Float::epsilon(wp), 10 * wp);
    result.re.setPrecision(precision);
    result.im.setPrecision(precision);
    return result;
}

[[nodiscard]] inline Complex<Float> jacobiTheta4(const Complex<Float>& z, const Float& q) {
    return jacobiTheta4(z, q, Float::defaultPrecision());
}

// ---- Complex<Float> z + Complex<Float> q (complex nome, arbitrary precision) ----
// The core of SPECIAL-X4 (TODO_BEYOND_MPFR). For any upper half-plane τ ∈ ℍ,
// q = exp(iπτ) (complex) can be accepted directly.

namespace detail {
inline Float complex_modSq_float(const Complex<Float>& q) {
    return q.re * q.re + q.im * q.im;
}
} // namespace detail

[[nodiscard]] inline Complex<Float> jacobiTheta1(const Complex<Float>& z, const Complex<Float>& q, int precision) {
    using C = Complex<Float>;
    int wp = precision + 20;
    Float::PrecisionScope _ps(wp);   // compute exact÷exact inside the q-series at wp digits
    C zw = z;
    zw.re.setResultPrecision(wp);
    zw.im.setResultPrecision(wp);
    C qw = q;
    qw.re.setResultPrecision(wp);
    qw.im.setResultPrecision(wp);
    C result = detail::jacobiTheta1_cqcz<Float>(zw, qw,
        Float::epsilon(wp), 10 * wp);
    result.re.setPrecision(precision);
    result.im.setPrecision(precision);
    return result;
}

[[nodiscard]] inline Complex<Float> jacobiTheta1(const Complex<Float>& z, const Complex<Float>& q) {
    return jacobiTheta1(z, q, Float::defaultPrecision());
}

[[nodiscard]] inline Complex<Float> jacobiTheta2(const Complex<Float>& z, const Complex<Float>& q, int precision) {
    using C = Complex<Float>;
    int wp = precision + 20;
    Float::PrecisionScope _ps(wp);   // compute exact÷exact inside the q-series at wp digits
    C zw = z;
    zw.re.setResultPrecision(wp);
    zw.im.setResultPrecision(wp);
    C qw = q;
    qw.re.setResultPrecision(wp);
    qw.im.setResultPrecision(wp);
    C result = detail::jacobiTheta2_cqcz<Float>(zw, qw,
        Float::epsilon(wp), 10 * wp);
    result.re.setPrecision(precision);
    result.im.setPrecision(precision);
    return result;
}

[[nodiscard]] inline Complex<Float> jacobiTheta2(const Complex<Float>& z, const Complex<Float>& q) {
    return jacobiTheta2(z, q, Float::defaultPrecision());
}

[[nodiscard]] inline Complex<Float> jacobiTheta3(const Complex<Float>& z, const Complex<Float>& q, int precision) {
    using C = Complex<Float>;
    int wp = precision + 20;
    Float::PrecisionScope _ps(wp);   // compute exact÷exact inside the q-series at wp digits
    C zw = z;
    zw.re.setResultPrecision(wp);
    zw.im.setResultPrecision(wp);
    C qw = q;
    qw.re.setResultPrecision(wp);
    qw.im.setResultPrecision(wp);
    C result = detail::jacobiTheta3_cqcz<Float>(zw, qw,
        Float::epsilon(wp), 10 * wp);
    result.re.setPrecision(precision);
    result.im.setPrecision(precision);
    return result;
}

[[nodiscard]] inline Complex<Float> jacobiTheta3(const Complex<Float>& z, const Complex<Float>& q) {
    return jacobiTheta3(z, q, Float::defaultPrecision());
}

[[nodiscard]] inline Complex<Float> jacobiTheta4(const Complex<Float>& z, const Complex<Float>& q, int precision) {
    using C = Complex<Float>;
    int wp = precision + 20;
    Float::PrecisionScope _ps(wp);   // compute exact÷exact inside the q-series at wp digits
    C zw = z;
    zw.re.setResultPrecision(wp);
    zw.im.setResultPrecision(wp);
    C qw = q;
    qw.re.setResultPrecision(wp);
    qw.im.setResultPrecision(wp);
    C result = detail::jacobiTheta4_cqcz<Float>(zw, qw,
        Float::epsilon(wp), 10 * wp);
    result.re.setPrecision(precision);
    result.im.setPrecision(precision);
    return result;
}

[[nodiscard]] inline Complex<Float> jacobiTheta4(const Complex<Float>& z, const Complex<Float>& q) {
    return jacobiTheta4(z, q, Float::defaultPrecision());
}

// ----------------------------------------------------------------
// AGM-accelerated version — high-precision computation of θ_2, θ_3, θ_4 at z=0
// ----------------------------------------------------------------
//
// Reference: J.M. Borwein & P.B. Borwein, "Pi and the AGM" (1987) Algorithm 2.1
// Idea: a_n = θ_3(0, q^{2^n})², b_n = θ_4(0, q^{2^n})² satisfy the forward AGM
//   a_{n+1} = (a_n + b_n)/2,  b_{n+1} = √(a_n b_n)
// Since q^{2^n} → 0 (n→∞), a_n, b_n → 1 with quadratic convergence.
//
// Computation strategy (inverse AGM):
//   1. Choose N and square until q^{2^N} is small enough (~2^{-(B/2)})
//   2. With q_N = q^{2^N}, compute a_N, b_N using the series (2-3 terms)
//   3. Apply the inverse AGM N times:
//      d = √(a_n² - b_n²),  a_{n-1} = a_n + d,  b_{n-1} = a_n - d
//   4. θ_3 = √a_0, θ_4 = √b_0, θ_2 = (θ_3⁴ - θ_4⁴)^{1/4} (Jacobi identity)
//
// Cost: at B digits, N ≈ log_2(B / |log_2 q|) iterations, each one sqrt.
//   For q=0.5, B=1000 digits, N≈10 (vs direct series ~32 terms + each q^{n²} computation).
// Applicability: only real nome 0 < q < 1. Complex q is future work (as is modular reduction).

namespace detail {

/// Compute θ_3(0, q_small)² and θ_4(0, q_small)² by direct series (only when q_small is small).
/// Since q_small = q^{2^N} is ≤ ~2^{-B/2}, 2-3 terms reach B-digit precision.
inline void thetaZeroSeriesSquares_(const Float& q_small, int wp,
                                    Float& a_out, Float& b_out) {
    // θ_3(0, q) = 1 + 2 (q + q^4 + q^9 + q^16 + ...)
    // θ_4(0, q) = 1 + 2 (-q + q^4 - q^9 + q^16 - ...)
    // Recurrence: q^{(n+1)²} = q^{n²} · q^{2n+1}
    //   ratio_n = q^{2n+1} is q³ at n=1, multiplied by q² each step
    Float t3(Float(1)); t3.setResultPrecision(wp);
    Float t4(Float(1)); t4.setResultPrecision(wp);
    Float qn = q_small; qn.setResultPrecision(wp);  // q^1 (= q^{1²})
    Float q_squared = q_small * q_small; q_squared.setResultPrecision(wp);  // q²
    Float ratio = q_small * q_squared; ratio.setResultPrecision(wp);  // q³ = ratio_1

    Float eps = Float::epsilon(wp);
    Float two(Float(2)); two.setResultPrecision(wp);
    int sign = -1;  // sign at n=1 (θ_4)

    for (int n = 1; n < 100; n++) {
        Float two_qn = two * qn;
        t3 = t3 + two_qn;
        if (sign > 0) t4 = t4 + two_qn; else t4 = t4 - two_qn;

        if (sangi::abs(two_qn) < eps) break;

        // qn → q^{(n+1)²}, ratio → q^{2(n+1)+1} = ratio · q²
        qn = qn * ratio;
        ratio = ratio * q_squared;
        sign = -sign;
    }

    a_out = t3 * t3;  // θ_3²
    b_out = t4 * t4;  // θ_4²
}

/// Compute θ_2, θ_3, θ_4 via AGM for a given q (0 < q < 1) and precision.
inline void jacobiThetaZeroAGM_impl_(const Float& q, int precision,
                                     Float& theta2_out, Float& theta3_out, Float& theta4_out) {
    // The old q≤0.2 dispatch (BUG_AGM_DYADIC_Q_PRECISION_20260509 workaround) became
    // unnecessary after the BUG_FLOAT_SPARSE_MANTISSA_PRECISION_20260509 H4 fix (commit 9e79cd5)
    // and has been removed (2026-05-09).

    // Choice of internal precision (important):
    //   At the deepest level a_N - b_N ≈ 4 q_N ≈ 2^{2 - B/2}. Since a_N, b_N are kept in wp bits,
    //   the relative error of the difference is ~ 2^{B/2 - wp}. Back-iteration scales the difference up
    //   to O(1) while the relative error propagates unchanged. A B-bit-precision final result needs wp ≥ 3B/2 + safety.
    int wp = precision + precision / 2 + 64;
    Float qw = q;
    qw.setPrecision(wp);
    qw.setResultPrecision(wp);

    // The old q-perturbation workaround (BUG_AGM_DYADIC_Q_PRECISION_20260509) became
    // unnecessary after the BUG_FLOAT_SPARSE_MANTISSA_PRECISION_20260509 H4 fix (commit 9e79cd5)
    // and has been removed (2026-05-09).

    // Compute N corresponding to q^{2^N} ≤ 2^{-B/2}
    double q_d = qw.toDouble();
    if (!(q_d > 0.0 && q_d < 1.0)) {
        Float nan_v = Float::nan();
        theta2_out = nan_v; theta3_out = nan_v; theta4_out = nan_v;
        return;
    }
    double log2_q_inv = -std::log2(q_d);  // > 0
    int N = 0;
    {
        double target = static_cast<double>(precision) / (2.0 * log2_q_inv);
        if (target > 1.0) N = static_cast<int>(std::ceil(std::log2(target)));
        if (N < 0) N = 0;
    }

    // q_N = q^{2^N}
    Float q_N = qw;
    for (int i = 0; i < N; i++) {
        q_N = q_N * q_N;
    }

    // The old q_N-perturbation workaround became unnecessary with the H4 fix and has been removed (2026-05-09).

    // a_N = θ_3(0, q_N)², b_N = θ_4(0, q_N)²
    Float a, b;
    a.setResultPrecision(wp); b.setResultPrecision(wp);
    thetaZeroSeriesSquares_(q_N, wp, a, b);

    // Inverse AGM: in each iteration
    //   d² = a² - b² = (a-b)(a+b)
    //   a_new = a + d, b_new = a - d
    // Note: since a > b > 0, d² > 0
    for (int n = 0; n < N; n++) {
        Float diff = a - b;
        Float sum = a + b;
        Float d_sq = diff * sum;
        Float d = sangi::sqrt(d_sq);
        Float a_prev = a + d;
        Float b_prev = a - d;
        a = a_prev;
        b = b_prev;
    }

    // a = θ_3(0,q)², b = θ_4(0,q)²
    Float t3 = sangi::sqrt(a);
    Float t4 = sangi::sqrt(b);

    // θ_2 = (θ_3⁴ - θ_4⁴)^{1/4} = (a² - b²)^{1/4}
    Float a_sq = a * a;
    Float b_sq = b * b;
    Float t2_4 = a_sq - b_sq;  // = θ_2⁴
    Float t2 = sangi::sqrt(sangi::sqrt(t2_4));

    t2.setPrecision(precision);
    t3.setPrecision(precision);
    t4.setPrecision(precision);
    theta2_out = t2;
    theta3_out = t3;
    theta4_out = t4;
}

} // namespace detail

/// Compute θ_2(0,q), θ_3(0,q), θ_4(0,q) simultaneously with AGM acceleration (real q ∈ (0,1)).
/// Faster than the series at high precision (above ~200 digits).
inline void jacobiThetaZeroAGM(const Float& q, int precision,
                               Float& theta2, Float& theta3, Float& theta4) {
    detail::jacobiThetaZeroAGM_impl_(q, precision, theta2, theta3, theta4);
}

inline void jacobiThetaZeroAGM(const Float& q,
                               Float& theta2, Float& theta3, Float& theta4) {
    detail::jacobiThetaZeroAGM_impl_(q, Float::defaultPrecision(), theta2, theta3, theta4);
}

/// θ_3(0, q) via AGM, standalone. Internally also computes θ_2, θ_4 but discards them (convenience function when only θ_3 is needed).
[[nodiscard]] inline Float jacobiTheta3ZeroAGM(const Float& q, int precision) {
    Float t2, t3, t4;
    detail::jacobiThetaZeroAGM_impl_(q, precision, t2, t3, t4);
    return t3;
}

[[nodiscard]] inline Float jacobiTheta4ZeroAGM(const Float& q, int precision) {
    Float t2, t3, t4;
    detail::jacobiThetaZeroAGM_impl_(q, precision, t2, t3, t4);
    return t4;
}

[[nodiscard]] inline Float jacobiTheta2ZeroAGM(const Float& q, int precision) {
    Float t2, t3, t4;
    detail::jacobiThetaZeroAGM_impl_(q, precision, t2, t3, t4);
    return t2;
}

// ---- qFromTau: τ → q = exp(iπτ) (arbitrary precision) ----
// Convert τ ∈ ℍ (upper half-plane, Im(τ) > 0) to the nome q. If Im(τ) > 0 then |q| < 1.
// The double version already exists in ThetaFunctions.hpp (sangi namespace). Here we provide
// the Float version in the sangi::special namespace.

[[nodiscard]] inline Complex<Float> qFromTau(const Complex<Float>& tau, int precision) {
    using C = Complex<Float>;
    int wp = precision + 20;
    Float::PrecisionScope _ps(wp);   // compute exact÷exact inside the function at wp digits
    C tw = tau;
    tw.re.setResultPrecision(wp);
    tw.im.setResultPrecision(wp);
    Float pi_val = Float::pi(wp);
    // q = exp(iπτ) = exp(C(0, π) * τ)
    C i_pi(Float(0), pi_val);
    C result = sangi::exp(i_pi * tw);
    result.re.setPrecision(precision);
    result.im.setPrecision(precision);
    return result;
}

[[nodiscard]] inline Complex<Float> qFromTau(const Complex<Float>& tau) {
    return qFromTau(tau, Float::defaultPrecision());
}

// ----------------------------------------------------------------
// τ-form API + Modular Reduction (DLMF 20.7)
// ----------------------------------------------------------------
//
// Compute θ_n(z | τ) for any τ in the upper half-plane (Im(τ) > 0).
// When |q| = exp(-π Im(τ)) is close to 1 (= small Im(τ)), the series converges slowly.
// Using T (τ → τ+1) and S (τ → -1/τ) of SL(2,Z), reduce τ into the fundamental domain
// (|τ| ≥ 1 and |Re(τ)| ≤ 1/2, i.e. Im(τ) ≥ √3/2 ≈ 0.866),
// then evaluate the q-series.
//
// Transformation rules (DLMF 20.7.5/20.7.30, etc.):
//   T: θ_1, θ_2 acquire an e^{iπ/4} phase, θ_3 ↔ θ_4 swap
//   S: θ_n(z|τ) = ζ_n · (-iτ)^{-1/2} · exp(iπ z²/τ) · θ_{S(n)}(z/τ | -1/τ)
//        ζ_1 = -i, S(1)=1; ζ_2=1, S(2)=4; ζ_3=1, S(3)=3; ζ_4=1, S(4)=2

namespace detail {

/// Round a Float to the nearest int64_t (for k = round(Re(τ)) in T-reduction).
inline int64_t roundFloatToInt_(const Float& x) {
    double d = x.toDouble();
    if (!std::isfinite(d)) return 0;
    if (d > 1e18) return static_cast<int64_t>(1e18);
    if (d < -1e18) return -static_cast<int64_t>(1e18);
    return static_cast<int64_t>(std::llround(d));
}

/// Reduce τ into the fundamental domain via SL(2,Z), returning the corresponding theta index
/// and the accumulated phase. Up to max_iter T+S iterations.
struct ThetaReduction_ {
    Complex<Float> z;
    Complex<Float> tau;
    int n;                  // index of the θ_n to evaluate finally (1..4)
    Complex<Float> phase;   // phase factor to multiply into the result
};

inline ThetaReduction_ reduceTau_(int n_in, const Complex<Float>& z_in,
                                  const Complex<Float>& tau_in, int wp,
                                  int max_iter = 20) {
    Complex<Float> z = z_in;
    Complex<Float> tau = tau_in;
    z.re.setResultPrecision(wp); z.im.setResultPrecision(wp);
    tau.re.setResultPrecision(wp); tau.im.setResultPrecision(wp);

    Complex<Float> phase(Float(1), Float(0));
    phase.re.setResultPrecision(wp); phase.im.setResultPrecision(wp);
    int tracked_n = n_in;

    Float pi_val = Float::pi(wp);
    Float one(Float(1)); one.setResultPrecision(wp);

    for (int iter = 0; iter < max_iter; iter++) {
        // T-reduction: k = round(Re(τ)), τ ← τ - k
        int64_t k = roundFloatToInt_(tau.re);
        if (k != 0) {
            Float kf(static_cast<double>(k)); kf.setResultPrecision(wp);
            tau.re = tau.re - kf;

            if (tracked_n == 1 || tracked_n == 2) {
                // phase *= exp(iπk/4)
                Float angle = pi_val * kf / Float(4);
                Complex<Float> rot(sangi::cos(angle), sangi::sin(angle));
                phase = phase * rot;
            } else {
                // 3 or 4: swap if |k| odd
                if (k % 2 != 0) tracked_n = 7 - tracked_n;  // 3↔4
            }
        }

        // fundamental domain if |τ|² ≥ 1
        Float tau_abs_sq = tau.re * tau.re + tau.im * tau.im;
        if (tau_abs_sq >= one) break;

        // S-transform: τ → -1/τ, z → z/τ
        // prefactor = (-iτ)^{-1/2} · exp(iπz²/τ)
        // -iτ = (Im(τ), -Re(τ))
        Complex<Float> minus_i_tau(tau.im, -tau.re);
        Complex<Float> sqrt_minus_i_tau = sangi::sqrt(minus_i_tau);
        Complex<Float> inv_sqrt = Complex<Float>(one, Float(0)) / sqrt_minus_i_tau;

        Complex<Float> i_pi(Float(0), pi_val);
        Complex<Float> z_sq = z * z;
        Complex<Float> exp_arg = (i_pi * z_sq) / tau;
        Complex<Float> exp_part = sangi::exp(exp_arg);

        Complex<Float> prefactor = inv_sqrt * exp_part;

        if (tracked_n == 1) {
            // ζ_1 = -i
            Complex<Float> minus_i(Float(0), Float(-1));
            phase = phase * minus_i * prefactor;
        } else if (tracked_n == 2) {
            phase = phase * prefactor;
            tracked_n = 4;
        } else if (tracked_n == 3) {
            phase = phase * prefactor;
        } else {  // 4
            phase = phase * prefactor;
            tracked_n = 2;
        }

        Complex<Float> new_z = z / tau;
        Complex<Float> new_tau = Complex<Float>(Float(-1), Float(0)) / tau;
        z = new_z;
        tau = new_tau;
    }

    return {z, tau, tracked_n, phase};
}

/// Based on the reduction result, evaluate the q-form at n_reduced and multiply by the phase.
inline Complex<Float> evaluateReduced_(const ThetaReduction_& red, int prec, int wp) {
    Complex<Float> q = sangi::exp(Complex<Float>(Float(0), Float::pi(wp)) * red.tau);
    Complex<Float> theta_val;
    switch (red.n) {
        case 1: theta_val = jacobiTheta1(red.z, q, prec); break;
        case 2: theta_val = jacobiTheta2(red.z, q, prec); break;
        case 3: theta_val = jacobiTheta3(red.z, q, prec); break;
        default: theta_val = jacobiTheta4(red.z, q, prec); break;  // 4
    }
    Complex<Float> result = red.phase * theta_val;
    result.re.setPrecision(prec);
    result.im.setPrecision(prec);
    return result;
}

} // namespace detail

/// τ-form API: compute θ_n(z | τ) including modular reduction (T+S).
/// τ may be any point in the upper half-plane (Im(τ) > 0). Faster when |q| is close to 1.

[[nodiscard]] inline Complex<Float> jacobiTheta1Tau(const Complex<Float>& z, const Complex<Float>& tau, int precision) {
    int wp = precision + 64;
    Float::PrecisionScope _ps(wp);   // compute exact÷exact inside reduceTau_/evaluateReduced_ at wp digits
    auto red = detail::reduceTau_(1, z, tau, wp);
    return detail::evaluateReduced_(red, precision, wp);
}
[[nodiscard]] inline Complex<Float> jacobiTheta1Tau(const Complex<Float>& z, const Complex<Float>& tau) {
    return jacobiTheta1Tau(z, tau, Float::defaultPrecision());
}

[[nodiscard]] inline Complex<Float> jacobiTheta2Tau(const Complex<Float>& z, const Complex<Float>& tau, int precision) {
    int wp = precision + 64;
    Float::PrecisionScope _ps(wp);   // compute exact÷exact inside reduceTau_/evaluateReduced_ at wp digits
    auto red = detail::reduceTau_(2, z, tau, wp);
    return detail::evaluateReduced_(red, precision, wp);
}
[[nodiscard]] inline Complex<Float> jacobiTheta2Tau(const Complex<Float>& z, const Complex<Float>& tau) {
    return jacobiTheta2Tau(z, tau, Float::defaultPrecision());
}

[[nodiscard]] inline Complex<Float> jacobiTheta3Tau(const Complex<Float>& z, const Complex<Float>& tau, int precision) {
    int wp = precision + 64;
    Float::PrecisionScope _ps(wp);   // compute exact÷exact inside reduceTau_/evaluateReduced_ at wp digits
    auto red = detail::reduceTau_(3, z, tau, wp);
    return detail::evaluateReduced_(red, precision, wp);
}
[[nodiscard]] inline Complex<Float> jacobiTheta3Tau(const Complex<Float>& z, const Complex<Float>& tau) {
    return jacobiTheta3Tau(z, tau, Float::defaultPrecision());
}

[[nodiscard]] inline Complex<Float> jacobiTheta4Tau(const Complex<Float>& z, const Complex<Float>& tau, int precision) {
    int wp = precision + 64;
    Float::PrecisionScope _ps(wp);   // compute exact÷exact inside reduceTau_/evaluateReduced_ at wp digits
    auto red = detail::reduceTau_(4, z, tau, wp);
    return detail::evaluateReduced_(red, precision, wp);
}
[[nodiscard]] inline Complex<Float> jacobiTheta4Tau(const Complex<Float>& z, const Complex<Float>& tau) {
    return jacobiTheta4Tau(z, tau, Float::defaultPrecision());
}

} // namespace special
} // namespace sangi

#endif // SANGI_SPECIAL_THETA_HPP