delaunay 0.8.3

//! Finite binary64 and exact rational matrix operations.
//!
//! This module is Delaunay's boundary around the stack-allocated linear algebra
//! functionality provided by `la-stack`.  Geometry code should depend on the
//! local [`Matrix`] alias, checked access helpers, determinant filters, and
//! error conversions here rather than reaching into `la-stack` internals
//! directly.
//!
//! Keeping that shim in one file preserves a narrow API boundary: `la-stack`
//! can evolve its dispatch macros, exact-arithmetic fallbacks, tolerance names,
//! and diagnostic variants while the rest of Delaunay keeps speaking in
//! geometry-level concepts such as predicate matrices, checked active blocks,
//! and public construction errors.
//!
//! Delaunay assembles geometric coefficients and interprets predicate signs;
//! `la-stack` owns interval arithmetic, certified vector reductions, and exact
//! matrix elimination. Already-derived rational inputs
//! stay exact through `RationalMatrix` and `RationalVector`, with explicit
//! rounding only at a public coordinate boundary.
//! Provisional intersection bases use backend LU, Gram construction, and
//! checked dot products; their tolerance and exact certificate checks remain
//! geometry policy. Rectangular rank selection, rational scalar expressions,
//! and symbolic perturbation polynomials remain local because the backend's
//! square-matrix and binary64 vector APIs do not represent those operations.
//!
//! Fixed-size distances use [`la_stack::Vector::norm`], which can return a finite
//! length even when its square overflows. Rounded squared quantities use
//! [`la_stack::Vector::norm_squared`]; predicate-specific exact lifting still forms squares
//! from rational source coordinates. The former `geometry::util::{hypot,
//! squared_norm}` helpers and their prelude exports have been removed. Downstream
//! callers replacing those imports must declare their own direct `la-stack`
//! dependency. Construct a vector with [`la_stack::Vector::try_new`], then call
//! [`norm`](la_stack::Vector::norm) or [`norm_squared`](la_stack::Vector::norm_squared)
//! and handle the typed backend errors.
//!
//! The former `determinant(&matrix)` helper has also been removed. Use
//! [`Matrix::det`] for a floating-point estimate, preserving numerical
//! singularity errors, or [`Matrix::det_sign_exact`] when an exact sign or
//! exact singularity classification is required.

#![forbid(unsafe_code)]

/// Typed errors from the stack-allocated linear algebra backend.
///
/// Delaunay re-exports this backend error type at the matrix boundary so public
/// geometry errors can preserve exact failure context without
/// exposing the rest of `la-stack` to downstream callers.
pub use la_stack::LaError;
use la_stack::Matrix as LaMatrix;
pub(crate) use la_stack::{
    BigRational, DEFAULT_SINGULAR_TOL, ExactF64Conversion, FromPrimitive, Interval,
    IntervalDeterminantSign, IntervalMatrix, RationalMatrix, RationalVector, Signed,
    SingularityReason, Tolerance, Vector as LaVector, gram_matrix,
};
use thiserror::Error;

/// Stack-matrix dispatch limit.
///
/// This is chosen so that common predicate matrices can be built as:
/// - orientation: (D+1)×(D+1)
/// - relative-coordinate insphere: (D+1)×(D+1)
///
/// With `MAX_STACK_MATRIX_DIM = 7`, the relative predicate supports through
/// `D = 6`. Absolute lifted matrices still require `D+2` rows and therefore
/// have a lower stack limit.
pub const MAX_STACK_MATRIX_DIM: usize = la_stack::MAX_STACK_MATRIX_DISPATCH_DIM;

/// Stack-allocated matrix type used by this crate for fixed-size linear algebra.
///
/// This alias is Delaunay's public matrix boundary around `la-stack`: callers
/// can build small matrices for diagnostics and helper APIs without depending
/// on backend module paths.
pub type Matrix<const D: usize> = LaMatrix<D>;

