geographdb-core 0.5.4

//! Number-theoretic algorithms for GeoGraphDB Core
//!
//! Inspired by PARI/GP functionality, focused on:
//! - Prime number generation and testing
//! - Integer factorization (including semiprime weak-family detection)
//! - Modular arithmetic
//! - Riemann zeta function approximations
//! - Continued fractions
//!
//! # Crop Circle Connection
//!
//! The Arecibo Reply (2001) uses semiprime dimensions 23×73=1679.
//! The ratio q/p ≈ π suggests a fast factoring algorithm for this family.

use std::collections::HashMap;

// ── Prime utilities ───────────────────────────────────────────────────────────

/// Test if n is prime using trial division.
/// For n < 2^64, deterministic Miller-Rabin would be better,
/// but trial division is sufficient for demo purposes.
pub fn is_prime(n: u64) -> bool {
    let mut result = true;
    if n < 2 {
        result = false;
    } else if n == 2 || n == 3 {
        result = true;
    } else if n.is_multiple_of(2) {
        result = false;
    } else {
        let mut i = 3;
        while i * i <= n {
            if n.is_multiple_of(i) {
                result = false;
                break;
            }
            i += 2;
        }
    }
    result
}

/// Deterministic Miller-Rabin primality test for n < 2^64.
///
/// Uses known witness sets that are sufficient for all 64-bit integers.
/// For n < 3,317,044,064,679,887,385,641, the witnesses {2, 3, 5, 7, 11, 13, 17}
/// are sufficient (Jaeschke 1993, extended by deterministic bounds).
///
/// # Example
/// ```
/// use geographdb_core::algorithms::number_theory::is_prime_miller_rabin;
/// assert!(is_prime_miller_rabin(409));
/// assert!(!is_prime_miller_rabin(1679));
/// ```
pub fn is_prime_miller_rabin(n: u64) -> bool {
    let mut result = true;
    if n < 2 {
        result = false;
    } else {
        // Small primes check
        let small_primes: [u64; 12] = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37];
        let mut divisible_by_small = false;
        for &p in &small_primes {
            if n == p {
                divisible_by_small = true;
                result = true;
                break;
            }
            if n.is_multiple_of(p) {
                divisible_by_small = true;
                result = false;
                break;
            }
        }

        if !divisible_by_small {
            // Write n-1 as d * 2^s
            let mut d = n - 1;
            let mut s = 0;
            while d.is_multiple_of(2) {
                d /= 2;
                s += 1;
            }

            // Witnesses sufficient for n < 2^64 (Jim Sinclair, 2011)
            let witnesses: [u64; 12] = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37];

            let mut composite = false;
            for &a in &witnesses {
                if a >= n {
                    continue;
                }
                let mut x = mod_pow(a, d, n);
                if x == 1 || x == n - 1 {
                    continue;
                }
                let mut witness_passed = false;
                for _ in 0..(s - 1) {
                    x = mul_mod(x, x, n);
                    if x == n - 1 {
                        witness_passed = true;
                        break;
                    }
                }
                if !witness_passed {
                    composite = true;
                    break;
                }
            }
            result = !composite;
        }
    }
    result
}

/// Modular multiplication (a * b) % m without overflow for u64.
/// Uses the fact that 128-bit intermediate is available on modern CPUs.
fn mul_mod(a: u64, b: u64, m: u64) -> u64 {
    ((a as u128) * (b as u128) % (m as u128)) as u64
}

// ── Advanced Factorization (PARI/GP-level) ────────────────────────────────────

/// Pollard's Rho factorization algorithm.
///
/// A probabilistic algorithm for finding a non-trivial factor of n.
/// Expected time: O(n^(1/4)) for finding a factor p where p ≈ √n.
/// Much faster than trial division for large semiprimes.
///
/// # Algorithm
/// Uses Floyd's cycle detection on the pseudo-random sequence:
/// x_{i+1} = f(x_i) mod n where f(x) = x² + c
///
/// # Example
/// ```
/// use geographdb_core::algorithms::number_theory::pollard_rho;
/// let factor = pollard_rho(1679, 1);
/// assert!(factor.is_some());
/// assert_eq!(1679 % factor.unwrap(), 0);
/// ```
pub fn pollard_rho(n: u64, seed: u64) -> Option<u64> {
    if n.is_multiple_of(2) {
        return Some(2);
    }
    if is_prime_miller_rabin(n) {
        return None;
    }

    let f = |x: u64| mul_mod(x, x, n).wrapping_add(seed) % n;
    let mut x = 2u64;
    let mut y = 2u64;
    let mut d = 1u64;

    while d == 1 {
        x = f(x);
        y = f(f(y));
        let diff = x.abs_diff(y);
        d = gcd(diff, n);
    }

    if d == n {
        None
    } else {
        Some(d)
    }
}

/// Greatest common divisor using Euclid's algorithm.
pub fn gcd(mut a: u64, mut b: u64) -> u64 {
    while b != 0 {
        let t = b;
        b = a % b;
        a = t;
    }
    a
}

/// Fully factor n using a combination of trial division and Pollard's Rho.
///
/// For small n: uses trial division.
/// For large n: uses Pollard's Rho with Miller-Rabin primality testing.
///
/// Returns a map of prime -> exponent.
pub fn factor_advanced(n: u64) -> HashMap<u64, u32> {
    let mut factors = HashMap::new();
    let mut n = n;

    // Trial division for small factors
    let small_primes = sieve_primes(10000);
    for &p in &small_primes {
        if p * p > n {
            break;
        }
        while n.is_multiple_of(p) {
            *factors.entry(p).or_insert(0) += 1;
            n /= p;
        }
    }

    // Pollard's Rho for remaining large factors
    let mut stack = vec![n];
    while let Some(m) = stack.pop() {
        if m == 1 {
            continue;
        }
        if is_prime_miller_rabin(m) {
            *factors.entry(m).or_insert(0) += 1;
            continue;
        }

        // Try Pollard's Rho with different seeds
        let mut found = false;
        for seed in 1..=10 {
            if let Some(d) = pollard_rho(m, seed) {
                stack.push(d);
                stack.push(m / d);
                found = true;
                break;
            }
        }

        if !found {
            // Fallback to trial division
            let tf = factor_trial(m);
            for (p, exp) in tf {
                *factors.entry(p).or_insert(0) += exp;
            }
        }
    }

    factors
}

