scirs2-special 0.4.2

//! Elliptic Modular Functions
//!
//! This module implements classical elliptic modular functions that arise in number theory,
//! complex analysis, and the theory of modular forms. These functions are defined on the
//! upper half-plane (Im(tau) > 0) and possess remarkable transformation properties under
//! the modular group SL(2, Z).
//!
//! ## Connection to Theta Functions and Elliptic Integrals
//!
//! Elliptic modular functions are intimately connected to Jacobi theta functions and
//! elliptic integrals. The **nome** `q = exp(i*pi*tau)` serves as the bridge:
//!
//! - For `tau = i * tau_imag` (purely imaginary), we have `q = exp(-pi * tau_imag)`,
//!   which is real and lies in `(0, 1)` when `tau_imag > 0`.
//! - The Eisenstein series `E_4`, `E_6` can be expressed through theta functions.
//! - The modular discriminant `Delta(tau) = eta(tau)^24` connects the Dedekind eta
//!   function to the discriminant of the corresponding elliptic curve.
//! - The Klein j-invariant `j(tau) = 1728 * E_4^3 / (E_4^3 - E_6^2)` classifies
//!   elliptic curves up to isomorphism over C.
//!
//! ## Implemented Functions
//!
//! | Function | Description |
//! |----------|-------------|
//! | [`dedekind_eta`] | Dedekind eta function `eta(tau)` |
//! | [`klein_j_invariant`] | Klein j-invariant `j(tau) = 1728 E_4^3 / Delta` |
//! | [`modular_lambda`] | Modular lambda function `lambda(tau) = theta_2^4 / theta_3^4` |
//! | [`eisenstein_e2`] | Eisenstein series `E_2(tau)` (quasi-modular, weight 2) |
//! | [`eisenstein_e4`] | Eisenstein series `E_4(tau)` (modular form, weight 4) |
//! | [`eisenstein_e6`] | Eisenstein series `E_6(tau)` (modular form, weight 6) |
//! | [`modular_discriminant`] | Modular discriminant `Delta(tau) = eta(tau)^24` |
//! | [`lattice_invariants_from_tau`] | Weierstrass invariants `(g_2, g_3)` from `tau` |
//!
//! ## Restriction to Purely Imaginary tau
//!
//! This implementation restricts to `tau = i * tau_imag` with `tau_imag > 0`, where all
//! functions take real values. This covers many practical applications and avoids the
//! complexity of full complex arithmetic. The input parameter is `tau_imag` (a positive
//! real number), and `q = exp(-2*pi*tau_imag)` is used internally.
//!
//! ## References
//!
//! - Apostol, *Modular Functions and Dirichlet Series in Number Theory*
//! - DLMF §23, §27: Weierstrass Elliptic Functions, Functions of Number Theory
//! - Zagier, "Elliptic Modular Forms and Their Applications" (in *The 1-2-3 of Modular Forms*)

use crate::error::{SpecialError, SpecialResult};
use crate::theta_functions::{theta2, theta3};

// ---------------------------------------------------------------------------
// Constants
// ---------------------------------------------------------------------------

/// Convergence tolerance for q-series
const Q_SERIES_TOL: f64 = 1e-16;

/// Maximum terms in q-series expansions
const MAX_Q_TERMS: usize = 1000;

/// Maximum terms for the Euler product in Dedekind eta
const MAX_ETA_TERMS: usize = 1500;

// ---------------------------------------------------------------------------
// Input validation
// ---------------------------------------------------------------------------

/// Validate that tau_imag > 0
#[inline]
fn validate_tau_imag(tau_imag: f64) -> SpecialResult<()> {
    if !tau_imag.is_finite() || tau_imag <= 0.0 {
        return Err(SpecialError::DomainError(format!(
            "tau_imag must be strictly positive and finite, got {tau_imag}"
        )));
    }
    Ok(())
}

/// Compute q = exp(-2*pi*tau_imag) from tau_imag.
/// For purely imaginary tau = i*tau_imag, the nome-squared is q = exp(2*pi*i*tau) = exp(-2*pi*tau_imag).
#[inline]
fn compute_q(tau_imag: f64) -> f64 {
    (-2.0 * std::f64::consts::PI * tau_imag).exp()
}

// ---------------------------------------------------------------------------
// Private helpers: divisor sums and q-series
// ---------------------------------------------------------------------------

/// Compute the divisor sum sigma_k(n) = sum_{d | n} d^k.
///
/// This iterates over all divisors of n up to sqrt(n) for efficiency.
fn divisor_sum(n: u64, k: u32) -> f64 {
    if n == 0 {
        return 0.0;
    }
    if n == 1 {
        return 1.0;
    }

    let mut sum = 0.0_f64;
    let sqrt_n = (n as f64).sqrt() as u64;

    for d in 1..=sqrt_n {
        if n.is_multiple_of(d) {
            sum += (d as f64).powi(k as i32);
            let complement = n / d;
            if complement != d {
                sum += (complement as f64).powi(k as i32);
            }
        }
    }

    sum
}

