implicit-predicates 0.0.2

// 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 on explicit points, and
//! [`staged_in`] on points of any type.

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

/// The degree `G_k` of the predicate in the differences of coordinates of each
/// axis `k`, for x, y and z.
pub const DEGREES: [u32; 3] = [1, 1, 1];

/// The coefficient of `d^G`, the product of `d_k^G_k` over the axes, in the
/// error bound of [`filter`], rounded up. On each axis `k`, `d_k` is the
/// largest magnitude of a difference of coordinates on that axis plus `e[k]`,
/// the uncertainty given to the filter, and the bound is `EPSILON*d^G` plus
/// `KAPPA[k]*e[k]*d^G/d_k` for each axis of positive degree, computed without
/// divisions.
pub const EPSILON: f64 = 5.10702591803201e-15;

/// The coefficient of `e[k]*d^G/d_k` in the error bound of [`filter`], for each
/// axis `k`, rounded up; zero on an axis of degree zero.
pub const KAPPA: [f64; 3] = [6.000000005587936, 6.000000005587936, 6.000000005587936];

/// The smallest `d_k`, on each axis of positive degree, 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_k`, on each axis of positive degree, 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 most components of a coordinate that [`expansions_translated_in`] and
/// [`expansions_untranslated_in`] accept.
pub const MAX_COORD_COMPONENTS: usize = 16;

/// `(min, max)`, `(2^-188, 2^234)`: [`expansions_translated_in`] decides when
/// every nonzero component of every coordinate has a magnitude in `min..=max`.
pub const SPANS_TRANSLATED_RANGE: (f64, f64) = (2.5489470578119236e-57, 2.7606985387162255e70);

/// `(min, max)`, `(2^-188, 2^235)`: [`expansions_untranslated_in`] decides when
/// every nonzero component of every coordinate has a magnitude in `min..=max`.
pub const SPANS_UNTRANSLATED_RANGE: (f64, f64) = (2.5489470578119236e-57, 5.521397077432451e70);

/// 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[k]` is the sum of
/// the two largest errors of the points on axis `k` (0, 1 and 2 for x, y and
/// z), where the error of a point on an axis is an upper bound on the absolute
/// error of its coordinate on that axis; each `e[k]` must be finite,
/// non-negative and rounded upward, and they are `0.0` for explicit points. A
/// NaN or negative `e[k]` never decides.
///
/// `Some(Sign::Zero)` means that some axis is flat: every difference of
/// coordinates on it is zero and so is its `e[k]`, so the predicate is exactly
/// zero. `None` means that this stage cannot decide: the value is too close to
/// zero for its error bound, the scale `d_k` of the differences of some axis is
/// outside `D_MIN..=D_MAX`, some `e[k]` is negative or NaN, or an input is not
/// finite.
#[allow(clippy::manual_range_contains)]
#[inline(always)]
pub fn filter(p: [&[f64; 3]; 4], e: [f64; 3]) -> 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;
    // The largest magnitude of the leaves of each axis. A NaN leaf can be
    // skipped by the comparisons, but it makes `det` NaN.
    let delta_x = t0.abs();
    let delta_x = max(delta_x, t3.abs());
    let delta_x = max(delta_x, t6.abs());
    let delta_y = t1.abs();
    let delta_y = max(delta_y, t4.abs());
    let delta_y = max(delta_y, t7.abs());
    let delta_z = t2.abs();
    let delta_z = max(delta_z, t5.abs());
    let delta_z = max(delta_z, t8.abs());
    let dx = delta_x + e[0];
    let dy = delta_y + e[1];
    let dz = delta_z + e[2];
    // Each term is a chain of products that starts with its constant. A term
    // with `e[k] == 0` is zero whenever the filter can decide, and is left out,
    // so that it costs nothing when `e` is a constant zero.
    let bound = EPSILON * dx * dy * dz;
    let bound = if e[0] > 0.0 {
        bound + KAPPA[0] * e[0] * dy * dz
    } else {
        bound
    };
    let bound = if e[1] > 0.0 {
        bound + KAPPA[1] * e[1] * dx * dz
    } else {
        bound
    };
    let bound = if e[2] > 0.0 {
        bound + KAPPA[2] * e[2] * dx * dy
    } else {
        bound
    };
    let decided = (e[0] >= 0.0)
        & (e[1] >= 0.0)
        & (e[2] >= 0.0)
        & (D_MIN <= dx)
        & (dx <= D_MAX)
        & (D_MIN <= dy)
        & (dy <= D_MAX)
        & (D_MIN <= dz)
        & (dz <= D_MAX)
        & (det.abs() > bound);
    let sign = if det > 0.0 {
        Sign::Positive
    } else {
        Sign::Negative
    };
    // A flat axis, with every leaf and its `e[k]` zero, makes the predicate
    // exactly zero; a NaN leaf skipped by the maxima makes `det` NaN, and each
    // `e[k]` must be finite. A flat axis is below `D_MIN`, so `decided` is
    // false then.
    let flat = (dx == 0.0) | (dy == 0.0) | (dz == 0.0);
    let zero = (e[0] >= 0.0)
        & (e[1] >= 0.0)
        & (e[2] >= 0.0)
        & (e[0] <= f64::MAX)
        & (e[1] <= f64::MAX)
        & (e[2] <= f64::MAX)
        & flat
        & (det == 0.0);
    let sign = if zero { Sign::Zero } else { sign };
    (decided | zero).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 exact sign of the predicate for points given by their exact coordinates:
