implicit-predicates 0.0.1

// Generated by implicit-predicates-gen from predicates/orient3d.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 stages of the predicate `orient3d`.
//!
//! Every stage takes the points `a`, `b`, `c` and `d` of the definition as `p`,
//! in that order, and returns the sign of the predicate, or `None` when it
//! cannot decide. [`staged`] runs them in order.

use implicit_predicates_core::{Sign, eft, expansion};

/// The degree `k` of the predicate in the differences of coordinates.
pub const DEGREE: u32 = 3;

/// The coefficient of `D^k` in the error bound of [`filter`], `(EPSILON * D +
/// KAPPA * E) * D^(k - 1)`, rounded up. `E` is the `e` given to the filter and
/// `D` is `E` plus the largest magnitude of a difference of coordinates.
pub const EPSILON: f64 = 5.10702591803201e-15;

/// The coefficient of `E * D^(k - 1)` in the error bound of [`filter`], rounded
/// up.
pub const KAPPA: f64 = 18.00000001676381;

/// The smallest `D` for which [`filter`] decides, `2^-324`: below it, underflow
/// could exceed the error bound.
pub const D_MIN: f64 = 2.926047721682624e-98;

/// The largest `D` for which [`filter`] decides, `2^339`: above it, a computed
/// value could overflow.
pub const D_MAX: f64 = 1.1198723710889021e102;

/// `(min, max)`, `(2^-188, 2^239)`: [`exact_differences`] decides when every
/// difference of coordinates is exact in `f64` and every nonzero one has a
/// magnitude in `min..=max`.
pub const DIFFERENCES_RANGE: (f64, f64) = (2.5489470578119236e-57, 8.834235323891922e71);

/// `(min, max)`, `(2^-188, 2^238)`: [`expansions_translated`] decides when
/// every nonzero coordinate has a magnitude in `min..=max`.
pub const TRANSLATED_RANGE: (f64, f64) = (2.5489470578119236e-57, 4.417117661945961e71);

/// `(min, max)`, `(2^-188, 2^239)`: [`expansions_untranslated`] decides when
/// every nonzero coordinate has a magnitude in `min..=max`.
pub const UNTRANSLATED_RANGE: (f64, f64) = (2.5489470578119236e-57, 8.834235323891922e71);

/// The sign of the predicate from an evaluation of its main form in `f64`, when
/// the rounding error cannot change it.
///
/// The points given are approximations of exact points. `e` is the sum of the
/// two largest per-point errors, where the error of a point is an upper bound
/// on the absolute error of each coordinate of its approximation; `e` must be
/// finite, non-negative and rounded upward, and it is `0.0` for explicit
/// points. A NaN or negative `e` never decides.
///
/// `None` means that this stage cannot decide: the value is too close to zero
/// for its error bound, the scale of the differences is outside
/// `D_MIN..=D_MAX`, `e` is negative or NaN, or an input is not finite. The
/// result is never `Some(Sign::Zero)`.
#[inline]
pub fn filter(p: [&[f64; 3]; 4], e: f64) -> Option<Sign> {
    let t0 = p[0][0] - p[3][0];
    let t1 = p[0][1] - p[3][1];
    let t2 = p[0][2] - p[3][2];
    let t3 = p[1][0] - p[3][0];
    let t4 = p[1][1] - p[3][1];
    let t5 = p[1][2] - p[3][2];
    let t6 = p[2][0] - p[3][0];
    let t7 = p[2][1] - p[3][1];
    let t8 = p[2][2] - p[3][2];
    let t9 = t4 * t8;
    let t10 = t5 * t7;
    let t11 = t9 - t10;
    let t12 = t0 * t11;
    let t13 = t3 * t8;
    let t14 = t5 * t6;
    let t15 = t13 - t14;
    let t16 = t1 * t15;
    let t17 = t12 - t16;
    let t18 = t3 * t7;
    let t19 = t4 * t6;
    let t20 = t18 - t19;
    let t21 = t2 * t20;
    let det = t17 + t21;
    // A NaN leaf can be skipped by the comparisons, but it makes `det` NaN.
    let delta = t0.abs();
    let delta = max(delta, t1.abs());
    let delta = max(delta, t2.abs());
    let delta = max(delta, t3.abs());
    let delta = max(delta, t4.abs());
    let delta = max(delta, t5.abs());
    let delta = max(delta, t6.abs());
    let delta = max(delta, t7.abs());
    let delta = max(delta, t8.abs());
    // The bound takes 6 f64 operations; with the rounding of `d` counted at
    // each of its 3 uses, at most 8 roundings reach it.
    let d = delta + e;
    let bound = (EPSILON * d + KAPPA * e) * d * d;
    let decided = (e >= 0.0) & (D_MIN..=D_MAX).contains(&d) & (det.abs() > bound);
    let sign = if det > 0.0 {
        Sign::Positive
    } else {
        Sign::Negative
    };
    decided.then_some(sign)
}

