1 /// This module is about comparison of floating point arithmetics.
2 /// Supported by this very informative article:
3 /// https://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/
4 module gfx.math.approx;
5
6 import std.traits : isFloatingPoint;
7
8 pure @safe @nogc nothrow:
9
10 /// Different methods to check if two floating point values are close to each other.
11 enum ApproxMethod
12 {
13 /// Unit of Least Precision.
14 ulp,
15 /// Relative error
16 rel,
17 /// ULP, with absolute error
18 ulpAndAbs,
19 /// Relative and absolute error
20 relAndAbs,
21 }
22
23 /// Check with method given at runtime with default parameters
24 /// a and b can be any data sets supported by the approx template
25 bool approx(T)(ApproxMethod method, in T a, in T b)
26 {
27 final switch(method)
28 {
29 case ApproxMethod.ulp:
30 return approxUlp(a, b);
31 case ApproxMethod.rel:
32 return approxRel(a, b);
33 case ApproxMethod.ulpAndAbs:
34 return approxUlpAndAbs(a, b);
35 case ApproxMethod.relAndAbs:
36 return approxRelAndAbs(a, b);
37 }
38 }
39
40 /// Compare two data sets with ApproxMethod.ulp.
41 /// That is call scalarApproxUlp on each couple of the data sets.
42 alias approxUlp = approx!(ApproxMethod.ulp);
43 /// Compare two data sets with ApproxMethod.rel.
44 /// That is call scalarApproxRel on each couple of the data sets.
45 alias approxRel = approx!(ApproxMethod.rel);
46 /// Compare two data sets with ApproxMethod.ulpAndAbs.
47 /// That is call scalarApproxUlpAndAbs on each couple of the data sets.
48 alias approxUlpAndAbs = approx!(ApproxMethod.ulpAndAbs);
49 /// Compare two data sets with ApproxMethod.relAndAbs.
50 /// That is call scalarApproxRelAndAbs on each couple of the data sets.
51 alias approxRelAndAbs = approx!(ApproxMethod.relAndAbs);
52
53
54 /// Compute the ULP difference between two floating point numbers
55 /// Negative result indicates that b has higher ULP value than a.
56 template ulpDiff(T) if (isFloatingPoint!T)
57 {
58 long ulpDiff(in T a, in T b)
59 {
60 immutable fnA = FloatNum!T(a);
61 immutable fnB = FloatNum!T(b);
62
63 return (fnA - fnB);
64 }
65 }
66
67 /// Determines if two floating point scalars are maxUlps close to each other.
68 template scalarApproxUlp(T) if (isFloatingPoint!T)
69 {
70 bool scalarApproxUlp(in T a, in T b, in ulong maxUlps=4)
71 {
72 import std.math : abs;
73
74 immutable fnA = FloatNum!T(a);
75 immutable fnB = FloatNum!T(b);
76
77 if (fnA.negative != fnB.negative)
78 {
79 return a == b; // check for +0 / -0
80 }
81
82 return (abs(fnA - fnB) <= maxUlps);
83 }
84 }
85
86 /// Check whether the relative error between a and b is smaller than maxEps
87 template scalarApproxRel(T) if (isFloatingPoint!T)
88 {
89 bool scalarApproxRel (in T a, in T b, in T maxEps=4*T.epsilon)
90 {
91 import std.algorithm : max;
92 import std.math : abs;
93 const diff = abs(b-a);
94 const absA = abs(a);
95 const absB = abs(b);
96 const largest = max(absA, absB);
97 return diff <= maxEps*largest;
98 }
99 }
100
101 /// Determines if two floating point scalars are maxUlps close to each other.
102 /// If the absolute error is less than maxAbs, the test succeeds however.
103 /// This is useful when comparing against zero the result of a subtraction.
104 template scalarApproxUlpAndAbs(T) if (isFloatingPoint!T)
105 {
106 bool scalarApproxUlpAndAbs(in T a, in T b, in T maxAbs=0.0001, in ulong maxUlps=4)
107 {
108 import std.math : abs;
109 if (abs(b-a) <= maxAbs) return true;
110 return scalarApproxUlp(a, b, maxUlps);
111 }
112 }
113
114 /// Check whether the relative error between a and b is smaller than maxEps.
115 /// If the absolute error is less than maxAbs, the test succeeds however.
116 /// This is useful when comparing against zero the result of a subtraction.
117 template scalarApproxRelAndAbs(T) if (isFloatingPoint!T)
118 {
119 bool scalarApproxRelAndAbs(in T a, in T b, in T maxAbs=0.0001, in T maxEps=4*T.epsilon)
120 {
121 import std.math : abs;
122 if (abs(b-a) <= maxAbs) return true;
123 return scalarApproxRel(a, b, maxEps);
124 }
125 }
126
127 /// Generic template to check approx method on different sets of data.
128 template approx(ApproxMethod method)
129 {
130 import gfx.math.mat : Mat;
131 import gfx.math.vec : Vec;
132
133 /// Apply method on scalar
134 bool approx(T, Args...)(in T a, in T b, Args args)
135 if (isFloatingPoint!T)
136 {
137 static if (method == ApproxMethod.ulp)
138 {
139 return scalarApproxUlp(a, b, args);
140 }
141 else static if (method == ApproxMethod.rel)
142 {
143 return scalarApproxRel(a, b, args);
144 }
145 else static if (method == ApproxMethod.ulpAndAbs)
146 {
147 return scalarApproxUlpAndAbs(a, b, args);
148 }
149 else static if (method == ApproxMethod.relAndAbs)
150 {
151 return scalarApproxRelAndAbs(a, b, args);
152 }
153 else
154 {
155 static assert(false, "unimplemented");
156 }
157 }
158 /// Apply method check on vectors
159 bool approx(T, size_t N, Args...)(in T[N] v1, in T[N] v2, Args args)
160 if (isFloatingPoint!T)
161 {
162 static foreach (i; 0 .. N) {
163 if (!approx(v1[i], v2[i], args)) return false;
164 }
165 return true;
166 }
167 /// ditto
168 bool approx(T, size_t N, Args...)(in Vec!(T, N) v1, in Vec!(T, N) v2, Args args)
169 if (isFloatingPoint!T)
170 {
171 return approx(v1.data, v2.data, args);
172 }
173 /// Apply method check on arrays
174 bool approx(T, Args...)(in T[] arr1, in T[] arr2, Args args)
175 if (isFloatingPoint!T)
176 {
177 if (arr1.length != arr2.length) return false;
178 foreach (i; 0 .. arr1.length) {
179 if (!approx(arr1[i], arr2[i], args)) return false;
180 }
181 return true;
182 }
183 /// Apply method check on matrices
184 bool approx(T, size_t R, size_t C, Args...)(in T[C][R] m1, in T[C][R] m2, Args args)
185 if (isFloatingPoint!T)
186 {
187 static foreach (r; 0 .. R) {
188 if (!approx(m1[r], m2[r], args)) return false;
189 }
190 return true;
191 }
192 // /// ditto
193 bool approx(T, size_t R, size_t C, Args...)(in Mat!(T, R, C) m1, in Mat!(T, R, C) m2, Args args)
194 if (isFloatingPoint!T)
195 {
196 static foreach (r; 0 .. R) {
197 if (!approx(m1[r], m2[r], args)) return false;
198 }
199 return true;
200 }
201 }
202
203
204 private:
205
206 template FloatTraits(T)
207 if (isFloatingPoint!T)
208 {
209 static if (is(T == float))
210 {
211 alias IntType = int;
212 enum exponentMask = 0x7f80_0000;
213 }
214 else static if (is(T == double) || (is(T == real) && T.sizeof == 8))
215 {
216 alias IntType = long;
217 enum exponentMask = 0x7ff0_0000_0000_0000;
218 }
219 else static if (is(T == real) && T.sizeof > 8)
220 {
221 static assert(false, "FloatTraits unsupported for real.");
222 }
223 }
224
225 union FloatNum(T)
226 if (is(T == float) || is(T == double) || is(T == real) && T.sizeof == 8)
227 {
228 alias F = FloatTraits!T;
229
230 alias FloatType = T;
231 alias IntType = F.IntType;
232
233 this (FloatType f)
234 {
235 this.f = f;
236 }
237
238 FloatType f;
239 IntType i;
240
241 long opBinary(string op)(FloatNum rhs) const
242 if (op == "-")
243 {
244 return i - rhs.i;
245 }
246
247 @property bool negative() const
248 {
249 return i < 0;
250 }
251
252 debug @property IntType mantissa() const
253 {
254 enum IntType one = 1;
255 return i & ((one << T.mant_dig) - one);
256 }
257
258 debug @property IntType exponent() const
259 {
260 return ((i & F.exponentMask) >> T.mant_dig);
261 }
262 }
263
264
265 union FloatNum(T)
266 if (is(T == real) && T.sizeof > 8)
267 {
268 static assert(T.sizeof == 16 || T.sizeof == 12 || T.sizeof == 10, "Unexpected real size.");
269
270 static if (T.sizeof == 10)
271 {
272 struct Int
273 {
274 long l;
275 ushort h;
276 }
277 }
278 static if (T.sizeof == 12)
279 {
280 struct Int
281 {
282 long l;
283 uint h;
284 }
285 }
286 static if (T.sizeof == 16)
287 {
288 struct Int
289 {
290 long l;
291 ulong h;
292 }
293 }
294
295 real f;
296 Int i;
297
298 this (real f)
299 {
300 this.f = f;
301 }
302
303 long opBinary(string op)(FloatNum rhs) const
304 if (op == "-")
305 {
306 immutable lo = i.l - rhs.i.l;
307 // do not check for overflow of the hi long
308 return lo;
309 }
310
311 @property bool negative() const
312 {
313 return f < 0;
314 }
315 }
316
317 unittest
318 {
319 immutable fn1 = FloatNum!real(1);
320 immutable fn2 = FloatNum!real(-1);
321
322 assert(! fn1.negative);
323 assert( fn2.negative);
324 }
325
326 unittest
327 {
328 import std.math : nextUp;
329 immutable real r1 = 0;
330 immutable real r2 = nextUp(r1);
331
332 immutable fn1 = FloatNum!real(r1);
333 immutable fn2 = FloatNum!real(r2);
334
335 assert(fn2 - fn1 == 1);
336 }