/// its main form evaluated with expansions in `scratch`, each difference of
/// coordinates computed exactly.
///
/// `Ok(None)` means that this stage cannot decide: some point has a weight
/// other than [`Weight::One`], or some coordinate has more than
/// [`MAX_COORD_COMPONENTS`] components or a nonzero component of magnitude
/// outside [`SPANS_TRANSLATED_RANGE`] or not finite.
///
/// Requirement: every coordinate is a valid expansion of `scratch`: nonzero
/// components of strictly increasing magnitude that do not overlap, as
/// [`Point::exact`](implicit_predicates_core::Point::exact) writes them. This
/// is not checked, and the result is unspecified otherwise.
///
/// The arena is left as it was, on success and on error.
///
/// # Errors
///
/// [`ScratchExhausted`] if `scratch` runs out of room.
#[cold]
#[inline(never)]
pub fn expansions_translated_in(
    scratch: &mut Scratch<'_>,
    p: [&ExactCoords; 4],
) -> Result<Option<Sign>, ScratchExhausted> {
    fn evaluate(scratch: &mut Scratch<'_>, c: [[Span; 3]; 4]) -> Result<Sign, ScratchExhausted> {
        let e0 = scratch.diff(c[0][0], c[3][0])?;
        let e1 = scratch.diff(c[0][1], c[3][1])?;
        let e2 = scratch.diff(c[0][2], c[3][2])?;
        let e3 = scratch.diff(c[1][0], c[3][0])?;
        let e4 = scratch.diff(c[1][1], c[3][1])?;
        let e5 = scratch.diff(c[1][2], c[3][2])?;
        let e6 = scratch.diff(c[2][0], c[3][0])?;
        let e7 = scratch.diff(c[2][1], c[3][1])?;
        let e8 = scratch.diff(c[2][2], c[3][2])?;
        let e9 = scratch.product(e4, e8)?;
        let e10 = scratch.product(e5, e7)?;
        let e11 = scratch.diff(e9, e10)?;
        let e12 = scratch.product(e0, e11)?;
        let e13 = scratch.product(e3, e8)?;
        let e14 = scratch.product(e5, e6)?;
        let e15 = scratch.diff(e13, e14)?;
        let e16 = scratch.product(e1, e15)?;
        let e17 = scratch.diff(e12, e16)?;
        let e18 = scratch.product(e3, e7)?;
        let e19 = scratch.product(e4, e6)?;
        let e20 = scratch.diff(e18, e19)?;
        let e21 = scratch.product(e2, e20)?;
        let e22 = scratch.sum(e17, e21)?;
        Ok(scratch.sign(e22))
    }
    // Every coordinate read must be short enough and within the range; the
    // arena is reset to its mark whatever the evaluation returns.
    if !p.iter().all(|q| q.w == Weight::One) {
        return Ok(None);
    }
    let c = p.map(|q| [q.x, q.y, q.z]);
    let coordinates = [
        c[0][0], c[0][1], c[0][2], c[1][0], c[1][1], c[1][2], c[2][0], c[2][1], c[2][2], c[3][0],
        c[3][1], c[3][2],
    ];
    if !coordinates
        .into_iter()
        .all(|s| span_in_range(scratch.get(s), SPANS_TRANSLATED_RANGE))
    {
        return Ok(None);
    }
    let mark = scratch.mark();
    let sign = evaluate(scratch, c);
    scratch.reset(mark);
    sign.map(Some)
}

