// Copyright (C) 2026 Kiyotsugu Arai // SPDX-License-Identifier: LGPL-3.0-or-later // weierstrass.hpp // Weierstrass elliptic functions ℘ / ζ_W / σ_W // // Functions provided (take the half-periods ω₁, ω₂ as input): // weierstrassP(z, omega1, omega2) — Weierstrass ℘ function // weierstrassZeta(z, omega1, omega2) — Weierstrass ζ_W function (≠ Riemann ζ) // weierstrassSigma(z, omega1, omega2) — Weierstrass σ_W function // // Properties: // ℘ is an even function with double periods (2ω₁, 2ω₂), with a pole of order 2 at lattice points // (℘')² = 4℘³ - g₂℘ - g₃ (Weierstrass differential equation) // ζ_W'(z) = -℘(z), σ_W'(z)/σ_W(z) = ζ_W(z) // // Implementation: // Uses the theta-function representation from DLMF 23.6 // ℘(z) = e₃ + (π/(2ω₁) · θ₂(0,q) θ₃(0,q) · θ₄(v,q)/θ₁(v,q))² // q is generally complex (the real part of τ = ω₂/ω₁ is nonzero); internally the theta series for complex q is expanded directly. // // Supported types (v0.1): // Complex — fully supported (any half-periods / any z) // double — only for real z, ω₁ + purely imaginary ω₂ (rectangular lattice). Real-valued if poles are avoided // // TODO: // - Float / Complex arbitrary-precision version (in a future release) // - Initialization from the invariants (g₂, g₃) (solve the cubic 4w³ - g₂ w - g₃ = 0) #ifndef SANGI_SPECIAL_WEIERSTRASS_HPP #define SANGI_SPECIAL_WEIERSTRASS_HPP #include