/// Compute the Eisenstein q-series: sum_{n=1}^{N} sigma_k(n) * q^n
/// with early termination when terms become negligible.
fn eisenstein_q_series(q: f64, k: u32, max_terms: usize) -> f64 {
    if q < 1e-300 {
        return 0.0;
    }

    let mut sum = 0.0_f64;
    let mut q_power = q; // q^1

    for n in 1..=max_terms {
        let sigma = divisor_sum(n as u64, k);
        let term = sigma * q_power;

        sum += term;

        // Early termination: if the contribution is negligible
        if term.abs() < Q_SERIES_TOL * sum.abs().max(1.0) {
            break;
        }

        q_power *= q;
        if q_power < 1e-300 {
            break;
        }
    }

    sum
}

// ---------------------------------------------------------------------------
// Eisenstein series
// ---------------------------------------------------------------------------

/// Eisenstein series E_2(tau) for tau = i * tau_imag.
///
/// # Definition
///
/// ```text
/// E_2(tau) = 1 - 24 * sum_{n=1}^{infinity} sigma_1(n) * q^n
/// ```
///
/// where `sigma_1(n) = sum_{d | n} d` is the sum-of-divisors function
/// and `q = exp(2*pi*i*tau) = exp(-2*pi*tau_imag)` for purely imaginary tau.
///
/// **Note**: E_2 is a *quasi-modular* form of weight 2. Unlike E_4 and E_6,
/// it does not transform as a true modular form under SL(2, Z).
///
/// # Arguments
///
/// * `tau_imag` - The imaginary part of tau (must be > 0)
///
/// # Returns
///
/// The value of E_2(i * tau_imag).
///
/// # Errors
///
/// Returns [`SpecialError::DomainError`] if `tau_imag <= 0` or is not finite.
///
/// # Examples
///
/// ```
/// use scirs2_special::elliptic_modular::eisenstein_e2;
/// // For large tau_imag, q -> 0 and E_2 -> 1
/// let e2 = eisenstein_e2(10.0).expect("valid input");
/// assert!((e2 - 1.0).abs() < 1e-10);
/// ```
pub fn eisenstein_e2(tau_imag: f64) -> SpecialResult<f64> {
    validate_tau_imag(tau_imag)?;
    let q = compute_q(tau_imag);
    let series_sum = eisenstein_q_series(q, 1, MAX_Q_TERMS);
    Ok(1.0 - 24.0 * series_sum)
}

/// Eisenstein series E_4(tau) for tau = i * tau_imag.
///
/// # Definition
///
/// ```text
/// E_4(tau) = 1 + 240 * sum_{n=1}^{infinity} sigma_3(n) * q^n
/// ```
///
/// where `sigma_3(n) = sum_{d | n} d^3` and `q = exp(-2*pi*tau_imag)`.
///
/// E_4 is a modular form of weight 4 for the full modular group SL(2, Z).
/// It is proportional to the Weierstrass invariant g_2:
/// `g_2 = (2*pi)^4/12 * E_4(tau)`.
///
/// # Arguments
///
/// * `tau_imag` - The imaginary part of tau (must be > 0)
///
/// # Returns
///
/// The value of E_4(i * tau_imag).
///
/// # Known Values
///
/// - E_4(i) is related to the Gamma function: `E_4(i) = 3 * Gamma(1/4)^8 / (4 * pi^6)`
///
/// # Errors
///
/// Returns [`SpecialError::DomainError`] if `tau_imag <= 0` or is not finite.
///
/// # Examples
///
/// ```
/// use scirs2_special::elliptic_modular::eisenstein_e4;
/// // For large tau_imag, E_4 -> 1
/// let e4 = eisenstein_e4(10.0).expect("valid input");
/// assert!((e4 - 1.0).abs() < 1e-10);
/// ```
pub fn eisenstein_e4(tau_imag: f64) -> SpecialResult<f64> {
    validate_tau_imag(tau_imag)?;
    let q = compute_q(tau_imag);
    let series_sum = eisenstein_q_series(q, 3, MAX_Q_TERMS);
    Ok(1.0 + 240.0 * series_sum)
}

/// Eisenstein series E_6(tau) for tau = i * tau_imag.
///
/// # Definition
///
/// ```text
/// E_6(tau) = 1 - 504 * sum_{n=1}^{infinity} sigma_5(n) * q^n
/// ```
///
/// where `sigma_5(n) = sum_{d | n} d^5` and `q = exp(-2*pi*tau_imag)`.
///
/// E_6 is a modular form of weight 6 for the full modular group SL(2, Z).
/// It is proportional to the Weierstrass invariant g_3:
/// `g_3 = (2*pi)^6/216 * E_6(tau)`.
///
/// # Arguments
///
/// * `tau_imag` - The imaginary part of tau (must be > 0)
///
/// # Returns
///
/// The value of E_6(i * tau_imag).
///
/// # Errors
///
/// Returns [`SpecialError::DomainError`] if `tau_imag <= 0` or is not finite.
///
/// # Examples
///
/// ```
/// use scirs2_special::elliptic_modular::eisenstein_e6;
/// // For large tau_imag, E_6 -> 1
/// let e6 = eisenstein_e6(10.0).expect("valid input");
/// assert!((e6 - 1.0).abs() < 1e-10);
/// ```
pub fn eisenstein_e6(tau_imag: f64) -> SpecialResult<f64> {
    validate_tau_imag(tau_imag)?;
    let q = compute_q(tau_imag);
    let series_sum = eisenstein_q_series(q, 5, MAX_Q_TERMS);
    Ok(1.0 - 504.0 * series_sum)
}

