axiolid_reference/
arithmetic.rs

1//! Arbitrary-length expansion arithmetic.
2//!
3//! An *expansion* is a list of non-overlapping f64 components whose exact sum
4//! is the value it represents. Two f64s can hold a product exactly; a list can
5//! hold an arbitrary determinant exactly. This is what lets `orient3d`,
6//! `incircle`, and `insphere` certify a sign rather than guess one.
7//!
8//! Components are ordered smallest to largest, so the sign of a non-zero
9//! expansion is the sign of its last component.
10//!
11//! # Cost and where it is paid
12//!
13//! These operations allocate. That is deliberate and confined to the *exact*
14//! path, which a filtered predicate reaches only when the floating-point
15//! determinant is too close to zero to trust. The benchmark harness measures
16//! the escalation rate precisely so this cost is a number, not a hope.
17
18use axiolid_contracts::Sign;
19
20use crate::expansion::{two_product, two_sum};
21
22/// Sum of `a` and `b` where `|a| >= |b|` is already known.
23///
24/// Cheaper than [`two_sum`] by two operations. The precondition is not checked
25/// in release builds; violating it silently produces a non-expansion, so every
26/// caller here derives the ordering structurally rather than assuming it.
27#[inline]
28#[must_use]
29fn fast_two_sum(a: f64, b: f64) -> (f64, f64) {
30    let sum = a + b;
31    let b_virtual = sum - a;
32    (sum, b - b_virtual)
33}
34
35/// Grow an expansion by one scalar: `e + b`, exactly.
36///
37/// Sweeps `b` through the components from smallest to largest, carrying the
38/// rounding error forward. Zero components are dropped: they carry no value
39/// and would break the non-overlapping invariant later.
40#[must_use]
41pub fn grow_expansion(e: &[f64], b: f64) -> Vec<f64> {
42    let mut out = Vec::with_capacity(e.len() + 1);
43    let mut carry = b;
44    for &component in e {
45        let (sum, error) = two_sum(carry, component);
46        if error != 0.0 {
47            out.push(error);
48        }
49        carry = sum;
50    }
51    if carry != 0.0 || out.is_empty() {
52        out.push(carry);
53    }
54    out
55}
56
57/// Sum of two expansions, exactly.
58///
59/// Merges the two component lists in magnitude order, then runs a single
60/// carry-propagating pass. Both inputs must be non-overlapping expansions in
61/// increasing magnitude order; the result is the same.
62#[must_use]
63pub fn expansion_sum(a: &[f64], b: &[f64]) -> Vec<f64> {
64    let mut result = a.to_vec();
65    for &component in b {
66        result = grow_expansion(&result, component);
67    }
68    result
69}
70
71/// Scale an expansion by a scalar, exactly.
72///
73/// Each component contributes a two-term product, so the result has at most
74/// twice as many components as the input.
75#[must_use]
76pub fn scale_expansion(e: &[f64], b: f64) -> Vec<f64> {
77    let mut out = Vec::with_capacity(e.len() * 2);
78    let mut carry = 0.0;
79    for &component in e {
80        let (product, product_error) = two_product(component, b);
81        let (sum, error) = two_sum(carry, product_error);
82        if error != 0.0 {
83            out.push(error);
84        }
85        let (new_carry, hi_error) = fast_two_sum(product, sum);
86        if hi_error != 0.0 {
87            out.push(hi_error);
88        }
89        carry = new_carry;
90    }
91    if carry != 0.0 || out.is_empty() {
92        out.push(carry);
93    }
94    out
95}
96
97/// Multiply two expansions exactly.
98///
99/// Each component of `b` scales `a` exactly; the partial expansions are then
100/// accumulated without rounding away any component.
101#[must_use]
102pub fn expansion_product(a: &[f64], b: &[f64]) -> Vec<f64> {
103    let mut out = Vec::from([0.0]);
104    for &component in b {
105        out = expansion_sum(&out, &scale_expansion(a, component));
106    }
107    out
108}
109
110/// Negate every component. Exact: negation is always representable.
111#[must_use]
112pub fn negate_expansion(e: &[f64]) -> Vec<f64> {
113    e.iter().map(|c| -c).collect()
114}
115
116/// Sign of an expansion.
117///
118/// The components are non-overlapping and ordered by increasing magnitude, so
119/// the largest non-zero component dominates the sum and decides the sign.
120#[must_use]
121pub fn expansion_sign(e: &[f64]) -> Sign {
122    for &component in e.iter().rev() {
123        if component > 0.0 {
124            return Sign::Positive;
125        }
126        if component < 0.0 {
127            return Sign::Negative;
128        }
129    }
130    Sign::Zero
131}
axiolid_reference/arithmetic.rs

axiolid_reference/
arithmetic.rs