/// Quadratic Sieve factorization.
///
/// The Quadratic Sieve is the second-fastest known factorization algorithm
/// (after the General Number Field Sieve). It factors integers up to ~100 digits.
///
/// # Algorithm
/// 1. Choose smoothness bound B
/// 2. Build factor base of primes ≤ B
/// 3. Find relations x² ≡ y (mod n) where y is B-smooth
/// 4. Solve linear system over GF(2) to find subset product that is a square
/// 5. Compute gcd(x ± √y, n) to find factors
///
/// # Reference
/// See Pomerance, "The Quadratic Sieve Factoring Algorithm" (1984).
pub fn factor_quadratic_sieve(n: u64) -> Option<HashMap<u64, u32>> {
    if n < 2 {
        return None;
    }
    if n.is_multiple_of(2) {
        let mut result = HashMap::new();
        result.insert(2, 1);
        if let Some(rest) = factor_quadratic_sieve(n / 2) {
            for (p, e) in rest {
                *result.entry(p).or_insert(0) += e;
            }
        }
        return Some(result);
    }

    // Step 1: Smoothness bound estimation
    // B ≈ exp(0.5 * sqrt(ln(n) * ln(ln(n))))
    let nf = n as f64;
    let ln_n = nf.ln();
    let ln_ln_n = ln_n.ln();
    let b_est = (0.5 * (ln_n * ln_ln_n).sqrt()).exp() as usize;
    let b = b_est.clamp(10, 1000); // Clamp for practical bounds

    // Step 2: Build factor base (primes p where legendre_symbol(n, p) = 1)
    let primes = sieve_primes(b);
    let mut factor_base: Vec<u64> = Vec::new();
    for &p in &primes {
        if p == 2 {
            continue;
        }
        if legendre_symbol(n, p) == 1 || n.is_multiple_of(p) {
            factor_base.push(p);
        }
    }

    if factor_base.is_empty() {
        // Fallback to trial division for very small n
        return Some(factor_trial(n));
    }

    // Step 3: Sieving for relations
    // Find x such that (x² mod n) is B-smooth
    let m = (nf.sqrt().floor() as u64) + 1;
    let sieve_range = 10000usize; // Number of x values to test
    let mut relations: Vec<(u64, Vec<u32>)> = Vec::new(); // (x, exponents vector)

    for i in 0..sieve_range {
        let x = m + i as u64;
        let y = if x * x > n { (x * x) % n } else { x * x };
        let _y_f = y as f64;

        // Check if y is B-smooth over factor base
        let mut exponents = vec![0u32; factor_base.len()];
        let mut remaining = y;
        let mut _is_smooth = true;

        for (j, &p) in factor_base.iter().enumerate() {
            while remaining % p == 0 {
                remaining /= p;
                exponents[j] += 1;
            }
        }

        if remaining == 1 {
            relations.push((x, exponents));
        }

        if relations.len() >= factor_base.len() + 5 {
            break;
        }
    }

    if relations.len() < factor_base.len() {
        // Not enough relations found — fallback to Pollard's Rho
        return Some(factor_advanced(n));
    }

    // Step 4: Gaussian elimination over GF(2) to find dependent rows
    // Simplified: try random subsets to find a square product
    let mut rng_seed = 1u64;
    for _ in 0..1000 {
        rng_seed = rng_seed.wrapping_mul(1103515245).wrapping_add(12345);
        let mask = rng_seed;
        let mut combined_exp = vec![0u32; factor_base.len()];
        let mut x_product = 1u64;
        let mut count = 0u32;

        for (i, (x, exp)) in relations.iter().enumerate() {
            if (mask >> (i % 64)) & 1 == 1 {
                for j in 0..factor_base.len() {
                    combined_exp[j] += exp[j];
                }
                x_product = mul_mod(x_product, *x, n);
                count += 1;
            }
        }

        // Check if all exponents are even (square)
        let all_even = combined_exp.iter().all(|&e| e % 2 == 0);
        if all_even && count > 0 {
            // Compute y = product of p^(e/2)
            let mut y_sqrt = 1u64;
            for (j, &p) in factor_base.iter().enumerate() {
                let half_exp = combined_exp[j] / 2;
                for _ in 0..half_exp {
                    y_sqrt = mul_mod(y_sqrt, p, n);
                }
            }

            // Step 5: Compute gcd(x ± y, n)
            let diff = x_product.abs_diff(y_sqrt);
            let factor1 = gcd(diff, n);

            if factor1 > 1 && factor1 < n {
                let mut result = HashMap::new();
                let factor2 = n / factor1;

                // Recursively factor if needed
                if is_prime_miller_rabin(factor1) {
                    result.insert(factor1, 1);
                } else if let Some(sub) = factor_quadratic_sieve(factor1) {
                    for (p, e) in sub {
                        *result.entry(p).or_insert(0) += e;
                    }
                }

                if is_prime_miller_rabin(factor2) {
                    result.insert(factor2, 1);
                } else if let Some(sub) = factor_quadratic_sieve(factor2) {
                    for (p, e) in sub {
                        *result.entry(p).or_insert(0) += e;
                    }
                }

                return Some(result);
            }
        }
    }

    // Fallback if QS fails to find factors
    Some(factor_advanced(n))
}

