implicit_predicates/

lib.rs

//! Exact geometric predicates over explicit and implicit points.
//!
//! An implicit point is defined by an unevaluated construction, such as a
//! point on an edge or inside a triangle. Predicates receive the construction
//! itself and return the exact sign, without rounding the point's
//! coordinates.
//!
//! **Status: early development.** Today the crate provides one predicate,
//! `orient3d`, over points of any type that implements [`Point`]: an
//! approximation of its coordinates, a bound on the error of each one, and
//! its exact coordinates. The point types of the crate are `[f64; 3]`, an
//! explicit point that is not validated; [`Explicit`], a validated explicit
//! point; [`Lnc`], a point on the segment between two explicit points; and
//! [`AnyPoint`], either of the last two. Like every predicate of the crate,
//! `orient3d` is generated by `implicit-predicates-gen` from a short
//! definition: a floating-point filter decides the common cases, and exact
//! stages built on the expansions of `implicit-predicates-core` decide the
//! rest. The point types are generated from their declarations in the same
//! way. A type of your own can implement [`Point`] too, but the predicates
//! rely on its contract without checking it: an implementation that breaks
//! it makes them return wrong signs. `insphere` and the other kinds of
//! implicit points do not exist yet.
//!
//! ```
//! # #[cfg(feature = "std")] {
//! use implicit_predicates::{Lnc, Sign, orient3d};
//!
//! let (a, b, c) = ([1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]);
//! // a, b and c appear counter-clockwise from (1, 1, 1), above the plane.
//! assert_eq!(orient3d(&a, &b, &c, &[1.0, 1.0, 1.0]), Some(Sign::Negative));
//! assert_eq!(orient3d(&a, &b, &c, &[0.0, 0.0, 0.0]), Some(Sign::Positive));
//! // A point on the plane.
//! assert_eq!(orient3d(&a, &b, &c, &[0.5, 0.25, 0.25]), Some(Sign::Zero));
//! // A point on the edge between b and c, c + 0.3 (b - c): on the plane,
//! // exactly.
//! let m = Lnc::new(b, c, 0.3)?;
//! assert_eq!(orient3d(&a, &b, &c, &m), Some(Sign::Zero));
//! # }
//! # Ok::<(), implicit_predicates::PointError>(())
//! ```
//!
//! A point of type [`Lnc`] is `q + t (p - q)`; its error bound is zero on
//! the coordinates that are exact, here those that `p` and `q` share:
//!
//! ```
//! use implicit_predicates::{Lnc, Point};
//!
//! let m = Lnc::new([0.3, 2.0, 5.0], [0.1, 2.0, 5.0], 0.5)?;
//! assert!(m.err()[0] > 0.0);
//! assert_eq!(&m.approx()[1..], &[2.0, 5.0]);
//! assert_eq!(&m.err()[1..], &[0.0, 0.0]);
//! # Ok::<(), implicit_predicates::PointError>(())
//! ```
//!
//! # Features
//!
//! - `std`, on by default: `orient3d`, which keeps the exact coordinates
//!   of implicit points in an arena of each thread. Without it the crate is
//!   `no_std`, and [`orient3d_in`], which takes the arena, is the way to
//!   call the predicate.
#![cfg_attr(not(feature = "std"), no_std)]
#![forbid(unsafe_code)]

#[cfg(feature = "std")]
mod arena;
mod generated;

pub use generated::points::{AnyPoint, Explicit, Lnc, PointError};
pub use implicit_predicates_core::{
    ExactCoords, ExactError, Point, Scratch, ScratchExhausted, Sign, Weight,
};

use generated::orient3d as orient3d_stages;

