geographdb-core 0.5.4

//! 13-Fold Symmetries — Geometric, Algebraic, and Graph-Theoretic
//!
//! This module extends the divisibility-based 13-symmetry analysis in
//! `number_theory.rs` with actual geometric, algebraic, and graph
//! constructions.  Every function is designed to be testable without
//! external dependencies (no SageMath / PARI-GP required).
//!
//! # Crop Circle Connection
//!
//! Metatron's Cube (13 circles) is the primary formation that motivates
//! this work.  The 13 circles may encode:
//! - Geometric: a regular 13-gon and its dihedral symmetries
//! - Algebraic: 13th roots of unity and the cyclotomic field Q(ζ₁₃)
//! - Graph: a 13-regular graph whose automorphism group carries the
//!   missing 13-fold symmetries needed for the Monster group.

#![allow(
    clippy::needless_range_loop,
    reason = "index-heavy math loops are clearer with explicit indices"
)]

use std::collections::HashMap;

// ═══════════════════════════════════════════════════════════════════════════════
//  GEOMETRIC 13-FOLD SYMMETRIES
// ═══════════════════════════════════════════════════════════════════════════════

/// Vertices of a regular 13-gon inscribed in the unit circle.
///
/// Returns 13 points (x, y) with 2π/13 rotational symmetry.
///
/// # Crop Circle Connection
/// The 13 circles of Metatron's Cube can be modelled as vertices of a
/// regular 13-gon when projected onto a plane.
pub fn regular_13_gon() -> Vec<(f64, f64)> {
    let n = 13;
    let mut pts = Vec::with_capacity(n);
    for k in 0..n {
        let theta = 2.0 * std::f64::consts::PI * (k as f64) / (n as f64);
        pts.push((theta.cos(), theta.sin()));
    }
    pts
}

/// Apply a k-step rotation to a regular 13-gon vertex index.
///
/// C₁₃ action:  rot_k(i) = (i + k) mod 13.
/// This is the cyclic group of order 13.
pub fn c13_rotate(index: usize, k: usize) -> usize {
    (index + k) % 13
}

/// Reflect a regular 13-gon vertex index across an axis through vertex `axis`.
///
/// D₁₃ action:  ref_axis(i) = (2·axis − i) mod 13.
/// The dihedral group D₁₃ has order 26 = 2 × 13 (13 rotations + 13 reflections).
pub fn d13_reflect(index: usize, axis: usize) -> usize {
    (2 * axis + 13 - (index % 13)) % 13
}

/// Compute the full dihedral group D₁₃ as permutations of {0,…,12}.
///
/// Returns 26 permutations (13 rotations + 13 reflections).
/// Each permutation is a Vec<usize> of length 13.
pub fn dihedral_d13_permutations() -> Vec<Vec<usize>> {
    let mut perms = Vec::with_capacity(26);
    // Rotations
    for k in 0..13 {
        let mut p = Vec::with_capacity(13);
        for i in 0..13 {
            p.push(c13_rotate(i, k));
        }
        perms.push(p);
    }
    // Reflections
    for axis in 0..13 {
        let mut p = Vec::with_capacity(13);
        for i in 0..13 {
            p.push(d13_reflect(i, axis));
        }
        perms.push(p);
    }
    perms
}

