1// Copyright 2014-2020 Optimal Computing (NZ) Ltd.
2// Licensed under the MIT license. See LICENSE for details.
3
4use super::Ulps;
5use core::{f32, f64};
6#[cfg(feature = "num-traits")]
7#[allow(unused_imports)]
8use num_traits::float::FloatCore;
9
10/// A margin specifying a maximum distance two floating point values can be while
11/// still being considered equal enough.
12pub trait FloatMargin: Copy + Default {
13 /// A floating-point type used for epsilon values
14 type F;
15
16 /// An integer type used for ulps values
17 type I;
18
19 /// Zero margin
20 fn zero() -> Self;
21
22 /// Set the epsilon value for this margin
23 fn epsilon(self, epsilon: Self::F) -> Self;
24
25 /// Set the ulps value for this margin
26 fn ulps(self, ulps: Self::I) -> Self;
27}
28
29/// A trait for approximate equality comparisons.
30pub trait ApproxEq: Sized {
31 /// This type type defines a margin within which two values are to be
32 /// considered approximately equal. It must implement `Default` so that
33 /// `approx_eq()` can be called on unknown types.
34 type Margin: FloatMargin;
35
36 /// This method tests that the `self` and `other` values are equal within `margin`
37 /// of each other.
38 fn approx_eq<M: Into<Self::Margin>>(self, other: Self, margin: M) -> bool;
39
40 /// This method tests that the `self` and `other` values are not within `margin`
41 /// of each other.
42 fn approx_ne<M: Into<Self::Margin>>(self, other: Self, margin: M) -> bool {
43 !self.approx_eq(other, margin)
44 }
45}
46
47/// This type defines a margin within two `f32` values might be considered equal,
48/// and is intended as the associated type for the `ApproxEq` trait.
49///
50/// Two tests are used to determine approximate equality.
51///
52/// The first test considers two values approximately equal if they differ by <=
53/// `epsilon`. This will only succeed for very small numbers. Note that it may
54/// succeed even if the parameters are of differing signs, straddling zero.
55///
56/// The second test considers how many ULPs (units of least precision, units in
57/// the last place, which is the integer number of floating-point representations
58/// that the parameters are separated by) different the parameters are and considers
59/// them approximately equal if this is <= `ulps`. For large floating-point numbers,
60/// an ULP can be a rather large gap, but this kind of comparison is necessary
61/// because floating-point operations must round to the nearest representable value
62/// and so larger floating-point values accumulate larger errors.
63#[repr(C)]
64#[derive(Debug, Clone, Copy)]
65pub struct F32Margin {
66 pub epsilon: f32,
67 pub ulps: i32,
68}
69impl Default for F32Margin {
70 #[inline]
71 fn default() -> F32Margin {
72 F32Margin {
73 epsilon: f32::EPSILON,
74 ulps: 4,
75 }
76 }
77}
78impl FloatMargin for F32Margin {
79 type F = f32;
80 type I = i32;
81
82 #[inline]
83 fn zero() -> F32Margin {
84 F32Margin {
85 epsilon: 0.0,
86 ulps: 0,
87 }
88 }
89 fn epsilon(self, epsilon: f32) -> Self {
90 F32Margin { epsilon, ..self }
91 }
92 fn ulps(self, ulps: i32) -> Self {
93 F32Margin { ulps, ..self }
94 }
95}
96impl From<(f32, i32)> for F32Margin {
97 fn from(m: (f32, i32)) -> F32Margin {
98 F32Margin {
99 epsilon: m.0,
100 ulps: m.1,
101 }
102 }
103}
104
105// no-std compatible abs function
106#[inline(always)]
107fn f32abs(x: f32) -> f32 {
108 f32::from_bits(x.to_bits() & !(1 << 31))
109}
110
111impl ApproxEq for f32 {
112 type Margin = F32Margin;
113
114 fn approx_eq<M: Into<Self::Margin>>(self, other: f32, margin: M) -> bool {
115 let margin = margin.into();
116
117 // Check for exact equality first. This is often true, and so we get the
118 // performance benefit of only doing one compare in most cases.
119 self == other || {
120 // Perform epsilon comparison next
121 let eps: f32 = f32abs(self - other);
122 (eps <= margin.epsilon) || {
123 // Perform ulps comparison last
124 let diff: i32 = self.ulps(&other);
125 saturating_abs_i32!(diff) <= margin.ulps
126 }
127 }
128 }
129}
130
131#[test]
132fn f32_approx_eq_test1() {
133 let f: f32 = 0.0_f32;
134 let g: f32 = -0.0000000000000005551115123125783_f32;
135 assert!(f != g); // Should not be directly equal
136 assert!(f.approx_eq(g, (f32::EPSILON, 0)) == true);
137}
138#[test]
139fn f32_approx_eq_test2() {
140 let f: f32 = 0.0_f32;
141 let g: f32 = -0.0_f32;
142 assert!(f.approx_eq(g, (f32::EPSILON, 0)) == true);
143}
144#[test]
145fn f32_approx_eq_test3() {
146 let f: f32 = 0.0_f32;
147 let g: f32 = 0.00000000000000001_f32;
148 assert!(f.approx_eq(g, (f32::EPSILON, 0)) == true);
149}
150#[test]
151fn f32_approx_eq_test4() {
152 let f: f32 = 0.00001_f32;
153 let g: f32 = 0.00000000000000001_f32;
154 assert!(f.approx_eq(g, (f32::EPSILON, 0)) == false);
155}
156#[test]
157fn f32_approx_eq_test5() {
158 let f: f32 = 0.1_f32;
159 let mut sum: f32 = 0.0_f32;
160 for _ in 0_isize..10_isize {
161 sum += f;
162 }
163 let product: f32 = f * 10.0_f32;
164 assert!(sum != product); // Should not be directly equal:
165 assert!(sum.approx_eq(product, (f32::EPSILON, 1)) == true);
166 assert!(sum.approx_eq(product, F32Margin::zero()) == false);
167}
168#[test]
169fn f32_approx_eq_test6() {
170 let x: f32 = 1000000_f32;
171 let y: f32 = 1000000.1_f32;
172 assert!(x != y); // Should not be directly equal
173 assert!(x.approx_eq(y, (0.0, 2)) == true); // 2 ulps does it
174 // epsilon method no good here:
175 assert!(x.approx_eq(y, (1000.0 * f32::EPSILON, 0)) == false);
176}
177
178/// This type defines a margin within two `f64` values might be considered equal,
179/// and is intended as the associated type for the `ApproxEq` trait.
180///
181/// Two tests are used to determine approximate equality.
182///
183/// The first test considers two values approximately equal if they differ by <=
184/// `epsilon`. This will only succeed for very small numbers. Note that it may
185/// succeed even if the parameters are of differing signs, straddling zero.
186///
187/// The second test considers how many ULPs (units of least precision, units in
188/// the last place, which is the integer number of floating-point representations
189/// that the parameters are separated by) different the parameters are and considers
190/// them approximately equal if this is <= `ulps`. For large floating-point numbers,
191/// an ULP can be a rather large gap, but this kind of comparison is necessary
192/// because floating-point operations must round to the nearest representable value
193/// and so larger floating-point values accumulate larger errors.
194#[derive(Debug, Clone, Copy)]
195pub struct F64Margin {
196 pub epsilon: f64,
197 pub ulps: i64,
198}
199impl Default for F64Margin {
200 #[inline]
201 fn default() -> F64Margin {
202 F64Margin {
203 epsilon: f64::EPSILON,
204 ulps: 4,
205 }
206 }
207}
208impl FloatMargin for F64Margin {
209 type F = f64;
210 type I = i64;
211
212 #[inline]
213 fn zero() -> F64Margin {
214 F64Margin {
215 epsilon: 0.0,
216 ulps: 0,
217 }
218 }
219 fn epsilon(self, epsilon: f64) -> Self {
220 F64Margin { epsilon, ..self }
221 }
222 fn ulps(self, ulps: i64) -> Self {
223 F64Margin { ulps, ..self }
224 }
225}
226impl From<(f64, i64)> for F64Margin {
227 fn from(m: (f64, i64)) -> F64Margin {
228 F64Margin {
229 epsilon: m.0,
230 ulps: m.1,
231 }
232 }
233}
234
235// no-std compatible abs function
236#[inline(always)]
237fn f64abs(x: f64) -> f64 {
238 f64::from_bits(x.to_bits() & !(1 << 63))
239}
240
241impl ApproxEq for f64 {
242 type Margin = F64Margin;
243
244 fn approx_eq<M: Into<Self::Margin>>(self, other: f64, margin: M) -> bool {
245 let margin = margin.into();
246
247 // Check for exact equality first. This is often true, and so we get the
248 // performance benefit of only doing one compare in most cases.
249 self == other || {
250 // Perform epsilon comparison next
251 let eps: f64 = f64abs(self - other);
252 (eps <= margin.epsilon) || {
253 // Perform ulps comparison last
254 let diff: i64 = self.ulps(&other);
255 saturating_abs_i64!(diff) <= margin.ulps
256 }
257 }
258 }
259}
260
261#[test]
262fn f64_approx_eq_test1() {
263 let f: f64 = 0.0_f64;
264 let g: f64 = -0.0000000000000005551115123125783_f64;
265 assert!(f != g); // Should not be precisely equal.
266 assert!(f.approx_eq(g, (3.0 * f64::EPSILON, 0)) == true); // 3e is enough.
267 // ULPs test won't ever call these equal.
268}
269#[test]
270fn f64_approx_eq_test2() {
271 let f: f64 = 0.0_f64;
272 let g: f64 = -0.0_f64;
273 assert!(f.approx_eq(g, (f64::EPSILON, 0)) == true);
274}
275#[test]
276fn f64_approx_eq_test3() {
277 let f: f64 = 0.0_f64;
278 let g: f64 = 1e-17_f64;
279 assert!(f.approx_eq(g, (f64::EPSILON, 0)) == true);
280}
281#[test]
282fn f64_approx_eq_test4() {
283 let f: f64 = 0.00001_f64;
284 let g: f64 = 0.00000000000000001_f64;
285 assert!(f.approx_eq(g, (f64::EPSILON, 0)) == false);
286}
287#[test]
288fn f64_approx_eq_test5() {
289 let f: f64 = 0.1_f64;
290 let mut sum: f64 = 0.0_f64;
291 for _ in 0_isize..10_isize {
292 sum += f;
293 }
294 let product: f64 = f * 10.0_f64;
295 assert!(sum != product); // Should not be precisely equally.
296 assert!(sum.approx_eq(product, (f64::EPSILON, 0)) == true);
297 assert!(sum.approx_eq(product, (0.0, 1)) == true);
298}
299#[test]
300fn f64_approx_eq_test6() {
301 let x: f64 = 1000000_f64;
302 let y: f64 = 1000000.0000000003_f64;
303 assert!(x != y); // Should not be precisely equal.
304 assert!(x.approx_eq(y, (0.0, 3)) == true);
305}
306#[test]
307fn f64_code_triggering_issue_20() {
308 assert_eq!((-25.0f64).approx_eq(25.0, (0.00390625, 1)), false);
309}
310
311impl<T> ApproxEq for &[T]
312where
313 T: Copy + ApproxEq,
314{
315 type Margin = <T as ApproxEq>::Margin;
316
317 fn approx_eq<M: Into<Self::Margin>>(self, other: Self, margin: M) -> bool {
318 let margin = margin.into();
319 if self.len() != other.len() {
320 return false;
321 }
322 self.iter()
323 .zip(other.iter())
324 .all(|(a: &T, b: &T)| a.approx_eq(*b, margin))
325 }
326}
327
328#[test]
329fn test_slices() {
330 assert!([1.33, 2.4, 2.5].approx_eq(&[1.33, 2.4, 2.5], (0.0, 0_i64)));
331 assert!(![1.33, 2.4, 2.6].approx_eq(&[1.33, 2.4, 2.5], (0.0, 0_i64)));
332 assert!(![1.33, 2.4].approx_eq(&[1.33, 2.4, 2.5], (0.0, 0_i64)));
333 assert!(![1.33, 2.4, 2.5].approx_eq(&[1.33, 2.4], (0.0, 0_i64)));
334}
335
336impl<T> ApproxEq for Option<T>
337where
338 T: Copy + ApproxEq,
339{
340 type Margin = <T as ApproxEq>::Margin;
341
342 fn approx_eq<M: Into<Self::Margin>>(self, other: Self, margin: M) -> bool {
343 let margin = margin.into();
344 match (self, other) {
345 (None, None) => true,
346 (Some(slf: T), Some(oth: T)) => slf.approx_eq(other:oth, margin),
347 _ => false,
348 }
349 }
350}
351
352#[test]
353fn test_option() {
354 let x: Option<f32> = None;
355 assert!(x.approx_eq(None, (0.0, 0_i32)));
356 assert!(Some(5.3_f32).approx_eq(Some(5.3), (0.0, 0_i32)));
357 assert!(Some(5.3_f32).approx_ne(Some(5.7), (0.0, 0_i32)));
358 assert!(Some(5.3_f32).approx_ne(None, (0.0, 0_i32)));
359 assert!(x.approx_ne(Some(5.3), (0.0, 0_i32)));
360}
361