implicit-predicates 0.0.2

// Generated by implicit-predicates-gen from predicates/points.pred. Do not edit.
// To regenerate, run from the directory that the paths are relative to:
// implicit-predicates-gen crates/implicit-predicates/predicates --stages crates/implicit-predicates/src/generated --checks crates/implicit-predicates/tests/generated

//! The point types [`Explicit`], [`Lnc`], and [`AnyPoint`], which holds a point
//! of any of them.

use core::fmt;

use implicit_predicates_core::{
    ExactCoords, ExactError, Mark, Point, Scratch, ScratchExhausted, Span, Weight, eft, expansion,
};

/// Why the constructor of a point type rejected its inputs.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub enum PointError {
    /// An input is NaN or infinite.
    NonFinite,
    /// The parameters do not satisfy the conditions of the type.
    ParamOutOfRange,
}

impl fmt::Display for PointError {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        f.write_str(match self {
            PointError::NonFinite => "an input of the point is not finite",
            PointError::ParamOutOfRange => "a parameter of the point is out of range",
        })
    }
}

impl core::error::Error for PointError {}

/// A point of type `Explicit`: the point with the coordinates `p`, given
/// exactly. Its approximation is its coordinates and its error is zero.
#[derive(Clone, Copy, Debug, PartialEq)]
#[repr(transparent)]
pub struct Explicit([f64; 3]);

impl Explicit {
    /// The point with the coordinates `p`.
    ///
    /// # Errors
    ///
    /// [`PointError::NonFinite`] if a coordinate is NaN or infinite.
    #[inline]
    pub fn new(p: [f64; 3]) -> Result<Self, PointError> {
        if !(p.iter().all(|x| x.is_finite())) {
            return Err(PointError::NonFinite);
        }
        Ok(Self::new_unchecked(p))
    }

    /// The point with the coordinates `p`, without checking them. A coordinate that
    /// is not finite makes [`Point::exact`] fail with [`ExactError::OutOfRange`].
    #[inline]
    pub fn new_unchecked(p: [f64; 3]) -> Self {
        Self(p)
    }
}

impl Point for Explicit {
    #[inline]
    fn approx(&self) -> &[f64; 3] {
        &self.0
    }

    #[inline]
    fn err(&self) -> &[f64; 3] {
        &[0.0; 3]
    }

    #[inline]
    fn is_explicit(&self) -> bool {
        true
    }

    #[inline]
    fn exact(&self, scratch: &mut Scratch<'_>) -> Result<ExactCoords, ExactError> {
        self.0.exact(scratch)
    }
}

/// A point of type `Lnc`, computed from `p`, `q` and `t`:
///
/// ```text
/// x = q.x + t * (p.x - q.x)
/// y = q.y + t * (p.y - q.y)
/// z = q.z + t * (p.z - q.z)
/// where 0 <= t, t <= 1
/// ```
///
/// The constructors compute the coordinates in `f64` with a bound on the error
/// of each one, [`Point::approx`] and [`Point::err`]; [`Point::exact`] computes
/// them exactly.
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct Lnc {
    p: [f64; 3],
    q: [f64; 3],
    t: f64,
    approx: [f64; 3],
    err: [f64; 3],
}

impl Lnc {
    /// The slack factor of the error bound: the smallest `f64` `S` with `S * (1 -
    /// 2^-53)^6 >= 1`, where 6 is the number of roundings of the computation of the
    /// error of a coordinate on its longest chain of operations.
    const ERR_SLACK: f64 = 1.0000000000000009;

    /// The point of type `Lnc` with the parameters `p`, `q` and `t`.
    ///
    /// # Errors
    ///
    /// [`PointError::NonFinite`] if an input is NaN or infinite, and
    /// [`PointError::ParamOutOfRange`] if the parameters do not satisfy the
    /// conditions of the type.
    pub fn new(p: [f64; 3], q: [f64; 3], t: f64) -> Result<Self, PointError> {
        if !(p.iter().all(|x| x.is_finite()) && q.iter().all(|x| x.is_finite()) && t.is_finite()) {
            return Err(PointError::NonFinite);
        }
        if t < 0.0 {
            return Err(PointError::ParamOutOfRange);
        }
        if t > 1.0 {
            return Err(PointError::ParamOutOfRange);
        }
        Ok(Self::new_unchecked(p, q, t))
    }