/// The exact sign of the predicate for points given by their exact coordinates:
/// its form for exact evaluation, marked `#[exact]` in the definition,
/// evaluated with expansions in `scratch` on the coordinates.
///
/// `Ok(None)` means that this stage cannot decide: some point has a weight
/// other than [`Weight::One`], or some coordinate has more than
/// [`MAX_COORD_COMPONENTS`] components or a nonzero component of magnitude
/// outside [`SPANS_UNTRANSLATED_RANGE`] or not finite.
///
/// Requirement: every coordinate is a valid expansion of `scratch`: nonzero
/// components of strictly increasing magnitude that do not overlap, as
/// [`Point::exact`](implicit_predicates_core::Point::exact) writes them. This
/// is not checked, and the result is unspecified otherwise.
///
/// The arena is left as it was, on success and on error.
///
/// # Errors
///
/// [`ScratchExhausted`] if `scratch` runs out of room.
#[cold]
#[inline(never)]
pub fn expansions_untranslated_in(
    scratch: &mut Scratch<'_>,
    p: [&ExactCoords; 4],
) -> Result<Option<Sign>, ScratchExhausted> {
    fn evaluate(scratch: &mut Scratch<'_>, c: [[Span; 3]; 4]) -> Result<Sign, ScratchExhausted> {
        let e3 = scratch.product(c[2][1], c[3][2])?;
        let e6 = scratch.product(c[2][2], c[3][1])?;
        let e7 = scratch.diff(e3, e6)?;
        let e8 = scratch.product(c[1][0], e7)?;
        let e11 = scratch.product(c[2][0], c[3][2])?;
        let e13 = scratch.product(c[2][2], c[3][0])?;
        let e14 = scratch.diff(e11, e13)?;
        let e15 = scratch.product(c[1][1], e14)?;
        let e16 = scratch.diff(e8, e15)?;
        let e18 = scratch.product(c[2][0], c[3][1])?;
        let e19 = scratch.product(c[2][1], c[3][0])?;
        let e20 = scratch.diff(e18, e19)?;
        let e21 = scratch.product(c[1][2], e20)?;
        let e22 = scratch.sum(e16, e21)?;
        let e24 = scratch.product(c[0][0], e7)?;
        let e26 = scratch.product(c[0][1], e14)?;
        let e27 = scratch.diff(e24, e26)?;
        let e29 = scratch.product(c[0][2], e20)?;
        let e30 = scratch.sum(e27, e29)?;
        let e31 = scratch.product(c[1][1], c[3][2])?;
        let e32 = scratch.product(c[1][2], c[3][1])?;
        let e33 = scratch.diff(e31, e32)?;
        let e34 = scratch.product(c[0][0], e33)?;
        let e35 = scratch.product(c[1][0], c[3][2])?;
        let e36 = scratch.product(c[1][2], c[3][0])?;
        let e37 = scratch.diff(e35, e36)?;
        let e38 = scratch.product(c[0][1], e37)?;
        let e39 = scratch.diff(e34, e38)?;
        let e40 = scratch.product(c[1][0], c[3][1])?;
        let e41 = scratch.product(c[1][1], c[3][0])?;
        let e42 = scratch.diff(e40, e41)?;
        let e43 = scratch.product(c[0][2], e42)?;
        let e44 = scratch.sum(e39, e43)?;
        let e45 = scratch.product(c[1][1], c[2][2])?;
        let e46 = scratch.product(c[1][2], c[2][1])?;
        let e47 = scratch.diff(e45, e46)?;
        let e48 = scratch.product(c[0][0], e47)?;
        let e49 = scratch.product(c[1][0], c[2][2])?;
        let e50 = scratch.product(c[1][2], c[2][0])?;
        let e51 = scratch.diff(e49, e50)?;
        let e52 = scratch.product(c[0][1], e51)?;
        let e53 = scratch.diff(e48, e52)?;
        let e54 = scratch.product(c[1][0], c[2][1])?;
        let e55 = scratch.product(c[1][1], c[2][0])?;
        let e56 = scratch.diff(e54, e55)?;
        let e57 = scratch.product(c[0][2], e56)?;
        let e58 = scratch.sum(e53, e57)?;
        let e59 = scratch.push(&[])?;
        let e60 = scratch.diff(e59, e22)?;
        let e61 = scratch.sum(e60, e30)?;
        let e62 = scratch.diff(e61, e44)?;
        let e63 = scratch.sum(e62, e58)?;
        Ok(scratch.sign(e63))
    }
    // Every coordinate read must be short enough and within the range; the
    // arena is reset to its mark whatever the evaluation returns.
    if !p.iter().all(|q| q.w == Weight::One) {
        return Ok(None);
    }
    let c = p.map(|q| [q.x, q.y, q.z]);
    let coordinates = [
        c[0][0], c[0][1], c[0][2], c[1][0], c[1][1], c[1][2], c[2][0], c[2][1], c[2][2], c[3][0],
        c[3][1], c[3][2],
    ];
    if !coordinates
        .into_iter()
        .all(|s| span_in_range(scratch.get(s), SPANS_UNTRANSLATED_RANGE))
    {
        return Ok(None);
    }
    let mark = scratch.mark();
    let sign = evaluate(scratch, c);
    scratch.reset(mark);
    sign.map(Some)
}