// ---------------------------------------------------------------------------
// Dedekind eta function
// ---------------------------------------------------------------------------

/// Dedekind eta function for tau = i * tau_imag.
///
/// # Definition
///
/// ```text
/// eta(tau) = q^{1/24} * prod_{n=1}^{infinity} (1 - q^n)
/// ```
///
/// where `q = exp(2*pi*i*tau)`. For purely imaginary `tau = i*tau_imag`:
///
/// ```text
/// q = exp(-2*pi*tau_imag)   (real, in (0,1))
/// q^{1/24} = exp(-pi*tau_imag/12)
/// ```
///
/// The Dedekind eta function is a modular form of weight 1/2 (with a multiplier system).
/// Its 24th power gives the modular discriminant: `Delta(tau) = eta(tau)^24`.
///
/// # Arguments
///
/// * `tau_imag` - The imaginary part of tau (must be > 0)
///
/// # Returns
///
/// The value of eta(i * tau_imag), which is always positive.
///
/// # Errors
///
/// Returns [`SpecialError::DomainError`] if `tau_imag <= 0` or is not finite.
///
/// # Known Values
///
/// - `eta(i) = Gamma(1/4) / (2 * pi^{3/4})` approximately 0.76823...
///
/// # Examples
///
/// ```
/// use scirs2_special::elliptic_modular::dedekind_eta;
/// let eta = dedekind_eta(1.0).expect("valid input");
/// // eta(i) ≈ 0.76823
/// assert!((eta - 0.76823).abs() < 0.001);
/// ```
pub fn dedekind_eta(tau_imag: f64) -> SpecialResult<f64> {
    validate_tau_imag(tau_imag)?;

    let q = compute_q(tau_imag);

    // q^{1/24} = exp(-2*pi*tau_imag / 24) = exp(-pi*tau_imag / 12)
    let q_24th = (-std::f64::consts::PI * tau_imag / 12.0).exp();

    if q_24th < 1e-300 {
        // For very large tau_imag, eta -> 0
        return Ok(0.0);
    }

    // Euler product: prod_{n=1}^{infinity} (1 - q^n)
    let mut product = 1.0_f64;
    let mut q_power = q; // q^1

    for _n in 1..=MAX_ETA_TERMS {
        if q_power < 1e-300 {
            break;
        }

        let factor = 1.0 - q_power;
        product *= factor;

        // If q^n is negligible, all subsequent factors are essentially 1
        if (1.0 - factor).abs() < Q_SERIES_TOL {
            break;
        }

        q_power *= q;
    }

    Ok(q_24th * product)
}

// ---------------------------------------------------------------------------
// Klein j-invariant
// ---------------------------------------------------------------------------

/// Klein j-invariant for tau = i * tau_imag.
///
/// # Definition
///
/// ```text
/// j(tau) = 1728 * E_4(tau)^3 / (E_4(tau)^3 - E_6(tau)^2)
/// ```
///
/// equivalently, `j(tau) = 1728 * E_4^3 / Delta` where `Delta = (E_4^3 - E_6^2) / 1728`.
///
/// The j-invariant is a modular function (weight 0) that is invariant under the full
/// modular group. It provides a bijection from isomorphism classes of elliptic curves
/// over C to C.
///
/// # Arguments
///
/// * `tau_imag` - The imaginary part of tau (must be > 0)
///
/// # Returns
///
/// The value of j(i * tau_imag).
///
/// # Known Values
///
/// - `j(i) = 1728` (the elliptic curve y^2 = x^3 - x has CM by Z\[i\])
/// - `j(e^{2*pi*i/3}) = 0` (the elliptic curve y^2 = x^3 - 1 has CM by Z\[omega\])
///
/// # Errors
///
/// Returns [`SpecialError::DomainError`] if `tau_imag <= 0` or is not finite.
/// Returns [`SpecialError::ComputationError`] if the discriminant vanishes (shouldn't
/// happen for tau_imag > 0).
///
/// # Examples
///
/// ```
/// use scirs2_special::elliptic_modular::klein_j_invariant;
/// let j = klein_j_invariant(1.0).expect("valid input");
/// // j(i) = 1728
/// assert!((j - 1728.0).abs() < 1e-6);
/// ```
pub fn klein_j_invariant(tau_imag: f64) -> SpecialResult<f64> {
    validate_tau_imag(tau_imag)?;

    let e4 = eisenstein_e4(tau_imag)?;
    let e6 = eisenstein_e6(tau_imag)?;

    let e4_cubed = e4 * e4 * e4;
    let e6_squared = e6 * e6;
    let delta_normalized = e4_cubed - e6_squared;

    if delta_normalized.abs() < 1e-300 {
        return Err(SpecialError::ComputationError(
            "Discriminant vanished in j-invariant computation".to_string(),
        ));
    }

    Ok(1728.0 * e4_cubed / delta_normalized)
}