/// The exact sign of the predicate: its main form evaluated with expansions on
/// the differences of coordinates, each computed in `f64`.
///
/// `None` means that this stage cannot decide: some difference is not exact in
/// `f64`, or some nonzero difference has a magnitude outside
/// [`DIFFERENCES_RANGE`], as when a coordinate is not finite.
#[cold]
#[inline(never)]
pub fn exact_differences(p: [&[f64; 3]; 4]) -> Option<Sign> {
    // Stack: 88 f64 (704 bytes) of expansions.
    let (v0, r0) = eft::two_diff(p[0][0], p[3][0]);
    let (v1, r1) = eft::two_diff(p[0][1], p[3][1]);
    let (v2, r2) = eft::two_diff(p[0][2], p[3][2]);
    let (v3, r3) = eft::two_diff(p[1][0], p[3][0]);
    let (v4, r4) = eft::two_diff(p[1][1], p[3][1]);
    let (v5, r5) = eft::two_diff(p[1][2], p[3][2]);
    let (v6, r6) = eft::two_diff(p[2][0], p[3][0]);
    let (v7, r7) = eft::two_diff(p[2][1], p[3][1]);
    let (v8, r8) = eft::two_diff(p[2][2], p[3][2]);
    // Each difference must be exact, and within the range.
    let errors = [r0, r1, r2, r3, r4, r5, r6, r7, r8];
    let differences = [v0, v1, v2, v3, v4, v5, v6, v7, v8];
    if !(errors.into_iter().all(|r| r == 0.0)
        && differences
            .into_iter()
            .all(|x| in_range(x, DIFFERENCES_RANGE)))
    {
        return None;
    }
    let mut e9 = [0.0; 2];
    let l9 = expansion::from_two_product(v4, v8, &mut e9);
    let mut e10 = [0.0; 2];
    let l10 = expansion::from_two_product(v5, v7, &mut e10);
    let mut e11 = [0.0; 4];
    let l11 = expansion::diff(&e9[..l9], &e10[..l10], &mut e11);
    let mut e12 = [0.0; 8];
    let l12 = expansion::scale(&e11[..l11], v0, &mut e12);
    let mut e13 = [0.0; 2];
    let l13 = expansion::from_two_product(v3, v8, &mut e13);
    let mut e14 = [0.0; 2];
    let l14 = expansion::from_two_product(v5, v6, &mut e14);
    let mut e15 = [0.0; 4];
    let l15 = expansion::diff(&e13[..l13], &e14[..l14], &mut e15);
    let mut e16 = [0.0; 8];
    let l16 = expansion::scale(&e15[..l15], v1, &mut e16);
    let mut e17 = [0.0; 16];
    let l17 = expansion::diff(&e12[..l12], &e16[..l16], &mut e17);
    let mut e18 = [0.0; 2];
    let l18 = expansion::from_two_product(v3, v7, &mut e18);
    let mut e19 = [0.0; 2];
    let l19 = expansion::from_two_product(v4, v6, &mut e19);
    let mut e20 = [0.0; 4];
    let l20 = expansion::diff(&e18[..l18], &e19[..l19], &mut e20);
    let mut e21 = [0.0; 8];
    let l21 = expansion::scale(&e20[..l20], v2, &mut e21);
    let mut e22 = [0.0; 24];
    let l22 = expansion::sum(&e17[..l17], &e21[..l21], &mut e22);
    Some(expansion::sign(&e22[..l22]))
}

