| // Copyright 2014-2020 Optimal Computing (NZ) Ltd. |
| // Licensed under the MIT license. See LICENSE for details. |
| |
| #[cfg(feature="num_traits")] |
| use num_traits::NumCast; |
| |
| /// A trait for floating point numbers which computes the number of representable |
| /// values or ULPs (Units of Least Precision) that separate the two given values. |
| #[cfg(feature="num_traits")] |
| pub trait Ulps { |
| type U: Copy + NumCast; |
| |
| /// The number of representable values or ULPs (Units of Least Precision) that |
| /// separate `self` and `other`. The result `U` is an integral value, and will |
| /// be zero if `self` and `other` are exactly equal. |
| fn ulps(&self, other: &Self) -> <Self as Ulps>::U; |
| |
| /// The next representable number above this one |
| fn next(&self) -> Self; |
| |
| /// The previous representable number below this one |
| fn prev(&self) -> Self; |
| } |
| |
| #[cfg(not(feature="num_traits"))] |
| pub trait Ulps { |
| type U: Copy; |
| |
| /// The number of representable values or ULPs (Units of Least Precision) that |
| /// separate `self` and `other`. The result `U` is an integral value, and will |
| /// be zero if `self` and `other` are exactly equal. |
| fn ulps(&self, other: &Self) -> <Self as Ulps>::U; |
| |
| /// The next representable number above this one |
| fn next(&self) -> Self; |
| |
| /// The previous representable number below this one |
| fn prev(&self) -> Self; |
| } |
| |
| impl Ulps for f32 { |
| type U = i32; |
| |
| fn ulps(&self, other: &f32) -> i32 { |
| |
| // IEEE754 defined floating point storage representation to |
| // maintain their order when their bit patterns are interpreted as |
| // integers. This is a huge boon to the task at hand, as we can |
| // reinterpret them as integers to find out how many ULPs apart any |
| // two floats are |
| |
| // Setup integer representations of the input |
| let ai32: i32 = self.to_bits() as i32; |
| let bi32: i32 = other.to_bits() as i32; |
| |
| ai32.wrapping_sub(bi32) |
| } |
| |
| fn next(&self) -> Self { |
| if self.is_infinite() && *self > 0.0 { |
| *self |
| } else if *self == -0.0 && self.is_sign_negative() { |
| 0.0 |
| } else { |
| let mut u = self.to_bits(); |
| if *self >= 0.0 { |
| u += 1; |
| } else { |
| u -= 1; |
| } |
| f32::from_bits(u) |
| } |
| } |
| |
| fn prev(&self) -> Self { |
| if self.is_infinite() && *self < 0.0 { |
| *self |
| } else if *self == 0.0 && self.is_sign_positive() { |
| -0.0 |
| } else { |
| let mut u = self.to_bits(); |
| if *self <= -0.0 { |
| u += 1; |
| } else { |
| u -= 1; |
| } |
| f32::from_bits(u) |
| } |
| } |
| } |
| |
| #[test] |
| fn f32_ulps_test1() { |
| let x: f32 = 1000000_f32; |
| let y: f32 = 1000000.1_f32; |
| assert!(x.ulps(&y) == -2); |
| } |
| |
| #[test] |
| fn f32_ulps_test2() { |
| let pzero: f32 = f32::from_bits(0x00000000_u32); |
| let nzero: f32 = f32::from_bits(0x80000000_u32); |
| assert!(pzero.ulps(&nzero) == -2147483648); |
| } |
| #[test] |
| fn f32_ulps_test3() { |
| let pinf: f32 = f32::from_bits(0x7f800000_u32); |
| let ninf: f32 = f32::from_bits(0xff800000_u32); |
| assert!(pinf.ulps(&ninf) == -2147483648); |
| } |
| |
| #[test] |
| fn f32_ulps_test4() { |
| let x: f32 = f32::from_bits(0x63a7f026_u32); |
| let y: f32 = f32::from_bits(0x63a7f023_u32); |
| assert!(x.ulps(&y) == 3); |
| } |
| |
| #[test] |
| fn f32_ulps_test5() { |
| let x: f32 = 2.0; |
| let ulps: i32 = x.to_bits() as i32; |
| let x2: f32 = <f32>::from_bits(ulps as u32); |
| assert_eq!(x, x2); |
| } |
| |
| #[test] |
| fn f32_ulps_test6() { |
| let negzero: f32 = -0.; |
| let zero: f32 = 0.; |
| assert_eq!(negzero.next(), zero); |
| assert_eq!(zero.prev(), negzero); |
| assert!(negzero.prev() < 0.0); |
| assert!(zero.next() > 0.0); |
| } |
| |
| impl Ulps for f64 { |
| type U = i64; |
| |
| fn ulps(&self, other: &f64) -> i64 { |
| |
| // IEEE754 defined floating point storage representation to |
| // maintain their order when their bit patterns are interpreted as |
| // integers. This is a huge boon to the task at hand, as we can |
| // reinterpret them as integers to find out how many ULPs apart any |
| // two floats are |
| |
| // Setup integer representations of the input |
| let ai64: i64 = self.to_bits() as i64; |
| let bi64: i64 = other.to_bits() as i64; |
| |
| ai64.wrapping_sub(bi64) |
| } |
| |
| fn next(&self) -> Self { |
| if self.is_infinite() && *self > 0.0 { |
| *self |
| } else if *self == -0.0 && self.is_sign_negative() { |
| 0.0 |
| } else { |
| let mut u = self.to_bits(); |
| if *self >= 0.0 { |
| u += 1; |
| } else { |
| u -= 1; |
| } |
| f64::from_bits(u) |
| } |
| } |
| |
| fn prev(&self) -> Self { |
| if self.is_infinite() && *self < 0.0 { |
| *self |
| } else if *self == 0.0 && self.is_sign_positive() { |
| -0.0 |
| } else { |
| let mut u = self.to_bits(); |
| if *self <= -0.0 { |
| u += 1; |
| } else { |
| u -= 1; |
| } |
| f64::from_bits(u) |
| } |
| } |
| } |
| |
| #[test] |
| fn f64_ulps_test1() { |
| let x: f64 = 1000000_f64; |
| let y: f64 = 1000000.00000001_f64; |
| assert!(x.ulps(&y) == -86); |
| } |
| |
| #[test] |
| fn f64_ulps_test2() { |
| let pzero: f64 = f64::from_bits(0x0000000000000000_u64); |
| let nzero: f64 = f64::from_bits(0x8000000000000000_u64); |
| assert!(pzero.ulps(&nzero) == -9223372036854775808i64); |
| } |
| #[test] |
| fn f64_ulps_test3() { |
| let pinf: f64 = f64::from_bits(0x7f80000000000000_u64); |
| let ninf: f64 = f64::from_bits(0xff80000000000000_u64); |
| assert!(pinf.ulps(&ninf) == -9223372036854775808i64); |
| } |
| |
| #[test] |
| fn f64_ulps_test4() { |
| let x: f64 = f64::from_bits(0xd017f6cc63a7f026_u64); |
| let y: f64 = f64::from_bits(0xd017f6cc63a7f023_u64); |
| assert!(x.ulps(&y) == 3); |
| } |
| |
| #[test] |
| fn f64_ulps_test5() { |
| let x: f64 = 2.0; |
| let ulps: i64 = x.to_bits() as i64; |
| let x2: f64 = <f64>::from_bits(ulps as u64); |
| assert_eq!(x, x2); |
| } |
| |
| #[test] |
| fn f64_ulps_test6() { |
| let negzero: f64 = -0.; |
| let zero: f64 = 0.; |
| assert_eq!(negzero.next(), zero); |
| assert_eq!(zero.prev(), negzero); |
| assert!(negzero.prev() < 0.0); |
| assert!(zero.next() > 0.0); |
| } |
| |