#include #include #include #include #include namespace sangi {
namespace special {

namespace detail_weier {

// ------------------------------------------------------------------
// Internal theta functions (arbitrary complex v, complex q)
// theta.hpp handles only real q (or the Float version of Complex),
// so this is a local implementation for the "general complex q" needed by the Weierstrass functions.
// ------------------------------------------------------------------

/// θ₁(v, q) = 2 Σ_{n=0}^∞ (-1)^n q^{(n+1/2)²} sin((2n+1)v)
template [[nodiscard]] Complex theta1(const Complex & v, const Complex & q,
                                int max_iter = 200) {
    using C = Complex;
    C sum(R(0));
    R sign = R(1);
    R eps = std::numeric_limits ::epsilon();

    for (int n = 0; n < max_iter; n++) {
        R k = R(n) + R(0.5);
        C qpow = pow(q, C(k * k));
        C term = qpow * sin(C(R(2 * n + 1)) * v);
        if (sign < R(0)) term = -term;
        sum = sum + term;
        if (n >= 4 && abs(term) < eps * (R(1) + abs(sum))) break;
        sign = -sign;
    }
    return C(R(2)) * sum;
}

/// θ₂(v, q) = 2 Σ_{n=0}^∞ q^{(n+1/2)²} cos((2n+1)v)
template [[nodiscard]] Complex theta2(const Complex & v, const Complex & q,
                                int max_iter = 200) {
    using C = Complex;
    C sum(R(0));
    R eps = std::numeric_limits ::epsilon();

    for (int n = 0; n < max_iter; n++) {
        R k = R(n) + R(0.5);
        C qpow = pow(q, C(k * k));
        C term = qpow * cos(C(R(2 * n + 1)) * v);
        sum = sum + term;
        if (n >= 4 && abs(term) < eps * (R(1) + abs(sum))) break;
    }
    return C(R(2)) * sum;
}

/// θ₃(v, q) = 1 + 2 Σ_{n=1}^∞ q^{n²} cos(2nv)
template [[nodiscard]] Complex theta3(const Complex & v, const Complex & q,
                                int max_iter = 200) {
    using C = Complex;
    C sum(R(0));
    R eps = std::numeric_limits ::epsilon();

    for (int n = 1; n < max_iter; n++) {
        C qpow = pow(q, C(R(n) * R(n)));
        C term = qpow * cos(C(R(2 * n)) * v);
        sum = sum + term;
        if (n >= 4 && abs(term) < eps * (R(1) + abs(sum))) break;
    }
    return C(R(1)) + C(R(2)) * sum;
}

/// θ₄(v, q) = 1 + 2 Σ_{n=1}^∞ (-1)^n q^{n²} cos(2nv)
template [[nodiscard]] Complex theta4(const Complex & v, const Complex & q,
                                int max_iter = 200) {
    using C = Complex;
    C sum(R(0));
    R sign = R(-1);
    R eps = std::numeric_limits ::epsilon();

    for (int n = 1; n < max_iter; n++) {
        C qpow = pow(q, C(R(n) * R(n)));
        C term = qpow * cos(C(R(2 * n)) * v);
        if (sign < R(0)) term = -term;
        sum = sum + term;
        if (n >= 4 && abs(term) < eps * (R(1) + abs(sum))) break;
        sign = -sign;
    }
    return C(R(1)) + C(R(2)) * sum;
}

// ∂θ₁/∂v at v=0  =  2 Σ (-1)^n q^{(n+1/2)²} (2n+1)
template [[nodiscard]] Complex theta1_prime_at_zero(const Complex & q,
                                              int max_iter = 200) {
    using C = Complex;
    C sum(R(0));
    R sign = R(1);
    R eps = std::numeric_limits ::epsilon();

    for (int n = 0; n < max_iter; n++) {
        R k = R(n) + R(0.5);
        C qpow = pow(q, C(k * k));
        C term = qpow * C(R(2 * n + 1));
        if (sign < R(0)) term = -term;
        sum = sum + term;
        if (n >= 4 && abs(term) < eps * (R(1) + abs(sum))) break;
        sign = -sign;
    }
    return C(R(2)) * sum;
}

// ∂³θ₁/∂v³ at v=0  =  -2 Σ (-1)^n q^{(n+1/2)²} (2n+1)³
template [[nodiscard]] Complex theta1_triple_prime_at_zero(const Complex & q,
                                                    int max_iter = 200) {
    using C = Complex;
    C sum(R(0));
    R sign = R(1);
    R eps = std::numeric_limits ::epsilon();

    for (int n = 0; n < max_iter; n++) {
        R k = R(n) + R(0.5);
        C qpow = pow(q, C(k * k));
        R cube = R(2 * n + 1) * R(2 * n + 1) * R(2 * n + 1);
        C term = qpow * C(cube);
        if (sign < R(0)) term = -term;
        sum = sum + term;
        if (n >= 4 && abs(term) < eps * (R(1) + abs(sum))) break;
        sign = -sign;
    }
    return C(R(-2)) * sum;
}

// ∂θ₁/∂v at general v  =  2 Σ (-1)^n q^{(n+1/2)²} (2n+1) cos((2n+1)v)
template [[nodiscard]] Complex theta1_prime(const Complex & v, const Complex & q,
                                     int max_iter = 200) {
    using C = Complex;
    C sum(R(0));
    R sign = R(1);
    R eps = std::numeric_limits ::epsilon();

    for (int n = 0; n < max_iter; n++) {
        R k = R(n) + R(0.5);
        C qpow = pow(q, C(k * k));
        C term = qpow * C(R(2 * n + 1)) * cos(C(R(2 * n + 1)) * v);
        if (sign < R(0)) term = -term;
        sum = sum + term;
        if (n >= 4 && abs(term) < eps * (R(1) + abs(sum))) break;
        sign = -sign;
    }
    return C(R(2)) * sum;
}

// Validity check and normalization of the half-periods (ω₁, ω₂):
// Invalid unless τ = ω₂/ω₁ is in the upper half-plane (Im(τ) > 0).
// When valid, returns q = exp(iπτ) (|q| < 1 is guaranteed).
struct LatticeData {
    Complex alpha;  // π/(2ω₁)
    Complex q;      // nome
    bool valid;
};

inline LatticeData prepare_lattice(const Complex & omega1,
                                   const Complex & omega2) {
    LatticeData out{};
    out.valid = false;
    if (omega1.re == 0.0 && omega1.im == 0.0) return out;
    Complex