/// Error type for matrix operations.
///
/// # Examples
///
/// ```rust
/// use delaunay::prelude::geometry::MatrixError;
///
/// let err = MatrixError::SingularMatrix;
/// std::assert_matches!(err, MatrixError::SingularMatrix);
/// ```
#[derive(Clone, Debug, Error, PartialEq, Eq)]
#[non_exhaustive]
pub enum MatrixError {
    /// Matrix is singular.
    #[error("Matrix is singular!")]
    SingularMatrix,
    /// Matrix row or column index is outside the concrete stack matrix.
    #[error("matrix index out of bounds: ({row}, {column}) for {dimension}x{dimension}")]
    OutOfBounds {
        /// Requested row index.
        row: usize,
        /// Requested column index.
        column: usize,
        /// Concrete matrix dimension.
        dimension: usize,
    },
}

/// Error type for stack-matrix dispatch and active-block access.
#[derive(Clone, Debug, Error, PartialEq)]
#[non_exhaustive]
pub(crate) enum StackMatrixDispatchError {
    /// The requested matrix size is not supported by the stack-matrix dispatcher.
    #[error("unsupported stack matrix size: {k} (max {max})")]
    UnsupportedDim {
        /// Requested matrix dimension.
        k: usize,
        /// Maximum supported matrix dimension.
        max: usize,
    },
    /// The requested active block size does not match the concrete matrix type.
    #[error("active matrix block size {k} does not match concrete matrix dimension {dim}")]
    ActiveBlockDimensionMismatch {
        /// Requested active matrix dimension.
        k: usize,
        /// Concrete matrix dimension.
        dim: usize,
    },
    /// A linear algebra error originating from `la-stack`.
    #[error(transparent)]
    La {
        /// Typed source error from the linear algebra backend.
        source: LaError,
    },
    /// A matrix access failed inside a dispatched stack-matrix operation.
    #[error(transparent)]
    Matrix {
        /// Typed source error from matrix operations.
        #[from]
        source: MatrixError,
    },
}

impl From<LaError> for StackMatrixDispatchError {
    fn from(source: LaError) -> Self {
        match source {
            LaError::UnsupportedDimension { requested, max, .. } => {
                Self::UnsupportedDim { k: requested, max }
            }
            LaError::IndexOutOfBounds { row, col, dim, .. } => Self::Matrix {
                source: MatrixError::OutOfBounds {
                    row,
                    column: col,
                    dimension: dim,
                },
            },
            source => Self::La { source },
        }
    }
}

/// Dispatch a runtime matrix dimension to a stack matrix, returning an error if unsupported.
///
/// Unsupported upstream dispatch dimensions are converted from [`LaError`], so callers
/// may return [`StackMatrixDispatchError`] directly or a public error type that implements
/// `From<LaError>` and `From<StackMatrixDispatchError>`.
macro_rules! try_with_la_stack_matrix {
    ($k:expr, |$m:ident| $body:block) => {{
        la_stack::try_with_stack_matrix!($k, |mut $m| -> _ $body)
    }};
}

/// Write one finite entry into a stack matrix, preserving backend diagnostics.
///
/// This wrapper is the boundary where predicate matrix construction rejects
/// non-finite values and out-of-bounds indices before later determinant stages
/// can accidentally classify invalid matrix state as geometric degeneracy.
#[inline]
pub(crate) fn matrix_set<const D: usize>(
    m: &mut Matrix<D>,
    row: usize,
    column: usize,
    value: f64,
) -> Result<(), StackMatrixDispatchError> {
    m.set(row, column, value).map_err(Into::into)
}

/// Solve a runtime-sized finite `f64` system with exact fraction-free elimination.
///
/// The backend converts every IEEE 754 input to its exact rational value before
/// applying Bareiss elimination, so this is an exact solve rather than a
/// floating-point approximation.
pub(crate) fn solve_exact_runtime_system(
    matrix: &[Vec<f64>],
    rhs: &[f64],
) -> Option<Result<Vec<BigRational>, StackMatrixDispatchError>> {
    let dimension = rhs.len();
    if matrix.len() != dimension || matrix.iter().any(|row| row.len() != dimension) {
        return None;
    }

    Some(try_with_la_stack_matrix!(dimension, |stack_matrix| {
        for (row, values) in matrix.iter().enumerate() {
            for (column, value) in values.iter().copied().enumerate() {
                matrix_set(&mut stack_matrix, row, column, value)?;
            }
        }
        let rhs_vector = LaVector::try_new(std::array::from_fn(|index| rhs[index]))?;
        stack_matrix
            .solve_exact(rhs_vector)
            .map(|solution| solution.into_array().into_iter().collect())
            .map_err(Into::into)
    }))
}

