Oregami
Repositories/oxedyne/fe2o3

oxedyne/fe2o3/fe2o3_num/src/float.rs

11.7 KiB, 31 runs

created by r1870400018:631, which is this file's identity for as long as the history lasts, whatever it is later renamed to

download · who wrote it · its history

1use oxedyne_fe2o3_core::prelude::*;
2
3use std::{
4 hash::{
5 Hash,
6 Hasher,
7 },
8 string,
9};
10
11pub trait PrimitiveFloat: Sized + string::ToString {}
12
13impl PrimitiveFloat for f32 {}
14impl PrimitiveFloat for f64 {}
15
16pub fn round_to_sf(n: f64, sf: u8) -> f64 {
17 if sf == 0 || n == 0.0 || !n.is_finite() {
18 return n;
19 }
20 // Decimal place of the leading digit: 0 for one up to ten, -1 for a tenth
21 // up to one, and so on. The logarithm of a negative number is not a number,
22 // so the magnitude is taken first; the earlier form omitted that, and the
23 // resulting cast of a not-a-number to zero put every negative value's
24 // leading digit in the units place.
25 let mut mag = n.abs().log10().floor() as i32;
26 // A logarithm is not exact, so the place is checked against the value it is
27 // meant to describe and nudged if it names the wrong decade.
28 let lead = mul_pow10(n.abs(), -mag);
29 if lead >= 10.0 {
30 mag += 1;
31 } else if lead < 1.0 {
32 mag -= 1;
33 }
34 // Places the decimal point moves right to leave sf digits before it.
35 let shift = (sf as i32) - 1 - mag;
36 let scaled = mul_pow10(n, shift);
37 if !scaled.is_finite() {
38 return n;
39 }
40 let rnd = scaled.round();
41 // Past the range in which a power of ten is exactly representable the
42 // scaling itself carries an error of a few units in the last place. A value
43 // already at the requested precision must not be nudged by that error.
44 if shift.abs() > 22 && (scaled - rnd).abs() <= scaled.abs() * 16.0 * f64::EPSILON {
45 return n;
46 }
47 let out = mul_pow10(rnd, -shift);
48 if out.is_finite() { out } else { n }
49}
50
51pub fn mul_pow10(n: f64, exp: i32) -> f64 {
52 if exp.abs() <= 300 {
53 let p = 10.0f64.powi(exp.abs());
54 if exp >= 0 { n * p } else { n / p }
55 } else {
56 let half = exp / 2;
57 mul_pow10(mul_pow10(n, half), exp - half)
58 }
59}
60
61new_type!(Float32, f32, Clone, Debug, Default, PartialOrd);
62
63impl Ord for Float32 {
64 // total_cmp function currently yet to make it to stable
65 fn cmp(&self, other: &Self) -> std::cmp::Ordering {
66 let mut left = self.to_bits() as i32;
67 let mut right = other.to_bits() as i32;
68 left ^= (((left >> 31) as u32) >> 1) as i32;
69 right ^= (((right >> 31) as u32) >> 1) as i32;
70 left.cmp(&right)
71 }
72}
73
74impl PartialEq for Float32 {
75 fn eq(&self, other: &Float32) -> bool {
76 self.cmp(other) == std::cmp::Ordering::Equal
77 }
78 fn ne(&self, other: &Float32) -> bool {
79 self.cmp(other) != std::cmp::Ordering::Equal
80 }
81}
82
83impl Eq for Float32 {}
84
85impl Hash for Float32 {
86 fn hash<H: Hasher>(&self, state: &mut H) {
87 let (m, e, s) = self.integer_decode();
88 m.hash(state);
89 e.hash(state);
90 s.hash(state);
91 }
92}
93
94impl Float32 {
95 // A function deprecated from the std library, modified to take a reference and use the inner type.
96 // https://github.com/rust-lang/rust/blob/5c674a11471ec0569f616854d715941757a48a0a/src/libcore/num/f32.rs
97 fn integer_decode(&self) -> (u64, i16, i8) {
98 let bits: u32 = self.0.to_bits();
99 let sign: i8 = if bits >> 31 == 0 { 1 } else { -1 };
100 let mut exponent: i16 = ((bits >> 23) & 0xff) as i16;
101 let mantissa = if exponent == 0 {
102 (bits & 0x7fffff) << 1
103 } else {
104 (bits & 0x7fffff) | 0x800000
105 };
106 // Exponent bias + mantissa shift
107 exponent -= 127 + 23;
108 (mantissa as u64, exponent, sign)
109 }
110
111 pub fn is_zero(&self) -> bool {
112 let (m, _, _) = self.integer_decode();
113 m == 0
114 }
115}
116
117new_type!(Float64, f64, Clone, Debug, Default, PartialOrd);
118
119impl Ord for Float64 {
120 // total_cmp function currently yet to make it to stable
121 fn cmp(&self, other: &Self) -> std::cmp::Ordering {
122 let mut left = self.to_bits() as i64;
123 let mut right = other.to_bits() as i64;
124 left ^= (((left >> 63) as u64) >> 1) as i64;
125 right ^= (((right >> 63) as u64) >> 1) as i64;
126 left.cmp(&right)
127 }
128}
129
130impl PartialEq for Float64 {
131 fn eq(&self, other: &Float64) -> bool {
132 self.cmp(other) == std::cmp::Ordering::Equal
133 }
134 fn ne(&self, other: &Float64) -> bool {
135 self.cmp(other) != std::cmp::Ordering::Equal
136 }
137}
138
139impl Eq for Float64 {}
140
141impl Hash for Float64 {
142 fn hash<H: Hasher>(&self, state: &mut H) {
143 let (m, e, s) = self.integer_decode();
144 m.hash(state);
145 e.hash(state);
146 s.hash(state);
147 }
148}
149
150impl Float64 {
151 // A function deprecated from the std library, modified to take a reference and use the inner type.
152 // https://github.com/rust-lang/rust/blob/5c674a11471ec0569f616854d715941757a48a0a/src/libcore/num/f64.rs
153 fn integer_decode(&self) -> (u64, i16, i8) {
154 let bits: u64 = self.0.to_bits();
155 let sign: i8 = if bits >> 63 == 0 { 1 } else { -1 };
156 let mut exponent: i16 = ((bits >> 52) & 0x7ff) as i16;
157 let mantissa = if exponent == 0 {
158 (bits & 0xfffffffffffff) << 1
159 } else {
160 (bits & 0xfffffffffffff) | 0x10000000000000
161 };
162 // Exponent bias + mantissa shift
163 exponent -= 1023 + 52;
164 (mantissa, exponent, sign)
165 }
166
167 pub fn is_zero(&self) -> bool {
168 let (m, _, _) = self.integer_decode();
169 m == 0
170 }
171}
172
173#[cfg(test)]
174mod round_to_sf_tests {
175 use super::*;
176
177
178 // Values of one and above, which the function has always handled. These
179 // pin the behaviour that must not change. //
180
181 #[test]
182 fn test_round_to_sf_above_one_01() {
183 // 1234 is 1.234 x 10^3; three figures keep 1.23, so 1230.
184 assert_eq!(round_to_sf(1234.0, 3), 1230.0);
185 // 1236 is 1.236 x 10^3; the fourth digit is 6, so the third rounds up.
186 assert_eq!(round_to_sf(1236.0, 3), 1240.0);
187 // 98765 is 9.8765 x 10^4; two figures keep 9.9, so 99000.
188 assert_eq!(round_to_sf(98765.0, 2), 99000.0);
189 // 9.99 to two figures carries into a new decade: 10.
190 assert_eq!(round_to_sf(9.99, 2), 10.0);
191 assert_eq!(round_to_sf(1.0, 3), 1.0);
192 assert_eq!(round_to_sf(1.0e6, 3), 1.0e6);
193 }
194
195 // The same values negated. A magnitude does not depend on a sign, so each
196 // expectation is the mirror of the one above. //
197
198 #[test]
199 fn test_round_to_sf_negative_above_one_01() {
200 assert_eq!(round_to_sf(-1234.0, 3), -1230.0);
201 assert_eq!(round_to_sf(-1236.0, 3), -1240.0);
202 assert_eq!(round_to_sf(-98765.0, 2), -99000.0);
203 assert_eq!(round_to_sf(-9.99, 2), -10.0);
204 assert_eq!(round_to_sf(-1.0, 3), -1.0);
205 assert_eq!(round_to_sf(-1.0e6, 3), -1.0e6);
206 }
207
208 // Values between zero and one, where the leading digit sits to the right
209 // of the point. //
210
211 #[test]
212 fn test_round_to_sf_below_one_01() {
213 // 0.4749 is 4.749 x 10^-1; two figures keep 4.7, so 0.47.
214 assert_eq!(round_to_sf(0.4749, 2), 0.47);
215 // 0.4751 is 4.751 x 10^-1; the third digit is 5 with more behind it.
216 assert_eq!(round_to_sf(0.4751, 2), 0.48);
217 // 0.05512 is 5.512 x 10^-2; two figures keep 5.5, so 0.055.
218 assert_eq!(round_to_sf(0.05512, 2), 0.055);
219 // 0.9994 is 9.994 x 10^-1; three figures keep 9.99, so 0.999.
220 assert_eq!(round_to_sf(0.9994, 3), 0.999);
221 // 0.9996 rounds up through the decade to 1.00.
222 assert_eq!(round_to_sf(0.9996, 3), 1.0);
223 // 0.0999 is 9.99 x 10^-2; two figures carry to 1.0 x 10^-1.
224 assert_eq!(round_to_sf(0.0999, 2), 0.1);
225 }
226
227 #[test]
228 fn test_round_to_sf_negative_below_one_01() {
229 assert_eq!(round_to_sf(-0.4749, 2), -0.47);
230 assert_eq!(round_to_sf(-0.4751, 2), -0.48);
231 assert_eq!(round_to_sf(-0.05512, 2), -0.055);
232 assert_eq!(round_to_sf(-0.9994, 3), -0.999);
233 assert_eq!(round_to_sf(-0.9996, 3), -1.0);
234 assert_eq!(round_to_sf(-0.0999, 2), -0.1);
235 }
236
237 #[test]
238 fn test_round_to_sf_reported_slopes_01() {
239 assert_eq!(round_to_sf(-0.475, 2), -0.48);
240 assert_eq!(round_to_sf(-0.5496, 2), -0.55);
241 assert_eq!(round_to_sf(-0.055, 2), -0.055);
242 }
243
244 #[test]
245 fn test_round_to_sf_ties_go_away_from_zero_01() {
246 assert_eq!(round_to_sf(0.125, 2), 0.13);
247 assert_eq!(round_to_sf(-0.125, 2), -0.13);
248 assert_eq!(round_to_sf(2.5, 1), 3.0);
249 assert_eq!(round_to_sf(-2.5, 1), -3.0);
250 assert_eq!(round_to_sf(1.5, 1), 2.0);
251 assert_eq!(round_to_sf(-1.5, 1), -2.0);
252 assert_eq!(round_to_sf(0.25, 1), 0.3);
253 assert_eq!(round_to_sf(-0.25, 1), -0.3);
254 }
255
256 #[test]
257 fn test_round_to_sf_powers_of_ten_01() {
258 for exp in -320i32..=308 {
259 let v = match format!("1e{}", exp).parse::<f64>() {
260 Ok(v) => v,
261 Err(_) => continue,
262 };
263 for sf in 1..=6u8 {
264 assert_eq!(round_to_sf(v, sf), v, "10^{} at {} sf", exp, sf);
265 assert_eq!(round_to_sf(-v, sf), -v, "-10^{} at {} sf", exp, sf);
266 }
267 }
268 }
269
270 #[test]
271 fn test_round_to_sf_sign_symmetry_01() {
272 let mut v = 3.0e-7;
273 for _ in 0..2000 {
274 for sf in 1..=6u8 {
275 assert_eq!(round_to_sf(-v, sf), -round_to_sf(v, sf), "{:e} at {} sf", v, sf);
276 }
277 v *= 1.017;
278 }
279 }
280
281 #[test]
282 fn test_round_to_sf_degenerate_01() {
283 assert_eq!(round_to_sf(0.0, 3), 0.0);
284 assert_eq!(round_to_sf(-0.0, 3), 0.0);
285 assert_eq!(round_to_sf(12.34, 0), 12.34);
286 assert!(round_to_sf(f64::NAN, 3).is_nan());
287 assert_eq!(round_to_sf(f64::INFINITY, 3), f64::INFINITY);
288 assert_eq!(round_to_sf(f64::NEG_INFINITY, 3), f64::NEG_INFINITY);
289 }
290
291 #[test]
292 fn test_round_to_sf_extremes_01() {
293 assert_eq!(round_to_sf(1.0e300, 3), 1.0e300);
294 assert_eq!(round_to_sf(-1.0e300, 3), -1.0e300);
295 assert_eq!(round_to_sf(1.0e-23, 3), 1.0e-23);
296 assert!(round_to_sf(f64::MAX, 3).is_finite());
297 assert!(round_to_sf(1.0e-320, 3) > 0.0);
298 }
299
300 #[test]
301 fn test_round_to_sf_across_decades_01() {
302 for exp in -300i32..=300 {
303 let v = match format!("1.2345e{}", exp).parse::<f64>() {
304 Ok(v) => v,
305 Err(_) => continue,
306 };
307 let want = match format!("1.23e{}", exp).parse::<f64>() {
308 Ok(w) => w,
309 Err(_) => continue,
310 };
311 let got = round_to_sf(v, 3);
312 let rel = ((got - want) / want).abs();
313 assert!(rel < 1.0e-15, "1.2345e{} gave {:e}, wanted {:e}", exp, got, want);
314 }
315 }
316}
317
318#[cfg(test)]
319mod mul_pow10_tests {
320 use super::*;
321
322 #[test]
323 fn test_mul_pow10_exact_range_01() {
324 for exp in -22i32..=22 {
325 let want = match format!("1e{}", exp).parse::<f64>() {
326 Ok(v) => v,
327 Err(_) => continue,
328 };
329 assert_eq!(mul_pow10(1.0, exp), want, "10^{}", exp);
330 assert_eq!(mul_pow10(-1.0, exp), -want, "-10^{}", exp);
331 }
332 assert_eq!(mul_pow10(1.234, 0), 1.234);
333 }
334
335 #[test]
336 fn test_mul_pow10_beyond_a_single_power_01() {
337 assert_eq!(1.0e-300 * 10.0f64.powi(320), f64::INFINITY);
338 let got = mul_pow10(1.0e-300, 320);
339 assert!(((got - 1.0e20) / 1.0e20).abs() < 1.0e-15, "got {:e}", got);
340 let got = mul_pow10(1.0e300, -320);
341 assert!(((got - 1.0e-20) / 1.0e-20).abs() < 1.0e-15, "got {:e}", got);
342 }
343
344 #[test]
345 fn test_mul_pow10_round_trip_01() {
346 let v = 1.234567;
347 for exp in -22i32..=22 {
348 assert_eq!(mul_pow10(mul_pow10(v, exp), -exp), v, "10^{}", exp);
349 }
350 for exp in -300i32..=300 {
351 let got = mul_pow10(mul_pow10(v, exp), -exp);
352 assert!(((got - v) / v).abs() < 4.0 * f64::EPSILON, "10^{} gave {}", exp, got);
353 }
354 }
355}