/// The sign of the predicate for explicit points, from the first stage that
/// decides.
///
/// Runs [`filter`] with `e = [0.0; 3]`, 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; 3])
        .or_else(|| exact_differences(p))
        .or_else(|| expansions_translated(p))
        .or_else(|| expansions_untranslated(p))
}

/// [`filter`] on the approximations of points of any type, with the uncertainty
/// of their errors: `[0.0; 3]` when every point is explicit, and otherwise
/// `uncertainty` of their errors, from `implicit_predicates_core::filter`.
///
/// The errors are not read when every point is explicit. When that is known
/// from the types, as with `[f64; 3]`, whose `is_explicit` is the constant
/// `true`, this is `filter` with `e = [0.0; 3]`.
#[inline]
pub fn filter_points<P0: Point, P1: Point, P2: Point, P3: Point>(
    p0: &P0,
    p1: &P1,
    p2: &P2,
    p3: &P3,
) -> Option<Sign> {
    let p = [p0.approx(), p1.approx(), p2.approx(), p3.approx()];
    let e = if p0.is_explicit() & p1.is_explicit() & p2.is_explicit() & p3.is_explicit() {
        [0.0; 3]
    } else {
        implicit_predicates_core::filter::uncertainty([p0.err(), p1.err(), p2.err(), p3.err()])
    };
    filter(p, e)
}

/// The exact stages of [`staged`], after its filter: [`exact_differences`],
/// then [`expansions_translated`], then [`expansions_untranslated`]. `None`
/// means that no stage can decide.
#[cold]
#[inline(never)]
pub fn exact_stages(p: [&[f64; 3]; 4]) -> Option<Sign> {
    exact_differences(p)
        .or_else(|| expansions_translated(p))
        .or_else(|| expansions_untranslated(p))
}

