tilezz 0.2.0

//! Float-lattice toolbox: the shared plane-lattice primitives behind the
//! periodic detectors and the cert verifier -- norms and cross products,
//! basis picking/inversion/reduction, exact lattice membership, and the
//! (orientation, position mod L) coset partition of a placement set.
//!
//! The float discipline throughout is PROPOSAL-SIDE ONLY (the pipeline's
//! "float proposes, exact check disposes" rule): floats rank candidates, size
//! blocks, and read off near-integer coordinates, but every acceptance-relevant
//! decision reduces to an exact ring-element comparison ([`in_lattice`] rounds
//! then tests exact equality; [`gauss_reduce`] applies only unimodular integer
//! ops, so soundness never depends on its float guidance). A wrong float value
//! yields a missed candidate, never a wrong accept.

use rustc_hash::{FxHashMap as HashMap, FxHashSet as HashSet};

use crate::cyclotomic::IsRing;
use crate::cyclotomic::geometry::cmp_xy;
use crate::cyclotomic::geometry::float::{cross_f, norm2_f};
use crate::geom::iso::Iso;

/// Cross-product magnitude below which two float vectors are treated as
/// parallel (and a float determinant as degenerate). Coordinates here are
/// exact ring elements evaluated to f64, so genuine non-parallelism is far
/// above this; only true degeneracy lands below it.
pub(crate) const PARALLEL_EPS: f64 = 1e-9;

/// First pair of linearly independent vectors in `vs`: the first nonzero
/// vector, then the first later one whose cross product against it is nonzero
/// (float test, 1e-9). `None` if the sequence does not span the plane.
/// Order-sensitive on purpose -- a caller that wants the SHORTEST basis sorts
/// by norm first. The single implementation behind every "pick two lattice
/// generators" site (cert verify, Conway isometries, orbit lattices).
pub(crate) fn independent_pair<T: IsRing>(vs: impl IntoIterator<Item = T>) -> Option<(T, T)> {
    let mut first: Option<T> = None;
    for v in vs {
        if v.xy() == (0.0, 0.0) {
            continue;
        }
        match first {
            None => first = Some(v),
            Some(a) => {
                if cross_f(&a, &v).abs() > PARALLEL_EPS {
                    return Some((a, v));
                }
            }
        }
    }
    None
}

/// The shortest independent translation-lattice basis read off a placement
/// orbit: take the pure translations (`rot == 0`), drop the (near-)zero one,
/// sort shortest-first, and pick the first independent pair via
/// [`independent_pair`]. `None` if the orbit's translations do not span the
/// plane.
///
/// The single home for "orbit -> lattice basis", shared by the isohedral
/// detector and [the cert verifier](crate::classify::cert). Shortest-first is
/// what makes [`independent_pair`] most likely to return an actual basis rather
/// than a sublattice pair; and the choice of basis among bases of the same
/// lattice does not affect the cert's exact multiplicity gates (covolume is
/// basis-invariant).
pub(crate) fn lattice_from_orbit<T: IsRing>(placements: &[Iso<T>]) -> Option<(T, T)> {
    let mut trans: Vec<(f64, T)> = placements
        .iter()
        .filter(|p| p.rot == 0)
        .map(|p| (norm2_f(&p.shift), p.shift))
        .filter(|(d, _)| *d > PARALLEL_EPS)
        .collect();
    trans.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap());
    independent_pair(trans.into_iter().map(|(_, v)| v))
}

/// Float inverse of the basis matrix `[v1 v2]` (columns), used to read off
/// integer lattice coordinates of a ring vector. Returns `None` if degenerate.
pub(crate) fn basis_inverse<T: IsRing>(v1: &T, v2: &T) -> Option<[[f64; 2]; 2]> {
    let (v1x, v1y) = v1.xy();
    let (v2x, v2y) = v2.xy();
    let det = v1x * v2y - v1y * v2x;
    if det.abs() < PARALLEL_EPS {
        return None;
    }
    Some([[v2y / det, -v2x / det], [-v1y / det, v1x / det]])
}

/// Whether the ring vector `d` lies in the lattice `Z v1 + Z v2`. The float
/// inverse gives the (near-integer) coordinates; rounding then an **exact**
/// ring-equality test decides it (so float noise can never admit a wrong `d`).
pub(crate) fn in_lattice<T: IsRing>(d: T, v1: T, v2: T, inv: &[[f64; 2]; 2]) -> bool {
    let (dx, dy) = d.xy();
    let m1 = (inv[0][0] * dx + inv[0][1] * dy).round() as i64;
    let m2 = (inv[1][0] * dx + inv[1][1] * dy).round() as i64;
    d == v1.scale(m1) + v2.scale(m2)
}

/// Gauss-reduce a rank-2 basis: repeated exact integer subtractions
/// (`b -= mu * a`), guided by float norm/dot ratios. Unimodular column ops
/// only, so the returned pair spans the SAME lattice -- soundness never
/// depends on the float guidance or on reaching perfect reducedness (the
/// iteration cap only affects how small the gold-check block gets), it only
/// shrinks the vectors so the coverage reaches stay small. `None` if the
/// pair is (near-)parallel.
pub(crate) fn gauss_reduce<T: IsRing>(mut a: T, mut b: T) -> Option<(T, T)> {
    if cross_f(&a, &b).abs() < PARALLEL_EPS {
        return None;
    }
    let dot = |u: &T, w: &T| {
        let (ux, uy) = u.xy();
        let (wx, wy) = w.xy();
        ux * wx + uy * wy
    };
    for _ in 0..64 {
        if norm2_f(&a) > norm2_f(&b) {
            std::mem::swap(&mut a, &mut b);
        }
        let mu = (dot(&a, &b) / norm2_f(&a)).round() as i64;
        if mu == 0 {
            break;
        }
        b = b - a.scale(mu);
    }
    Some((a, b))
}