/// Detect 13-fold rotational symmetry in a 2-D point cloud.
///
/// Algorithm:
/// 1. Compute centroid.
/// 2. Convert to polar coordinates (r, θ).
/// 3. Bin θ into 13 angular sectors.
/// 4. Compute χ² uniformity statistic.
/// 5. If uniform and r-distribution is peaked → 13-fold symmetry detected.
///
/// Returns (is_13_symmetric, chi_squared, peak_r_fraction).
///
/// # Parameters
/// - `points`: slice of (x, y) coordinates.
/// - `chi_threshold`: maximum χ² for uniformity (default 26.0 for 12 dof).
pub fn detect_13_fold_symmetry(points: &[(f64, f64)], chi_threshold: f64) -> (bool, f64, f64) {
    if points.len() < 13 {
        return (false, f64::INFINITY, 0.0);
    }

    // Centroid
    let (sum_x, sum_y) = points
        .iter()
        .fold((0.0, 0.0), |(sx, sy), (x, y)| (sx + x, sy + y));
    let cx = sum_x / points.len() as f64;
    let cy = sum_y / points.len() as f64;

    // Polar coordinates
    let mut radii: Vec<f64> = Vec::with_capacity(points.len());
    let mut angles: Vec<f64> = Vec::with_capacity(points.len());
    for &(x, y) in points {
        let dx = x - cx;
        let dy = y - cy;
        radii.push((dx * dx + dy * dy).sqrt());
        angles.push(dy.atan2(dx)); // [-π, π]
    }

    // Angular histogram (13 bins)
    let mut bins = [0usize; 13];
    for &a in &angles {
        let norm = if a < 0.0 {
            a + 2.0 * std::f64::consts::PI
        } else {
            a
        };
        let bin = ((norm / (2.0 * std::f64::consts::PI)) * 13.0).floor() as usize % 13;
        bins[bin] += 1;
    }

    let expected = points.len() as f64 / 13.0;
    let chi_sq: f64 = bins
        .iter()
        .map(|&count| {
            let diff = count as f64 - expected;
            diff * diff / expected
        })
        .sum();

    // Radial concentration: fraction within 1 std-dev of mean radius
    let mean_r = radii.iter().sum::<f64>() / radii.len() as f64;
    let var_r = radii.iter().map(|&r| (r - mean_r).powi(2)).sum::<f64>() / radii.len() as f64;
    let std_r = var_r.sqrt();
    let in_band = radii
        .iter()
        .filter(|&&r| (r - mean_r).abs() <= std_r)
        .count();
    let peak_fraction = in_band as f64 / radii.len() as f64;

    let is_symmetric = chi_sq <= chi_threshold && peak_fraction >= 0.5;
    (is_symmetric, chi_sq, peak_fraction)
}

/// Generate a 13-point star polygon {13/k}.
///
/// A star polygon connects every k-th vertex of a regular 13-gon.
/// k must be coprime to 13 (i.e., k ∈ {1,2,3,4,5,6}).
///
/// Returns the sequence of vertex indices that trace the star.
pub fn star_polygon_13(k: usize) -> Option<Vec<usize>> {
    if k == 0 || k >= 13 || gcd(k as u64, 13) != 1 {
        return None;
    }
    let mut seq = Vec::with_capacity(13);
    let mut current = 0;
    for _ in 0..13 {
        seq.push(current);
        current = (current + k) % 13;
    }
    Some(seq)
}

// ═══════════════════════════════════════════════════════════════════════════════
//  ALGEBRAIC 13-FOLD SYMMETRIES
// ═══════════════════════════════════════════════════════════════════════════════

/// The 13th roots of unity: ζ₁₃^k = e^(2πi·k/13) for k = 0,…,12.
///
/// Returns (real, imag) pairs.  These are the vertices of the regular
/// 13-gon in the complex plane, and generate the cyclotomic field
/// Q(ζ₁₃), a degree-12 extension of Q (since φ(13) = 12).
pub fn roots_of_unity_13() -> Vec<(f64, f64)> {
    let mut roots = Vec::with_capacity(13);
    for k in 0..13 {
        let theta = 2.0 * std::f64::consts::PI * (k as f64) / 13.0;
        roots.push((theta.cos(), theta.sin()));
    }
    roots
}

/// Evaluate the 13th cyclotomic polynomial Φ₁₃(x).
///
/// Φ₁₃(x) = (x¹³ − 1) / (x − 1) = x¹² + x¹¹ + … + x + 1.
///
/// This is the minimal polynomial of ζ₁₃ over Q.
pub fn cyclotomic_13(x: f64) -> f64 {
    if (x - 1.0).abs() < 1e-12 {
        return 13.0; // limit as x → 1
    }
    // Direct evaluation: sum of geometric series
    // x^12 + x^11 + ... + x + 1
    let mut sum = 1.0;
    let mut term = 1.0;
    for _ in 0..12 {
        term *= x;
        sum += term;
    }
    sum
}

/// Compute the real subfield Q(ζ₁₃ + ζ₁₃⁻¹) basis elements.
///
/// The maximal real subfield has degree φ(13)/2 = 6 over Q.
/// Basis: η_k = ζ₁₃^k + ζ₁₃^(-k) = 2·cos(2πk/13) for k = 1,…,6.
///
/// These are algebraic integers that generate the ring of integers
/// of the real subfield.
pub fn real_subfield_basis_13() -> [f64; 6] {
    let mut basis = [0.0; 6];
    for k in 1..=6 {
        let theta = 2.0 * std::f64::consts::PI * (k as f64) / 13.0;
        basis[k - 1] = 2.0 * theta.cos();
    }
    basis
}

