1 /*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17 package org.apache.commons.numbers.gamma;
18
19 import java.math.BigDecimal;
20 import java.math.MathContext;
21 import java.util.function.DoubleConsumer;
22 import java.util.function.LongConsumer;
23 import java.util.function.Supplier;
24 import org.junit.jupiter.api.Assertions;
25
26 /**
27 * Test utilities.
28 */
29 final class TestUtils {
30 /** Positive zero bits. */
31 private static final long POSITIVE_ZERO_DOUBLE_BITS = Double.doubleToRawLongBits(+0.0);
32 /** Negative zero bits. */
33 private static final long NEGATIVE_ZERO_DOUBLE_BITS = Double.doubleToRawLongBits(-0.0);
34 /** Set this to true to report all deviations to System out when the maximum ULPs is negative. */
35 private static boolean reportAllDeviations = false;
36
37 /**
38 * Class to compute the error statistics.
39 *
40 * <p>This class can be used to summary errors if used as the DoubleConsumer
41 * argument to {@link TestUtils#assertEquals(BigDecimal, double, double)}.
42 *
43 * @see <a href="https://en.wikipedia.org/wiki/Root_mean_square">Wikipedia: RMS</a>
44 */
45 static class ErrorStatistics {
46 /** Sum of squared error. */
47 private double ss;
48 /** Maximum absolute error. */
49 private double maxAbs;
50 /** Number of terms. */
51 private int n;
52 /** Positive sum. */
53 private double ps;
54 /** Positive sum round-off compensation. */
55 private double psc;
56 /** Negative sum. */
57 private double ns;
58 /** Negative sum round-off compensation. */
59 private double nsc;
60
61 /**
62 * @param x Value
63 */
64 void add(double x) {
65 // Overflow is not supported.
66 // Assume the expected and actual are quite close when measuring the RMS.
67 ss += x * x;
68 n++;
69 // Summing terms of the same sign avoids cancellation in the working sums.
70 // There may be cancellation in the final sum so the sums are totalled
71 // to 106-bit precision. This is done by adding the term through the lower
72 // then higher bits of the split quad length number.
73 if (x < 0) {
74 final double s = nsc + x;
75 nsc = twoSumLow(nsc, x, s);
76 final double t = ns + s;
77 nsc += twoSumLow(ns, s, t);
78 ns = t;
79
80 maxAbs = maxAbs < -x ? -x : maxAbs;
81 } else {
82 final double s = psc + x;
83 psc = twoSumLow(psc, x, s);
84 final double t = ps + s;
85 psc += twoSumLow(ps, s, t);
86 ps = t;
87
88 maxAbs = maxAbs < x ? x : maxAbs;
89 }
90 }
91
92 /**
93 * Gets the count of recorded values.
94 *
95 * @return the size
96 */
97 int size() {
98 return n;
99 }
100
101 /**
102 * Gets the maximum absolute error.
103 *
104 * <p>This can be used to set maximum ULP thresholds for test data if the
105 * TestUtils.assertEquals method is used with a large maxUlps to measure the ulp
106 * (and effectively ignore failures) and the maximum reported as the end of
107 * testing.
108 *
109 * @return maximum absolute error
110 */
111 double getMaxAbs() {
112 return maxAbs;
113 }
114
115 /**
116 * Gets the root mean squared error (RMS).
117 *
118 * <p> Note: If no data has been added this will return 0/0 = nan.
119 * This prevents using in assertions without adding data.
120 *
121 * @return root mean squared error (RMS)
122 */
123 double getRMS() {
124 return Math.sqrt(ss / n);
125 }
126
127 /**
128 * Gets the mean error.
129 *
130 * <p>The mean can be used to determine if the error is consistently above or
131 * below zero.
132 *
133 * @return mean error
134 */
135 double getMean() {
136 // Sum the negative parts into the positive.
137 // (see Shewchuk (1997) Grow-Expansion and Expansion-Sum [1]).
138 // [1] Shewchuk (1997): Arbitrary Precision Floating-Point Arithmetic
139 // http://www-2.cs.cmu.edu/afs/cs/project/quake/public/papers/robust-arithmetic.ps
140
141 // This creates a 3 part number (c,b,a) in descending order of magnitude.
142 double s;
143 double t;
144 s = nsc + psc;
145 double a = twoSumLow(nsc, psc, s);
146 t = s + ps;
147 double b = twoSumLow(s, ps, t);
148 double c = t;
149
150 // Sum the remaining part to create a 4 part number (d,c,b,a).
151 // Note: If ns+nsc is a non-overlapping 2 part number we can skip
152 // adding the round-off from nsc and psc (a) as it would be smaller
153 // in magnitude than nsc and hence ns. But nsc is a compensation
154 // term that may overlap ns.
155 s = ns + a;
156 a = twoSumLow(ns, a, s);
157 t = s + b;
158 b = twoSumLow(s, b, t);
159 s = t + c;
160 c = twoSumLow(t, c, s);
161 final double d = s;
162
163 // Sum the parts in order of magnitude for 1 ULP error.
164 // Reducing the error requires a two-sum rebalancing of the terms
165 // iterated through the parts.
166 return (a + b + c + d) / n;
167 }
168
169 /**
170 * Compute the round-off from the sum of two numbers {@code a} and {@code b} using
171 * Knuth's two-sum algorithm. The values are not required to be ordered by magnitude.
172 * The standard precision sum must be provided.
173 *
174 * @param a First part of sum.
175 * @param b Second part of sum.
176 * @param sum Sum of the parts (a + b).
177 * @return <code>(b - (sum - (sum - b))) + (a - (sum - b))</code>
178 * @see <a href="http://www-2.cs.cmu.edu/afs/cs/project/quake/public/papers/robust-arithmetic.ps">
179 * Shewchuk (1997) Theorum 7</a>
180 */
181 static double twoSumLow(double a, double b, double sum) {
182 final double bVirtual = sum - a;
183 // sum - bVirtual == aVirtual.
184 // a - aVirtual == a round-off
185 // b - bVirtual == b round-off
186 return (a - (sum - bVirtual)) + (b - bVirtual);
187 }
188 }
189
190 /** Private constructor. */
191 private TestUtils() {
192 // intentionally empty.
193 }
194
195 /**
196 * Assert the two numbers are equal within the provided units of least precision.
197 * The maximum count of numbers allowed between the two values is {@code maxUlps - 1}.
198 *
199 * <p>The values -0.0 and 0.0 are considered equal.
200 *
201 * <p>Set {@code maxUlps} to negative to report the ulps to the stdout and ignore
202 * failures.
203 *
204 * <p>The ulp difference is signed. It may be truncated to +/-Long.MAX_VALUE. Use of
205 * {@link Math#abs(long)} on the value will always be positive. The sign of the error
206 * is the same as that returned from Double.compare(actual, expected).
207 *
208 * @param expected expected value
209 * @param actual actual value
210 * @param maxUlps maximum units of least precision between the two values
211 * @return ulp difference between the values (signed; may be truncated to +/-Long.MAX_VALUE)
212 */
213 static long assertEquals(double expected, double actual, long maxUlps) {
214 return assertEquals(expected, actual, maxUlps, null, (Supplier<String>) null);
215 }
216
217 /**
218 * Assert the two numbers are equal within the provided units of least precision.
219 * The maximum count of numbers allowed between the two values is {@code maxUlps - 1}.
220 *
221 * <p>The values -0.0 and 0.0 are considered equal.
222 *
223 * <p>Set {@code maxUlps} to negative to report the ulps to the stdout and ignore
224 * failures.
225 *
226 * <p>The ulp difference is signed. It may be truncated to +/-Long.MAX_VALUE. Use of
227 * {@link Math#abs(long)} on the value will always be positive. The sign of the error
228 * is the same as that returned from Double.compare(actual, expected).
229 *
230 * @param expected expected value
231 * @param actual actual value
232 * @param maxUlps maximum units of least precision between the two values
233 * @param msg failure message
234 * @return ulp difference between the values (signed; may be truncated to +/-Long.MAX_VALUE)
235 */
236 static long assertEquals(double expected, double actual, long maxUlps, String msg) {
237 return assertEquals(expected, actual, maxUlps, null, () -> msg);
238 }
239
240 /**
241 * Assert the two numbers are equal within the provided units of least
242 * precision. The maximum count of numbers allowed between the two values is
243 * {@code maxUlps - 1}.
244 *
245 * <p>The values -0.0 and 0.0 are considered equal.
246 *
247 * <p>Set {@code maxUlps} to negative to report the ulps to the stdout and
248 * ignore failures.
249 *
250 * <p>The ulp difference is signed. It may be truncated to +/-Long.MAX_VALUE. Use of
251 * {@link Math#abs(long)} on the value will always be positive. The sign of the error
252 * is the same as that returned from Double.compare(actual, expected).
253 *
254 * @param expected expected value
255 * @param actual actual value
256 * @param maxUlps maximum units of least precision between the two values
257 * @param error Consumer for the ulp difference between the values (signed)
258 * @param msg failure message
259 * @return ulp difference between the values (signed; may be truncated to +/-Long.MAX_VALUE)
260 */
261 static long assertEquals(double expected, double actual, long maxUlps,
262 LongConsumer error, Supplier<String> msg) {
263 final long e = Double.doubleToLongBits(expected);
264 final long a = Double.doubleToLongBits(actual);
265
266 // Code adapted from Precision#equals(double, double, int) so we maintain the delta
267 // for the message and return it for reporting. The sign is maintained separately
268 // to allow reporting errors above Long.MAX_VALUE.
269
270 int sign;
271 long delta;
272 boolean equal;
273 if (e == a) {
274 // Binary equal
275 equal = true;
276 sign = 0;
277 delta = 0;
278 } else if ((a ^ e) < 0L) {
279 // The difference is the count of numbers between each and zero.
280 // This makes -0.0 and 0.0 equal.
281 long d1;
282 long d2;
283 if (a < e) {
284 sign = -1;
285 d1 = e - POSITIVE_ZERO_DOUBLE_BITS;
286 d2 = a - NEGATIVE_ZERO_DOUBLE_BITS;
287 } else {
288 sign = 1;
289 d1 = a - POSITIVE_ZERO_DOUBLE_BITS;
290 d2 = e - NEGATIVE_ZERO_DOUBLE_BITS;
291 }
292 // This may overflow but we report it using an unsigned formatter.
293 delta = d1 + d2;
294 if (delta < 0) {
295 // Overflow
296 equal = false;
297 } else {
298 // Allow input of a negative maximum ULPs
299 equal = delta <= ((maxUlps < 0) ? (-maxUlps - 1) : maxUlps);
300 }
301 } else {
302 if (a < e) {
303 sign = -1;
304 delta = e - a;
305 } else {
306 sign = 1;
307 delta = a - e;
308 }
309 // The sign must be negated for negative doubles since the magnitude
310 // comparison (a < e) included the sign bit.
311 sign = a < 0 ? -sign : sign;
312
313 // Allow input of a negative maximum ULPs
314 equal = delta <= ((maxUlps < 0) ? (-maxUlps - 1) : maxUlps);
315 }
316
317 Assertions.assertEquals(sign, Double.compare(actual, expected));
318
319 // DEBUG:
320 if (maxUlps < 0) {
321 // CHECKSTYLE: stop Regex
322 if (!equal || reportAllDeviations) {
323 System.out.printf("%sexpected <%s> != actual <%s> (ulps=%c%s)%n",
324 prefix(msg), expected, actual, sign < 0 ? '-' : ' ', Long.toUnsignedString(delta));
325 }
326 // CHECKSTYLE: resume Regex
327 } else if (!equal) {
328 Assertions.fail(String.format("%sexpected <%s> != actual <%s> (ulps=%c%s)",
329 prefix(msg), expected, actual, sign < 0 ? '-' : ' ', Long.toUnsignedString(delta)));
330 }
331
332 // This may have overflowed.
333 delta = delta < 0 ? Long.MAX_VALUE : delta;
334 delta *= sign;
335 if (error != null) {
336 error.accept(delta);
337 }
338 return delta;
339 }
340
341 /**
342 * Assert the two numbers are equal within the provided units of least precision.
343 *
344 * <p>This method is for values that can be computed to arbitrary precision.
345 * It raises an exception when the actual value is not finite and the expected value
346 * has a non-infinite representation; or the actual value is finite and the expected
347 * value has a infinite representation. In this case the computed ulp difference
348 * is infinite.
349 *
350 * <p>This method expresses the error relative the units in the last place of the
351 * expected value when converted to a {@code double} type
352 * (see {@link #assertEquals(BigDecimal, double, double, DoubleConsumer, Supplier)} for details).
353 *
354 * <p>Set {@code maxUlps} to negative to report the ulps to the stdout and ignore
355 * failures.
356 *
357 * <p>The ulp difference is signed. The sign of the error
358 * is the same as that returned from Double.compare(actual, expected); it is
359 * computed using {@code actual - expected}.
360 *
361 * @param expected expected value
362 * @param actual actual value
363 * @param maxUlps maximum units of least precision between the two values
364 * @return ulp difference between the values (signed)
365 */
366 static double assertEquals(BigDecimal expected, double actual, double maxUlps) {
367 return assertEquals(expected, actual, maxUlps, null, (Supplier<String>) null);
368 }
369
370 /**
371 * Assert the two numbers are equal within the provided units of least precision.
372 *
373 * <p>This method is for values that can be computed to arbitrary precision.
374 * It raises an exception when the actual value is not finite and the expected value
375 * has a non-infinite representation; or the actual value is finite and the expected
376 * value has a infinite representation. In this case the computed ulp difference
377 * is infinite.
378 *
379 * <p>This method expresses the error relative the units in the last place of the
380 * expected value when converted to a {@code double} type
381 * (see {@link #assertEquals(BigDecimal, double, double, DoubleConsumer, Supplier)} for details).
382 *
383 * <p>Set {@code maxUlps} to negative to report the ulps to the stdout and ignore
384 * failures.
385 *
386 * <p>The ulp difference is signed. The sign of the error
387 * is the same as that returned from Double.compare(actual, expected); it is
388 * computed using {@code actual - expected}.
389 *
390 * @param expected expected value
391 * @param actual actual value
392 * @param maxUlps maximum units of least precision between the two values
393 * @param msg failure message
394 * @return ulp difference between the values (signed)
395 */
396 static double assertEquals(BigDecimal expected, double actual, double maxUlps, String msg) {
397 return assertEquals(expected, actual, maxUlps, null, () -> msg);
398 }
399
400 /**
401 * Assert the two numbers are equal within the provided units of least
402 * precision.
403 *
404 * <p>This method is for values that can be computed to arbitrary precision. It
405 * raises an exception when the actual value is not finite and the expected
406 * value has a non-infinite representation; or the actual value is finite and
407 * the expected value has a infinite representation. In this case the computed
408 * ulp difference is infinite.
409 *
410 * <p>This method expresses the error relative the units in the last place (ulp)
411 * of the expected value when converted to a {@code double} type. If the actual
412 * value equals the expected value the error is 0. Otherwise the error is
413 * computed relative to the ulp of the expected value. The minimum non-zero
414 * value for the error is 0.5. A |ulp| of < 1.0 indicates the value is the closest
415 * value to the result that is not exact.
416 *
417 * <pre>
418 * ulp -1 0 1 2
419 * --------------|---------------|---------------|---------------|--------
420 * ^
421 * |
422 * expected
423 *
424 * a b c
425 *
426 * a = -0.75 ulp
427 * b = 0 ulp
428 * c = 1.25 ulp
429 * </pre>
430 *
431 * <p>Set {@code maxUlps} to negative to report the ulps to the stdout and
432 * ignore failures.
433 *
434 * <p>The ulp difference is signed. The sign of the error
435 * is the same as that returned from Double.compare(actual, expected); it is
436 * computed using {@code actual - expected}.
437 *
438 * @param expected expected value
439 * @param actual actual value
440 * @param maxUlps maximum units of least precision between the two values
441 * @param error Consumer for the ulp difference between the values (always positive)
442 * @param msg failure message
443 * @return ulp difference between the values (signed)
444 */
445 static double assertEquals(BigDecimal expected, double actual, double maxUlps,
446 DoubleConsumer error, Supplier<String> msg) {
447 final double e = expected.doubleValue();
448
449 double delta;
450 boolean equal;
451 if (e == actual) {
452 // Binary equal. This will match infinity if expected is very large.
453 equal = true;
454 delta = 0;
455 } else if (!Double.isFinite(e) || !Double.isFinite(actual)) {
456 // No representable delta between infinite and non-infinite values
457 equal = false;
458 delta = Double.compare(actual, e) * Double.POSITIVE_INFINITY;
459 } else {
460 // Two finite numbers
461 delta = new BigDecimal(actual).subtract(expected)
462 .divide(new BigDecimal(Math.ulp(e)), MathContext.DECIMAL64).doubleValue();
463 // Allow input of a negative maximum ULPs
464 equal = Math.abs(delta) <= Math.abs(maxUlps);
465 }
466
467 if (error != null) {
468 error.accept(delta);
469 }
470
471 // DEBUG:
472 if (maxUlps < 0) {
473 // CHECKSTYLE: stop Regex
474 if (!equal || reportAllDeviations) {
475 System.out.printf("%sexpected <%s> != actual <%s> (ulps=%s)%n",
476 prefix(msg), expected, actual, delta);
477 }
478 // CHECKSTYLE: resume Regex
479 } else if (!equal) {
480 Assertions.fail(String.format("%sexpected <%s> != actual <%s> (ulps=%s)",
481 prefix(msg), expected, actual, delta));
482 }
483
484 return delta;
485 }
486
487 /**
488 * Get the prefix for the message.
489 *
490 * @param msg Message supplier
491 * @return the prefix
492 */
493 private static String prefix(Supplier<String> msg) {
494 return msg == null ? "" : msg.get() + ": ";
495 }
496 }