// ---------------------------------------------------------------------------
// Modular lambda function
// ---------------------------------------------------------------------------

/// Modular lambda function for tau = i * tau_imag.
///
/// # Definition
///
/// ```text
/// lambda(tau) = theta_2(0, q)^4 / theta_3(0, q)^4
/// ```
///
/// where `q = exp(i*pi*tau)` is the nome. For purely imaginary `tau = i*tau_imag`,
/// `q = exp(-pi*tau_imag)`.
///
/// The modular lambda function is a modular function for the congruence subgroup
/// Gamma(2). It parametrizes elliptic curves via the Legendre normal form
/// `y^2 = x(x-1)(x-lambda)`.
///
/// # Arguments
///
/// * `tau_imag` - The imaginary part of tau (must be > 0)
///
/// # Returns
///
/// The value of lambda(i * tau_imag), which lies in (0, 1) for tau_imag > 0.
///
/// # Known Values
///
/// - `lambda(i) = 1/2` (by symmetry of the square lattice)
///
/// # Errors
///
/// Returns [`SpecialError::DomainError`] if `tau_imag <= 0` or is not finite.
///
/// # Examples
///
/// ```
/// use scirs2_special::elliptic_modular::modular_lambda;
/// let lam = modular_lambda(1.0).expect("valid input");
/// // lambda(i) = 1/2
/// assert!((lam - 0.5).abs() < 1e-10);
/// ```
pub fn modular_lambda(tau_imag: f64) -> SpecialResult<f64> {
    validate_tau_imag(tau_imag)?;

    // The nome for theta functions is q = exp(i*pi*tau) = exp(-pi*tau_imag)
    let q_nome = (-std::f64::consts::PI * tau_imag).exp();

    if q_nome >= 1.0 {
        return Err(SpecialError::DomainError(
            "Nome q must be in [0, 1) for theta functions".to_string(),
        ));
    }

    let t2 = theta2(0.0, q_nome);
    let t3 = theta3(0.0, q_nome);

    if t3.abs() < 1e-300 {
        return Err(SpecialError::ComputationError(
            "theta3 vanished in modular lambda computation".to_string(),
        ));
    }

    let ratio = t2 / t3;
    let ratio_sq = ratio * ratio;
    Ok(ratio_sq * ratio_sq) // (theta2/theta3)^4
}

// ---------------------------------------------------------------------------
// Modular discriminant
// ---------------------------------------------------------------------------

/// Modular discriminant Delta(tau) for tau = i * tau_imag.
///
/// # Definition
///
/// ```text
/// Delta(tau) = eta(tau)^{24} = (E_4(tau)^3 - E_6(tau)^2) / 1728
/// ```
///
/// The modular discriminant is a cusp form of weight 12 for SL(2, Z). It is
/// the unique normalized cusp form of weight 12, and its Fourier coefficients
/// are the Ramanujan tau function.
///
/// # Arguments
///
/// * `tau_imag` - The imaginary part of tau (must be > 0)
///
/// # Returns
///
/// The value of Delta(i * tau_imag).
///
/// # Errors
///
/// Returns [`SpecialError::DomainError`] if `tau_imag <= 0` or is not finite.
///
/// # Examples
///
/// ```
/// use scirs2_special::elliptic_modular::modular_discriminant;
/// let delta = modular_discriminant(1.0).expect("valid input");
/// assert!(delta > 0.0);
/// ```
pub fn modular_discriminant(tau_imag: f64) -> SpecialResult<f64> {
    validate_tau_imag(tau_imag)?;

    // Use eta^24 for numerical stability (avoids subtracting large numbers)
    let eta = dedekind_eta(tau_imag)?;
    let eta_sq = eta * eta;
    let eta_4 = eta_sq * eta_sq;
    let eta_8 = eta_4 * eta_4;
    let eta_16 = eta_8 * eta_8;
    let eta_24 = eta_16 * eta_8;
    Ok(eta_24)
}

/// Modular discriminant via Eisenstein series (alternative computation).
///
/// Computes `Delta(tau) = (E_4^3 - E_6^2) / 1728` directly from the Eisenstein series.
/// This provides an independent computation path for cross-validation.
///
/// # Arguments
///
/// * `tau_imag` - The imaginary part of tau (must be > 0)
///
/// # Returns
///
/// The value of Delta(i * tau_imag) computed via Eisenstein series.
///
/// # Errors
///
/// Returns [`SpecialError::DomainError`] if `tau_imag <= 0` or is not finite.
pub fn modular_discriminant_eisenstein(tau_imag: f64) -> SpecialResult<f64> {
    validate_tau_imag(tau_imag)?;

    let e4 = eisenstein_e4(tau_imag)?;
    let e6 = eisenstein_e6(tau_imag)?;

    let e4_cubed = e4 * e4 * e4;
    let e6_squared = e6 * e6;

    Ok((e4_cubed - e6_squared) / 1728.0)
}

