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

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

3 | |

4 | use core::{f32, f64}; |

5 | #[cfg(feature = "num-traits")] |

6 | #[allow(unused_imports)] |

7 | use num_traits::float::FloatCore; |

8 | use super::Ulps; |

9 | |

10 | /// A trait for approximate equality comparisons. |

11 | pub trait ApproxEq: Sized { |

12 | /// This type type defines a margin within which two values are to be |

13 | /// considered approximately equal. It must implement `Default` so that |

14 | /// `approx_eq()` can be called on unknown types. |

15 | type Margin: Copy + Default; |

16 | |

17 | /// This method tests that the `self` and `other` values are equal within `margin` |

18 | /// of each other. |

19 | fn approx_eq<M: Into<Self::Margin>>(self, other: Self, margin: M) -> bool; |

20 | |

21 | /// This method tests that the `self` and `other` values are not within `margin` |

22 | /// of each other. |

23 | fn approx_ne<M: Into<Self::Margin>>(self, other: Self, margin: M) -> bool { |

24 | !self.approx_eq(other, margin) |

25 | } |

26 | } |

27 | |

28 | /// This type defines a margin within two `f32` values might be considered equal, |

29 | /// and is intended as the associated type for the `ApproxEq` trait. |

30 | /// |

31 | /// Two tests are used to determine approximate equality. |

32 | /// |

33 | /// The first test considers two values approximately equal if they differ by <= |

34 | /// `epsilon`. This will only succeed for very small numbers. Note that it may |

35 | /// succeed even if the parameters are of differing signs, straddling zero. |

36 | /// |

37 | /// The second test considers how many ULPs (units of least precision, units in |

38 | /// the last place, which is the integer number of floating-point representations |

39 | /// that the parameters are separated by) different the parameters are and considers |

40 | /// them approximately equal if this is <= `ulps`. For large floating-point numbers, |

41 | /// an ULP can be a rather large gap, but this kind of comparison is necessary |

42 | /// because floating-point operations must round to the nearest representable value |

43 | /// and so larger floating-point values accumulate larger errors. |

44 | #[repr(C)] |

45 | #[derive(Debug, Clone, Copy)] |

46 | pub struct F32Margin { |

47 | pub epsilon: f32, |

48 | pub ulps: i32 |

49 | } |

50 | impl Default for F32Margin { |

51 | #[inline] |

52 | fn default() -> F32Margin { |

53 | F32Margin { |

54 | epsilon: f32::EPSILON, |

55 | ulps: 4 |

56 | } |

57 | } |

58 | } |

59 | impl F32Margin { |

60 | #[inline] |

61 | pub fn zero() -> F32Margin { |

62 | F32Margin { |

63 | epsilon: 0.0, |

64 | ulps: 0 |

65 | } |

66 | } |

67 | pub fn epsilon(self, epsilon: f32) -> Self { |

68 | F32Margin { |

69 | epsilon, |

70 | ..self |

71 | } |

72 | } |

73 | pub fn ulps(self, ulps: i32) -> Self { |

74 | F32Margin { |

75 | ulps, |

76 | ..self |

77 | } |

78 | } |

79 | } |

80 | impl From<(f32, i32)> for F32Margin { |

81 | fn from(m: (f32, i32)) -> F32Margin { |

82 | F32Margin { |

83 | epsilon: m.0, |

84 | ulps: m.1 |

85 | } |

86 | } |

87 | } |

88 | |

89 | impl ApproxEq for f32 { |

90 | type Margin = F32Margin; |

91 | |

92 | fn approx_eq<M: Into<Self::Margin>>(self, other: f32, margin: M) -> bool { |

93 | let margin = margin.into(); |

94 | |

95 | // Check for exact equality first. This is often true, and so we get the |

96 | // performance benefit of only doing one compare in most cases. |

97 | self==other || |

98 | |

99 | // Perform epsilon comparison next |

100 | ((self - other).abs() <= margin.epsilon) || |

101 | |

102 | { |

103 | // Perform ulps comparion last |

104 | let diff: i32 = self.ulps(&other); |

105 | saturating_abs_i32!(diff) <= margin.ulps |

106 | } |

107 | } |

108 | } |

109 | |

110 | #[test] |

111 | fn f32_approx_eq_test1() { |

112 | let f: f32 = 0.0_f32; |

113 | let g: f32 = -0.0000000000000005551115123125783_f32; |

114 | assert!(f != g); // Should not be directly equal |

115 | assert!(f.approx_eq(g, (f32::EPSILON, 0)) == true); |

116 | } |