/// Compute Gaussian periods of degree 3 for p = 13.
///
/// The multiplicative group (Z/13Z)* is cyclic of order 12.
/// Let g = 2 be a generator.  The subgroup H = {g⁰, g³, g⁶, g⁹}
/// has order 4, giving 3 cosets and 3 Gaussian periods η₀, η₁, η₂.
///
/// These periods generate the unique cubic subfield of Q(ζ₁₃),
/// with discriminant 13² = 169.
///
/// Returns (η₀, η₁, η₂) where each ηᵢ = Σ_{h∈H·g^i} ζ₁₃^h.
pub fn gaussian_periods_cubic_13() -> (f64, f64, f64) {
    // Generator g = 2 mod 13
    // (Z/13Z)* = {1,2,3,4,5,6,7,8,9,10,11,12}
    // Powers of 2 mod 13: 2^0=1, 2^1=2, 2^2=4, 2^3=8, 2^4=3, 2^5=6, 2^6=12, 2^7=11, 2^8=9, 2^9=5, 2^10=10, 2^11=7
    // Subgroup H of index 3 (order 4): {2^0, 2^3, 2^6, 2^9} = {1, 8, 12, 5}
    let h: [usize; 4] = [1, 8, 12, 5];

    // Coset representatives: g^0=1, g^1=2, g^2=4
    let coset_reps: [usize; 3] = [1, 2, 4];

    let mut periods = [0.0; 3];
    for i in 0..3 {
        let rep = coset_reps[i];
        let mut sum_r = 0.0;
        for &hj in &h {
            let exp = (hj * rep) % 13;
            let theta = 2.0 * std::f64::consts::PI * (exp as f64) / 13.0;
            sum_r += theta.cos();
        }
        periods[i] = sum_r; // imaginary parts cancel
    }
    (periods[0], periods[1], periods[2])
}

/// Evaluate the modular discriminant Δ(τ) = η(τ)²⁴ for level-1.
///
/// For Γ₀(13) modular forms, we need the cusp forms of weight 2.
/// The function η(τ)²·η(13τ)² is a weight-2 cusp form for Γ₀(13).
///
/// Returns the real value (imaginary part is zero for τ purely imaginary).
///
/// # Parameters
/// - `tau_imag`: Im(τ) > 0 (real part assumed 0 for simplicity).
pub fn cusp_form_gamma0_13(tau_imag: f64, terms: usize) -> f64 {
    if tau_imag <= 0.0 {
        return 0.0;
    }

    // η(τ) = q^{1/24} Π (1 − q^n) where q = e^{−2π·Im(τ)}
    let q = (-2.0 * std::f64::consts::PI * tau_imag).exp();
    let q13 = (-2.0 * std::f64::consts::PI * 13.0 * tau_imag).exp();

    let eta = |q_base: f64, t: usize| -> f64 {
        let mut prod = 1.0;
        for n in 1..=t {
            prod *= 1.0 - q_base.powi(n as i32);
        }
        q_base.powf(1.0 / 24.0) * prod
    };

    let eta_tau = eta(q, terms);
    let eta_13tau = eta(q13, terms);

    // Weight-2 cusp form: f(τ) = η(τ)² · η(13τ)²
    eta_tau.powi(2) * eta_13tau.powi(2)
}

