1 | // Copyright 2014-2020 Optimal Computing (NZ) Ltd. |
---|---|

2 | // Licensed under the MIT license. See LICENSE for details. |

3 | |

4 | #[cfg(feature= "num_traits")] |

5 | use 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")] |

10 | pub 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"))] |

26 | pub 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 | |

41 | impl 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] |

93 | fn 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] |

100 | fn 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] |

106 | fn 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] |

113 | fn 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] |

120 | fn 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] |

128 | fn 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 | |

137 | impl 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] |

189 | fn 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] |

196 | fn 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] |

202 | fn 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] |

209 | fn 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] |

216 | fn 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] |

224 | fn 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 |