117 | #[test] |

118 | fn f32_approx_eq_test2() { |

119 | let f: f32 = 0.0_f32; |

120 | let g: f32 = -0.0_f32; |

121 | assert!(f.approx_eq(g, (f32::EPSILON, 0)) == true); |

122 | } |

123 | #[test] |

124 | fn f32_approx_eq_test3() { |

125 | let f: f32 = 0.0_f32; |

126 | let g: f32 = 0.00000000000000001_f32; |

127 | assert!(f.approx_eq(g, (f32::EPSILON, 0)) == true); |

128 | } |

129 | #[test] |

130 | fn f32_approx_eq_test4() { |

131 | let f: f32 = 0.00001_f32; |

132 | let g: f32 = 0.00000000000000001_f32; |

133 | assert!(f.approx_eq(g, (f32::EPSILON, 0)) == false); |

134 | } |

135 | #[test] |

136 | fn f32_approx_eq_test5() { |

137 | let f: f32 = 0.1_f32; |

138 | let mut sum: f32 = 0.0_f32; |

139 | for _ in 0_isize..10_isize { sum += f; } |

140 | let product: f32 = f * 10.0_f32; |

141 | assert!(sum != product); // Should not be directly equal: |

142 | assert!(sum.approx_eq(product, (f32::EPSILON, 1)) == true); |

143 | assert!(sum.approx_eq(product, F32Margin::zero()) == false); |

144 | } |

145 | #[test] |

146 | fn f32_approx_eq_test6() { |

147 | let x: f32 = 1000000_f32; |

148 | let y: f32 = 1000000.1_f32; |

149 | assert!(x != y); // Should not be directly equal |

150 | assert!(x.approx_eq(y, (0.0, 2)) == true); // 2 ulps does it |

151 | // epsilon method no good here: |

152 | assert!(x.approx_eq(y, (1000.0 * f32::EPSILON, 0)) == false); |

153 | } |

154 | |

155 | /// This type defines a margin within two `f64` values might be considered equal, |

156 | /// and is intended as the associated type for the `ApproxEq` trait. |

157 | /// |

158 | /// Two tests are used to determine approximate equality. |

159 | /// |

160 | /// The first test considers two values approximately equal if they differ by <= |

161 | /// `epsilon`. This will only succeed for very small numbers. Note that it may |

162 | /// succeed even if the parameters are of differing signs, straddling zero. |

163 | /// |

164 | /// The second test considers how many ULPs (units of least precision, units in |

165 | /// the last place, which is the integer number of floating-point representations |

166 | /// that the parameters are separated by) different the parameters are and considers |

167 | /// them approximately equal if this is <= `ulps`. For large floating-point numbers, |

168 | /// an ULP can be a rather large gap, but this kind of comparison is necessary |

169 | /// because floating-point operations must round to the nearest representable value |

170 | /// and so larger floating-point values accumulate larger errors. |

171 | #[derive(Debug, Clone, Copy)] |

172 | pub struct F64Margin { |

173 | pub epsilon: f64, |

174 | pub ulps: i64 |

175 | } |

176 | impl Default for F64Margin { |

177 | #[inline] |

178 | fn default() -> F64Margin { |

179 | F64Margin { |

180 | epsilon: f64::EPSILON, |

181 | ulps: 4 |

182 | } |

183 | } |

184 | } |

185 | impl F64Margin { |

186 | #[inline] |

187 | pub fn zero() -> F64Margin { |

188 | F64Margin { |

189 | epsilon: 0.0, |

190 | ulps: 0 |

191 | } |

192 | } |

193 | pub fn epsilon(self, epsilon: f64) -> Self { |

194 | F64Margin { |

195 | epsilon, |

196 | ..self |

197 | } |

198 | } |

199 | pub fn ulps(self, ulps: i64) -> Self { |

200 | F64Margin { |

201 | ulps, |

202 | ..self |

203 | } |

204 | } |

205 | } |

206 | impl From<(f64, i64)> for F64Margin { |

207 | fn from(m: (f64, i64)) -> F64Margin { |

208 | F64Margin { |

209 | epsilon: m.0, |

210 | ulps: m.1 |

211 | } |

212 | } |

213 | } |

214 | |