/// The exact sign of the predicate: its main form evaluated with expansions,
/// each difference of coordinates as an exact expansion of two components.
///
/// `None` means that this stage cannot decide: some nonzero coordinate has a
/// magnitude outside [`TRANSLATED_RANGE`] or is not finite.
#[cold]
#[inline(never)]
pub fn expansions_translated(p: [&[f64; 3]; 4]) -> Option<Sign> {
    // Stack: 722 f64 (5776 bytes) of expansions, 96 of them the scratch of the
    // products of two expansions.
    let coordinates = [
        p[0][0], p[0][1], p[0][2], p[1][0], p[1][1], p[1][2], p[2][0], p[2][1], p[2][2], p[3][0],
        p[3][1], p[3][2],
    ];
    if !coordinates
        .into_iter()
        .all(|x| in_range(x, TRANSLATED_RANGE))
    {
        return None;
    }
    let mut scaled = [0.0; 32];
    let mut spare = [0.0; 64];
    let mut e0 = [0.0; 2];
    let l0 = expansion::from_two_diff(p[0][0], p[3][0], &mut e0);
    let mut e1 = [0.0; 2];
    let l1 = expansion::from_two_diff(p[0][1], p[3][1], &mut e1);
    let mut e2 = [0.0; 2];
    let l2 = expansion::from_two_diff(p[0][2], p[3][2], &mut e2);
    let mut e3 = [0.0; 2];
    let l3 = expansion::from_two_diff(p[1][0], p[3][0], &mut e3);
    let mut e4 = [0.0; 2];
    let l4 = expansion::from_two_diff(p[1][1], p[3][1], &mut e4);
    let mut e5 = [0.0; 2];
    let l5 = expansion::from_two_diff(p[1][2], p[3][2], &mut e5);
    let mut e6 = [0.0; 2];
    let l6 = expansion::from_two_diff(p[2][0], p[3][0], &mut e6);
    let mut e7 = [0.0; 2];
    let l7 = expansion::from_two_diff(p[2][1], p[3][1], &mut e7);
    let mut e8 = [0.0; 2];
    let l8 = expansion::from_two_diff(p[2][2], p[3][2], &mut e8);
    let mut e9 = [0.0; 8];
    let l9 = product(&e4[..l4], &e8[..l8], &mut scaled, &mut spare, &mut e9);
    let mut e10 = [0.0; 8];
    let l10 = product(&e5[..l5], &e7[..l7], &mut scaled, &mut spare, &mut e10);
    let mut e11 = [0.0; 16];
    let l11 = expansion::diff(&e9[..l9], &e10[..l10], &mut e11);
    let mut e12 = [0.0; 64];
    let l12 = product(&e0[..l0], &e11[..l11], &mut scaled, &mut spare, &mut e12);
    let mut e13 = [0.0; 8];
    let l13 = product(&e3[..l3], &e8[..l8], &mut scaled, &mut spare, &mut e13);
    let mut e14 = [0.0; 8];
    let l14 = product(&e5[..l5], &e6[..l6], &mut scaled, &mut spare, &mut e14);
    let mut e15 = [0.0; 16];
    let l15 = expansion::diff(&e13[..l13], &e14[..l14], &mut e15);
    let mut e16 = [0.0; 64];
    let l16 = product(&e1[..l1], &e15[..l15], &mut scaled, &mut spare, &mut e16);
    let mut e17 = [0.0; 128];
    let l17 = expansion::diff(&e12[..l12], &e16[..l16], &mut e17);
    let mut e18 = [0.0; 8];
    let l18 = product(&e3[..l3], &e7[..l7], &mut scaled, &mut spare, &mut e18);
    let mut e19 = [0.0; 8];
    let l19 = product(&e4[..l4], &e6[..l6], &mut scaled, &mut spare, &mut e19);
    let mut e20 = [0.0; 16];
    let l20 = expansion::diff(&e18[..l18], &e19[..l19], &mut e20);
    let mut e21 = [0.0; 64];
    let l21 = product(&e2[..l2], &e20[..l20], &mut scaled, &mut spare, &mut e21);
    let mut e22 = [0.0; 192];
    let l22 = expansion::sum(&e17[..l17], &e21[..l21], &mut e22);
    Some(expansion::sign(&e22[..l22]))
}