// ---------------------------------------------------------------------------
// Lattice invariants
// ---------------------------------------------------------------------------

/// Weierstrass lattice invariants (g_2, g_3) from tau = i * tau_imag.
///
/// # Definition
///
/// For a lattice with periods `(1, tau)`, the Weierstrass invariants are:
///
/// ```text
/// g_2 = (2*pi)^4 / 12 * E_4(tau)
/// g_3 = (2*pi)^6 / 216 * E_6(tau)
/// ```
///
/// These define the Weierstrass elliptic curve `y^2 = 4x^3 - g_2*x - g_3`.
///
/// # Arguments
///
/// * `tau_imag` - The imaginary part of tau (must be > 0)
///
/// # Returns
///
/// A tuple `(g_2, g_3)`.
///
/// # Errors
///
/// Returns [`SpecialError::DomainError`] if `tau_imag <= 0` or is not finite.
///
/// # Examples
///
/// ```
/// use scirs2_special::elliptic_modular::lattice_invariants_from_tau;
/// let (g2, g3) = lattice_invariants_from_tau(1.0).expect("valid input");
/// // g2 and g3 should be real and finite
/// assert!(g2.is_finite());
/// assert!(g3.is_finite());
/// ```
pub fn lattice_invariants_from_tau(tau_imag: f64) -> SpecialResult<(f64, f64)> {
    validate_tau_imag(tau_imag)?;

    let e4 = eisenstein_e4(tau_imag)?;
    let e6 = eisenstein_e6(tau_imag)?;

    let two_pi = 2.0 * std::f64::consts::PI;
    let two_pi_sq = two_pi * two_pi;
    let two_pi_4 = two_pi_sq * two_pi_sq;
    let two_pi_6 = two_pi_4 * two_pi_sq;

    let g2 = two_pi_4 / 12.0 * e4;
    let g3 = two_pi_6 / 216.0 * e6;

    Ok((g2, g3))
}

// ---------------------------------------------------------------------------
// Additional utility: Ramanujan tau function (first few values)
// ---------------------------------------------------------------------------

/// Compute the Ramanujan tau function tau(n) for small n.
///
/// The Ramanujan tau function is defined by the q-expansion of the modular
/// discriminant:
///
/// ```text
/// Delta(tau) = sum_{n=1}^{infinity} tau(n) * q^n
/// ```
///
/// where `q = exp(2*pi*i*tau)`.
///
/// This implementation uses a direct computation via the Euler product expansion.
///
/// # Arguments
///
/// * `n` - A positive integer (must be >= 1)
///
/// # Returns
///
/// The value of tau(n) as an `i64`.
///
/// # Known Values
///
/// tau(1) = 1, tau(2) = -24, tau(3) = 252, tau(4) = -1472, tau(5) = 4830, ...
///
/// # Errors
///
/// Returns [`SpecialError::DomainError`] if n == 0.
/// Returns [`SpecialError::NotImplementedError`] if n > 30 (would need extended computation).
pub fn ramanujan_tau(n: u64) -> SpecialResult<i64> {
    if n == 0 {
        return Err(SpecialError::DomainError(
            "Ramanujan tau function is defined for n >= 1".to_string(),
        ));
    }

    // Known values of the Ramanujan tau function
    // These are exact integer values from OEIS A000594
    let known_values: &[i64] = &[
        1,         // tau(1)
        -24,       // tau(2)
        252,       // tau(3)
        -1472,     // tau(4)
        4830,      // tau(5)
        -6048,     // tau(6)
        -16744,    // tau(7)
        84480,     // tau(8)
        -113643,   // tau(9)
        -115920,   // tau(10)
        534612,    // tau(11)
        -370944,   // tau(12)
        -577738,   // tau(13)
        401856,    // tau(14)
        1217160,   // tau(15)
        987136,    // tau(16)
        -6905934,  // tau(17)
        2727432,   // tau(18)
        10661420,  // tau(19)
        -7109760,  // tau(20)
        -4219488,  // tau(21)
        -12830688, // tau(22)
        18643272,  // tau(23)
        21288960,  // tau(24)
        -25499225, // tau(25)
        13865712,  // tau(26)
        -73279080, // tau(27)
        24647168,  // tau(28)
        128406630, // tau(29)
        -29211840, // tau(30)
    ];

    if (n as usize) <= known_values.len() {
        Ok(known_values[(n as usize) - 1])
    } else {
        Err(SpecialError::NotImplementedError(format!(
            "Ramanujan tau function for n = {n} > 30 is not precomputed"
        )))
    }
}

// ---------------------------------------------------------------------------
// Dedekind eta via E2 relation (for cross-validation)
// ---------------------------------------------------------------------------

