1// Copyright 2014-2020 Optimal Computing (NZ) Ltd.
2// Licensed under the MIT license. See LICENSE for details.
34use super::Ulps;
56/// ApproxEqUlps is a trait for approximate equality comparisons.
7/// The associated type Flt is a floating point type which implements Ulps, and is
8/// required so that this trait can be implemented for compound types (e.g. vectors),
9/// not just for the floats themselves.
10pub trait ApproxEqUlps {
11type Flt: Ulps;
1213/// This method tests for `self` and `other` values to be approximately equal
14 /// within ULPs (Units of Least Precision) floating point representations.
15 /// Differing signs are always unequal with this method, and zeroes are only
16 /// equal to zeroes. Use approx_eq() from the ApproxEq trait if that is more
17 /// appropriate.
18fn approx_eq_ulps(&self, other: &Self, ulps: <Self::Flt as Ulps>::U) -> bool;
1920/// This method tests for `self` and `other` values to be not approximately
21 /// equal within ULPs (Units of Least Precision) floating point representations.
22 /// Differing signs are always unequal with this method, and zeroes are only
23 /// equal to zeroes. Use approx_eq() from the ApproxEq trait if that is more
24 /// appropriate.
25#[inline]
26fn approx_ne_ulps(&self, other: &Self, ulps: <Self::Flt as Ulps>::U) -> bool {
27 !self.approx_eq_ulps(other, ulps)
28 }
29}
3031impl ApproxEqUlps for f32 {
32type Flt = f32;
3334fn approx_eq_ulps(&self, other: &f32, ulps: i32) -> bool {
35// -0 and +0 are drastically far in ulps terms, so
36 // we need a special case for that.
37if *self==*other { return true; }
3839// Handle differing signs as a special case, even if
40 // they are very close, most people consider them
41 // unequal.
42if self.is_sign_positive() != other.is_sign_positive() { return false; }
4344let diff: i32 = self.ulps(other);
45 diff >= -ulps && diff <= ulps
46 }
47}
4849#[test]
50fn f32_approx_eq_ulps_test1() {
51let f: f32 = 0.1_f32;
52let mut sum: f32 = 0.0_f32;
53for _ in 0_isize..10_isize { sum += f; }
54let product: f32 = f * 10.0_f32;
55assert!(sum != product); // Should not be directly equal:
56assert!(sum.approx_eq_ulps(&product,1) == true); // But should be close
57assert!(sum.approx_eq_ulps(&product,0) == false);
58}
59#[test]
60fn f32_approx_eq_ulps_test2() {
61let x: f32 = 1000000_f32;
62let y: f32 = 1000000.1_f32;
63assert!(x != y); // Should not be directly equal
64assert!(x.approx_eq_ulps(&y,2) == true);
65assert!(x.approx_eq_ulps(&y,1) == false);
66}
67#[test]
68fn f32_approx_eq_ulps_test_zeroes() {
69let x: f32 = 0.0_f32;
70let y: f32 = -0.0_f32;
71assert!(x.approx_eq_ulps(&y,0) == true);
72}
7374impl ApproxEqUlps for f64 {
75type Flt = f64;
7677fn approx_eq_ulps(&self, other: &f64, ulps: i64) -> bool {
78// -0 and +0 are drastically far in ulps terms, so
79 // we need a special case for that.
80if *self==*other { return true; }
8182// Handle differing signs as a special case, even if
83 // they are very close, most people consider them
84 // unequal.
85if self.is_sign_positive() != other.is_sign_positive() { return false; }
8687let diff: i64 = self.ulps(other);
88 diff >= -ulps && diff <= ulps
89 }
90}
9192#[test]
93fn f64_approx_eq_ulps_test1() {
94let f: f64 = 0.1_f64;
95let mut sum: f64 = 0.0_f64;
96for _ in 0_isize..10_isize { sum += f; }
97let product: f64 = f * 10.0_f64;
98assert!(sum != product); // Should not be directly equal:
99assert!(sum.approx_eq_ulps(&product,1) == true); // But should be close
100assert!(sum.approx_eq_ulps(&product,0) == false);
101}
102#[test]
103fn f64_approx_eq_ulps_test2() {
104let x: f64 = 1000000_f64;
105let y: f64 = 1000000.0000000003_f64;
106assert!(x != y); // Should not be directly equal
107assert!(x.approx_eq_ulps(&y,3) == true);
108assert!(x.approx_eq_ulps(&y,2) == false);
109}
110#[test]
111fn f64_approx_eq_ulps_test_zeroes() {
112let x: f64 = 0.0_f64;
113let y: f64 = -0.0_f64;
114assert!(x.approx_eq_ulps(&y,0) == true);
115}