/// The orientation of `d` with respect to the plane through `a`, `b` and
/// `c`: the exact sign of the determinant
///
/// ```text
/// | a.x - d.x   a.y - d.y   a.z - d.z |
/// | b.x - d.x   b.y - d.y   b.z - d.z |
/// | c.x - d.x   c.y - d.y   c.z - d.z |
/// ```
///
/// at the exact coordinates of the points, which may be of any type that
/// implements [`Point`]: `[f64; 3]`, [`Explicit`], [`Lnc`], [`AnyPoint`],
/// or a type of your own.
///
/// The sign is [`Sign::Positive`] when `d` lies below the plane, taking
/// "above" to be the side from which `a`, `b` and `c` appear in
/// counter-clockwise order; [`Sign::Negative`] when `d` lies above it; and
/// [`Sign::Zero`] when the four points are coplanar. This is the convention
/// of Shewchuk's `orient3d`. Swapping any two points flips the sign.
///
/// A floating-point filter on the approximations of the points decides most
/// calls. When it cannot, explicit points go through exact stages whose
/// temporaries are on the stack; with an implicit point, the exact
/// coordinates of the points are written to an arena of the calling thread,
/// which starts at 16,384 `f64` and grows on the heap when a call needs
/// more. [`orient3d_in`] takes the arena instead, and needs no `std`.
///
/// Requires the `std` feature, which is on by default.
///
/// # When the result is `None`
///
/// The result is `None` whenever an input of a point is not finite. For
/// explicit points it is `Some` whenever every nonzero coordinate has a
/// magnitude between `2^-188` and `2^239` (about `2.5e-57` and `8.8e71`):
/// the last two stages evaluate with expansions, one on coordinates up to
/// `2^238` and the other on coordinates up to `2^239`, and together they
/// decide every such input. With an implicit point it is `Some` whenever
/// every point can compute its exact coordinates (for an [`Lnc`], its
/// inputs lie within the range of its exact arithmetic), with a weight of
/// [`Weight::One`], and every nonzero component of those coordinates has a
/// magnitude between `2^-188` and `2^235`, with at most 16 components per
/// coordinate. A point of your own type whose exact coordinates have
/// another weight gives `None`: the predicates do not support weights yet.
/// Outside those ranges the sign is still exact when it is returned, but
/// `None` is returned when no stage can decide. This is provisional: a
/// later version will decide every finite input.
///
/// An [`Lnc`] computes each exact coordinate `q_k + t (p_k - q_k)` as the
/// expansion of `p_k - q_k`, scaled by `t`, plus `q_k`: at most five
/// components. It fails, and the result is `None`, when a coordinate of `p`
/// or `q`, or a component of `t (p_k - q_k)`, reaches `2^1020` in magnitude,
/// or when `t` or a nonzero component of `p_k - q_k` has a magnitude below
/// `2^-480` or above `2^480`. Its components can be far smaller than its
/// inputs, though, so a condition on the inputs is more useful: the result
/// is `Some` whenever every point is
///
/// - explicit, with every nonzero coordinate of a magnitude between
///   `2^-188` and `2^235`, or
/// - an [`Lnc`] built with [`Lnc::new`] for which, for some integer `k`
///   with `0 <= k <= 317`, `t` is zero or at least `2^-k` and every
///   nonzero coordinate of `p` and `q` has a magnitude between `2^(k - 84)`
///   and `2^233`: for example, coordinates between `2^-30` and `2^233` and
///   `t` zero or at least `2^-54`. So `t` is always zero or at least
///   `2^-317`, even when every coordinate of `p` and `q` is zero.
///
/// For such an [`Lnc`], every coordinate of `p` and `q` is a multiple of
/// `2^(k - 136)` and `t` one of `2^(-k - 52)`, so every value its
/// computation handles is a multiple of `2^-188`: so are the exact sums and
/// products of such values, their roundings and their rounding errors. And,
/// as `0 <= t <= 1`, the components add up in magnitude to at most
/// `|t| (|p_k| + |q_k|) + |q_k|`, at most `3 * 2^233`, up to a relative
/// error far below `1/3`, so each is below `2^235`. And `t`, being zero or
/// between `2^-317` and 1, is in the product domain. The condition is
/// sufficient, not necessary. For `k >= 1` its lower ends are reached:
/// with `p_k = 2^(k - 84) (1 + 2^-52)` below the last bit of `q_k` and
/// `t = 2^-k (1 + 2^-52)`, a component is exactly `2^-188`, and with half
/// that `p_k` or half that `t` one is below `2^-188`, so a call that only
/// an exact stage can decide gives `None`. Three inputs outside the
/// condition that give `None`:
///
/// ```
/// use implicit_predicates::{Lnc, orient3d};
///
/// let (a, c) = ([1.0, 0.25, 0.5], [0.3, 0.7, 5.0]);
/// // `t = 0.1` is at least `2^-4`, but `q.x` is below `2^-80`.
/// let b = [2f64.powi(-140), 3.0, 0.75];
/// assert_eq!(orient3d(&a, &b, &c, &Lnc::new(a, b, 0.1)?), None);
/// // The coordinates are at least `2^-4`, so `k` can be up to 80, but `t`
/// // is below `2^-80`.
/// let b = [0.1, 3.0, 0.75];
/// assert_eq!(orient3d(&a, &b, &c, &Lnc::new(a, b, 2f64.powi(-200))?), None);
/// // Every coordinate is zero, but `t` is below `2^-317`, and below the
/// // product domain too.
/// let (a, b, c) = ([1.0, -1.0, 0.0], [0.0, 1.0, -1.0], [1.0, 0.0, -1.0]);
/// let m = Lnc::new([0.0; 3], [0.0; 3], 2f64.powi(-600))?;
/// assert_eq!(orient3d(&a, &b, &c, &m), None);
/// # Ok::<(), implicit_predicates::PointError>(())
/// ```
///
/// # Examples
///
/// ```
/// use implicit_predicates::{Lnc, Sign, orient3d};
///
/// let (a, b, c) = ([1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]);
/// // a, b and c appear counter-clockwise from (1, 1, 1), above the plane.
/// assert_eq!(orient3d(&a, &b, &c, &[1.0, 1.0, 1.0]), Some(Sign::Negative));
/// assert_eq!(orient3d(&a, &b, &c, &[0.0, 0.0, 0.0]), Some(Sign::Positive));
/// // A point on the plane, and one that is not a point.
/// assert_eq!(orient3d(&a, &b, &c, &[0.5, 0.25, 0.25]), Some(Sign::Zero));
/// assert_eq!(orient3d(&a, &b, &c, &[f64::NAN, 0.0, 0.0]), None);
///
/// // A point on the edge between a and b, b + 0.1 (a - b), is on the
/// // plane, exactly, although its coordinates are not representable.
/// let m = Lnc::new(a, b, 0.1)?;
/// assert_eq!(orient3d(&a, &b, &c, &m), Some(Sign::Zero));
/// assert_eq!(orient3d(&m, &b, &c, &[1.0, 1.0, 1.0]), Some(Sign::Negative));
/// # Ok::<(), implicit_predicates::PointError>(())
/// ```
#[cfg(feature = "std")]
#[inline]
pub fn orient3d(a: &impl Point, b: &impl Point, c: &impl Point, d: &impl Point) -> Option<Sign> {
    if let Some(sign) = orient3d_stages::filter_points(a, b, c, d) {
        return Some(sign);
    }
    if a.is_explicit() & b.is_explicit() & c.is_explicit() & d.is_explicit() {
        return orient3d_stages::exact_stages([a.approx(), b.approx(), c.approx(), d.approx()]);
    }
    orient3d_in_thread_arena(a, b, c, d)
}