/// The exact sign of the predicate: its form for exact evaluation, marked
/// `#[exact]` in the definition, evaluated with expansions on the coordinates.
///
/// `None` means that this stage cannot decide: some nonzero coordinate has a
/// magnitude outside [`UNTRANSLATED_RANGE`] or is not finite.
#[cold]
#[inline(never)]
pub fn expansions_untranslated(p: [&[f64; 3]; 4]) -> Option<Sign> {
    // Stack: 568 f64 (4544 bytes) of expansions.
    let coordinates = [
        p[0][0], p[0][1], p[0][2], p[1][0], p[1][1], p[1][2], p[2][0], p[2][1], p[2][2], p[3][0],
        p[3][1], p[3][2],
    ];
    if !coordinates
        .into_iter()
        .all(|x| in_range(x, UNTRANSLATED_RANGE))
    {
        return None;
    }
    let mut e3 = [0.0; 2];
    let l3 = expansion::from_two_product(p[2][1], p[3][2], &mut e3);
    let mut e6 = [0.0; 2];
    let l6 = expansion::from_two_product(p[2][2], p[3][1], &mut e6);
    let mut e7 = [0.0; 4];
    let l7 = expansion::diff(&e3[..l3], &e6[..l6], &mut e7);
    let mut e8 = [0.0; 8];
    let l8 = expansion::scale(&e7[..l7], p[1][0], &mut e8);
    let mut e11 = [0.0; 2];
    let l11 = expansion::from_two_product(p[2][0], p[3][2], &mut e11);
    let mut e13 = [0.0; 2];
    let l13 = expansion::from_two_product(p[2][2], p[3][0], &mut e13);
    let mut e14 = [0.0; 4];
    let l14 = expansion::diff(&e11[..l11], &e13[..l13], &mut e14);
    let mut e15 = [0.0; 8];
    let l15 = expansion::scale(&e14[..l14], p[1][1], &mut e15);
    let mut e16 = [0.0; 16];
    let l16 = expansion::diff(&e8[..l8], &e15[..l15], &mut e16);
    let mut e18 = [0.0; 2];
    let l18 = expansion::from_two_product(p[2][0], p[3][1], &mut e18);
    let mut e19 = [0.0; 2];
    let l19 = expansion::from_two_product(p[2][1], p[3][0], &mut e19);
    let mut e20 = [0.0; 4];
    let l20 = expansion::diff(&e18[..l18], &e19[..l19], &mut e20);
    let mut e21 = [0.0; 8];
    let l21 = expansion::scale(&e20[..l20], p[1][2], &mut e21);
    let mut e22 = [0.0; 24];
    let l22 = expansion::sum(&e16[..l16], &e21[..l21], &mut e22);
    let mut e24 = [0.0; 8];
    let l24 = expansion::scale(&e7[..l7], p[0][0], &mut e24);
    let mut e26 = [0.0; 8];
    let l26 = expansion::scale(&e14[..l14], p[0][1], &mut e26);
    let mut e27 = [0.0; 16];
    let l27 = expansion::diff(&e24[..l24], &e26[..l26], &mut e27);
    let mut e29 = [0.0; 8];
    let l29 = expansion::scale(&e20[..l20], p[0][2], &mut e29);
    let mut e30 = [0.0; 24];
    let l30 = expansion::sum(&e27[..l27], &e29[..l29], &mut e30);
    let mut e31 = [0.0; 2];
    let l31 = expansion::from_two_product(p[1][1], p[3][2], &mut e31);
    let mut e32 = [0.0; 2];
    let l32 = expansion::from_two_product(p[1][2], p[3][1], &mut e32);
    let mut e33 = [0.0; 4];
    let l33 = expansion::diff(&e31[..l31], &e32[..l32], &mut e33);
    let mut e34 = [0.0; 8];
    let l34 = expansion::scale(&e33[..l33], p[0][0], &mut e34);
    let mut e35 = [0.0; 2];
    let l35 = expansion::from_two_product(p[1][0], p[3][2], &mut e35);
    let mut e36 = [0.0; 2];
    let l36 = expansion::from_two_product(p[1][2], p[3][0], &mut e36);
    let mut e37 = [0.0; 4];
    let l37 = expansion::diff(&e35[..l35], &e36[..l36], &mut e37);
    let mut e38 = [0.0; 8];
    let l38 = expansion::scale(&e37[..l37], p[0][1], &mut e38);
    let mut e39 = [0.0; 16];
    let l39 = expansion::diff(&e34[..l34], &e38[..l38], &mut e39);
    let mut e40 = [0.0; 2];
    let l40 = expansion::from_two_product(p[1][0], p[3][1], &mut e40);
    let mut e41 = [0.0; 2];
    let l41 = expansion::from_two_product(p[1][1], p[3][0], &mut e41);
    let mut e42 = [0.0; 4];
    let l42 = expansion::diff(&e40[..l40], &e41[..l41], &mut e42);
    let mut e43 = [0.0; 8];
    let l43 = expansion::scale(&e42[..l42], p[0][2], &mut e43);
    let mut e44 = [0.0; 24];
    let l44 = expansion::sum(&e39[..l39], &e43[..l43], &mut e44);
    let mut e45 = [0.0; 2];
    let l45 = expansion::from_two_product(p[1][1], p[2][2], &mut e45);
    let mut e46 = [0.0; 2];
    let l46 = expansion::from_two_product(p[1][2], p[2][1], &mut e46);
    let mut e47 = [0.0; 4];
    let l47 = expansion::diff(&e45[..l45], &e46[..l46], &mut e47);
    let mut e48 = [0.0; 8];
    let l48 = expansion::scale(&e47[..l47], p[0][0], &mut e48);
    let mut e49 = [0.0; 2];
    let l49 = expansion::from_two_product(p[1][0], p[2][2], &mut e49);
    let mut e50 = [0.0; 2];
    let l50 = expansion::from_two_product(p[1][2], p[2][0], &mut e50);
    let mut e51 = [0.0; 4];
    let l51 = expansion::diff(&e49[..l49], &e50[..l50], &mut e51);
    let mut e52 = [0.0; 8];
    let l52 = expansion::scale(&e51[..l51], p[0][1], &mut e52);
    let mut e53 = [0.0; 16];
    let l53 = expansion::diff(&e48[..l48], &e52[..l52], &mut e53);
    let mut e54 = [0.0; 2];
    let l54 = expansion::from_two_product(p[1][0], p[2][1], &mut e54);
    let mut e55 = [0.0; 2];
    let l55 = expansion::from_two_product(p[1][1], p[2][0], &mut e55);
    let mut e56 = [0.0; 4];
    let l56 = expansion::diff(&e54[..l54], &e55[..l55], &mut e56);
    let mut e57 = [0.0; 8];
    let l57 = expansion::scale(&e56[..l56], p[0][2], &mut e57);
    let mut e58 = [0.0; 24];
    let l58 = expansion::sum(&e53[..l53], &e57[..l57], &mut e58);
    let mut e60 = [0.0; 24];
    e60[..l22].copy_from_slice(&e22[..l22]);
    expansion::negate(&mut e60[..l22]);
    let l60 = l22;
    let mut e61 = [0.0; 48];
    let l61 = expansion::sum(&e60[..l60], &e30[..l30], &mut e61);
    let mut e62 = [0.0; 72];
    let l62 = expansion::diff(&e61[..l61], &e44[..l44], &mut e62);
    let mut e63 = [0.0; 96];
    let l63 = expansion::sum(&e62[..l62], &e58[..l58], &mut e63);
    Some(expansion::sign(&e63[..l63]))
}