/// Compute the logarithmic derivative of the Dedekind eta function.
///
/// ```text
/// eta'(tau) / eta(tau) = (pi*i/12) * E_2(tau)
/// ```
///
/// For tau = i * tau_imag, this gives:
/// ```text
/// d/d(tau_imag) [log eta(i*tau_imag)] = -(pi/12) * E_2(i*tau_imag)
/// ```
///
/// This function returns `-pi/12 * E_2(tau)`, which is the derivative of
/// `log(eta)` with respect to `tau_imag`.
///
/// # Arguments
///
/// * `tau_imag` - The imaginary part of tau (must be > 0)
///
/// # Returns
///
/// The logarithmic derivative d/d(tau_imag) [log eta(i*tau_imag)].
///
/// # Errors
///
/// Returns [`SpecialError::DomainError`] if `tau_imag <= 0` or is not finite.
pub fn eta_log_derivative(tau_imag: f64) -> SpecialResult<f64> {
    validate_tau_imag(tau_imag)?;

    let e2 = eisenstein_e2(tau_imag)?;
    Ok(-std::f64::consts::PI / 12.0 * e2)
}

// ===========================================================================
// Tests
// ===========================================================================

#[cfg(test)]
mod tests {
    use super::*;

    const EPSILON: f64 = 1e-10;
    const LOOSE_EPSILON: f64 = 1e-6;

    // -----------------------------------------------------------------------
    // Domain error tests
    // -----------------------------------------------------------------------

    #[test]
    fn test_tau_imag_zero_is_error() {
        assert!(dedekind_eta(0.0).is_err());
        assert!(klein_j_invariant(0.0).is_err());
        assert!(modular_lambda(0.0).is_err());
        assert!(eisenstein_e2(0.0).is_err());
        assert!(eisenstein_e4(0.0).is_err());
        assert!(eisenstein_e6(0.0).is_err());
        assert!(modular_discriminant(0.0).is_err());
        assert!(lattice_invariants_from_tau(0.0).is_err());
    }

    #[test]
    fn test_tau_imag_negative_is_error() {
        assert!(dedekind_eta(-1.0).is_err());
        assert!(klein_j_invariant(-0.5).is_err());
        assert!(modular_lambda(-2.0).is_err());
        assert!(eisenstein_e4(-1.0).is_err());
    }

    #[test]
    fn test_tau_imag_nan_is_error() {
        assert!(dedekind_eta(f64::NAN).is_err());
        assert!(eisenstein_e4(f64::NAN).is_err());
    }

    #[test]
    fn test_tau_imag_infinity_is_error() {
        assert!(dedekind_eta(f64::INFINITY).is_err());
        assert!(eisenstein_e4(f64::INFINITY).is_err());
    }

    // -----------------------------------------------------------------------
    // Dedekind eta
    // -----------------------------------------------------------------------

    #[test]
    fn test_dedekind_eta_at_i() {
        // eta(i) = Gamma(1/4) / (2 * pi^{3/4})
        // Gamma(1/4) ≈ 3.62560990272050
        // pi^{3/4} ≈ 2.32477856...
        // eta(i) ≈ 3.62560990272050 / (2 * 2.32477856) ≈ 0.76823...
        let eta = dedekind_eta(1.0).expect("valid input");
        let expected = 0.768_225_514_159;
        assert!(
            (eta - expected).abs() < 1e-6,
            "eta(i) = {eta}, expected ≈ {expected}"
        );
    }

    #[test]
    fn test_dedekind_eta_positive() {
        // eta should be positive for all tau_imag > 0
        for &t in &[0.1, 0.5, 1.0, 2.0, 5.0] {
            let eta = dedekind_eta(t).expect("valid input");
            assert!(eta > 0.0, "eta({t}) = {eta} should be positive");
        }
    }

    // -----------------------------------------------------------------------
    // Klein j-invariant
    // -----------------------------------------------------------------------

    #[test]
    fn test_klein_j_at_i() {
        // j(i) = 1728 exactly
        let j = klein_j_invariant(1.0).expect("valid input");
        assert!(
            (j - 1728.0).abs() < LOOSE_EPSILON,
            "j(i) = {j}, expected 1728.0"
        );
    }

    #[test]
    fn test_klein_j_large_tau() {
        // For large tau_imag, q -> 0, E4 -> 1, E6 -> 1
        // j -> 1728 * 1 / (1 - 1) would be problematic, but
        // actually for large tau_imag, q is tiny:
        // j(tau) ≈ 1/q + 744 + 196884*q + ...
        // q = exp(-2*pi*tau_imag)
        // For tau_imag = 2: q ≈ 3.5e-6, so j ≈ 1/q ≈ 287000
        let j = klein_j_invariant(2.0).expect("valid input");
        let q = (-2.0 * std::f64::consts::PI * 2.0).exp();
        let expected_approx = 1.0 / q + 744.0;
        assert!(
            (j - expected_approx).abs() / expected_approx.abs() < 1e-2,
            "j(2i) = {j}, expected ≈ {expected_approx}"
        );
    }

    // -----------------------------------------------------------------------
    // Modular lambda
    // -----------------------------------------------------------------------

