monge-rs 0.1.0

//! Foundational types for Monge geometry — circles in ℝ², spheres in ℝⁿ,
//! 3D lifting (paraboloid map), and n-dimensional Lassak generalization.
//!
//! ## Overview
//! The central objects in Monge's theorem are **circles** (or n-spheres).
//! For three circles in the plane, the three external homothetic centers
//! are collinear — that's the classic theorem. This module provides the
//! type system to extend that insight to arbitrary dimensions.

use nalgebra::{SVector, Vector2, Vector3, Vector4};

// ---------------------------------------------------------------------------
// Core types
// ---------------------------------------------------------------------------

/// A circle in ℝ² with center `c` and radius `r`.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct Circle {
    pub center: Vector2<f64>,
    pub radius: f64,
}

impl Circle {
    pub fn new(center: Vector2<f64>, radius: f64) -> Self {
        assert!(radius >= 0.0, "radius must be non-negative");
        Self { center, radius }
    }

    /// Power of point `x` w.r.t. this circle: |x - c|² - r².
    pub fn power(&self, x: &Vector2<f64>) -> f64 {
        (x - self.center).norm_squared() - self.radius.powi(2)
    }
}

/// An n-sphere in ℝⁿ with center `c` and radius `r`.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct Sphere<const D: usize> {
    pub center: SVector<f64, D>,
    pub radius: f64,
}

impl<const D: usize> Sphere<D> {
    pub fn new(center: SVector<f64, D>, radius: f64) -> Self {
        assert!(radius >= 0.0, "radius must be non-negative");
        Self { center, radius }
    }

    /// Power of point `x` w.r.t. this sphere: |x - c|² - r².
    pub fn power(&self, x: &SVector<f64, D>) -> f64 {
        (x - self.center).norm_squared() - self.radius.powi(2)
    }
}

// ---------------------------------------------------------------------------
// 3D lifting (paraboloid map)
// ---------------------------------------------------------------------------

/// Lift a circle in ℝ² to a 3D point on the paraboloid z = x² + y² - r².
///
/// The lift map takes (c_x, c_y, r) ↦ (c_x, c_y, c_x² + c_y² - r²).
/// Under this lift, the radical axis of two circles becomes the intersection
/// of their lifted planes, and Monge collinearity becomes coplanarity.
pub fn lift_circle(c: &Circle) -> Vector3<f64> {
    let x = c.center.x;
    let y = c.center.y;
    Vector3::new(x, y, x * x + y * y - c.radius.powi(2))
}

/// Lift an n-sphere in ℝⁿ to an (n+1)-dimensional point on the paraboloid.
///
/// Concrete implementations for n = 2,3,4 as proof of concept.
/// The general case requires const-generic arithmetic (`feature(generic_const_exprs)`).
pub fn lift_sphere_2d(s: &Sphere<2>) -> Vector3<f64> {
    let x = s.center[0];
    let y = s.center[1];
    Vector3::new(x, y, s.center.norm_squared() - s.radius.powi(2))
}

/// Lift a 3-sphere to a 4D point on the paraboloid.
pub fn lift_sphere_3d(s: &Sphere<3>) -> Vector4<f64> {
    let mut v = Vector4::zeros();
    for i in 0..3 {
        v[i] = s.center[i];
    }
    v[3] = s.center.norm_squared() - s.radius.powi(2);
    v
}

// ---------------------------------------------------------------------------
// Radial axis helper
// ---------------------------------------------------------------------------

/// Compute the radical line (1-cocycle) between two circles.
/// Returns (a, b, c) such that the radical axis is a·x + b·y + c = 0.
///
/// The radical axis is the set of points with equal power w.r.t. both circles.
/// Algebraically: 2·(c₂ - c₁)·x = |c₂|² - r₂² - (|c₁|² - r₁²).
pub fn radical_line_coeffs(c1: &Circle, c2: &Circle) -> (f64, f64, f64) {
    let a = 2.0 * (c2.center.x - c1.center.x);
    let b = 2.0 * (c2.center.y - c1.center.y);
    let c = c1.center.norm_squared()
        - c1.radius.powi(2)
        - c2.center.norm_squared()
        + c2.radius.powi(2);
    (a, b, c)
}

// ---------------------------------------------------------------------------
// nD Lassak generalization
// ---------------------------------------------------------------------------

/// For n+1 spheres in ℝⁿ, their external homothetic centers lie in an
/// (n-1)-flat.  This function returns the normal vector of that flat
/// for n = 2 (the line of collinearity) as proof of concept.
///
/// For n = 2 (3 circles in ℝ²): returns the normal of the line containing
/// the three external homothetic centers.
///
/// For the general nD case, this is an (n-1)-flat whose normal is the
/// (n+1) × (n+1) determinant of the lifted points.
/// See: Lassak, "n-dimensional Monge theorem" (1988).
pub fn lassak_normal(circles: &[Circle; 3]) -> Vector2<f64> {
    // Lift all three circles
    let p1 = lift_circle(&circles[0]);
    let p2 = lift_circle(&circles[1]);
    let p3 = lift_circle(&circles[2]);

    // The normal is the cross product of (p2 - p1) and (p3 - p1),
    // which gives the normal of the plane containing the lifted points.
    // The projection onto the xy-plane gives the Monge line direction.
    let v1 = p2 - p1;
    let v2 = p3 - p1;
    let normal = v1.cross(&v2); // 3D cross product

    // The Monge line in ℝ² is the intersection of this plane with z=0,
    // which has normal = (normal.x, normal.y) — the first two components
    Vector2::new(normal.x, normal.y)
}