    /// The point of type `Lnc` with the parameters `p`, `q` and `t`, without
    /// checking them: any `f64` is accepted. A coordinate whose approximation or
    /// error bound is not finite gets an error of `+∞`, and [`Point::exact`]
    /// fails with [`ExactError::OutOfRange`] when an input is not finite.
    pub fn new_unchecked(p: [f64; 3], q: [f64; 3], t: f64) -> Self {
        let (approx, err) = Self::evaluate([p, q], [t]);
        Self {
            p,
            q,
            t,
            approx,
            err,
        }
    }

    /// The parameter `p`.
    #[inline]
    pub fn p(&self) -> &[f64; 3] {
        &self.p
    }

    /// The parameter `q`.
    #[inline]
    pub fn q(&self) -> &[f64; 3] {
        &self.q
    }

    /// The parameter `t`.
    #[inline]
    pub fn t(&self) -> f64 {
        self.t
    }

    /// The coordinates in `f64` and the bounds on their errors, from the `F64x3`
    /// parameters `c` and the `f64` parameters `s`, in the order of the
    /// declaration.
    fn evaluate(c: [[f64; 3]; 2], s: [f64; 1]) -> ([f64; 3], [f64; 3]) {
        let v0 = c[1][0];
        let v1 = c[1][1];
        let v2 = c[1][2];
        let v3 = s[0];
        let v4 = c[0][0];
        let v5 = c[0][1];
        let v6 = c[0][2];
        let v7 = v4 - v0;
        let e7 = sum_rounding(v4, v0, v7);
        let v8 = v5 - v1;
        let e8 = sum_rounding(v5, v1, v8);
        let v9 = v6 - v2;
        let e9 = sum_rounding(v6, v2, v9);
        let v10 = v3 * v7;
        let e10 = covered_product(v3.abs(), e7) + product_rounding(v3, v7, v10);
        let v11 = v3 * v8;
        let e11 = covered_product(v3.abs(), e8) + product_rounding(v3, v8, v11);
        let v12 = v3 * v9;
        let e12 = covered_product(v3.abs(), e9) + product_rounding(v3, v9, v12);
        let v13 = v0 + v10;
        let e13 = e10 + sum_rounding(v0, v10, v13);
        let v14 = v1 + v11;
        let e14 = e11 + sum_rounding(v1, v11, v14);
        let v15 = v2 + v12;
        let e15 = e12 + sum_rounding(v2, v12, v15);
        let approx = [v13, v14, v15];
        let err = [
            finish(approx[0], e13, Self::ERR_SLACK),
            finish(approx[1], e14, Self::ERR_SLACK),
            finish(approx[2], e15, Self::ERR_SLACK),
        ];
        (approx, err)
    }