/// Translation-invariant signature of a placed tile, plus its anchor (the
/// minimal vertex under `cmp_xy`). Two placements have the *same orientation*
/// -- one is a pure translate of the other -- exactly when their signatures are
/// equal; the translation between them is the difference of their anchors.
///
/// Keying on the signature (rather than the `Iso` rotation) is what makes the
/// fold robust to a tile's own rotational symmetry: a symmetric tile can carry
/// several `(rot, shift)` representations for one physical placement, but its
/// vertex *set* is unambiguous.
pub(crate) fn signature<T: IsRing>(tile_verts: &[T]) -> (Vec<T>, T) {
    let anchor = tile_verts
        .iter()
        .copied()
        .min_by(|a, b| cmp_xy(a, b))
        .expect("tile has vertices");
    let mut norm: Vec<T> = tile_verts.iter().map(|&x| x - anchor).collect();
    norm.sort_by(|a, b| cmp_xy(a, b));
    (norm, anchor)
}

/// Group placements by orientation ([`signature`] of the placed tile), each
/// member keyed by its anchor. This is the shared first step of the coset
/// partition, split out so a caller testing MANY candidate lattices over one
/// placement set (the torus cover search) builds it once.
pub(crate) fn signature_groups<T: IsRing>(
    placements: &[Iso<T>],
    verts: &[T],
) -> HashMap<Vec<T>, Vec<(T, Iso<T>)>> {
    let mut groups: HashMap<Vec<T>, Vec<(T, Iso<T>)>> = HashMap::default();
    for iso in placements {
        let (sig, anchor) = signature(&iso.tile(verts));
        groups.entry(sig).or_default().push((anchor, *iso));
    }
    groups
}

/// Partition a signature-grouped placement set into its `(orientation,
/// position mod L)` classes: within each orientation group, anchors that
/// differ by a lattice vector share a class (classes never merge across
/// orientations). Returns one `(member count, representative)` per class; the
/// representative is the class's FIRST member in the group's insertion order,
/// so for an orbit/patch listed near-to-far it sits near the origin. The
/// single coset-partition implementation behind the torus fold
/// (`domain_reps`) and the cert verifier (`coset_reps`).
pub(crate) fn lattice_classes<T: IsRing>(
    groups: &HashMap<Vec<T>, Vec<(T, Iso<T>)>>,
    v1: T,
    v2: T,
    inv: &[[f64; 2]; 2],
) -> Vec<(usize, Iso<T>)> {
    let mut classes: Vec<(usize, Iso<T>)> = Vec::new();
    for members in groups.values() {
        // (anchor, count, representative iso) per class within this group.
        let mut group: Vec<(T, usize, Iso<T>)> = Vec::new();
        for &(m, iso) in members {
            if let Some(slot) = group
                .iter_mut()
                .find(|(a, _, _)| in_lattice(m - *a, v1, v2, inv))
            {
                slot.1 += 1;
            } else {
                group.push((m, 1, iso));
            }
        }
        classes.extend(group.into_iter().map(|(_, c, iso)| (c, iso)));
    }
    classes
}

/// Lay `domain` across the block of lattice translates `i*v1 + j*v2` with
/// `|i| <= reach1`, `|j| <= reach2` (deduplicated, first occurrence kept).
/// The shared instantiation step behind the gold check's proven-coverage
/// block and the torus minter's laid cover.
pub(crate) fn lay_lattice_block<T: IsRing>(
    domain: &[Iso<T>],
    v1: T,
    v2: T,
    reach1: i64,
    reach2: i64,
) -> Vec<Iso<T>> {
    let mut seen: HashSet<Iso<T>> = HashSet::default();
    let mut orbit: Vec<Iso<T>> = Vec::new();
    for d in domain {
        for i in -reach1..=reach1 {
            for j in -reach2..=reach2 {
                let iso = Iso {
                    rot: d.rot,
                    shift: d.shift + v1.scale(i) + v2.scale(j),
                };
                if seen.insert(iso) {
                    orbit.push(iso);
                }
            }
        }
    }
    orbit
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::cyclotomic::ZZ12;
    use crate::cyclotomic::traits::{SymNum, Units};

    /// gauss_reduce spans the SAME lattice (covolume preserved up to sign),
    /// shrinks a skewed basis, and rejects a (near-)parallel pair.
    #[test]
    fn gauss_reduce_shrinks_and_preserves_lattice() {
        let u0 = ZZ12::unit(0);
        let u3 = ZZ12::unit(3);
        // Skewed: v2 = 1000*v1 + u3 reduces back to (u0, u3).
        let (a, b) = (u0, u0.scale(1000) + u3);
        let covol_before = cross_f(&a, &b).abs();
        let (ra, rb) = gauss_reduce(a, b).expect("independent pair reduces");
        assert!(
            (cross_f(&ra, &rb).abs() - covol_before).abs() < 1e-9,
            "same covolume"
        );
        assert!(
            norm2_f(&ra).max(norm2_f(&rb)) < 2.0,
            "skew removed: both vectors short"
        );
        // Parallel pair: no basis.
        assert!(gauss_reduce(u0, u0.scale(7)).is_none());
    }
}