mathverse_core/
precision.rs

1//! Floating-point comparison and rounding utilities.
2//!
3//! This module provides tools for safe floating-point arithmetic, including
4//! multiple comparison strategies (absolute, relative, ULP-based) and
5//! rounding utilities.
6
7use crate::traits::Real;
8
9/// Machine epsilon for `f64` (2^-52).
10pub const EPS: f64 = 2.220446049250313e-16;
11/// Machine epsilon for `f32` (2^-23).
12pub const F32_EPS: f32 = 1.1920929e-7;
13
14/// Absolute-tolerance equality: `|a - b| <= eps`.
15///
16/// Prefer [`almost_eq_rel`] when operands span many orders of magnitude.
17///
18/// # Examples
19///
20/// ```
21/// use mathverse_core::precision::almost_eq;
22///
23/// assert!(almost_eq(0.1 + 0.2, 0.3, 1e-15));
24/// assert!(!almost_eq(0.1 + 0.2, 0.3, 1e-17));
25/// ```
26#[must_use]
27#[inline]
28pub fn almost_eq<T: Real>(a: T, b: T, eps: T) -> bool {
29    (a - b).abs() <= eps
30}
31
32/// Relative-tolerance equality: `|a - b| <= tol * max(|a|, |b|)`.
33/// Equal when both are zero.
34///
35/// # Examples
36///
37/// ```
38/// use mathverse_core::precision::almost_eq_rel;
39///
40/// assert!(almost_eq_rel(1e15, 1e15 + 1.0, 1e-12));
41/// assert!(almost_eq_rel(0.0, 0.0, 1e-10));
42/// ```
43#[must_use]
44#[inline]
45pub fn almost_eq_rel<T: Real>(a: T, b: T, tol: T) -> bool {
46    let m = a.abs().max(b.abs());
47    m == T::zero() || (a - b).abs() <= tol * m
48}
49
50/// Combined absolute + relative tolerance comparison.
51///
52/// Returns `true` when `|a - b| <= max(abs_tol, rel_tol * max(|a|, |b|))`.
53/// This is the recommended general-purpose float comparison.
54///
55/// # Examples
56///
57/// ```
58/// use mathverse_core::precision::is_close;
59///
60/// assert!(is_close(1.0, 1.0 + 1e-16, 1e-12, 1e-15));
61/// assert!(!is_close(1.0, 2.0, 1e-12, 1e-15));
62/// assert!(is_close(0.0, 0.0, 1e-12, 1e-15));
63/// ```
64#[must_use]
65#[inline]
66pub fn is_close<T: Real>(a: T, b: T, rel_tol: T, abs_tol: T) -> bool {
67    let diff = (a - b).abs();
68    let m = a.abs().max(b.abs());
69    diff <= abs_tol || diff <= rel_tol * m
70}
71
72/// ULP-based floating-point equality for `f64`.
73///
74/// Returns `true` if `a` and `b` differ by at most `max_ulp` representable
75/// values. Handles signed zeros and infinities correctly.
76///
77/// # Examples
78///
79/// ```
80/// use mathverse_core::precision::almost_eq_ulp;
81///
82/// assert!(almost_eq_ulp(1.0, 1.0, 0));
83/// assert!(almost_eq_ulp(f64::INFINITY, f64::INFINITY, 0));
84/// assert!(!almost_eq_ulp(f64::NAN, f64::NAN, 0));
85/// ```
86#[must_use]
87pub fn almost_eq_ulp(a: f64, b: f64, max_ulp: i64) -> bool {
88    if a == b {
89        return true;
90    }
91    if a.is_nan() || b.is_nan() {
92        return false;
93    }
94    if a.is_infinite() || b.is_infinite() {
95        return a == b;
96    }
97    let bits_a = a.to_bits();
98    let bits_b = b.to_bits();
99    let diff = (bits_a as i64 - bits_b as i64).abs();
100    diff <= max_ulp
101}
102
103/// Check if `f64` is subnormal (denormal).
104///
105/// # Examples
106///
107/// ```
108/// use mathverse_core::precision::is_subnormal;
109///
110/// assert!(is_subnormal(f64::MIN_POSITIVE / 2.0));
111/// assert!(!is_subnormal(1.0));
112/// ```
113#[must_use]
114#[inline]
115pub fn is_subnormal(x: f64) -> bool {
116    x.is_finite() && x != 0.0 && x.abs() < f64::MIN_POSITIVE
117}
118
119/// Check if `f64` is NaN.
120///
121/// # Examples
122///
123/// ```
124/// use mathverse_core::precision::is_nan;
125///
126/// assert!(is_nan(f64::NAN));
127/// assert!(!is_nan(1.0));
128/// ```
129#[must_use]
130#[inline]
131pub const fn is_nan(x: f64) -> bool {
132    x.is_nan()
133}
134
135/// Check if `f64` is infinite (positive or negative).
136///
137/// # Examples
138///
139/// ```
140/// use mathverse_core::precision::is_infinite;
141///
142/// assert!(is_infinite(f64::INFINITY));
143/// assert!(!is_infinite(1.0));
144/// ```
145#[must_use]
146#[inline]
147pub const fn is_infinite(x: f64) -> bool {
148    x.is_infinite()
149}
150
151/// Machine epsilon for a generic `Real` type.
152///
153/// Returns the difference between 1.0 and the next representable value.
154///
155/// # Examples
156///
157/// ```
158/// use mathverse_core::precision::{epsilon, EPS, almost_eq};
159///
160/// assert!(almost_eq(epsilon::<f64>(), EPS, 1e-20));
161/// ```
162#[must_use]
163#[inline]
164pub fn epsilon<T: Real>() -> T {
165    T::from_f64(EPS)
166}
167
168/// Next representable `f64` toward positive infinity.
169///
170/// # Examples
171///
172/// ```
173/// use mathverse_core::precision::next_float;
174///
175/// let next = next_float(1.0);
176/// assert!(next > 1.0);
177/// ```
178#[must_use]
179pub fn next_float(x: f64) -> f64 {
180    if x.is_nan() {
181        return x;
182    }
183    if x == f64::INFINITY {
184        return f64::INFINITY;
185    }
186    let bits = x.to_bits();
187    if x.is_sign_negative() && x == 0.0 {
188        f64::from_bits(1)
189    } else if x >= 0.0 {
190        f64::from_bits(bits + 1)
191    } else {
192        f64::from_bits(bits - 1)
193    }
194}
195
196/// Next representable `f64` toward negative infinity.
197///
198/// # Examples
199///
200/// ```
201/// use mathverse_core::precision::{next_float, prev_float};
202///
203/// let next = next_float(1.0);
204/// assert_eq!(prev_float(next), 1.0);
205/// ```
206#[must_use]
207pub fn prev_float(x: f64) -> f64 {
208    if x.is_nan() {
209        return x;
210    }
211    if x == f64::NEG_INFINITY {
212        return f64::NEG_INFINITY;
213    }
214    let bits = x.to_bits();
215    if x.is_sign_positive() && x == 0.0 {
216        f64::from_bits(0x8000_0000_0000_0001)
217    } else if x > 0.0 {
218        f64::from_bits(bits - 1)
219    } else {
220        f64::from_bits(bits + 1)
221    }
222}
223
224/// Copy the sign of `b` to the magnitude of `a`.
225///
226/// # Examples
227///
228/// ```
229/// use mathverse_core::precision::copysign;
230///
231/// assert_eq!(copysign(3.0, -1.0), -3.0);
232/// assert_eq!(copysign(3.0, 1.0), 3.0);
233/// ```
234#[must_use]
235#[inline]
236pub fn copysign<T: Real>(a: T, b: T) -> T {
237    if b.is_negative() {
238        a.abs().neg()
239    } else {
240        a.abs()
241    }
242}
243
244/// Absolute difference: `|a - b|`.
245///
246/// # Examples
247///
248/// ```
249/// use mathverse_core::precision::abs_diff;
250///
251/// assert_eq!(abs_diff(5.0, 3.0), 2.0);
252/// assert_eq!(abs_diff(3.0, 5.0), 2.0);
253/// ```
254#[must_use]
255#[inline]
256pub fn abs_diff<T: Real>(a: T, b: T) -> T {
257    (a - b).abs()
258}
259
260/// Relative difference: `|a - b| / max(|a|, |b|)`.
261///
262/// Returns 0 when both are zero.
263///
264/// # Examples
265///
266/// ```
267/// use mathverse_core::precision::relative_diff;
268///
269/// assert!((relative_diff(10.0_f64, 9.0) - 0.1).abs() < 1e-12);
270/// assert_eq!(relative_diff(0.0, 0.0), 0.0);
271/// ```
272#[must_use]
273#[inline]
274pub fn relative_diff<T: Real>(a: T, b: T) -> T {
275    let diff = (a - b).abs();
276    let m = a.abs().max(b.abs());
277    if m == T::zero() {
278        T::zero()
279    } else {
280        diff / m
281    }
282}
283
284/// Round `x` to `decimals` decimal places.
285///
286/// # Examples
287///
288/// ```
289/// use mathverse_core::precision::round_to;
290///
291/// assert_eq!(round_to(2.345678, 3), 2.346);
292/// assert_eq!(round_to(-1.234, 2), -1.23);
293/// ```
294#[must_use]
295pub fn round_to<T: Real>(x: T, decimals: i32) -> T {
296    let f = T::from_f64(10f64.powi(decimals));
297    (x * f).round() / f
298}
299
300/// Round `x` to `n` significant figures.
301///
302/// # Examples
303///
304/// ```
305/// use mathverse_core::precision::significant_figures;
306///
307/// assert_eq!(significant_figures(0.00234567, 3), 0.00235);
308/// assert_eq!(significant_figures(12345.0, 2), 12000.0);
309/// ```
310#[must_use]
311pub fn significant_figures<T: Real>(x: T, n: i32) -> T {
312    if x == T::zero() {
313        return T::zero();
314    }
315    let shift = T::from_f64(f64::from(n - 1)) - x.abs().log10().floor();
316    round_to(x, shift.to_f64() as i32)
317}
318
319/// Distance from `x` to the next representable `f64` toward +infinity.
320///
321/// # Examples
322///
323/// ```
324/// use mathverse_core::precision::ulp;
325///
326/// assert!(ulp(1.0) < 1e-15);
327/// assert!(ulp(-1.0) < 1e-15);
328/// ```
329#[must_use]
330pub fn ulp(x: f64) -> f64 {
331    if !x.is_finite() || x == 0.0 {
332        return f64::MIN_POSITIVE;
333    }
334    let b = x.to_bits();
335    let b = if x.is_sign_negative() {
336        b.wrapping_sub(1)
337    } else {
338        b.wrapping_add(1)
339    };
340    f64::from_bits(b) - x
341}
342
343/// Round `x` up to the nearest multiple of `m`.
344///
345/// Returns `x` if `m` is zero.
346///
347/// # Examples
348///
349/// ```
350/// use mathverse_core::precision::ceil_to_multiple;
351///
352/// assert_eq!(ceil_to_multiple(7.0, 5.0), 10.0);
353/// assert_eq!(ceil_to_multiple(10.0, 5.0), 10.0);
354/// assert_eq!(ceil_to_multiple(0.0, 3.0), 0.0);
355/// ```
356#[must_use]
357pub fn ceil_to_multiple<T: Real>(x: T, m: T) -> T {
358    if m == T::zero() {
359        return x;
360    }
361    (x / m).ceil() * m
362}
363
364/// Round `x` down to the nearest multiple of `m`.
365///
366/// Returns `x` if `m` is zero.
367///
368/// # Examples
369///
370/// ```
371/// use mathverse_core::precision::floor_to_multiple;
372///
373/// assert_eq!(floor_to_multiple(7.0, 5.0), 5.0);
374/// assert_eq!(floor_to_multiple(10.0, 5.0), 10.0);
375/// assert_eq!(floor_to_multiple(0.0, 3.0), 0.0);
376/// ```
377#[must_use]
378pub fn floor_to_multiple<T: Real>(x: T, m: T) -> T {
379    if m == T::zero() {
380        return x;
381    }
382    (x / m).floor() * m
383}
384
385/// Round `x` to the nearest multiple of `m`.
386///
387/// Returns `x` if `m` is zero.
388///
389/// # Examples
390///
391/// ```
392/// use mathverse_core::precision::round_to_multiple;
393///
394/// assert_eq!(round_to_multiple(7.0, 5.0), 5.0);
395/// assert_eq!(round_to_multiple(8.0, 5.0), 10.0);
396/// assert_eq!(round_to_multiple(0.0, 3.0), 0.0);
397/// ```
398#[must_use]
399pub fn round_to_multiple<T: Real>(x: T, m: T) -> T {
400    if m == T::zero() {
401        return x;
402    }
403    (x / m).round() * m
404}
405
406#[cfg(test)]
407mod tests {
408    use super::*;
409
410    #[test]
411    fn comparisons() {
412        assert!(almost_eq(0.1 + 0.2, 0.3, 1e-15));
413        assert!(!almost_eq(0.1 + 0.2, 0.3, 1e-17));
414        assert!(almost_eq_rel(1e15, 1e15 + 1.0, 1e-12));
415        assert!(almost_eq_rel(0.0, 0.0, 1e-10));
416    }
417
418    #[test]
419    fn rounding() {
420        assert_eq!(round_to(3.14159, 2), 3.14);
421        assert_eq!(significant_figures(12345.0, 2), 12000.0);
422        assert_eq!(significant_figures(0.0, 3), 0.0);
423    }
424
425    #[test]
426    fn ulp_is_small() {
427        assert!(ulp(1.0) < 1e-15);
428        assert!(ulp(-1.0) < 1e-15);
429    }
430
431    #[test]
432    fn is_close_combined() {
433        assert!(is_close(1.0, 1.0 + 1e-16, 1e-12, 1e-15));
434        assert!(!is_close(1.0, 2.0, 1e-12, 1e-15));
435        assert!(is_close(0.0, 0.0, 1e-12, 1e-15));
436        assert!(is_close(0.0, 1e-20, 1e-12, 1e-15));
437    }
438
439    #[test]
440    fn ulp_comparison() {
441        assert!(almost_eq_ulp(1.0, 1.0, 0));
442        assert!(almost_eq_ulp(1.0, 1.0 + 1e-15, 5));
443        assert!(!almost_eq_ulp(1.0, 2.0, 100));
444        assert!(almost_eq_ulp(f64::INFINITY, f64::INFINITY, 0));
445        assert!(!almost_eq_ulp(f64::NAN, f64::NAN, 0));
446    }
447
448    #[test]
449    fn float_helpers() {
450        assert!(is_subnormal(f64::MIN_POSITIVE / 2.0));
451        assert!(!is_subnormal(1.0));
452        assert!(is_nan(f64::NAN));
453        assert!(!is_nan(1.0));
454        assert!(is_infinite(f64::INFINITY));
455        assert!(!is_infinite(1.0));
456        assert!(almost_eq(epsilon::<f64>(), EPS, 1e-20));
457    }
458
459    #[test]
460    fn next_prev_float() {
461        let next = next_float(1.0);
462        assert!(next > 1.0);
463        assert_eq!(prev_float(next), 1.0);
464        assert_eq!(next_float(f64::INFINITY), f64::INFINITY);
465    }
466
467    #[test]
468    fn copysign_abs_diff() {
469        assert_eq!(copysign(3.0, -1.0), -3.0);
470        assert_eq!(copysign(3.0, 1.0), 3.0);
471        assert_eq!(abs_diff(5.0, 3.0), 2.0);
472        assert_eq!(relative_diff(10.0, 9.0), 0.1);
473        assert_eq!(relative_diff(0.0, 0.0), 0.0);
474    }
475
476    #[test]
477    fn multiple_rounding() {
478        assert_eq!(ceil_to_multiple(7.0, 5.0), 10.0);
479        assert_eq!(ceil_to_multiple(10.0, 5.0), 10.0);
480        assert_eq!(ceil_to_multiple(0.0, 3.0), 0.0);
481        assert_eq!(floor_to_multiple(7.0, 5.0), 5.0);
482        assert_eq!(floor_to_multiple(10.0, 5.0), 10.0);
483        assert_eq!(floor_to_multiple(0.0, 3.0), 0.0);
484        assert_eq!(round_to_multiple(7.0, 5.0), 5.0);
485        assert_eq!(round_to_multiple(8.0, 5.0), 10.0);
486        assert_eq!(round_to_multiple(0.0, 3.0), 0.0);
487    }
488}
mathverse_core/precision.rs

mathverse_core/
precision.rs