    /// Writes the exact coordinates in `scratch`, from the `F64x3` parameters `c`
    /// and the `f64` parameters `s`, in the order of the declaration. Fails with
    /// [`ExactError::OutOfRange`] if an input is not finite, or an operand of a sum
    /// of expansions reaches `2^1020`, or an operand of an exact product leaves the
    /// product domain; on an error the arena is as it was.
    fn write_exact(
        scratch: &mut Scratch<'_>,
        c: [[f64; 3]; 2],
        s: [f64; 1],
    ) -> Result<ExactCoords, ExactError> {
        if !(c.iter().flatten().all(|x| x.is_finite()) && s.iter().all(|x| x.is_finite())) {
            return Err(ExactError::OutOfRange);
        }
        let mark = scratch.mark();
        let mut in_range = true;
        let x = scratch.write(5, |out| {
            let v0 = c[1][0];
            let v3 = s[0];
            let v4 = c[0][0];
            if !(summable(&[v4]) && summable(&[v0])) {
                in_range = false;
                return 0;
            }
            let mut x7 = [0.0; 2];
            let n7 = expansion::from_two_diff(v4, v0, &mut x7);
            if !(multipliable(&x7[..n7]) && multipliable(&[v3])) {
                in_range = false;
                return 0;
            }
            let mut x10 = [0.0; 4];
            let n10 = expansion::scale(&x7[..n7], v3, &mut x10);
            if !(summable(&x10[..n10]) && summable(&[v0])) {
                in_range = false;
                return 0;
            }
            expansion::grow(&x10[..n10], v0, out)
        });
        let x = written(scratch, mark, x, in_range)?;
        let y = scratch.write(5, |out| {
            let v1 = c[1][1];
            let v3 = s[0];
            let v5 = c[0][1];
            if !(summable(&[v5]) && summable(&[v1])) {
                in_range = false;
                return 0;
            }
            let mut x8 = [0.0; 2];
            let n8 = expansion::from_two_diff(v5, v1, &mut x8);
            if !(multipliable(&x8[..n8]) && multipliable(&[v3])) {
                in_range = false;
                return 0;
            }
            let mut x11 = [0.0; 4];
            let n11 = expansion::scale(&x8[..n8], v3, &mut x11);
            if !(summable(&x11[..n11]) && summable(&[v1])) {
                in_range = false;
                return 0;
            }
            expansion::grow(&x11[..n11], v1, out)
        });
        let y = written(scratch, mark, y, in_range)?;
        let z = scratch.write(5, |out| {
            let v2 = c[1][2];
            let v3 = s[0];
            let v6 = c[0][2];
            if !(summable(&[v6]) && summable(&[v2])) {
                in_range = false;
                return 0;
            }
            let mut x9 = [0.0; 2];
            let n9 = expansion::from_two_diff(v6, v2, &mut x9);
            if !(multipliable(&x9[..n9]) && multipliable(&[v3])) {
                in_range = false;
                return 0;
            }
            let mut x12 = [0.0; 4];
            let n12 = expansion::scale(&x9[..n9], v3, &mut x12);
            if !(summable(&x12[..n12]) && summable(&[v2])) {
                in_range = false;
                return 0;
            }
            expansion::grow(&x12[..n12], v2, out)
        });
        let z = written(scratch, mark, z, in_range)?;
        Ok(ExactCoords {
            x,
            y,
            z,
            w: Weight::One,
        })
    }
}

impl Point for Lnc {
    #[inline]
    fn approx(&self) -> &[f64; 3] {
        &self.approx
    }

    #[inline]
    fn err(&self) -> &[f64; 3] {
        &self.err
    }

    #[inline]
    fn is_explicit(&self) -> bool {
        false
    }

    fn exact(&self, scratch: &mut Scratch<'_>) -> Result<ExactCoords, ExactError> {
        Self::write_exact(scratch, [self.p, self.q], [self.t])
    }
}

/// A point of any of the types of this module.
#[derive(Clone, Copy, Debug, PartialEq)]
pub enum AnyPoint {
    /// A point of type [`Explicit`].
    Explicit(Explicit),
    /// A point of type [`Lnc`].
    Lnc(Lnc),
}

impl Point for AnyPoint {
    #[inline]
    fn approx(&self) -> &[f64; 3] {
        match self {
            AnyPoint::Explicit(p) => p.approx(),
            AnyPoint::Lnc(p) => p.approx(),
        }
    }

    #[inline]
    fn err(&self) -> &[f64; 3] {
        match self {
            AnyPoint::Explicit(p) => p.err(),
            AnyPoint::Lnc(p) => p.err(),
        }
    }

    #[inline]
    fn is_explicit(&self) -> bool {
        match self {
            AnyPoint::Explicit(p) => p.is_explicit(),
            AnyPoint::Lnc(p) => p.is_explicit(),
        }
    }

    #[inline]
    fn exact(&self, scratch: &mut Scratch<'_>) -> Result<ExactCoords, ExactError> {
        match self {
            AnyPoint::Explicit(p) => p.exact(scratch),
            AnyPoint::Lnc(p) => p.exact(scratch),
        }
    }
}