/// Elliptic Curve Method (ECM) factorization.
///
/// ECM is particularly effective for finding medium-sized factors (20-30 digits)
/// of large numbers. It uses the group structure of elliptic curves over finite fields.
///
/// # Algorithm
/// 1. Choose random elliptic curve E: y² = x³ + ax + b mod n and point P on E
/// 2. Compute kP where k is product of small primes raised to powers
/// 3. If the computation fails (division by non-invertible element), gcd gives a factor
/// 4. Otherwise retry with different curve
///
/// # Reference
/// See Lenstra, "Factoring Integers with Elliptic Curves" (1987).
pub fn factor_ecm(n: u64) -> Option<HashMap<u64, u32>> {
    if n < 2 {
        return None;
    }
    if n.is_multiple_of(2) {
        let mut result = HashMap::new();
        result.insert(2, 1);
        if let Some(rest) = factor_ecm(n / 2) {
            for (p, e) in rest {
                *result.entry(p).or_insert(0) += e;
            }
        }
        return Some(result);
    }

    // Stage 1: Try multiple curves
    for curve_seed in 1..=50 {
        // Random curve: y² = x³ + ax + b, with point P = (x, y)
        let a = (curve_seed * 1234567 + 1) % n;
        let x = (curve_seed * 7654321 + 2) % n;
        let y = (curve_seed * 9876543 + 3) % n;

        // Verify point is on curve: y² ≡ x³ + ax + b (mod n)
        let _b = if let Some(b_val) = mod_sub(
            mul_mod(y, y, n),
            add_mod(mul_mod(mul_mod(x, x, n), x, n), mul_mod(a, x, n), n),
            n,
        ) {
            b_val
        } else {
            continue;
        };

        // Bounds for stage 1 (B1 determines which primes we multiply)
        let b1 = if n < 10000 {
            100
        } else if n < 1000000 {
            1000
        } else {
            10000
        };

        // Compute k = lcm(1, 2, ..., B1) as product of prime powers
        let primes = sieve_primes(b1 as usize);
        let mut k = 1u64;
        for &p in &primes {
            if p == 0 {
                continue;
            }
            let mut pp = p;
            while pp * p <= b1 as u64 {
                pp *= p;
            }
            k = mul_mod(k, pp, n);
        }

        // Compute kP using elliptic curve point multiplication
        let mut px = x;
        let mut py = y;
        let mut qx = 0u64; // Point at infinity
        let mut qy = 0u64;

        let mut kk = k;
        while kk > 0 {
            if kk % 2 == 1 {
                // Q = Q + P
                if qx == 0 && qy == 0 {
                    qx = px;
                    qy = py;
                } else {
                    let (rx, ry, success) = ec_add(qx, qy, px, py, a, n);
                    if !success {
                        // Division failed — found a factor!
                        let factor = gcd(rx, n);
                        if factor > 1 && factor < n {
                            let mut result = HashMap::new();
                            let other = n / factor;

                            if is_prime_miller_rabin(factor) {
                                result.insert(factor, 1);
                            } else if let Some(sub) = factor_ecm(factor) {
                                for (p, e) in sub {
                                    *result.entry(p).or_insert(0) += e;
                                }
                            }

                            if is_prime_miller_rabin(other) {
                                result.insert(other, 1);
                            } else if let Some(sub) = factor_ecm(other) {
                                for (p, e) in sub {
                                    *result.entry(p).or_insert(0) += e;
                                }
                            }

                            return Some(result);
                        }
                        break;
                    }
                    qx = rx;
                    qy = ry;
                }
            }

            // P = 2P
            let (rx, ry, success) = ec_double(px, py, a, n);
            if !success {
                let factor = gcd(rx, n);
                if factor > 1 && factor < n {
                    let mut result = HashMap::new();
                    let other = n / factor;

                    if is_prime_miller_rabin(factor) {
                        result.insert(factor, 1);
                    } else if let Some(sub) = factor_ecm(factor) {
                        for (p, e) in sub {
                            *result.entry(p).or_insert(0) += e;
                        }
                    }

                    if is_prime_miller_rabin(other) {
                        result.insert(other, 1);
                    } else if let Some(sub) = factor_ecm(other) {
                        for (p, e) in sub {
                            *result.entry(p).or_insert(0) += e;
                        }
                    }

                    return Some(result);
                }
                break;
            }
            px = rx;
            py = ry;
            kk /= 2;
        }
    }

    // Fallback if ECM fails
    Some(factor_advanced(n))
}

// ── Elliptic Curve Arithmetic ─────────────────────────────────────────────────

/// Add two points on an elliptic curve: R = P + Q.
///
/// Curve: y² = x³ + ax + b (mod n)
/// Returns (rx, ry, success). If success is false, rx contains the
/// non-invertible element (potential factor).
fn ec_add(px: u64, py: u64, qx: u64, qy: u64, a: u64, n: u64) -> (u64, u64, bool) {
    let _ = a; // a not needed for point addition formula

    if px == qx && py == qy {
        return ec_double(px, py, a, n);
    }

    if px == qx {
        // P = -Q, result is point at infinity
        return (0, 0, true);
    }

    // Slope λ = (qy - py) / (qx - px) mod n
    let num = if qy >= py { qy - py } else { qy + n - py };
    let den = if qx >= px { qx - px } else { qx + n - px };

    let inv_den = match mod_inverse(den, n) {
        Some(inv) => inv,
        None => return (den, 0, false), // Non-invertible: potential factor
    };

    let lambda = mul_mod(num, inv_den, n);

    // rx = λ² - px - qx
    let lambda_sq = mul_mod(lambda, lambda, n);
    let rx = if lambda_sq >= px + qx {
        lambda_sq - px - qx
    } else {
        lambda_sq + n - ((px + qx) % n)
    };
    let rx = rx % n;

    // ry = λ(px - rx) - py
    let diff = if px >= rx { px - rx } else { px + n - rx };
    let ry_term = mul_mod(lambda, diff, n);
    let ry = if ry_term >= py {
        ry_term - py
    } else {
        ry_term + n - py
    };
    let ry = ry % n;

    (rx, ry, true)
}