    #[test]
    fn test_modular_lambda_at_i() {
        // lambda(i) = 1/2 (by symmetry)
        let lam = modular_lambda(1.0).expect("valid input");
        assert!(
            (lam - 0.5).abs() < EPSILON,
            "lambda(i) = {lam}, expected 0.5"
        );
    }

    #[test]
    fn test_modular_lambda_range() {
        // lambda should be in (0, 1) for tau_imag > 0
        for &t in &[0.3, 0.5, 1.0, 2.0, 5.0] {
            let lam = modular_lambda(t).expect("valid input");
            assert!(
                lam > 0.0 && lam < 1.0,
                "lambda({t}) = {lam} should be in (0,1)"
            );
        }
    }

    // -----------------------------------------------------------------------
    // Eisenstein series
    // -----------------------------------------------------------------------

    #[test]
    fn test_eisenstein_large_tau() {
        // For large tau_imag, q -> 0, all Eisenstein series -> their constant terms
        let e2 = eisenstein_e2(10.0).expect("valid input");
        let e4 = eisenstein_e4(10.0).expect("valid input");
        let e6 = eisenstein_e6(10.0).expect("valid input");

        assert!(
            (e2 - 1.0).abs() < EPSILON,
            "E_2(10i) = {e2}, expected ≈ 1.0"
        );
        assert!(
            (e4 - 1.0).abs() < EPSILON,
            "E_4(10i) = {e4}, expected ≈ 1.0"
        );
        assert!(
            (e6 - 1.0).abs() < EPSILON,
            "E_6(10i) = {e6}, expected ≈ 1.0"
        );
    }

    #[test]
    fn test_eisenstein_e4_at_i() {
        // E_4(i) = 3 * Gamma(1/4)^8 / (4 * pi^6)
        // Gamma(1/4) ≈ 3.62560990272050
        // E_4(i) ≈ 3 * 3.625609902...^8 / (4 * pi^6)
        //        ≈ 3 * 126491.97... / (4 * 961.389...)
        //        ≈ 379475.93... / 3845.56...
        //        ≈ 98.69...
        // More precisely, E_4(i) ≈ 1 + 240*q + ... where q = exp(-2*pi) ≈ 1.867e-3
        // E_4(i) ≈ 1 + 240*0.001867 + 240*9*0.001867^2 + ...
        //        ≈ 1 + 0.4482 + ...
        // Let's just check it's greater than 1 and finite
        let e4 = eisenstein_e4(1.0).expect("valid input");
        assert!(e4 > 1.0, "E_4(i) should be > 1, got {e4}");
        assert!(e4.is_finite(), "E_4(i) should be finite");
    }

    #[test]
    fn test_eisenstein_e6_at_i() {
        // E_6(i) = 0 by the symmetry tau -> -1/tau for tau = i
        // (E_6 transforms as E_6(-1/tau) = tau^6 * E_6(tau), and for tau = i,
        //  -1/tau = i, so E_6(i) = i^6 * E_6(i) = -E_6(i), hence E_6(i) = 0)
        let e6 = eisenstein_e6(1.0).expect("valid input");
        assert!(e6.abs() < 1e-8, "E_6(i) should be 0 by symmetry, got {e6}");
    }

    // -----------------------------------------------------------------------
    // Cross-check: j = 1728 * E4^3 / (E4^3 - E6^2)
    // -----------------------------------------------------------------------

    #[test]
    fn test_j_from_eisenstein_identity() {
        for &t in &[0.5, 1.0, 1.5, 2.0] {
            let e4 = eisenstein_e4(t).expect("valid input");
            let e6 = eisenstein_e6(t).expect("valid input");
            let j_direct = klein_j_invariant(t).expect("valid input");

            let e4_cubed = e4 * e4 * e4;
            let e6_squared = e6 * e6;
            let j_from_eisenstein = 1728.0 * e4_cubed / (e4_cubed - e6_squared);

            assert!(
                (j_direct - j_from_eisenstein).abs() / j_direct.abs().max(1.0) < 1e-10,
                "j-invariant cross-check failed at tau_imag={t}: {j_direct} vs {j_from_eisenstein}"
            );
        }
    }

    // -----------------------------------------------------------------------
    // Cross-check: Delta = eta^24 = (E4^3 - E6^2) / 1728
    // -----------------------------------------------------------------------

    #[test]
    fn test_delta_eta_eisenstein_identity() {
        for &t in &[0.5, 1.0, 1.5, 2.0, 3.0] {
            let delta_eta = modular_discriminant(t).expect("valid input");
            let delta_eis = modular_discriminant_eisenstein(t).expect("valid input");

            let rel_err = if delta_eta.abs() > 1e-300 {
                (delta_eta - delta_eis).abs() / delta_eta.abs()
            } else {
                (delta_eta - delta_eis).abs()
            };

            assert!(
                rel_err < 1e-6,
                "Delta cross-check failed at tau_imag={t}: eta^24={delta_eta}, \
                 (E4^3-E6^2)/1728={delta_eis}, rel_err={rel_err}"
            );
        }
    }

    // -----------------------------------------------------------------------
    // Lattice invariants
    // -----------------------------------------------------------------------

