float_cmp/
ulps.rs

Help
1// Copyright 2014-2020 Optimal Computing (NZ) Ltd.
2// Licensed under the MIT license.  See LICENSE for details.
3
4#[cfg(feature = "num-traits")]
5use num_traits::NumCast;
6
7#[inline]
8fn f32_ordered_bits(f: f32) -> u32 {
9    const SIGN_BIT: u32 = 1 << 31;
10    let bits = f.to_bits();
11    if bits & SIGN_BIT != 0 {
12        !bits
13    } else {
14        bits ^ SIGN_BIT
15    }
16}
17
18#[inline]
19fn f64_ordered_bits(f: f64) -> u64 {
20    const SIGN_BIT: u64 = 1 << 63;
21    let bits = f.to_bits();
22    if bits & SIGN_BIT != 0 {
23        !bits
24    } else {
25        bits ^ SIGN_BIT
26    }
27}
28
29/// A trait for floating point numbers which computes the number of representable
30/// values or ULPs (Units of Least Precision) that separate the two given values.
31#[cfg(feature = "num-traits")]
32pub trait Ulps {
33    type U: Copy + NumCast;
34
35    /// The number of representable values or ULPs (Units of Least Precision) that
36    /// separate `self` and `other`.  The result `U` is an integral value, and will
37    /// be zero if `self` and `other` are exactly equal.
38    fn ulps(&self, other: &Self) -> <Self as Ulps>::U;
39
40    /// The next representable number above this one
41    fn next(&self) -> Self;
42
43    /// The previous representable number below this one
44    fn prev(&self) -> Self;
45}
46
47#[cfg(not(feature = "num-traits"))]
48pub trait Ulps {
49    type U: Copy;
50
51    /// The number of representable values or ULPs (Units of Least Precision) that
52    /// separate `self` and `other`.  The result `U` is an integral value, and will
53    /// be zero if `self` and `other` are exactly equal.
54    fn ulps(&self, other: &Self) -> <Self as Ulps>::U;
55
56    /// The next representable number above this one
57    fn next(&self) -> Self;
58
59    /// The previous representable number below this one
60    fn prev(&self) -> Self;
61}
62
63impl Ulps for f32 {
64    type U = i32;
65
66    fn ulps(&self, other: &f32) -> i32 {
67        // IEEE754 defined floating point storage representation to
68        // maintain their order when their bit patterns are interpreted as
69        // integers.  This is a huge boon to the task at hand, as we can
70        // reinterpret them as integers to find out how many ULPs apart any
71        // two floats are
72
73        // Setup integer representations of the input
74        let ai32: i32 = f32_ordered_bits(*self) as i32;
75        let bi32: i32 = f32_ordered_bits(*other) as i32;
76
77        ai32.wrapping_sub(bi32)
78    }
79
80    fn next(&self) -> Self {
81        if self.is_infinite() && *self > 0.0 {
82            *self
83        } else if *self == -0.0 && self.is_sign_negative() {
84            0.0
85        } else {
86            let mut u = self.to_bits();
87            if *self >= 0.0 {
88                u += 1;
89            } else {
90                u -= 1;
91            }
92            f32::from_bits(u)
93        }
94    }
95
96    fn prev(&self) -> Self {
97        if self.is_infinite() && *self < 0.0 {
98            *self
99        } else if *self == 0.0 && self.is_sign_positive() {
100            -0.0
101        } else {
102            let mut u = self.to_bits();
103            if *self <= -0.0 {
104                u += 1;
105            } else {
106                u -= 1;
107            }
108            f32::from_bits(u)
109        }
110    }
111}
112
113#[test]
114fn f32_ulps_test1() {
115    let x: f32 = 1000000_f32;
116    let y: f32 = 1000000.1_f32;
117    assert!(x.ulps(&y) == -2);
118}
119
120#[test]
121fn f32_ulps_test2() {
122    let pzero: f32 = f32::from_bits(0x00000000_u32);
123    let nzero: f32 = f32::from_bits(0x80000000_u32);
124    assert_eq!(pzero.ulps(&nzero), 1);
125}
126#[test]
127fn f32_ulps_test3() {
128    let pinf: f32 = f32::from_bits(0x7f800000_u32);
129    let ninf: f32 = f32::from_bits(0xff800000_u32);
130    assert_eq!(pinf.ulps(&ninf), -16777215);
131}
132
133#[test]
134fn f32_ulps_test4() {
135    let x: f32 = f32::from_bits(0x63a7f026_u32);
136    let y: f32 = f32::from_bits(0x63a7f023_u32);
137    assert!(x.ulps(&y) == 3);
138}
139
140#[test]
141fn f32_ulps_test5() {
142    let x: f32 = 2.0;
143    let ulps: i32 = x.to_bits() as i32;
144    let x2: f32 = <f32>::from_bits(ulps as u32);
145    assert_eq!(x, x2);
146}
147
148#[test]
149fn f32_ulps_test6() {
150    let negzero: f32 = -0.;
151    let zero: f32 = 0.;
152    assert_eq!(negzero.next(), zero);
153    assert_eq!(zero.prev(), negzero);
154    assert!(negzero.prev() < 0.0);
155    assert!(zero.next() > 0.0);
156}
157
158impl Ulps for f64 {
159    type U = i64;
160
161    fn ulps(&self, other: &f64) -> i64 {
162        // IEEE754 defined floating point storage representation to
163        // maintain their order when their bit patterns are interpreted as
164        // integers.  This is a huge boon to the task at hand, as we can
165        // reinterpret them as integers to find out how many ULPs apart any
166        // two floats are
167
168        // Setup integer representations of the input
169        let ai64: i64 = f64_ordered_bits(*self) as i64;
170        let bi64: i64 = f64_ordered_bits(*other) as i64;
171
172        ai64.wrapping_sub(bi64)
173    }
174
175    fn next(&self) -> Self {
176        if self.is_infinite() && *self > 0.0 {
177            *self
178        } else if *self == -0.0 && self.is_sign_negative() {
179            0.0
180        } else {
181            let mut u = self.to_bits();
182            if *self >= 0.0 {
183                u += 1;
184            } else {
185                u -= 1;
186            }
187            f64::from_bits(u)
188        }
189    }
190
191    fn prev(&self) -> Self {
192        if self.is_infinite() && *self < 0.0 {
193            *self
194        } else if *self == 0.0 && self.is_sign_positive() {
195            -0.0
196        } else {
197            let mut u = self.to_bits();
198            if *self <= -0.0 {
199                u += 1;
200            } else {
201                u -= 1;
202            }
203            f64::from_bits(u)
204        }
205    }
206}
207
208#[test]
209fn f64_ulps_test1() {
210    let x: f64 = 1000000_f64;
211    let y: f64 = 1000000.00000001_f64;
212    assert!(x.ulps(&y) == -86);
213}
214
215#[test]
216fn f64_ulps_test2() {
217    let pzero: f64 = f64::from_bits(0x0000000000000000_u64);
218    let nzero: f64 = f64::from_bits(0x8000000000000000_u64);
219    assert_eq!(pzero.ulps(&nzero), 1);
220}
221#[test]
222fn f64_ulps_test3() {
223    let pinf: f64 = f64::from_bits(0x7f80000000000000_u64);
224    let ninf: f64 = f64::from_bits(0xff80000000000000_u64);
225    assert_eq!(pinf.ulps(&ninf), -72057594037927935);
226}
227
228#[test]
229fn f64_ulps_test4() {
230    let x: f64 = f64::from_bits(0xd017f6cc63a7f026_u64);
231    let y: f64 = f64::from_bits(0xd017f6cc63a7f023_u64);
232    assert!(x.ulps(&y) == -3);
233}
234
235#[test]
236fn f64_ulps_test5() {
237    let x: f64 = 2.0;
238    let ulps: i64 = x.to_bits() as i64;
239    let x2: f64 = <f64>::from_bits(ulps as u64);
240    assert_eq!(x, x2);
241}
242
243#[test]
244fn f64_ulps_test6() {
245    let negzero: f64 = -0.;
246    let zero: f64 = 0.;
247    assert_eq!(negzero.next(), zero);
248    assert_eq!(zero.prev(), negzero);
249    assert!(negzero.prev() < 0.0);
250    assert!(zero.next() > 0.0);
251}
float_cmp/ulps.rs

float_cmp/
ulps.rs