/// Double a point on an elliptic curve: R = 2P.
///
/// Curve: y² = x³ + ax + b (mod n)
/// Returns (rx, ry, success). If success is false, rx contains the
/// non-invertible element (potential factor).
fn ec_double(px: u64, py: u64, a: u64, n: u64) -> (u64, u64, bool) {
    if py == 0 {
        return (0, 0, true); // Point at infinity
    }

    // Slope λ = (3x² + a) / (2y) mod n
    let x_sq = mul_mod(px, px, n);
    let three_x_sq = mul_mod(3, x_sq, n);
    let num = add_mod(three_x_sq, a, n);
    let den = mul_mod(2, py, n);

    let inv_den = match mod_inverse(den, n) {
        Some(inv) => inv,
        None => return (den, 0, false), // Non-invertible: potential factor
    };

    let lambda = mul_mod(num, inv_den, n);

    // rx = λ² - 2px
    let lambda_sq = mul_mod(lambda, lambda, n);
    let two_px = mul_mod(2, px, n);
    let rx = if lambda_sq >= two_px {
        lambda_sq - two_px
    } else {
        lambda_sq + n - two_px
    };
    let rx = rx % n;

    // ry = λ(px - rx) - py
    let diff = if px >= rx { px - rx } else { px + n - rx };
    let ry_term = mul_mod(lambda, diff, n);
    let ry = if ry_term >= py {
        ry_term - py
    } else {
        ry_term + n - py
    };
    let ry = ry % n;

    (rx, ry, true)
}

/// Safe modular addition: (a + b) % m.
fn add_mod(a: u64, b: u64, m: u64) -> u64 {
    ((a as u128 + b as u128) % (m as u128)) as u64
}

/// Safe modular subtraction: (a - b) % m.
fn mod_sub(a: u64, b: u64, m: u64) -> Option<u64> {
    if m == 0 {
        return None;
    }
    Some(((a as u128 + m as u128 - (b % m) as u128) % (m as u128)) as u64)
}

// ── Riemann Zeta Function ─────────────────────────────────────────────────────
/// Generate all primes up to n using the Sieve of Eratosthenes.
pub fn sieve_primes(n: usize) -> Vec<u64> {
    let mut is_prime = vec![true; n + 1];
    is_prime[0] = false;
    is_prime[1] = false;
    for i in 2..=((n as f64).sqrt() as usize) {
        if is_prime[i] {
            let mut j = i * i;
            while j <= n {
                is_prime[j] = false;
                j += i;
            }
        }
    }
    is_prime
        .iter()
        .enumerate()
        .filter(|(_, &p)| p)
        .map(|(i, _)| i as u64)
        .collect()
}

/// Return the n-th prime (1-indexed: nth_prime(1) = 2).
pub fn nth_prime(n: usize) -> u64 {
    if n == 0 {
        return 0;
    }
    // Overestimate using p_n ~ n * (ln n + ln ln n)
    let estimate = if n < 6 {
        15
    } else {
        let nf = n as f64;
        (nf * (nf.ln() + nf.ln().ln()) * 1.5) as usize + 100
    };
    let primes = sieve_primes(estimate);
    primes[n - 1]
}

// ── Factorization ─────────────────────────────────────────────────────────────

/// Trial division factorization.
/// Returns a map of prime -> exponent.
pub fn factor_trial(n: u64) -> HashMap<u64, u32> {
    let mut factors = HashMap::new();
    let mut n = n;
    let mut p = 2;
    while p * p <= n {
        while n.is_multiple_of(p) {
            *factors.entry(p).or_insert(0) += 1;
            n /= p;
        }
        p += if p == 2 { 1 } else { 2 };
    }
    if n > 1 {
        *factors.entry(n).or_insert(0) += 1;
    }
    factors
}

/// Check if n is a semiprime (product of exactly two primes).
pub fn is_semiprime(n: u64) -> bool {
    let factors = factor_trial(n);
    factors.values().sum::<u32>() == 2
}

/// Find the semiprime factors of n.
/// Returns (p, q) such that n = p * q, or None if not semiprime.
pub fn semiprime_factors(n: u64) -> Option<(u64, u64)> {
    let factors = factor_trial(n);
    if factors.values().sum::<u32>() != 2 {
        return None;
    }
    let mut result = Vec::new();
    for (p, exp) in factors {
        for _ in 0..exp {
            result.push(p);
        }
    }
    if result.len() == 2 {
        Some((result[0], result[1]))
    } else {
        None
    }
}

// ── Crop Circle: Arecibo Reply Factoring ──────────────────────────────────────

/// Fast factoring for semiprimes where q/p ≈ π.
///
/// # Algorithm
/// 1. Compute p_est = √(n/π)
/// 2. Test integers in [p_est - window, p_est + window]
/// 3. Return (p, q) if found
///
/// # Crop Circle Connection
/// The Arecibo Reply uses 23×73 = 1679, where 73/23 ≈ 3.174 ≈ π.
/// This algorithm factors such semiprimes in O(window) time.
pub fn factor_semiprime_pi_family(n: u64, window: u64) -> Option<(u64, u64)> {
    let pi = std::f64::consts::PI;
    let p_est = ((n as f64) / pi).sqrt();
    let p_floor = p_est as u64;

    let start = if p_floor > window {
        p_floor - window
    } else {
        2
    };
    let end = p_floor + window;

    for p in start..=end {
        if n.is_multiple_of(p) {
            let q = n / p;
            // Verify it's actually a semiprime factorization
            if p * q == n && p > 1 && q > 1 {
                return Some((p, q));
            }
        }
    }
    None
}

/// Compute the ratio q/p for a semiprime and compare to π.
pub fn semiprime_pi_ratio(n: u64) -> Option<f64> {
    let (p, q) = semiprime_factors(n)?;
    let (p, q) = if p < q { (p, q) } else { (q, p) };
    Some(q as f64 / p as f64)
}

// ── Riemann Zeta Function ─────────────────────────────────────────────────────

/// Approximate the Riemann zeta function ζ(s) for real s > 1.
/// Uses the series: ζ(s) = Σ_{n=1}^∞ 1/n^s
pub fn zeta_real(s: f64, terms: usize) -> f64 {
    if s <= 1.0 {
        return f64::NAN;
    }
    let mut sum = 0.0;
    for n in 1..=terms {
        sum += 1.0 / (n as f64).powf(s);
    }
    sum
}