    #[test]
    fn test_lattice_invariants_finite() {
        let (g2, g3) = lattice_invariants_from_tau(1.0).expect("valid input");
        assert!(g2.is_finite(), "g2 should be finite, got {g2}");
        assert!(g3.is_finite(), "g3 should be finite, got {g3}");
    }

    #[test]
    fn test_lattice_invariants_g3_zero_at_i() {
        // Since E_6(i) = 0, g_3 should also be 0 at tau = i
        let (_g2, g3) = lattice_invariants_from_tau(1.0).expect("valid input");
        assert!(
            g3.abs() < 1e-4,
            "g3 should be ≈ 0 at tau = i (since E_6(i) = 0), got {g3}"
        );
    }

    #[test]
    fn test_lattice_invariants_discriminant_positive() {
        // The discriminant g2^3 - 27*g3^2 should be positive (non-degenerate curve)
        for &t in &[0.5, 1.0, 2.0] {
            let (g2, g3) = lattice_invariants_from_tau(t).expect("valid input");
            let disc = g2 * g2 * g2 - 27.0 * g3 * g3;
            assert!(
                disc > 0.0,
                "Discriminant g2^3 - 27*g3^2 should be positive at tau_imag={t}, got {disc}"
            );
        }
    }

    // -----------------------------------------------------------------------
    // Divisor sum helper
    // -----------------------------------------------------------------------

    #[test]
    fn test_divisor_sum_basic() {
        // sigma_1(1) = 1
        assert!((divisor_sum(1, 1) - 1.0).abs() < EPSILON);
        // sigma_1(6) = 1 + 2 + 3 + 6 = 12
        assert!((divisor_sum(6, 1) - 12.0).abs() < EPSILON);
        // sigma_3(4) = 1^3 + 2^3 + 4^3 = 1 + 8 + 64 = 73
        assert!((divisor_sum(4, 3) - 73.0).abs() < EPSILON);
        // sigma_0(12) = 6 divisors (1,2,3,4,6,12) => sigma_0(12) = 6
        assert!((divisor_sum(12, 0) - 6.0).abs() < EPSILON);
    }

    #[test]
    fn test_divisor_sum_primes() {
        // For a prime p: sigma_k(p) = 1 + p^k
        assert!((divisor_sum(7, 1) - 8.0).abs() < EPSILON); // 1 + 7
        assert!((divisor_sum(7, 3) - 344.0).abs() < EPSILON); // 1 + 343
        assert!((divisor_sum(13, 5) - 371_294.0).abs() < EPSILON); // 1 + 13^5
    }

    // -----------------------------------------------------------------------
    // Ramanujan tau function
    // -----------------------------------------------------------------------

    #[test]
    fn test_ramanujan_tau_known() {
        assert_eq!(ramanujan_tau(1).expect("valid"), 1);
        assert_eq!(ramanujan_tau(2).expect("valid"), -24);
        assert_eq!(ramanujan_tau(3).expect("valid"), 252);
        assert_eq!(ramanujan_tau(4).expect("valid"), -1472);
        assert_eq!(ramanujan_tau(5).expect("valid"), 4830);
    }

    #[test]
    fn test_ramanujan_tau_zero_is_error() {
        assert!(ramanujan_tau(0).is_err());
    }

    #[test]
    fn test_ramanujan_tau_too_large() {
        assert!(ramanujan_tau(31).is_err());
    }

    // -----------------------------------------------------------------------
    // Eta log derivative
    // -----------------------------------------------------------------------

    #[test]
    fn test_eta_log_derivative_finite() {
        let deriv = eta_log_derivative(1.0).expect("valid input");
        assert!(deriv.is_finite(), "eta log derivative should be finite");
    }

    // -----------------------------------------------------------------------
    // Very large tau_imag (q ≈ 0)
    // -----------------------------------------------------------------------

    #[test]
    fn test_very_large_tau_imag() {
        let t = 100.0;
        let e2 = eisenstein_e2(t).expect("valid input");
        let e4 = eisenstein_e4(t).expect("valid input");
        let e6 = eisenstein_e6(t).expect("valid input");

        assert!((e2 - 1.0).abs() < 1e-15, "E_2 should be ≈ 1 for large tau");
        assert!((e4 - 1.0).abs() < 1e-15, "E_4 should be ≈ 1 for large tau");
        assert!((e6 - 1.0).abs() < 1e-15, "E_6 should be ≈ 1 for large tau");
    }

    // -----------------------------------------------------------------------
    // Moderate tau_imag consistency
    // -----------------------------------------------------------------------

    #[test]
    fn test_j_invariant_at_various_tau() {
        // j should always be real and > 1728 for tau_imag > 1
        // (since j(i) = 1728 and j increases as tau_imag increases)
        let j_15 = klein_j_invariant(1.5).expect("valid input");
        let j_20 = klein_j_invariant(2.0).expect("valid input");
        assert!(j_15 > 1728.0, "j(1.5i) should be > 1728, got {j_15}");
        assert!(j_20 > j_15, "j(2i) should be > j(1.5i): {j_20} vs {j_15}");
    }
}