/// For 4 spheres in ℝ³, their external homothetic centers lie in a 2-flat (plane).
/// Returns the normal vector of that plane.
///
/// Generic nD: for n+1 spheres in ℝⁿ, the (n-1)-flat normal is given
/// by the (n+1)×(n+1) Cayley-Menger-style determinant of the lifted points.
/// Here we compute it concretely for n=3 using the 4D cross product.
pub fn lassak_normal_3d(spheres: &[Sphere<3>; 4]) -> Vector3<f64> {
    let p1 = lift_sphere_3d(&spheres[0]);
    let p2 = lift_sphere_3d(&spheres[1]);
    let p3 = lift_sphere_3d(&spheres[2]);
    let p4 = lift_sphere_3d(&spheres[3]);

    // 4D cross product = determinant of basis vectors with (p2-p1), (p3-p1), (p4-p1)
    // Result is a 4-vector normal to the hyperplane containing the lifted points.
    let v1 = p2 - p1;
    let v2 = p3 - p1;
    let v3 = p4 - p1;

    // Compute 4D cross product
    let normal_4d = Vector4::new(
        det3(
            (v1[1], v1[2], v1[3]),
            (v2[1], v2[2], v2[3]),
            (v3[1], v3[2], v3[3]),
        ),
        -det3(
            (v1[0], v1[2], v1[3]),
            (v2[0], v2[2], v2[3]),
            (v3[0], v3[2], v3[3]),
        ),
        det3(
            (v1[0], v1[1], v1[3]),
            (v2[0], v2[1], v2[3]),
            (v3[0], v3[1], v3[3]),
        ),
        -det3(
            (v1[0], v1[1], v1[2]),
            (v2[0], v2[1], v2[2]),
            (v3[0], v3[1], v3[2]),
        ),
    );

    // The (n-1)-flat in ℝ³ is given by the first 3 components of the normal
    // (the intersection of the hyperplane with the w=0 slice)
    Vector3::new(normal_4d[0], normal_4d[1], normal_4d[2])
}

// ---------------------------------------------------------------------------
// Internal helpers
// ---------------------------------------------------------------------------

fn det3(a: (f64, f64, f64), b: (f64, f64, f64), c: (f64, f64, f64)) -> f64 {
    a.0 * (b.1 * c.2 - b.2 * c.1) - a.1 * (b.0 * c.2 - b.2 * c.0) + a.2 * (b.0 * c.1 - b.1 * c.0)
}

// ---------------------------------------------------------------------------
// Tests
// ---------------------------------------------------------------------------

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

    #[test]
    fn test_circle_power() {
        let c = Circle::new(Vector2::new(0.0, 0.0), 1.0);
        // On the circle: |(1,0)|² - 1 = 0
        assert!((c.power(&Vector2::new(1.0, 0.0))).abs() < 1e-12);
        // Inside: |(0,0)|² - 1 = -1
        assert!((c.power(&Vector2::new(0.0, 0.0)) + 1.0).abs() < 1e-12);
        // Outside: |(2,0)|² - 1 = 3
        assert!((c.power(&Vector2::new(2.0, 0.0)) - 3.0).abs() < 1e-12);
    }

    #[test]
    fn test_lift_identity() {
        // Lifting preserves the center coordinates
        let c = Circle::new(Vector2::new(3.0, 4.0), 2.0);
        let p = lift_circle(&c);
        assert!((p.x - 3.0).abs() < 1e-12);
        assert!((p.y - 4.0).abs() < 1e-12);
        // z = x² + y² - r² = 9 + 16 - 4 = 21
        assert!((p.z - 21.0).abs() < 1e-12);
    }

    #[test]
    fn test_sphere_generic() {
        let s: Sphere<3> = Sphere::new(Vector3::new(1.0, 2.0, 3.0), 5.0);
        assert!((s.center.norm() - (14.0_f64).sqrt()).abs() < 1e-12);
        assert!((s.radius - 5.0).abs() < 1e-12);
    }

    #[test]
    fn test_lassak_normal_is_nonzero() {
        let circles = [
            Circle::new(Vector2::new(0.0, 0.0), 1.0),
            Circle::new(Vector2::new(4.0, 0.0), 2.0),
            Circle::new(Vector2::new(2.0, 3.0), 1.5),
        ];
        let n = lassak_normal(&circles);
        // The normal should be non-zero for non-degenerate circles
        assert!(n.norm() > 1e-10, "Lassak normal should be non-zero");
    }

    #[test]
    fn test_radical_line_coeffs_simple() {
        // Two circles centered on x-axis
        let c1 = Circle::new(Vector2::new(0.0, 0.0), 1.0);
        let c2 = Circle::new(Vector2::new(4.0, 0.0), 2.0);
        let (a, b, c) = radical_line_coeffs(&c1, &c2);
        // a=8, b=0, c = 0-1 - 16+4 = -13
        // Line: 8x - 13 = 0 => x = 13/8 = 1.625
        assert!((a - 8.0).abs() < 1e-12);
        assert!((b).abs() < 1e-12);
        assert!((c + 13.0).abs() < 1e-12);
    }
}