/// Approximate the prime-counting function π(x) using the logarithmic integral.
/// li(x) = ∫_0^x dt / ln(t)
pub fn logarithmic_integral(x: f64) -> f64 {
    if x <= 1.0 {
        return 0.0;
    }
    // Simple trapezoidal integration
    let n = 10000;
    let a = 1.0 + 1e-10; // Avoid singularity at 1
    let b = x;
    let h = (b - a) / n as f64;
    let mut sum = 0.5 * (1.0 / a.ln() + 1.0 / b.ln());
    for i in 1..n {
        let t = a + i as f64 * h;
        sum += 1.0 / t.ln();
    }
    sum * h
}

/// Approximate the n-th zero of the Riemann zeta function.
/// For large n: t_n ≈ 2πn / ln(n)
pub fn zeta_zero_approx(n: u64) -> f64 {
    if n == 0 {
        return 0.0;
    }
    let nf = n as f64;
    2.0 * std::f64::consts::PI * nf / nf.ln()
}

/// Approximate the number of zeros up to height T.
/// N(T) ≈ (T/2π) · ln(T/2π) - T/2π
pub fn zeta_zero_count_approx(t: f64) -> f64 {
    if t <= 0.0 {
        return 0.0;
    }
    let t_norm = t / (2.0 * std::f64::consts::PI);
    t_norm * t_norm.ln() - t_norm
}

// ── Continued Fractions ───────────────────────────────────────────────────────

/// Compute the continued fraction expansion of a real number.
/// Returns the convergents as (numerator, denominator) pairs.
pub fn continued_fraction(x: f64, max_terms: usize) -> Vec<(u64, u64)> {
    let mut convergents = Vec::new();
    let mut x = x;
    let mut a0 = x as u64;
    let mut p0 = a0;
    let mut q0 = 1;
    convergents.push((p0, q0));

    if (x - a0 as f64).abs() < 1e-15 {
        return convergents;
    }

    let mut p1 = 1;
    let mut q1 = 0;

    for _ in 1..max_terms {
        x = 1.0 / (x - a0 as f64);
        if !x.is_finite() {
            break;
        }
        a0 = x as u64;

        let p2 = a0 * p0 + p1;
        let q2 = a0 * q0 + q1;

        convergents.push((p2, q2));

        p1 = p0;
        q1 = q0;
        p0 = p2;
        q0 = q2;
    }

    convergents
}

/// Find the best rational approximation to x with denominator ≤ max_den.
pub fn best_rational_approx(x: f64, max_den: u64) -> (u64, u64) {
    let conv = continued_fraction(x, 50);
    let mut best = conv[0];
    for (p, q) in conv {
        if q > max_den {
            break;
        }
        best = (p, q);
    }
    best
}

// ── Modular Arithmetic ────────────────────────────────────────────────────────

/// Compute a^b mod m using binary exponentiation.
pub fn mod_pow(mut a: u64, mut b: u64, m: u64) -> u64 {
    if m == 1 {
        return 0;
    }
    let mut result = 1;
    a %= m;
    while b > 0 {
        if b % 2 == 1 {
            result = (result * a) % m;
        }
        a = (a * a) % m;
        b /= 2;
    }
    result
}

/// Compute the modular multiplicative inverse of a mod m.
/// Returns None if inverse doesn't exist.
pub fn mod_inverse(a: u64, m: u64) -> Option<u64> {
    // Extended Euclidean algorithm
    let (mut t, mut new_t) = (0i64, 1i64);
    let (mut r, mut new_r) = (m as i64, a as i64);

    while new_r != 0 {
        let quotient = r / new_r;
        let tmp_t = t - quotient * new_t;
        t = new_t;
        new_t = tmp_t;

        let tmp_r = r - quotient * new_r;
        r = new_r;
        new_r = tmp_r;
    }

    if r > 1 {
        return None;
    }
    if t < 0 {
        t += m as i64;
    }
    Some(t as u64)
}

/// Compute the Legendre symbol (a/p) using Euler's criterion.
/// Returns 1 if a is quadratic residue mod p, -1 if non-residue, 0 if a ≡ 0 mod p.
pub fn legendre_symbol(a: u64, p: u64) -> i32 {
    if p == 2 {
        return 1;
    }
    let a = a % p;
    if a == 0 {
        return 0;
    }
    let result = mod_pow(a, (p - 1) / 2, p);
    if result == 1 {
        1
    } else {
        -1
    }
}

// ── Riemann-Siegel Formula ────────────────────────────────────────────────────

/// Riemann-Siegel theta function.
///
/// θ(t) = arg Γ(1/4 + it/2) - (t/2) ln π
///
/// This is the phase function used in the Riemann-Siegel formula
/// for computing zeta zeros on the critical line.
///
/// # Approximation
/// Uses Stirling's approximation for the Gamma function.
pub fn riemann_siegel_theta(t: f64) -> f64 {
    if t <= 0.0 {
        return 0.0;
    }

    // Stirling approximation: arg Γ(1/4 + it/2)
    // ≈ (t/2) ln(t/2π) - t/2 - π/8 + 1/(6t) + ...
    let t2 = t / 2.0;
    t2 * (t2 / std::f64::consts::PI).ln() - t2 - std::f64::consts::PI / 8.0 + 1.0 / (48.0 * t)
}

/// Riemann-Siegel Z function.
///
/// Z(t) = e^{iθ(t)} ζ(1/2 + it)
///
/// Z(t) is real-valued for real t. Zeros of Z(t) correspond to
/// zeros of ζ(s) on the critical line Re(s) = 1/2.
///
/// # Approximation
/// Uses the first few terms of the Riemann-Siegel formula.
pub fn riemann_siegel_z(t: f64, _terms: usize) -> f64 {
    if t <= 0.0 {
        return 0.0;
    }

    let n = ((t / (2.0 * std::f64::consts::PI)).sqrt().floor() as usize).max(1);
    let theta = riemann_siegel_theta(t);

    // Main sum: 2 Σ_{k=1}^n cos(θ(t) - t ln k) / √k
    let mut sum = 0.0;
    for k in 1..=n {
        sum += (theta - t * (k as f64).ln()).cos() / (k as f64).sqrt();
    }

    2.0 * sum
}