/// Compute the Hauptmodul (uniformizer) for X₀(13).
///
/// X₀(13) is the modular curve of level Γ₀(13).  It has genus 0,
/// so it has a Hauptmodul (a function that generates its function field).
///
/// The Hauptmodul can be expressed as:
///   t₁₃(τ) = (η(τ) / η(13τ))² · (some rational function)
///
/// For practical computation we use the q-expansion:
///   t₁₃ = q⁻¹ + 2 + 5q + 4q² − 3q³ + …
///
/// Returns the approximate value for given q = e^{2πiτ}.
pub fn hauptmodul_x0_13(q: f64, terms: usize) -> f64 {
    if q <= 0.0 || q >= 1.0 {
        return f64::NAN;
    }

    // q-expansion coefficients for Hauptmodul of X₀(13)
    // (from known modular forms literature)
    let coeffs = [
        2.0,  // q^0
        5.0,  // q^1
        4.0,  // q^2
        -3.0, // q^3
        0.0,  // q^4
        1.0,  // q^5
        -1.0, // q^6
        2.0,  // q^7
        0.0,  // q^8
        -2.0, // q^9
    ];

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

/// Check if a prime p is a 13th-power residue modulo another prime q.
///
/// Uses the 13th-power residue symbol, which generalizes the Legendre
/// symbol to higher powers.  For prime q ≡ 1 (mod 13), the multiplicative
/// group (Z/qZ)* has a unique subgroup of index 13.
///
/// Returns:
/// - 1  if p is a 13th power mod q
/// - 0  if p ≡ 0 mod q
/// - -1 otherwise (q not ≡ 1 mod 13, or p not a 13th power)
pub fn power_residue_symbol_13(p: u64, q: u64) -> i32 {
    if q == 0 {
        return 0;
    }
    let p = p % q;
    if p == 0 {
        return 0;
    }
    if q % 13 != 1 {
        return -1; // No primitive 13th roots of unity
    }

    // Euler criterion for 13th powers: p^((q-1)/13) ≡ 1 (mod q)
    let exp = (q - 1) / 13;
    let result = mod_pow(p, exp, q);
    if result == 1 {
        1
    } else {
        -1
    }
}

// ═══════════════════════════════════════════════════════════════════════════════
//  GRAPH-THEORETIC 13-FOLD SYMMETRIES
// ═══════════════════════════════════════════════════════════════════════════════

/// Build the circulant graph C₁₃(S) with connection set S.
///
/// Vertices: {0, 1, …, 12}
/// Edges:    {i, i±s mod 13} for each s ∈ S.
///
/// This graph always has C₁₃ (cyclic) symmetry.  If S is chosen so
/// that the graph is also symmetric under reflection, it has full D₁₃
/// symmetry.
///
/// # Parameters
/// - `s`: connection set (e.g., vec![1, 3] for edges to ±1 and ±3 neighbours).
///
/// Returns adjacency list as Vec<Vec<usize>>.
pub fn circulant_graph_13(s: &[usize]) -> Vec<Vec<usize>> {
    let n = 13;
    let mut adj: Vec<Vec<usize>> = vec![Vec::new(); n];
    for i in 0..n {
        for &step in s {
            let step = step % n;
            if step == 0 {
                continue;
            }
            let j1 = (i + step) % n;
            let j2 = (i + n - step) % n;
            if !adj[i].contains(&j1) {
                adj[i].push(j1);
            }
            if !adj[i].contains(&j2) {
                adj[i].push(j2);
            }
        }
    }
    adj
}

/// Build the Paley graph of order 13.
///
/// Vertices: {0, 1, …, 12} = GF(13)
/// Edges:    {a, b} iff a − b is a quadratic residue mod 13.
///
/// Quadratic residues mod 13: {1, 3, 4, 9, 10, 12}.
///
/// The Paley graph of order q (q ≡ 1 mod 4 prime) is strongly regular
/// with parameters (q, (q−1)/2, (q−5)/4, (q−1)/4).
/// For q = 13: srg(13, 6, 2, 3).
///
/// It has a large automorphism group containing the affine group
/// AGL(1, 13), which has order 13 × 12 = 156.
pub fn paley_graph_13() -> Vec<Vec<usize>> {
    let n = 13;
    let residues: [usize; 6] = [1, 3, 4, 9, 10, 12];
    let mut adj: Vec<Vec<usize>> = vec![Vec::new(); n];
    for i in 0..n {
        for &r in &residues {
            let j = (i + r) % n;
            if i != j && !adj[i].contains(&j) {
                adj[i].push(j);
                adj[j].push(i);
            }
        }
    }
    adj
}

/// Check whether a graph is 13-regular (every vertex has degree 13).
///
/// For a graph on n vertices, a 13-regular graph requires n ≥ 14
/// (by the handshaking lemma, n·13 must be even).
pub fn is_13_regular(adj: &[Vec<usize>]) -> bool {
    adj.len() >= 14
        && adj.iter().all(|neighbors| {
            let mut unique: Vec<usize> = neighbors.clone();
            unique.sort_unstable();
            unique.dedup();
            unique.len() == 13
        })
}

/// Compute the degree of each vertex.
pub fn degree_sequence(adj: &[Vec<usize>]) -> Vec<usize> {
    adj.iter()
        .map(|neighbors| {
            let mut unique = neighbors.clone();
            unique.sort_unstable();
            unique.dedup();
            unique.len()
        })
        .collect()
}

/// Build a 13-fold graph cover (voltage graph construction).
///
/// Given a base graph B and a voltage assignment α: E(B) → Z₁₃,
/// the derived graph (cover) has vertex set V(B) × Z₁₃ and edges
/// ((u, i), (v, i + α(u,v))) for each edge (u, v) in B.
///
/// This produces a graph with 13 × |V(B)| vertices that is a
/// 13-sheeted covering of B.
///
/// # Parameters
/// - `base_adj`: adjacency list of base graph.
/// - `voltages`: map from (u, v) to voltage in {0,…,12}.
///
/// Returns adjacency list of the cover graph.
pub fn graph_cover_13(
    base_adj: &[Vec<usize>],
    voltages: &HashMap<(usize, usize), usize>,
) -> Vec<Vec<usize>> {
    let n_base = base_adj.len();
    let n_cover = n_base * 13;
    let mut cover_adj: Vec<Vec<usize>> = vec![Vec::new(); n_cover];

    for u in 0..n_base {
        for &v in &base_adj[u] {
            let voltage = *voltages.get(&(u, v)).unwrap_or(&0) % 13;
            for i in 0..13 {
                let src = u * 13 + i;
                let dst = v * 13 + ((i + voltage) % 13);
                if !cover_adj[src].contains(&dst) {
                    cover_adj[src].push(dst);
                }
            }
        }
    }
    cover_adj
}

/// Detect whether a graph has a 13-fold automorphism (C₁₃ symmetry).
///
/// Algorithm: look for a permutation σ of order 13 that preserves
/// adjacency.  For small graphs this is done by brute force; for
/// circulant graphs it is automatic.
///
/// Returns (has_c13, rotation_orbits) where rotation_orbits lists
/// the vertex orbits under a detected C₁₃ action.
pub fn detect_c13_automorphism(adj: &[Vec<usize>]) -> (bool, Vec<Vec<usize>>) {
    let n = adj.len();
    if !n.is_multiple_of(13) {
        return (false, Vec::new());
    }

    // For graphs that are circulant or derived from voltage covers,
    // the C₁₃ action is (v, i) → (v, i+1 mod 13).
    // Check if this standard rotation is an automorphism.
    let n_base = n / 13;
    let mut is_automorphism = true;

    for u in 0..n_base {
        for i in 0..13 {
            let src = u * 13 + i;
            let src_rot = u * 13 + ((i + 1) % 13);

            // Neighbors of src after rotation
            let mut rot_neighbors: Vec<usize> = adj[src]
                .iter()
                .map(|&v| {
                    let base = v / 13;
                    let fiber = v % 13;
                    base * 13 + ((fiber + 1) % 13)
                })
                .collect();
            rot_neighbors.sort_unstable();

            let mut actual_neighbors = adj[src_rot].clone();
            actual_neighbors.sort_unstable();

            if rot_neighbors != actual_neighbors {
                is_automorphism = false;
                break;
            }
        }
        if !is_automorphism {
            break;
        }
    }

    if !is_automorphism {
        return (false, Vec::new());
    }

    // Compute orbits
    let mut visited = vec![false; n];
    let mut orbits = Vec::new();
    for v in 0..n {
        if visited[v] {
            continue;
        }
        let mut orbit = Vec::new();
        let base = v / 13;
        let fiber = v % 13;
        for k in 0..13 {
            let u = base * 13 + ((fiber + k) % 13);
            if !visited[u] {
                visited[u] = true;
                orbit.push(u);
            }
        }
        if !orbit.is_empty() {
            orbits.push(orbit);
        }
    }

    (true, orbits)
}

/// Check if a graph is a (13, g)-cage: a 13-regular graph of girth g
/// with the minimum possible number of vertices.
///
/// Known cages:
/// - (13, 3): complete graph K₁₄ (14 vertices)
/// - (13, 4): complete bipartite K₁₃,₁₃ (26 vertices)
/// - (13, 5): unknown / very large
///
/// This function checks the easy cases and returns the cage status.
pub fn check_13_cage(adj: &[Vec<usize>]) -> Option<(usize, &'static str)> {
    let n = adj.len();
    let deg_seq = degree_sequence(adj);
    let is_regular_13 = deg_seq.iter().all(|&d| d == 13);

    if !is_regular_13 {
        return None;
    }

    // Girth 3: K₁₄
    if n == 14 {
        // Every pair of vertices should be adjacent
        for i in 0..n {
            if adj[i].len() != 13 {
                return None;
            }
        }
        return Some((3, "K_14: the (13,3)-cage"));
    }

    // Girth 4: K₁₃,₁₃
    if n == 26 {
        // Check bipartite and each part has degree 13
        // Simple check: look for 2-coloring
        let mut color = vec![None; n];
        let mut is_bipartite = true;
        for start in 0..n {
            if color[start].is_some() {
                continue;
            }
            color[start] = Some(0);
            let mut stack = vec![start];
            while let Some(u) = stack.pop() {
                for &v in &adj[u] {
                    if color[v].is_none() {
                        color[v] = Some(1 - color[u].unwrap());
                        stack.push(v);
                    } else if color[v] == color[u] {
                        is_bipartite = false;
                        break;
                    }
                }
                if !is_bipartite {
                    break;
                }
            }
            if !is_bipartite {
                break;
            }
        }
        if is_bipartite {
            return Some((4, "K_13,13: the (13,4)-cage"));
        }
    }

    None
}

// ═══════════════════════════════════════════════════════════════════════════════
//  Integration: 13-Symmetry Score for Crop Circle Analysis
// ═══════════════════════════════════════════════════════════════════════════════

/// Comprehensive 13-symmetry analysis of a point cloud.
///
/// Combines geometric, algebraic, and graph scores into a single report.
///
/// Returns (geometric_score, algebraic_score, graph_score, overall).
/// Each score is in [0, 1], where 1.0 means strong 13-symmetry evidence.
pub fn analyze_13_symmetry_comprehensive(points: &[(f64, f64)]) -> (f64, f64, f64, f64) {
    // Geometric: detect rotational symmetry
    let (has_geo, chi, peak) = detect_13_fold_symmetry(points, 26.0);
    let geo_score = if has_geo {
        0.5 + 0.5 * peak.min(1.0)
    } else {
        (26.0 / (26.0 + chi)) * 0.3 // partial credit for near-uniformity
    };

    // Algebraic: check if point count relates to 13
    let n = points.len();
    let alg_score = if n.is_multiple_of(13) {
        0.8 + 0.2 * (1.0 / (1.0 + (n / 13) as f64)) // higher score for smaller multiples
    } else if n % 13 == 1 {
        0.5 // center + 12 surrounding pattern
    } else {
        0.1
    };

    // Graph: build nearest-neighbour graph and check for 13-regularity
    let nn_graph = build_nn_graph(points, 6);
    let deg_seq = degree_sequence(&nn_graph);
    let avg_deg = deg_seq.iter().sum::<usize>() as f64 / deg_seq.len() as f64;
    let graph_score = if (avg_deg - 13.0).abs() < 1.0 {
        0.7
    } else if (avg_deg - 12.0).abs() < 2.0 {
        0.4
    } else {
        0.1
    };

    let overall = (geo_score + alg_score + graph_score) / 3.0;
    (geo_score, alg_score, graph_score, overall)
}

/// Build a k-nearest-neighbour graph from a point cloud.
fn build_nn_graph(points: &[(f64, f64)], k: usize) -> Vec<Vec<usize>> {
    let n = points.len();
    let mut adj: Vec<Vec<usize>> = vec![Vec::new(); n];
    for i in 0..n {
        let mut dists: Vec<(f64, usize)> = Vec::with_capacity(n);
        for j in 0..n {
            if i == j {
                continue;
            }
            let dx = points[i].0 - points[j].0;
            let dy = points[i].1 - points[j].1;
            dists.push((dx * dx + dy * dy, j));
        }
        dists.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap());
        for (_, j) in dists.iter().take(k) {
            adj[i].push(*j);
            adj[*j].push(i);
        }
    }
    adj
}