/// The exact stages of [`orient3d`] for points that are not all explicit,
/// in the arena of the thread.
#[cfg(feature = "std")]
#[cold]
#[inline(never)]
fn orient3d_in_thread_arena(
    a: &impl Point,
    b: &impl Point,
    c: &impl Point,
    d: &impl Point,
) -> Option<Sign> {
    // The arena grows as needed: it is exhausted only past `u32::MAX`
    // components, far more than any call of this predicate takes.
    arena::with(|scratch| orient3d_stages::exact_stages_in(scratch, a, b, c, d)).unwrap_or(None)
}

/// `orient3d` with the arena `scratch` for the exact coordinates of
/// implicit points, instead of the arena of the thread: the same result,
/// or an error if `scratch` runs out of room. Available without `std`.
///
/// `scratch` is left as it was, on success and on error. It is not used
/// when the filter decides, nor when every point is explicit. With implicit
/// points the exact stages may need several thousand `f64`; an arena made
/// with `Scratch::growable` (with the `alloc` feature of
/// `implicit-predicates-core`) never runs out.
///
/// # Errors
///
/// [`ScratchExhausted`] if `scratch` runs out of room; a larger arena then
/// gives the result.
///
/// # Examples
///
/// ```
/// use implicit_predicates::{Lnc, Scratch, Sign, orient3d_in};
///
/// let (a, b, c) = ([1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]);
/// let m = Lnc::new(a, b, 0.1)?;
/// let mut buf = vec![0.0; 4096];
/// let mut scratch = Scratch::borrowed(&mut buf);
/// assert_eq!(orient3d_in(&mut scratch, &a, &b, &c, &m), Ok(Some(Sign::Zero)));
/// assert_eq!(scratch.used(), 0);
///
/// // An arena too small is reported.
/// let mut small = [0.0; 8];
/// assert!(orient3d_in(&mut Scratch::borrowed(&mut small), &a, &b, &c, &m).is_err());
/// # Ok::<(), implicit_predicates::PointError>(())
/// ```
#[inline]
pub fn orient3d_in(
    scratch: &mut Scratch<'_>,
    a: &impl Point,
    b: &impl Point,
    c: &impl Point,
    d: &impl Point,
) -> Result<Option<Sign>, ScratchExhausted> {
    orient3d_stages::staged_in(scratch, a, b, c, d)
}

/// The individual stages of each predicate, for tests and benchmarks.
#[doc(hidden)]
pub mod stages {
    pub use crate::generated::orient3d;
}