/// The exact sign of the predicate for points of any type:
/// [`expansions_translated_in`], then [`expansions_untranslated_in`], on the
/// exact coordinates of the points written to `scratch`, which is reset to its
/// mark whatever they return. `Ok(None)` means that no stage can decide: the
/// inputs of some point are outside the range of its exact arithmetic
/// ([`ExactError::OutOfRange`]), or none of the stages decides.
///
/// # Errors
///
/// [`ScratchExhausted`] if `scratch` runs out of room.
#[cold]
#[inline(never)]
pub fn exact_stages_in<P0: Point, P1: Point, P2: Point, P3: Point>(
    scratch: &mut Scratch<'_>,
    p0: &P0,
    p1: &P1,
    p2: &P2,
    p3: &P3,
) -> Result<Option<Sign>, ScratchExhausted> {
    fn run<P0: Point, P1: Point, P2: Point, P3: Point>(
        scratch: &mut Scratch<'_>,
        p0: &P0,
        p1: &P1,
        p2: &P2,
        p3: &P3,
    ) -> Result<Option<Sign>, ScratchExhausted> {
        let Some(q0) = exact(scratch, p0)? else {
            return Ok(None);
        };
        let Some(q1) = exact(scratch, p1)? else {
            return Ok(None);
        };
        let Some(q2) = exact(scratch, p2)? else {
            return Ok(None);
        };
        let Some(q3) = exact(scratch, p3)? else {
            return Ok(None);
        };
        let q = [&q0, &q1, &q2, &q3];
        if let Some(sign) = expansions_translated_in(scratch, q)? {
            return Ok(Some(sign));
        }
        expansions_untranslated_in(scratch, q)
    }
    let mark = scratch.mark();
    let result = run(scratch, p0, p1, p2, p3);
    scratch.reset(mark);
    result
}

/// The sign of the predicate for points of any type, from the first stage that
/// decides.
///
/// Runs [`filter_points`]; then, when every point is explicit, [`exact_stages`]
/// on their coordinates, which does not use `scratch`; otherwise
/// [`exact_stages_in`], on their exact coordinates written to `scratch`.
/// `scratch` is left as it was. `Ok(None)` means that no stage can decide.
///
/// # Errors
///
/// [`ScratchExhausted`] if `scratch` runs out of room.
#[inline]
pub fn staged_in<P0: Point, P1: Point, P2: Point, P3: Point>(
    scratch: &mut Scratch<'_>,
    p0: &P0,
    p1: &P1,
    p2: &P2,
    p3: &P3,
) -> Result<Option<Sign>, ScratchExhausted> {
    if let Some(sign) = filter_points(p0, p1, p2, p3) {
        return Ok(Some(sign));
    }
    if p0.is_explicit() & p1.is_explicit() & p2.is_explicit() & p3.is_explicit() {
        return Ok(exact_stages([
            p0.approx(),
            p1.approx(),
            p2.approx(),
            p3.approx(),
        ]));
    }
    exact_stages_in(scratch, p0, p1, p2, p3)
}

/// The exact coordinates of the point `p`, written to `scratch`, or `None` if
/// its inputs are outside the range of its exact arithmetic.
fn exact(
    scratch: &mut Scratch<'_>,
    p: &impl Point,
) -> Result<Option<ExactCoords>, ScratchExhausted> {
    match p.exact(scratch) {
        Ok(c) => Ok(Some(c)),
        Err(ExactError::OutOfRange) => Ok(None),
        Err(ExactError::Exhausted(e)) => Err(e),
    }
}

/// 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())
}

/// Whether the expansion `e` has at most [`MAX_COORD_COMPONENTS`] components,
/// each zero or of a magnitude in `min..=max`.
fn span_in_range(e: &[f64], range: (f64, f64)) -> bool {
    e.len() <= MAX_COORD_COMPONENTS && e.iter().all(|&x| in_range(x, range))
}

/// 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
}