/// Converts one finite IEEE-754 value to the exact rational it represents.
#[inline]
pub(crate) fn rational_from_f64(value: f64) -> Option<BigRational> {
    value
        .is_finite()
        .then(|| BigRational::from_f64(value))
        .flatten()
}

/// Return a determinant and its certified error bound when the f64 fast filter supports the matrix size.
///
/// `Ok(None)` means the closed-form direct determinant path is unavailable or
/// inconclusive for this matrix size, including arithmetic overflow or
/// underflow-sensitive evaluation. Callers should continue to exact arithmetic;
/// `la-stack` matrices are finite by construction.
#[inline]
pub(crate) fn matrix_fast_filter<const D: usize>(
    m: &Matrix<D>,
) -> Result<Option<(f64, f64)>, StackMatrixDispatchError> {
    match m.det_direct_with_errbound() {
        Ok(Some(estimate)) => Ok(Some((
            estimate.determinant(),
            estimate.absolute_error_bound(),
        ))),
        Ok(None) | Err(LaError::NonFinite { .. }) => Ok(None),
        Err(source) => Err(source.into()),
    }
}

#[cfg(test)]
pub(crate) mod test_support {
    use super::{Matrix, StackMatrixDispatchError};

    /// Keeps checked matrix reads available to predicate diagnostics in unit tests.
    pub fn matrix_get<const D: usize>(
        matrix: &Matrix<D>,
        row: usize,
        column: usize,
    ) -> Result<f64, StackMatrixDispatchError> {
        matrix.try_get(row, column).map_err(Into::into)
    }

    /// Dispatch a runtime `k` to a stack-allocated matrix for concise unit tests.
    macro_rules! with_la_stack_matrix {
        ($k:expr, |$m:ident| $body:block) => {{
            la_stack::try_with_stack_matrix!($k, |mut $m| -> Result<_, la_stack::LaError> {
                Ok($body)
            })
            .expect("test requested an unsupported stack matrix size")
        }};
    }

    pub(crate) use with_la_stack_matrix;
}

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

    #[test]
    fn try_with_la_stack_matrix_returns_err_on_unsupported_dim() {
        let k = MAX_STACK_MATRIX_DIM + 1;
        let res: Result<(), StackMatrixDispatchError> =
            try_with_la_stack_matrix!(k, |_m| { Ok(()) });
        assert_matches!(
            res,
            Err(StackMatrixDispatchError::UnsupportedDim {
                k: requested,
                max
            }) if requested == k && max == MAX_STACK_MATRIX_DIM
        );
    }

    #[test]
    fn solve_exact_runtime_system_rejects_malformed_shapes() {
        assert_eq!(
            solve_exact_runtime_system(&[vec![1.0, 0.0]], &[1.0, 0.0]),
            None
        );
        assert_eq!(
            solve_exact_runtime_system(&[vec![1.0], vec![0.0, 1.0]], &[1.0, 0.0]),
            None
        );
    }

    #[test]
    fn la_index_error_maps_to_matrix_error_with_context() {
        let err = StackMatrixDispatchError::from(LaError::index_out_of_bounds(3, 4, 2));

        assert_eq!(
            err,
            StackMatrixDispatchError::Matrix {
                source: MatrixError::OutOfBounds {
                    row: 3,
                    column: 4,
                    dimension: 2,
                },
            }
        );
    }

    #[test]
    fn stack_matrix_dispatch_error_clones_and_displays_la_error_source() {
        let source = LaError::singular_exact(3);
        let error = StackMatrixDispatchError::La { source };

        assert_eq!(error.clone(), error);
        assert_eq!(error.to_string(), source.to_string());
    }

    #[test]
    fn matrix_set_returns_error_on_out_of_bounds_index() {
        let mut matrix = Matrix::<2>::zero();
        let err = matrix_set(&mut matrix, 0, 2, 1.0).unwrap_err();
        assert_eq!(
            err,
            StackMatrixDispatchError::Matrix {
                source: MatrixError::OutOfBounds {
                    row: 0,
                    column: 2,
                    dimension: 2,
                },
            }
        );
    }
}