215 | impl ApproxEq for f64 { |

216 | type Margin = F64Margin; |

217 | |

218 | fn approx_eq<M: Into<Self::Margin>>(self, other: f64, margin: M) -> bool { |

219 | let margin = margin.into(); |

220 | |

221 | // Check for exact equality first. This is often true, and so we get the |

222 | // performance benefit of only doing one compare in most cases. |

223 | self == other || |

224 | |

225 | // Perform epsilon comparison next |

226 | ((self - other).abs() <= margin.epsilon) || |

227 | |

228 | { |

229 | // Perform ulps comparion last |

230 | let diff: i64 = self.ulps(&other); |

231 | saturating_abs_i64!(diff) <= margin.ulps |

232 | } |

233 | } |

234 | } |

235 | |

236 | #[test] |

237 | fn f64_approx_eq_test1() { |

238 | let f: f64 = 0.0_f64; |

239 | let g: f64 = -0.0000000000000005551115123125783_f64; |

240 | assert!(f != g); // Should not be precisely equal. |

241 | assert!(f.approx_eq(g, (3.0 * f64::EPSILON, 0)) == true); // 3e is enough. |

242 | // ULPs test won't ever call these equal. |

243 | } |

244 | #[test] |

245 | fn f64_approx_eq_test2() { |

246 | let f: f64 = 0.0_f64; |

247 | let g: f64 = -0.0_f64; |

248 | assert!(f.approx_eq(g, (f64::EPSILON, 0)) == true); |

249 | } |

250 | #[test] |

251 | fn f64_approx_eq_test3() { |

252 | let f: f64 = 0.0_f64; |

253 | let g: f64 = 1e-17_f64; |

254 | assert!(f.approx_eq(g, (f64::EPSILON, 0)) == true); |

255 | } |

256 | #[test] |

257 | fn f64_approx_eq_test4() { |

258 | let f: f64 = 0.00001_f64; |

259 | let g: f64 = 0.00000000000000001_f64; |

260 | assert!(f.approx_eq(g, (f64::EPSILON, 0)) == false); |

261 | } |

262 | #[test] |

263 | fn f64_approx_eq_test5() { |

264 | let f: f64 = 0.1_f64; |

265 | let mut sum: f64 = 0.0_f64; |

266 | for _ in 0_isize..10_isize { sum += f; } |

267 | let product: f64 = f * 10.0_f64; |

268 | assert!(sum != product); // Should not be precisely equaly. |

269 | assert!(sum.approx_eq(product, (f64::EPSILON, 0)) == true); |

270 | assert!(sum.approx_eq(product, (0.0, 1)) == true); |

271 | } |

272 | #[test] |

273 | fn f64_approx_eq_test6() { |

274 | let x: f64 = 1000000_f64; |

275 | let y: f64 = 1000000.0000000003_f64; |

276 | assert!(x != y); // Should not be precisely equal. |

277 | assert!(x.approx_eq(y, (0.0, 3)) == true); |

278 | } |

279 | #[test] |

280 | fn f64_code_triggering_issue_20() { |

281 | assert_eq!((-25.0f64).approx_eq(25.0, (0.00390625, 1)), false); |

282 | } |

283 | |

284 | impl<T> ApproxEq for &[T] |

285 | where T: Copy + ApproxEq { |

286 | type Margin = <T as ApproxEq>::Margin; |

287 | |

288 | fn approx_eq<M: Into<Self::Margin>>(self, other: Self, margin: M) -> bool { |

289 | let margin = margin.into(); |

290 | if self.len() != other.len() { return false; } |

291 | self.iter().zip(other.iter()).all(|(a,b)| { |

292 | a.approx_eq(*b, margin) |

293 | }) |

294 | } |

295 | } |

296 | |

297 | #[test] |

298 | fn test_slices() { |

299 | assert!( [1.33, 2.4, 2.5].approx_eq(&[1.33, 2.4, 2.5], (0.0, 0_i64)) ); |

300 | assert!( ! [1.33, 2.4, 2.6].approx_eq(&[1.33, 2.4, 2.5], (0.0, 0_i64)) ); |

301 | assert!( ! [1.33, 2.4].approx_eq(&[1.33, 2.4, 2.5], (0.0, 0_i64)) ); |

302 | assert!( ! [1.33, 2.4, 2.5].approx_eq(&[1.33, 2.4], (0.0, 0_i64)) ); |

303 | } |

304 | |

305 |