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: u32 = 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: u64 = 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}
252