// ── Helper: GCD (re-exported from number_theory for standalone use) ───────────

fn gcd(mut a: u64, mut b: u64) -> u64 {
    while b != 0 {
        let t = b;
        b = a % b;
        a = t;
    }
    a
}

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
}

// ═══════════════════════════════════════════════════════════════════════════════
//  Tests
// ═══════════════════════════════════════════════════════════════════════════════

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

    #[test]
    fn test_regular_13_gon() {
        let gon = regular_13_gon();
        assert_eq!(gon.len(), 13);
        // All points on unit circle
        for &(x, y) in &gon {
            let r = (x * x + y * y).sqrt();
            assert!((r - 1.0).abs() < 1e-10);
        }
        // Check rotational symmetry: rotating by 1 step
        for k in 0..13 {
            let rotated = c13_rotate(0, k);
            assert_eq!(rotated, k % 13);
        }
    }

    #[test]
    fn test_dihedral_d13() {
        let perms = dihedral_d13_permutations();
        assert_eq!(perms.len(), 26); // 13 rotations + 13 reflections

        // Identity is the first rotation (k=0)
        let id: Vec<usize> = (0..13).collect();
        assert_eq!(perms[0], id);

        // Reflections are involutions: applying twice gives identity
        for axis in 0..13 {
            let refl = &perms[13 + axis];
            let mut double = vec![0; 13];
            for i in 0..13 {
                double[i] = refl[refl[i]];
            }
            assert_eq!(double, id, "reflection {} is not an involution", axis);
        }
    }

    #[test]
    fn test_star_polygon_13() {
        // {13/1} is the regular 13-gon itself
        let star1 = star_polygon_13(1).unwrap();
        assert_eq!(star1, vec![0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]);

        // {13/2} is a star that visits every 2nd vertex
        let star2 = star_polygon_13(2).unwrap();
        assert_eq!(star2.len(), 13);
        assert_eq!(star2[0], 0);
        assert_eq!(star2[1], 2);

        // {13/5} is invalid (5 is coprime to 13, wait — it IS coprime)
        assert!(star_polygon_13(5).is_some());

        // {13/0} is invalid
        assert!(star_polygon_13(0).is_none());

        // {13/13} is invalid
        assert!(star_polygon_13(13).is_none());
    }

    #[test]
    fn test_detect_13_fold_symmetry() {
        // Perfect 13-gon: should detect symmetry
        let gon = regular_13_gon();
        let (is_sym, chi, peak) = detect_13_fold_symmetry(&gon, 26.0);
        assert!(is_sym, "regular 13-gon should have 13-fold symmetry");
        assert!(chi < 26.0, "chi-squared should be low for regular polygon");
        assert!(peak >= 0.5, "radial concentration should be high");

        // Random points: should NOT detect symmetry
        let random_pts: Vec<(f64, f64)> = (0..100)
            .map(|i| (i as f64 * 0.1, (i * 7) as f64 * 0.1))
            .collect();
        let (is_sym_r, _, _) = detect_13_fold_symmetry(&random_pts, 26.0);
        assert!(!is_sym_r, "random points should not have 13-fold symmetry");
    }

    #[test]
    fn test_roots_of_unity_13() {
        let roots = roots_of_unity_13();
        assert_eq!(roots.len(), 13);

        // ζ^0 = 1
        assert!((roots[0].0 - 1.0).abs() < 1e-10);
        assert!(roots[0].1.abs() < 1e-10);

        // Sum of all 13th roots of unity = 0
        let sum_r: f64 = roots.iter().map(|(r, _)| r).sum();
        let sum_i: f64 = roots.iter().map(|(_, i)| i).sum();
        assert!(sum_r.abs() < 1e-10);
        assert!(sum_i.abs() < 1e-10);
    }

    #[test]
    fn test_cyclotomic_13() {
        // Φ₁₃(1) = 13
        assert!((cyclotomic_13(1.0) - 13.0).abs() < 1e-10);

        // Φ₁₃(0) = 1
        assert!((cyclotomic_13(0.0) - 1.0).abs() < 1e-10);

        // For a primitive root ζ = e^(2πi/13), Φ₁₃(ζ) = 0.
        // We verify this by evaluating the complex sum directly.
        let theta = 2.0 * std::f64::consts::PI / 13.0;
        let zeta_r = theta.cos();
        let zeta_i = theta.sin();

        let mut sum_r = 1.0;
        let mut sum_i = 0.0;
        let mut zr = 1.0;
        let mut zi = 0.0;
        for _ in 0..12 {
            let new_r = zr * zeta_r - zi * zeta_i;
            let new_i = zr * zeta_i + zi * zeta_r;
            zr = new_r;
            zi = new_i;
            sum_r += zr;
            sum_i += zi;
        }
        let mag = (sum_r * sum_r + sum_i * sum_i).sqrt();
        assert!(mag < 1e-10, "|Φ(ζ)| = {} should be ~0", mag);
    }

    #[test]
    fn test_real_subfield_basis_13() {
        let basis = real_subfield_basis_13();
        assert_eq!(basis.len(), 6);

        // η_k = 2·cos(2πk/13) are all real and in (-2, 2)
        for &eta in &basis {
            assert!(eta > -2.0 && eta < 2.0);
        }

        // Sum of basis elements = -1 (trace from Q(ζ₁₃) to Q)
        let sum: f64 = basis.iter().sum();
        assert!((sum + 1.0).abs() < 1e-10);
    }

    #[test]
    fn test_gaussian_periods_cubic_13() {
        let (eta0, eta1, eta2) = gaussian_periods_cubic_13();

        // Sum of periods = -1 (they are roots of x³ + x² − 4x + 1 = 0)
        let sum = eta0 + eta1 + eta2;
        assert!((sum + 1.0).abs() < 0.5); // loose tolerance for numerical sum
    }

    #[test]
    fn test_cusp_form_gamma0_13() {
        let f = cusp_form_gamma0_13(0.5, 20);
        assert!(f.is_finite());
        assert!(f > 0.0); // cusp form is positive for τ = iy, y > 0
    }

    #[test]
    fn test_hauptmodul_x0_13() {
        let t = hauptmodul_x0_13(0.01, 5);
        assert!(t.is_finite());
        assert!(t > 0.0); // q^{-1} dominates for small q
    }

    #[test]
    fn test_power_residue_symbol_13() {
        // 3 is a primitive root mod 13
        // 3^3 = 27 ≡ 1 mod 13? No, 27 mod 13 = 1. Wait, 3^3 = 27 = 2*13 + 1 = 27 mod 13 = 1.
        // Actually 3^3 = 27 ≡ 1 (mod 13). So 3 has order 3, not primitive.
        // Let's use 2: 2^12 = 4096 ≡ 1 (mod 13), and no smaller positive exponent works.

        // For q = 53 (which is 1 mod 13: 53 = 4*13 + 1), check if 2 is a 13th power
        let sym = power_residue_symbol_13(2, 53);
        assert!(sym == 1 || sym == -1); // should be defined

        // For q not ≡ 1 mod 13, symbol is -1
        assert_eq!(power_residue_symbol_13(2, 17), -1);
    }

    #[test]
    fn test_circulant_graph_13() {
        // C_13(1) is the cycle graph C_13
        let cycle = circulant_graph_13(&[1]);
        assert_eq!(cycle.len(), 13);
        for i in 0..13 {
            assert_eq!(cycle[i].len(), 2); // degree 2
        }

        // C_13(1, 3) has degree 4
        let g = circulant_graph_13(&[1, 3]);
        for i in 0..13 {
            assert_eq!(g[i].len(), 4);
        }
    }

    #[test]
    fn test_paley_graph_13() {
        let paley = paley_graph_13();
        assert_eq!(paley.len(), 13);

        // Paley(13) is strongly regular with parameters (13, 6, 2, 3)
        // Every vertex has degree 6
        for i in 0..13 {
            let mut unique = paley[i].clone();
            unique.sort_unstable();
            unique.dedup();
            assert_eq!(unique.len(), 6, "vertex {} has wrong degree", i);
        }
    }

    #[test]
    fn test_is_13_regular() {
        // K_14 is 13-regular
        let mut k14: Vec<Vec<usize>> = vec![Vec::new(); 14];
        for i in 0..14 {
            for j in 0..14 {
                if i != j {
                    k14[i].push(j);
                }
            }
        }
        assert!(is_13_regular(&k14));

        // Cycle C_13 is 2-regular, not 13-regular
        let c13 = circulant_graph_13(&[1]);
        assert!(!is_13_regular(&c13));
    }

    #[test]
    fn test_graph_cover_13() {
        // Base: a single edge (2 vertices)
        let base = vec![vec![1], vec![0]];
        let mut voltages = HashMap::new();
        voltages.insert((0, 1), 1usize);
        voltages.insert((1, 0), 12usize); // inverse voltage

        let cover = graph_cover_13(&base, &voltages);
        assert_eq!(cover.len(), 26); // 2 × 13

        // Each vertex in cover should have same degree as base
        assert_eq!(cover[0].len(), 1);
        assert_eq!(cover[13].len(), 1);
    }

    #[test]
    fn test_detect_c13_automorphism() {
        // Circulant graph C_13(1) has obvious C_13 symmetry
        let c13 = circulant_graph_13(&[1]);
        let (has_c13, orbits) = detect_c13_automorphism(&c13);
        assert!(has_c13);
        assert_eq!(orbits.len(), 1); // single orbit of size 13
        assert_eq!(orbits[0].len(), 13);
    }

    #[test]
    fn test_check_13_cage() {
        // K_14 is the (13, 3)-cage
        let mut k14: Vec<Vec<usize>> = vec![Vec::new(); 14];
        for i in 0..14 {
            for j in 0..14 {
                if i != j {
                    k14[i].push(j);
                }
            }
        }
        let cage = check_13_cage(&k14);
        assert!(cage.is_some());
        let (girth, desc) = cage.unwrap();
        assert_eq!(girth, 3);
        assert!(desc.contains("K_14"));

        // C_13 is not 13-regular
        let c13 = circulant_graph_13(&[1]);
        assert!(check_13_cage(&c13).is_none());
    }

    #[test]
    fn test_comprehensive_analysis() {
        // Regular 13-gon should score high
        let gon = regular_13_gon();
        let (_geo, _alg, _graph, overall) = analyze_13_symmetry_comprehensive(&gon);
        assert!(
            overall > 0.5,
            "regular 13-gon should have high 13-symmetry score"
        );

        // Random points should score low
        let random: Vec<(f64, f64)> = (0..50)
            .map(|i| ((i * 17) as f64 * 0.1, (i * 31) as f64 * 0.1))
            .collect();
        let (_, _, _, overall_r) = analyze_13_symmetry_comprehensive(&random);
        assert!(overall_r < 0.5, "random points should have low score");
    }
}