Function orient3d 

pub fn orient3d(
    a: &impl Point,
    b: &impl Point,
    c: &impl Point,
    d: &impl Point,
) -> Option<Sign>

The orientation of d with respect to the plane through a, b and c: the exact sign of the determinant

| 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);

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));