/// The sign of the predicate for explicit points, from the first stage that
/// decides.
///
/// Runs [`filter`] with `e = 0.0`, then [`exact_differences`], then
/// [`expansions_translated`], then [`expansions_untranslated`]. `None` means
/// that no stage can decide.
#[inline]
pub fn staged(p: [&[f64; 3]; 4]) -> Option<Sign> {
    filter(p, 0.0)
        .or_else(|| exact_differences(p))
        .or_else(|| expansions_translated(p))
        .or_else(|| expansions_untranslated(p))
}

/// The larger of `a` and `b`, or `b` if they are unordered.
#[inline]
fn max(a: f64, b: f64) -> f64 {
    if a > b { a } else { b }
}

/// Whether `x` is zero or has a magnitude in `min..=max`; false for NaN and
/// infinities.
fn in_range(x: f64, (min, max): (f64, f64)) -> bool {
    x == 0.0 || (min..=max).contains(&x.abs())
}

/// Writes the expansion of `x * y` to `out` and returns its length: the sum of
/// the scalings of `y` by each component of `x`, each scaling written to
/// `scaled`, the partial sums alternating between `spare` and `out`.
fn product(x: &[f64], y: &[f64], scaled: &mut [f64], spare: &mut [f64], out: &mut [f64]) -> usize {
    let Some((&first, rest)) = x.split_first() else {
        return 0;
    };
    // Start so that the last partial sum is written to `out`.
    let (mut partial, mut next) = if rest.len() % 2 == 0 {
        (out, spare)
    } else {
        (spare, out)
    };
    let mut length = expansion::scale(y, first, partial);
    for &component in rest {
        let scaled_length = expansion::scale(y, component, scaled);
        length = expansion::sum(&partial[..length], &scaled[..scaled_length], next);
        core::mem::swap(&mut partial, &mut next);
    }
    length
}