/// Find approximate location of the n-th zero of ζ(s) on the critical line.
///
/// Uses the Riemann-Siegel Z function and Newton's method refinement.
///
/// # Algorithm
/// 1. Start with Gram point approximation: t_n ≈ 2πn / ln(n)
/// 2. Evaluate Z(t) near the approximation
/// 3. Use bisection to find sign changes
///
/// # Accuracy
/// Good to ~0.1 for small n, degrades for large n.
pub fn zeta_zero_riemann_siegel(n: u64, iterations: usize) -> f64 {
    if n == 0 {
        return 0.0;
    }

    // Initial approximation
    let mut t = zeta_zero_approx(n);
    let dt = 0.5;

    // Refine using bisection on Z(t)
    for _ in 0..iterations {
        let z_left = riemann_siegel_z(t - dt, 3);
        let z_mid = riemann_siegel_z(t, 3);
        let z_right = riemann_siegel_z(t + dt, 3);

        if z_left.signum() != z_mid.signum() {
            t -= dt / 2.0;
        } else if z_mid.signum() != z_right.signum() {
            t += dt / 2.0;
        }
    }

    t
}

// ── Modular Forms / Monster Moonshine ─────────────────────────────────────────

/// Compute the Dedekind eta function η(τ) for τ in the upper half-plane.
///
/// η(τ) = q^{1/24} Π_{n=1}^∞ (1 - q^n) where q = e^{2πiτ}
///
/// The eta function is a modular form of weight 1/2.
/// It appears in the Monster moonshine connection via the j-invariant.
///
/// # Crop Circle Connection
/// The 13 circles of Metatron's Cube may encode the 13³ factor
/// in the Monster group order. Modular forms connect the Monster
/// to the j-function, which has coefficients related to the
/// dimensions of Monster irreducible representations.
pub fn dedekind_eta(tau_real: f64, tau_imag: f64, terms: usize) -> (f64, f64) {
    if tau_imag <= 0.0 {
        return (0.0, 0.0);
    }

    let q = (-2.0 * std::f64::consts::PI * tau_imag).exp()
        * (2.0 * std::f64::consts::PI * tau_real).cos();
    let _q_imag = (-2.0 * std::f64::consts::PI * tau_imag).exp()
        * (2.0 * std::f64::consts::PI * tau_real).sin();

    // q^{1/24}
    let q24_real = (std::f64::consts::PI * tau_imag / 12.0).exp()
        * (std::f64::consts::PI * tau_real / 12.0).cos();
    let q24_imag = (std::f64::consts::PI * tau_imag / 12.0).exp()
        * (std::f64::consts::PI * tau_real / 12.0).sin();

    // Product Π (1 - q^n)
    let mut prod_real = 1.0;
    let mut prod_imag = 0.0;

    for n in 1..=terms {
        let qn = q.powi(n as i32);
        let qn_imag = 0.0; // Simplified: assumes q is real

        // (1 - q^n) as complex
        let factor_real = 1.0 - qn;
        let factor_imag = -qn_imag;

        // Multiply: (a+bi)(c+di) = (ac-bd) + (ad+bc)i
        let new_real = prod_real * factor_real - prod_imag * factor_imag;
        let new_imag = prod_real * factor_imag + prod_imag * factor_real;
        prod_real = new_real;
        prod_imag = new_imag;
    }

    // η(τ) = q^{1/24} * product
    let eta_real = q24_real * prod_real - q24_imag * prod_imag;
    let eta_imag = q24_real * prod_imag + q24_imag * prod_real;

    (eta_real, eta_imag)
}

/// Compute the j-invariant from the Dedekind eta function.
///
/// j(τ) = (θ_2(τ)^8 + θ_3(τ)^8 + θ_4(τ)^8)^3 / (8 η(τ)^8 θ_2(τ)^4 θ_3(τ)^4 θ_4(τ)^4)
///
/// The j-invariant has a q-expansion:
/// j(τ) = q^{-1} + 744 + 196884q + 21493760q^2 + ...
///
/// # Monster Moonshine Connection
/// The coefficients 196884, 21493760, ... are dimensions of
/// irreducible representations of the Monster group.
/// This is the famous Monstrous Moonshine conjecture (proved by Borcherds).
///
/// # Crop Circle Connection
/// Metatron's Cube (13 circles) may encode the 13³ factor in |Monster|.
/// The j-invariant coefficients connect to the Monster via moonshine.
pub fn j_invariant_q_expansion(q: f64, terms: usize) -> f64 {
    if q <= 0.0 || q >= 1.0 {
        return f64::NAN;
    }

    // j(τ) = q^{-1} + 744 + Σ c_n q^n
    // where c_n are Monster dimensions
    let mut sum = q.powi(-1) + 744.0;

    // First few Monster dimensions (from OEIS A014708)
    let coeffs = [
        196884.0,
        21493760.0,
        864299970.0,
        20245856256.0,
        333202640600.0,
    ];

    for (i, &c) in coeffs.iter().take(terms).enumerate() {
        sum += c * q.powi((i + 1) as i32);
    }

    sum
}

/// Check if n divides the order of the Monster group.
///
/// The Monster group has order:
/// |M| = 2^46 · 3^20 · 5^9 · 7^6 · 11^2 · 13^3 · 17 · 19 · 23 · 29 · 31 · 41 · 47 · 59 · 71
///
/// # Crop Circle Connection
/// Metatron's Cube has 13 circles. 13³ divides |Monster|.
/// This function checks if a given number is a factor of |Monster|.
pub fn divides_monster_order(n: u64) -> bool {
    n != 0 && {
        // Prime factorization of |Monster|
        let monster_factors: [(u64, u32); 15] = [
            (2, 46),
            (3, 20),
            (5, 9),
            (7, 6),
            (11, 2),
            (13, 3),
            (17, 1),
            (19, 1),
            (23, 1),
            (29, 1),
            (31, 1),
            (41, 1),
            (47, 1),
            (59, 1),
            (71, 1),
        ];

        let mut remaining = n;
        for (p, max_exp) in &monster_factors {
            let mut exp = 0;
            while remaining.is_multiple_of(*p) && exp < *max_exp {
                remaining /= p;
                exp += 1;
            }
        }

        remaining == 1
    }
}

