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, Transcendental};
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::epsilon()
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()
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(pow10(decimals));
297    (x * f).round() / f
298}
299
300fn pow10(n: i32) -> f64 {
301    let mut acc = 1.0;
302    for _ in 0..n.unsigned_abs() {
303        acc *= 10.0;
304    }
305    if n < 0 { 1.0 / acc } else { acc }
306}
307
308/// Round `x` to `n` significant figures.
309///
310/// # Examples
311///
312/// ```
313/// use mathverse_core::precision::significant_figures;
314///
315/// assert_eq!(significant_figures(0.00234567, 3), 0.00235);
316/// assert_eq!(significant_figures(12345.0, 2), 12000.0);
317/// ```
318#[must_use]
319#[allow(clippy::cast_lossless)] // From<i32> for f64 is unavailable in no_std
320pub fn significant_figures<T: Real + Transcendental>(x: T, n: i32) -> T {
321    if x == T::zero() {
322        return T::zero();
323    }
324    let shift = T::from_f64((n - 1) as f64) - x.abs().log10().floor();
325    round_to(x, shift.to_f64() as i32)
326}
327
328/// Distance from `x` to the next representable `f64` toward +infinity.
329///
330/// # Examples
331///
332/// ```
333/// use mathverse_core::precision::ulp;
334///
335/// assert!(ulp(1.0) < 1e-15);
336/// assert!(ulp(-1.0) < 1e-15);
337/// ```
338#[must_use]
339pub fn ulp(x: f64) -> f64 {
340    if !x.is_finite() || x == 0.0 {
341        return f64::MIN_POSITIVE;
342    }
343    let b = x.to_bits();
344    let b = if x.is_sign_negative() {
345        b.wrapping_sub(1)
346    } else {
347        b.wrapping_add(1)
348    };
349    f64::from_bits(b) - x
350}
351
352/// Round `x` up to the nearest multiple of `m`.
353///
354/// Returns `x` if `m` is zero.
355///
356/// # Examples
357///
358/// ```
359/// use mathverse_core::precision::ceil_to_multiple;
360///
361/// assert_eq!(ceil_to_multiple(7.0, 5.0), 10.0);
362/// assert_eq!(ceil_to_multiple(10.0, 5.0), 10.0);
363/// assert_eq!(ceil_to_multiple(0.0, 3.0), 0.0);
364/// ```
365#[must_use]
366pub fn ceil_to_multiple<T: Real>(x: T, m: T) -> T {
367    if m == T::zero() {
368        return x;
369    }
370    (x / m).ceil() * m
371}
372
373/// Round `x` down to the nearest multiple of `m`.
374///
375/// Returns `x` if `m` is zero.
376///
377/// # Examples
378///
379/// ```
380/// use mathverse_core::precision::floor_to_multiple;
381///
382/// assert_eq!(floor_to_multiple(7.0, 5.0), 5.0);
383/// assert_eq!(floor_to_multiple(10.0, 5.0), 10.0);
384/// assert_eq!(floor_to_multiple(0.0, 3.0), 0.0);
385/// ```
386#[must_use]
387pub fn floor_to_multiple<T: Real>(x: T, m: T) -> T {
388    if m == T::zero() {
389        return x;
390    }
391    (x / m).floor() * m
392}
393
394/// Round `x` to the nearest multiple of `m`.
395///
396/// Returns `x` if `m` is zero.
397///
398/// # Examples
399///
400/// ```
401/// use mathverse_core::precision::round_to_multiple;
402///
403/// assert_eq!(round_to_multiple(7.0, 5.0), 5.0);
404/// assert_eq!(round_to_multiple(8.0, 5.0), 10.0);
405/// assert_eq!(round_to_multiple(0.0, 3.0), 0.0);
406/// ```
407#[must_use]
408pub fn round_to_multiple<T: Real>(x: T, m: T) -> T {
409    if m == T::zero() {
410        return x;
411    }
412    (x / m).round() * m
413}
414
415#[cfg(test)]
416mod tests {
417    use super::*;
418
419    #[test]
420    fn comparisons() {
421        assert!(almost_eq(0.1 + 0.2, 0.3, 1e-15));
422        assert!(!almost_eq(0.1 + 0.2, 0.3, 1e-17));
423        assert!(almost_eq_rel(1e15, 1e15 + 1.0, 1e-12));
424        assert!(almost_eq_rel(0.0, 0.0, 1e-10));
425    }
426
427    #[test]
428    fn rounding() {
429        assert_eq!(round_to(3.14159, 2), 3.14);
430        assert_eq!(significant_figures(12345.0, 2), 12000.0);
431        assert_eq!(significant_figures(0.0, 3), 0.0);
432    }
433
434    #[test]
435    fn ulp_is_small() {
436        assert!(ulp(1.0) < 1e-15);
437        assert!(ulp(-1.0) < 1e-15);
438    }
439
440    #[test]
441    fn is_close_combined() {
442        assert!(is_close(1.0, 1.0 + 1e-16, 1e-12, 1e-15));
443        assert!(!is_close(1.0, 2.0, 1e-12, 1e-15));
444        assert!(is_close(0.0, 0.0, 1e-12, 1e-15));
445        assert!(is_close(0.0, 1e-20, 1e-12, 1e-15));
446    }
447
448    #[test]
449    fn ulp_comparison() {
450        assert!(almost_eq_ulp(1.0, 1.0, 0));
451        assert!(almost_eq_ulp(1.0, 1.0 + 1e-15, 5));
452        assert!(!almost_eq_ulp(1.0, 2.0, 100));
453        assert!(almost_eq_ulp(f64::INFINITY, f64::INFINITY, 0));
454        assert!(!almost_eq_ulp(f64::NAN, f64::NAN, 0));
455    }
456
457    #[test]
458    fn float_helpers() {
459        assert!(is_subnormal(f64::MIN_POSITIVE / 2.0));
460        assert!(!is_subnormal(1.0));
461        assert!(is_nan(f64::NAN));
462        assert!(!is_nan(1.0));
463        assert!(is_infinite(f64::INFINITY));
464        assert!(!is_infinite(1.0));
465        assert!(almost_eq(epsilon::<f64>(), EPS, 1e-20));
466    }
467
468    #[test]
469    fn next_prev_float() {
470        let next = next_float(1.0);
471        assert!(next > 1.0);
472        assert_eq!(prev_float(next), 1.0);
473        assert_eq!(next_float(f64::INFINITY), f64::INFINITY);
474    }
475
476    #[test]
477    fn copysign_abs_diff() {
478        assert_eq!(copysign(3.0, -1.0), -3.0);
479        assert_eq!(copysign(3.0, 1.0), 3.0);
480        assert_eq!(abs_diff(5.0, 3.0), 2.0);
481        assert_eq!(relative_diff(10.0, 9.0), 0.1);
482        assert_eq!(relative_diff(0.0, 0.0), 0.0);
483    }
484
485    #[test]
486    fn multiple_rounding() {
487        assert_eq!(ceil_to_multiple(7.0, 5.0), 10.0);
488        assert_eq!(ceil_to_multiple(10.0, 5.0), 10.0);
489        assert_eq!(ceil_to_multiple(0.0, 3.0), 0.0);
490        assert_eq!(floor_to_multiple(7.0, 5.0), 5.0);
491        assert_eq!(floor_to_multiple(10.0, 5.0), 10.0);
492        assert_eq!(floor_to_multiple(0.0, 3.0), 0.0);
493        assert_eq!(round_to_multiple(7.0, 5.0), 5.0);
494        assert_eq!(round_to_multiple(8.0, 5.0), 10.0);
495        assert_eq!(round_to_multiple(0.0, 3.0), 0.0);
496    }
497}
mathverse_core/precision.rs

mathverse_core/
precision.rs
