geop_core_math/scalars/
scal_in_f64.rs

1use core::fmt::Display;
2
3use super::{Field, Ring, Scalar};
4use crate::geop_error::{GeopError, GeopResult};
5
6// ── IEEE 754 outward-rounding helpers ─────────────────────────────────────────
7
8#[inline]
9fn next_up(x: f64) -> f64 {
10    if x.is_nan() {
11        return x;
12    }
13    if x == f64::INFINITY {
14        return f64::INFINITY;
15    }
16    if x == 0.0 {
17        return f64::MIN_POSITIVE;
18    }
19    let bits = x.to_bits();
20    let bits = if x > 0.0 { bits + 1 } else { bits - 1 };
21    f64::from_bits(bits)
22}
23
24#[inline]
25fn next_down(x: f64) -> f64 {
26    if x.is_nan() {
27        return x;
28    }
29    if x == f64::NEG_INFINITY {
30        return f64::NEG_INFINITY;
31    }
32    if x == 0.0 {
33        return -f64::MIN_POSITIVE;
34    }
35    let bits = x.to_bits();
36    let bits = if x < 0.0 { bits + 1 } else { bits - 1 };
37    f64::from_bits(bits)
38}
39
40// ── Interval trigonometry ─────────────────────────────────────────────────────
41
42/// True iff `[lo, hi]` contains some `target + k * period` for an integer `k`
43/// — i.e. whether the interval passes through one of `sin`/`cos`'s
44/// extrema, which sampling only the endpoints could miss entirely.
45#[inline]
46fn contains_periodic(lo: f64, hi: f64, target: f64, period: f64) -> bool {
47    let k = ((lo - target) / period).ceil();
48    target + k * period <= hi
49}
50
51/// Outward-rounded enclosure of `sin`/`cos([lo, hi])`, correct even though
52/// `f64::sin`/`f64::cos` are not guaranteed correctly rounded: the true
53/// extrema at the given `max_at`/`min_at` phases are recognized structurally
54/// (`could_be_periodic`) rather than hunted for numerically, and the sampled
55/// endpoint values are widened outward by a further ULP against roundoff in
56/// the libm call itself.
57fn interval_trig(lo: f64, hi: f64, f: impl Fn(f64) -> f64, max_at: f64, min_at: f64) -> (f64, f64) {
58    use std::f64::consts::TAU;
59    if !lo.is_finite() || !hi.is_finite() || hi - lo >= TAU {
60        return (-1.0, 1.0);
61    }
62    let (a, b) = (f(lo), f(hi));
63    let mut out_lo = next_down(a.min(b));
64    let mut out_hi = next_up(a.max(b));
65    if contains_periodic(lo, hi, max_at, TAU) {
66        out_hi = 1.0;
67    }
68    if contains_periodic(lo, hi, min_at, TAU) {
69        out_lo = -1.0;
70    }
71    (out_lo.max(-1.0), out_hi.min(1.0))
72}
73
74// ── Type ──────────────────────────────────────────────────────────────────────
75
76/// Interval f64 scalar: `[lo, hi]` with outward-rounded arithmetic.
77#[derive(Copy, Clone, Debug, PartialEq)]
78pub struct ScalInF64 {
79    pub lo: f64,
80    pub hi: f64,
81}
82
83impl ScalInF64 {
84    /// Create a proper interval. Panics in debug if lo > hi.
85    #[inline]
86    pub fn new(lo: f64, hi: f64) -> Self {
87        debug_assert!(lo <= hi, "ScalInF64::new: lo ({lo}) > hi ({hi})");
88        ScalInF64 { lo, hi }
89    }
90
91    /// Degenerate (point) interval.
92    #[inline]
93    pub fn degenerate(v: f64) -> Self {
94        ScalInF64 { lo: v, hi: v }
95    }
96}
97
98impl Display for ScalInF64 {
99    fn fmt(&self, f: &mut core::fmt::Formatter<'_>) -> core::fmt::Result {
100        write!(f, "{:.3}", self.to_f64())
101    }
102}
103
104impl Ring for ScalInF64 {
105    fn add(self, other: Self) -> Self {
106        ScalInF64::new(next_down(self.lo + other.lo), next_up(self.hi + other.hi))
107    }
108
109    fn sub(self, other: Self) -> Self {
110        ScalInF64::new(next_down(self.lo - other.hi), next_up(self.hi - other.lo))
111    }
112
113    fn mul(self, other: Self) -> Self {
114        let products = [
115            self.lo * other.lo,
116            self.lo * other.hi,
117            self.hi * other.lo,
118            self.hi * other.hi,
119        ];
120        let lo = products.iter().cloned().fold(f64::INFINITY, f64::min);
121        let hi = products.iter().cloned().fold(f64::NEG_INFINITY, f64::max);
122        ScalInF64::new(next_down(lo), next_up(hi))
123    }
124
125    fn neg(self) -> Self {
126        ScalInF64::new(-self.hi, -self.lo)
127    }
128}
129
130impl Field for ScalInF64 {
131    fn div(self, other: Self) -> GeopResult<Self> {
132        if other.lo <= 0.0 && other.hi >= 0.0 {
133            return Err(GeopError::new(
134                "ScalInF64::div: divisor interval contains zero",
135            ));
136        }
137        // [lo,hi] * [1/hi', 1/lo']
138        let inv_lo = next_down(1.0 / other.hi);
139        let inv_hi = next_up(1.0 / other.lo);
140        let inv = ScalInF64::new(inv_lo, inv_hi);
141        Ok(self.mul(inv))
142    }
143}
144
145// ── Scalar impl ───────────────────────────────────────────────────────────────
146
147impl Scalar for ScalInF64 {
148    const ZERO: Self = ScalInF64 { lo: 0.0, hi: 0.0 };
149    const ONE: Self = ScalInF64 { lo: 1.0, hi: 1.0 };
150    const TWO: Self = ScalInF64 { lo: 2.0, hi: 2.0 };
151
152    // π ≈ 3.141592653589793   (next_down/up computed at compile time as literals)
153    const PI: Self = ScalInF64 {
154        lo: std::f64::consts::PI,   // rounded value is a lower bound for π
155        hi: 3.1415926535897936_f64, // next_up(π)
156    };
157    const E: Self = ScalInF64 {
158        lo: std::f64::consts::E, // rounded value is a lower bound for e
159        hi: 2.7182818284590455_f64,
160    };
161    const INFINITY: Self = ScalInF64 {
162        lo: f64::INFINITY,
163        hi: f64::INFINITY,
164    };
165    const ENTIRE: Self = ScalInF64 {
166        lo: f64::NEG_INFINITY,
167        hi: f64::INFINITY,
168    };
169
170    fn from_i64(v: i64) -> Self {
171        ScalInF64::degenerate(v as f64)
172    }
173    fn from_f64(v: f64) -> Self {
174        ScalInF64::degenerate(v)
175    }
176
177    fn from_ratio(num: i64, den: i64) -> GeopResult<Self> {
178        if den == 0 {
179            return Err(GeopError::new("ScalInF64::from_ratio: denominator is zero"));
180        }
181        let exact = num as f64 / den as f64;
182        Ok(ScalInF64::new(next_down(exact), next_up(exact)))
183    }
184
185    fn abs(self) -> Self {
186        if self.lo >= 0.0 {
187            self
188        } else if self.hi <= 0.0 {
189            ScalInF64::new(-self.hi, -self.lo)
190        } else {
191            ScalInF64::new(0.0, self.lo.abs().max(self.hi.abs()))
192        }
193    }
194
195    fn sqrt(self) -> GeopResult<Self> {
196        if self.hi < 0.0 {
197            return Err(GeopError::new(
198                "ScalInF64::sqrt: interval is definitely negative",
199            ));
200        }
201        let lo_clamped = if self.lo < 0.0 { 0.0 } else { self.lo };
202        Ok(ScalInF64::new(
203            next_down(lo_clamped.sqrt()),
204            next_up(self.hi.sqrt()),
205        ))
206    }
207
208    fn sin(self) -> Self {
209        use std::f64::consts::FRAC_PI_2;
210        let (lo, hi) = interval_trig(self.lo, self.hi, f64::sin, FRAC_PI_2, -FRAC_PI_2);
211        ScalInF64::new(lo, hi)
212    }
213
214    fn cos(self) -> Self {
215        use std::f64::consts::PI;
216        let (lo, hi) = interval_trig(self.lo, self.hi, f64::cos, 0.0, PI);
217        ScalInF64::new(lo, hi)
218    }
219
220    fn could_be_equal(self, other: Self) -> bool {
221        self.lo <= other.hi && other.lo <= self.hi
222    }
223
224    fn definitely_not_equal(self, other: Self) -> bool {
225        self.hi < other.lo || self.lo > other.hi
226    }
227
228    fn could_be_greater(self, other: Self) -> bool {
229        self.hi > other.lo
230    }
231
232    fn definitely_greater(self, other: Self) -> bool {
233        self.lo > other.hi
234    }
235
236    fn could_be_less(self, other: Self) -> bool {
237        self.lo < other.hi
238    }
239
240    fn definitely_less(self, other: Self) -> bool {
241        self.hi < other.lo
242    }
243
244    fn is_infinite(self) -> bool {
245        self.lo == f64::NEG_INFINITY || self.hi == f64::INFINITY
246    }
247
248    fn is_finite(self) -> bool {
249        self.lo.is_finite() && self.hi.is_finite()
250    }
251
252    fn midpoint(self) -> Self {
253        let m = (self.lo + self.hi) / 2.0;
254        ScalInF64::degenerate(m)
255    }
256
257    fn is_sharp(self) -> bool {
258        self.lo == self.hi
259    }
260
261    fn lower(self) -> Self {
262        ScalInF64::degenerate(self.lo)
263    }
264
265    fn upper(self) -> Self {
266        ScalInF64::degenerate(self.hi)
267    }
268
269    fn width(self) -> Self {
270        // Non-negative and sharp, so it can serve as a decidable threshold.
271        let w = (self.hi - self.lo).max(0.0);
272        ScalInF64::new(w, w)
273    }
274
275    fn intersect(self, other: Self) -> Self {
276        let lo = self.lo.max(other.lo);
277        let hi = self.hi.min(other.hi);
278        if lo <= hi {
279            ScalInF64::new(lo, hi)
280        } else if self.hi - self.lo <= other.hi - other.lo {
281            self
282        } else {
283            other
284        }
285    }
286
287    fn to_f64(self) -> f64 {
288        (self.lo + self.hi) / 2.0
289    }
290
291    fn union(self, other: Self) -> Self {
292        ScalInF64::new(self.lo.min(other.lo), self.hi.max(other.hi))
293    }
294
295    fn is_subset_of(self, other: Self) -> bool {
296        other.lo <= self.lo && self.hi <= other.hi
297    }
298}
299
300impl core::ops::Add for ScalInF64 {
301    type Output = Self;
302    fn add(self, rhs: Self) -> Self {
303        Ring::add(self, rhs)
304    }
305}
306impl core::ops::Sub for ScalInF64 {
307    type Output = Self;
308    fn sub(self, rhs: Self) -> Self {
309        Ring::sub(self, rhs)
310    }
311}
312impl core::ops::Mul for ScalInF64 {
313    type Output = Self;
314    fn mul(self, rhs: Self) -> Self {
315        Ring::mul(self, rhs)
316    }
317}
318impl core::ops::Neg for ScalInF64 {
319    type Output = Self;
320    fn neg(self) -> Self {
321        Ring::neg(self)
322    }
323}
324
325impl From<i64> for ScalInF64 {
326    fn from(v: i64) -> Self {
327        ScalInF64::from_i64(v)
328    }
329}
330
331impl From<f64> for ScalInF64 {
332    fn from(v: f64) -> Self {
333        ScalInF64::degenerate(v)
334    }
335}
336
337impl Default for ScalInF64 {
338    fn default() -> Self {
339        ScalInF64::ZERO
340    }
341}
342
343// ── Tests ─────────────────────────────────────────────────────────────────────
344
345#[cfg(test)]
346mod tests {
347    use super::*;
348
349    fn pt(v: f64) -> ScalInF64 {
350        ScalInF64::degenerate(v)
351    }
352
353    fn iv(lo: f64, hi: f64) -> ScalInF64 {
354        ScalInF64::new(lo, hi)
355    }
356
357    #[test]
358    fn arithmetic_add() {
359        let r = pt(1.0).add(pt(2.0));
360        assert!(r.could_be_equal(pt(3.0)));
361    }
362
363    #[test]
364    fn arithmetic_sub() {
365        let r = pt(5.0).sub(pt(3.0));
366        assert!(r.could_be_equal(pt(2.0)));
367    }
368
369    #[test]
370    fn arithmetic_mul() {
371        let r = pt(3.0).mul(pt(4.0));
372        assert!(r.could_be_equal(pt(12.0)));
373    }
374
375    #[test]
376    fn arithmetic_div() {
377        let r = pt(10.0).div(pt(2.0)).unwrap();
378        assert!(r.could_be_equal(pt(5.0)));
379    }
380
381    #[test]
382    fn arithmetic_sqrt_4() {
383        let r = pt(4.0).sqrt().unwrap();
384        assert!(r.could_be_equal(pt(2.0)));
385    }
386
387    #[test]
388    fn arithmetic_sqrt_9() {
389        let r = pt(9.0).sqrt().unwrap();
390        assert!(r.could_be_equal(pt(3.0)));
391    }
392
393    #[test]
394    fn arithmetic_abs() {
395        assert!(pt(-3.0).abs().could_be_equal(pt(3.0)));
396    }
397
398    #[test]
399    fn arithmetic_neg() {
400        assert!(pt(-5.0).neg().could_be_equal(pt(5.0)));
401    }
402
403    #[test]
404    fn sin_at_a_point() {
405        assert!(pt(0.0).sin().could_be_equal(pt(0.0)));
406        assert!(
407            pt(std::f64::consts::FRAC_PI_2)
408                .sin()
409                .could_be_equal(pt(1.0))
410        );
411    }
412
413    #[test]
414    fn cos_at_a_point() {
415        assert!(pt(0.0).cos().could_be_equal(pt(1.0)));
416    }
417
418    #[test]
419    fn sin_over_a_peak_reaches_exactly_one() {
420        // [0, pi] straddles the sin maximum at pi/2; the lower bound stays
421        // outward-rounded (sin(0) = 0 exactly, but sin(pi) is a hair above
422        // zero in f64, and the enclosure must not round that away).
423        let r = iv(0.0, std::f64::consts::PI).sin();
424        assert!(r.hi == 1.0 && r.lo >= -1e-15);
425    }
426
427    #[test]
428    fn cos_over_a_full_turn_is_entire_range() {
429        let r = iv(0.0, std::f64::consts::TAU).cos();
430        assert!(r.lo == -1.0 && r.hi == 1.0);
431    }
432
433    #[test]
434    fn sin_narrow_interval_encloses_true_value() {
435        // A narrow enclosure of 0.5 radians; sin(0.5) should lie in the result.
436        let r = iv(0.5 - 1e-9, 0.5 + 1e-9).sin();
437        let truth = 0.5_f64.sin();
438        assert!(r.lo <= truth && truth <= r.hi);
439    }
440
441    #[test]
442    fn err_div_by_zero() {
443        assert!(pt(1.0).div(ScalInF64::ZERO).is_err());
444    }
445
446    #[test]
447    fn err_sqrt_negative() {
448        assert!(pt(-1.0).sqrt().is_err());
449    }
450
451    #[test]
452    fn sqrt_straddles_zero_ok() {
453        // Interval [-1, 4] straddles zero — sqrt should succeed and contain 2
454        let r = iv(-1.0, 4.0).sqrt().unwrap();
455        assert!(r.could_be_equal(pt(2.0)));
456    }
457
458    #[test]
459    fn overflow_infinity() {
460        let big = pt(f64::MAX);
461        let result = big.mul(big);
462        assert!(result.is_infinite());
463    }
464
465    #[test]
466    fn cmp_definitely_greater() {
467        assert!(pt(3.0).definitely_greater(pt(2.0)));
468    }
469
470    #[test]
471    fn cmp_could_be_less_false() {
472        assert!(!pt(3.0).definitely_less(pt(2.0)));
473    }
474
475    #[test]
476    fn cmp_overlapping_intervals_could_be_equal() {
477        assert!(iv(1.0, 3.0).could_be_equal(iv(2.0, 4.0)));
478    }
479
480    #[test]
481    fn cmp_disjoint_intervals_definitely_not_equal() {
482        assert!(iv(1.0, 2.0).definitely_not_equal(iv(3.0, 4.0)));
483    }
484
485    #[test]
486    fn entire_could_be_equal_point() {
487        assert!(ScalInF64::ENTIRE.could_be_equal(pt(42.0)));
488        assert!(pt(-1e300).could_be_equal(ScalInF64::ENTIRE));
489    }
490
491    #[test]
492    fn entire_could_be_equal_interval() {
493        assert!(ScalInF64::ENTIRE.could_be_equal(iv(-5.0, 5.0)));
494    }
495
496    #[test]
497    fn entire_could_be_greater_and_less() {
498        assert!(ScalInF64::ENTIRE.could_be_greater(pt(1e300)));
499        assert!(ScalInF64::ENTIRE.could_be_less(pt(-1e300)));
500    }
501
502    #[test]
503    fn entire_never_definitely() {
504        assert!(!ScalInF64::ENTIRE.definitely_not_equal(pt(0.0)));
505        assert!(!ScalInF64::ENTIRE.definitely_greater(pt(1e300)));
506        assert!(!ScalInF64::ENTIRE.definitely_less(pt(-1e300)));
507    }
508
509    #[test]
510    fn fp_enclosure_01_plus_02() {
511        // 0.1 + 0.2 is famously not exactly 0.3; the interval should still enclose the true sum
512        let a = iv(next_down(0.1), next_up(0.1));
513        let b = iv(next_down(0.2), next_up(0.2));
514        let c = iv(next_down(0.3), next_up(0.3));
515        let sum = a.add(b);
516        assert!(sum.could_be_equal(c));
517    }
518}
geop_core_math/scalars/scal_in_f64.rs

geop_core_math/scalars/
scal_in_f64.rs