/// Compute the Monster group order.
///
/// Returns the exact order as a string (too large for u64).
pub fn monster_order_string() -> &'static str {
    "808017424794512875886459904961710757005754368000000000"
}

/// Check if a number has a 13-fold symmetry connection to the Monster.
///
/// Returns true if n is divisible by 13 and the quotient also has
/// 13-fold symmetry (suggesting 13³ = 2197 structure).
pub fn has_13_cube_structure(n: u64) -> bool {
    if n.is_multiple_of(13) {
        let q = n / 13;
        q.is_multiple_of(13) && (q / 13).is_multiple_of(13)
    } else {
        false
    }
}

/// Enumerate all ways to represent n as a product of 13-fold symmetries.
///
/// The Monster group requires 3 independent 13-fold symmetries (13³).
/// Metatron's Cube shows 1 visible 13-fold symmetry (13 circles).
/// This function finds how many additional 13-fold symmetries are encoded
/// in the divisibility structure of n.
///
/// Returns (visible_count, hidden_count, total_13_power) where:
/// - visible_count: number of 13-factors in n (e.g., 13² → 2)
/// - hidden_count: estimated higher-dimensional 13-symmetries needed for Monster
/// - total_13_power: total 13-exponent in n
///
/// # Example
/// ```
/// use geographdb_core::algorithms::number_theory::analyze_13_symmetries;
/// let (visible, hidden, total) = analyze_13_symmetries(2197); // 13³
/// assert_eq!(visible, 3);
/// assert_eq!(hidden, 0);
/// assert_eq!(total, 3);
/// ```
pub fn analyze_13_symmetries(n: u64) -> (u32, u32, u32) {
    if n == 0 {
        return (0, 0, 0);
    }

    // Count visible 13-factors
    let mut visible_count = 0u32;
    let mut remaining = n;
    while remaining.is_multiple_of(13) {
        remaining /= 13;
        visible_count += 1;
    }

    // Monster requires 13³ = 3 independent symmetries
    let monster_required = 3u32;
    let hidden_count = monster_required.saturating_sub(visible_count);

    (visible_count, hidden_count, visible_count)
}

/// Find the minimal number that has k independent 13-fold symmetries
/// and also encodes other Monster divisibility properties.
///
/// This searches for numbers of the form 13^k × m where m encodes
/// additional Monster-compatible structure (divisible by other Monster primes).
pub fn find_monster_compatible_13_symmetry(k: u32, max_search: u64) -> Vec<u64> {
    let mut results = Vec::new();
    let base = 13u64.pow(k);

    // Monster primes and their exponents
    let monster_primes: [(u64, u32); 14] = [
        (2, 46),
        (3, 20),
        (5, 9),
        (7, 6),
        (11, 2),
        (17, 1),
        (19, 1),
        (23, 1),
        (29, 1),
        (31, 1),
        (41, 1),
        (47, 1),
        (59, 1),
        (71, 1),
    ];

    for m in 1..=max_search / base {
        let n = base * m;
        if n > max_search {
            break;
        }

        // Check if m is divisible by at least 3 other Monster primes
        let mut monster_prime_count = 0;
        for (p, _max_exp) in &monster_primes {
            if m % p == 0 {
                monster_prime_count += 1;
            }
        }

        if monster_prime_count >= 3 {
            results.push(n);
        }
    }

    results
}

/// Compute the 13-adic valuation of n (highest power of 13 dividing n).
pub fn valuation_13(n: u64) -> u32 {
    if n == 0 {
        return 0;
    }
    let mut count = 0;
    let mut remaining = n;
    while remaining.is_multiple_of(13) {
        remaining /= 13;
        count += 1;
    }
    count
}

/// Check if n encodes a "13-symmetric" structure compatible with Metatron's Cube.
///
/// Metatron's Cube has:
/// - 1 center circle
/// - 12 surrounding circles (2 rings of 6, or 3 rings of 4)
/// - Total: 13 circles
///
/// This function checks if n has divisibility properties that suggest
/// a multi-level 13-symmetric structure (center + rings).
pub fn has_metatron_13_structure(n: u64) -> bool {
    if n.is_multiple_of(13) {
        // Level 1: Single 13 (center + 12 surrounding)
        let level1 = n / 13;

        // Level 2: Check if level1 has further 13-structure
        // (e.g., 13×12 = 156 could encode 13 centers each with 12 satellites)
        let has_level2 =
            level1.is_multiple_of(12) || level1.is_multiple_of(6) || level1.is_multiple_of(4);

        // Level 3: Check for Monster-compatible remaining factor
        let remaining = if has_level2 {
            if level1.is_multiple_of(12) {
                level1 / 12
            } else if level1.is_multiple_of(6) {
                level1 / 6
            } else {
                level1 / 4
            }
        } else {
            level1
        };

        // The remaining factor should be small (like the center point)
        has_level2 && remaining <= 100
    } else {
        false
    }
}

