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