impl From<Explicit> for AnyPoint {
    #[inline]
    fn from(p: Explicit) -> Self {
        AnyPoint::Explicit(p)
    }
}

impl From<Lnc> for AnyPoint {
    #[inline]
    fn from(p: Lnc) -> Self {
        AnyPoint::Lnc(p)
    }
}

/// `2^-53`, the unit roundoff of `f64`.
const UNIT: f64 = 1.1102230246251565e-16;

/// `2^-1074`, the smallest positive `f64`.
const TINY: f64 = 5e-324;

/// `fl(fl(a * b) + 2^-1074)` for non-negative `a` and `b`, or zero if either is
/// zero: at least `(1 - 2^-53)^2 * a * b`, also when `a * b` underflows, since
/// adding `2^-1074` to a value below `2^-1022` is exact.
#[inline]
fn covered_product(a: f64, b: f64) -> f64 {
    if a == 0.0 || b == 0.0 {
        0.0
    } else {
        a * b + TINY
    }
}

/// The bound on the rounding error of the sum or difference `s` of `x` and `y`:
/// the covered product of `2^-53` and `|s|`, or zero if an operand is zero,
/// when the operation is exact.
#[inline]
fn sum_rounding(x: f64, y: f64, s: f64) -> f64 {
    if x != 0.0 && y != 0.0 {
        covered_product(UNIT, s.abs())
    } else {
        0.0
    }
}

/// The bound on the rounding error of the product `p` of `x` and `y`: the
/// covered product of `2^-53` and `|p|`, plus `2^-1074` for underflow, or zero
/// if an operand is zero, when the operation is exact. The test is on the
/// operands: a product that underflows is zero without being exact.
#[inline]
fn product_rounding(x: f64, y: f64, p: f64) -> f64 {
    if x != 0.0 && y != 0.0 {
        covered_product(UNIT, p.abs()) + TINY
    } else {
        0.0
    }
}

/// The error bound stored for a coordinate whose value is `approx` and whose
/// computed error is `err`: `err * slack` rounded up, which is at least the
/// error since the slack factor covers the roundings of its computation; zero
/// when `err` is zero; and `+∞` when `approx` or the bound is not finite.
#[inline]
fn finish(approx: f64, err: f64, slack: f64) -> f64 {
    let bound = err * slack;
    if !(approx.is_finite() && bound.is_finite()) {
        f64::INFINITY
    } else if err == 0.0 {
        0.0
    } else {
        bound.next_up()
    }
}

/// `2^1020`: the expansions of the core require every operand of a sum to be
/// below it in magnitude.
const SUM_OPERAND_LIMIT: f64 = 1.1235582092889474e307;

/// Whether every component of `e` is below `SUM_OPERAND_LIMIT` in magnitude, as
/// the components of an operand of a sum of expansions must be.
#[inline]
fn summable(e: &[f64]) -> bool {
    e.iter().all(|x| x.abs() < SUM_OPERAND_LIMIT)
}

/// Whether every component of `e` lies in the product domain of the core, zero
/// or a magnitude from `2^-480` to `2^480`, as the components of an operand of
/// an exact product must.
#[inline]
fn multipliable(e: &[f64]) -> bool {
    e.iter().all(|&x| eft::in_product_domain(x))
}

/// The span of a coordinate just written, or the error of `write_exact` once
/// the arena is reset to `mark`: [`ExactError::OutOfRange`] when `in_range` is
/// false.
#[inline]
fn written(
    scratch: &mut Scratch<'_>,
    mark: Mark,
    span: Result<Span, ScratchExhausted>,
    in_range: bool,
) -> Result<Span, ExactError> {
    match span {
        Ok(span) if in_range => Ok(span),
        Ok(_) => {
            scratch.reset(mark);
            Err(ExactError::OutOfRange)
        }
        Err(e) => {
            scratch.reset(mark);
            Err(e.into())
        }
    }
}