// ── Tests ─────────────────────────────────────────────────────────────────────

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

    #[test]
    fn test_is_prime() {
        assert!(!is_prime(0));
        assert!(!is_prime(1));
        assert!(is_prime(2));
        assert!(is_prime(3));
        assert!(!is_prime(4));
        assert!(is_prime(409));
        assert!(!is_prime(1679)); // 23 * 73
    }

    #[test]
    fn test_1679_is_not_prime() {
        assert!(!is_prime(1679)); // 23 * 73
    }

    #[test]
    fn test_sieve() {
        let primes = sieve_primes(30);
        assert_eq!(primes, vec![2, 3, 5, 7, 11, 13, 17, 19, 23, 29]);
    }

    #[test]
    fn test_factor_trial() {
        let factors = factor_trial(1679);
        assert_eq!(factors.get(&23), Some(&1));
        assert_eq!(factors.get(&73), Some(&1));
    }

    #[test]
    fn test_semiprime_factors() {
        let result = semiprime_factors(1679);
        assert!(result.is_some());
        let (p, q) = result.unwrap();
        assert_eq!(p * q, 1679);
        assert_eq!(semiprime_factors(100), None); // 2^2 * 5^2
    }

    #[test]
    fn test_factor_semiprime_pi_family() {
        // Arecibo Reply: 23 * 73 = 1679
        let result = factor_semiprime_pi_family(1679, 10);
        assert!(result.is_some());
        let (p, q) = result.unwrap();
        assert_eq!(p * q, 1679);

        // Other test cases
        assert!(factor_semiprime_pi_family(154, 10).is_some()); // 7 * 22
        assert!(factor_semiprime_pi_family(385, 10).is_some()); // 11 * 35
    }

    #[test]
    fn test_zeta() {
        let z2 = zeta_real(2.0, 100000);
        assert!((z2 - std::f64::consts::PI * std::f64::consts::PI / 6.0).abs() < 0.01);
    }

    #[test]
    fn test_continued_fraction() {
        let conv = continued_fraction(std::f64::consts::PI, 10);
        let (p, q) = conv[3]; // 355/113 approximation
        assert_eq!(p, 355);
        assert_eq!(q, 113);
    }

    #[test]
    fn test_mod_pow() {
        assert_eq!(mod_pow(2, 10, 1000), 24);
        assert_eq!(mod_pow(3, 5, 7), 5);
    }

    #[test]
    fn test_pollard_rho() {
        // Arecibo Reply
        let factor = pollard_rho(1679, 1);
        assert!(factor.is_some());
        let f = factor.unwrap();
        assert_eq!(1679 % f, 0);
        assert!(f > 1 && f < 1679);

        // Large semiprime
        let n = 10403; // 101 * 103
        let factor = pollard_rho(n, 1);
        assert!(factor.is_some());
        let f = factor.unwrap();
        assert_eq!(n % f, 0);
    }

    #[test]
    fn test_factor_advanced() {
        // Arecibo Reply
        let factors = factor_advanced(1679);
        assert_eq!(factors.get(&23), Some(&1));
        assert_eq!(factors.get(&73), Some(&1));

        // Large semiprime
        let factors = factor_advanced(10403);
        assert_eq!(factors.get(&101), Some(&1));
        assert_eq!(factors.get(&103), Some(&1));
    }

    #[test]
    fn test_gcd() {
        assert_eq!(gcd(48, 18), 6);
        assert_eq!(gcd(17, 13), 1);
        assert_eq!(gcd(100, 25), 25);
    }

    #[test]
    fn test_miller_rabin() {
        // Known primes
        assert!(is_prime_miller_rabin(2));
        assert!(is_prime_miller_rabin(3));
        assert!(is_prime_miller_rabin(409));
        assert!(is_prime_miller_rabin(2803)); // 409th prime

        // Known composites
        assert!(!is_prime_miller_rabin(0));
        assert!(!is_prime_miller_rabin(1));
        assert!(!is_prime_miller_rabin(4));
        assert!(!is_prime_miller_rabin(1679)); // 23 * 73
        assert!(!is_prime_miller_rabin(100));

        // Carmichael number (hard case for trial division)
        assert!(!is_prime_miller_rabin(561)); // 3 * 11 * 17
        assert!(!is_prime_miller_rabin(1105)); // 5 * 13 * 17
        assert!(!is_prime_miller_rabin(1729)); // 7 * 13 * 19

        // Large prime
        assert!(is_prime_miller_rabin(104729)); // 10000th prime
    }

    #[test]
    fn test_mod_inverse() {
        assert_eq!(mod_inverse(3, 11), Some(4)); // 3 * 4 = 12 ≡ 1 mod 11
        assert_eq!(mod_inverse(2, 4), None);
    }

    #[test]
    fn test_riemann_siegel() {
        // Test theta function at t = 14.13 (near first zero)
        let theta = riemann_siegel_theta(14.13);
        assert!(theta.is_finite());

        // Test Z function (simplified, no correction terms)
        let z = riemann_siegel_z(14.13, 3);
        assert!(z.is_finite());
    }

    #[test]
    fn test_monster_moonshine() {
        // 13³ = 2197 divides Monster order
        assert!(divides_monster_order(2197));
        assert!(has_13_cube_structure(2197));

        // 13² does not have cube structure
        assert!(!has_13_cube_structure(169));

        // Monster order string
        assert_eq!(
            monster_order_string(),
            "808017424794512875886459904961710757005754368000000000"
        );

        // j-invariant first coefficient
        let j = j_invariant_q_expansion(0.001, 1);
        assert!(j > 1000.0); // Should be large

        // Dedekind eta
        let (eta_r, eta_i) = dedekind_eta(0.0, 1.0, 10);
        assert!(eta_r.is_finite());
        assert!(eta_i.is_finite());

        // Test 13-symmetry analysis
        let (visible, hidden, total) = analyze_13_symmetries(2197);
        assert_eq!(visible, 3);
        assert_eq!(hidden, 0);
        assert_eq!(total, 3);

        let (visible, hidden, total) = analyze_13_symmetries(13);
        assert_eq!(visible, 1);
        assert_eq!(hidden, 2); // Need 2 more for Monster
        assert_eq!(total, 1);

        // Test 13-adic valuation
        assert_eq!(valuation_13(2197), 3);
        assert_eq!(valuation_13(13), 1);
        assert_eq!(valuation_13(1), 0);

        // Test Metatron structure detection
        assert!(has_metatron_13_structure(156)); // 13 × 12
        assert!(!has_metatron_13_structure(13)); // Just 13, no ring structure

        // Test Monster-compatible search
        let compatible = find_monster_compatible_13_symmetry(1, 10000);
        assert!(!compatible.is_empty());
        // All results should be divisible by 13 and at least 3 other Monster primes
        for &n in &compatible {
            assert_eq!(n % 13, 0);
        }
    }
}