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
18 package org.apache.commons.numbers.core;
19
20 import java.math.BigDecimal;
21 import java.math.RoundingMode;
22
23 /**
24 * Utilities for comparing numbers.
25 */
26 public final class Precision {
27 /**
28 * <p>
29 * Largest double-precision floating-point number such that
30 * {@code 1 + EPSILON} is numerically equal to 1. This value is an upper
31 * bound on the relative error due to rounding real numbers to double
32 * precision floating-point numbers.
33 * </p>
34 * <p>
35 * In IEEE 754 arithmetic, this is 2<sup>-53</sup>.
36 * </p>
37 *
38 * @see <a href="http://en.wikipedia.org/wiki/Machine_epsilon">Machine epsilon</a>
39 */
40 public static final double EPSILON;
41
42 /**
43 * Safe minimum, such that {@code 1 / SAFE_MIN} does not overflow.
44 * In IEEE 754 arithmetic, this is also the smallest normalized
45 * number 2<sup>-1022</sup>.
46 *
47 * @see Double#MIN_NORMAL
48 */
49 public static final double SAFE_MIN = Double.MIN_NORMAL;
50
51 /** Exponent offset in IEEE754 representation. */
52 private static final long EXPONENT_OFFSET = 1023L;
53
54 /** Positive zero. */
55 private static final double POSITIVE_ZERO = 0d;
56
57 static {
58 /*
59 * This was previously expressed as = 0x1.0p-53
60 * However, OpenJDK (Sparc Solaris) cannot handle such small
61 * constants: MATH-721
62 */
63 EPSILON = Double.longBitsToDouble((EXPONENT_OFFSET - 53L) << 52);
64 }
65
66 /**
67 * Private constructor.
68 */
69 private Precision() {}
70
71 /**
72 * Compares two numbers given some amount of allowed error.
73 * The returned value is:
74 * <ul>
75 * <li>zero if considered equal using {@link #equals(double,double,double) equals(x, y, eps)}
76 * <li>negative if not equal and {@code x < y}
77 * <li>positive if not equal and {@code x > y}
78 * </ul>
79 *
80 * <p>NaN values are handled as if using {@link Double#compare(double, double)} where the
81 * returned value is:
82 * <ul>
83 * <li>zero if {@code NaN, NaN}
84 * <li>negative if {@code !NaN, NaN}
85 * <li>positive if {@code NaN, !NaN}
86 * </ul>
87 *
88 * @param x First value.
89 * @param y Second value.
90 * @param eps Allowed error when checking for equality.
91 * @return 0 if the value are considered equal, -1 if the first is smaller than
92 * the second, 1 if the first is larger than the second.
93 * @see #equals(double, double, double)
94 */
95 public static int compareTo(double x, double y, double eps) {
96 if (equals(x, y, eps)) {
97 return 0;
98 } else if (x < y) {
99 return -1;
100 } else if (x > y) {
101 return 1;
102 }
103 // NaN input.
104 return Double.compare(x, y);
105 }
106
107 /**
108 * Compares two numbers given some amount of allowed error.
109 * The returned value is:
110 * <ul>
111 * <li>zero if considered equal using {@link #equals(double,double,int) equals(x, y, maxUlps)}
112 * <li>negative if not equal and {@code x < y}
113 * <li>positive if not equal and {@code x > y}
114 * </ul>
115 *
116 * <p>NaN values are handled as if using {@link Double#compare(double, double)} where the
117 * returned value is:
118 * <ul>
119 * <li>zero if {@code NaN, NaN}
120 * <li>negative if {@code !NaN, NaN}
121 * <li>positive if {@code NaN, !NaN}
122 * </ul>
123 *
124 * @param x First value.
125 * @param y Second value.
126 * @param maxUlps {@code (maxUlps - 1)} is the number of floating point
127 * values between {@code x} and {@code y}.
128 * @return 0 if the value are considered equal, -1 if the first is smaller than
129 * the second, 1 if the first is larger than the second.
130 * @see #equals(double, double, int)
131 */
132 public static int compareTo(final double x, final double y, final int maxUlps) {
133 if (equals(x, y, maxUlps)) {
134 return 0;
135 } else if (x < y) {
136 return -1;
137 } else if (x > y) {
138 return 1;
139 }
140 // NaN input.
141 return Double.compare(x, y);
142 }
143
144 /**
145 * Returns true iff they are equal as defined by
146 * {@link #equals(float,float,int) equals(x, y, 1)}.
147 *
148 * @param x first value
149 * @param y second value
150 * @return {@code true} if the values are equal.
151 */
152 public static boolean equals(float x, float y) {
153 return equals(x, y, 1);
154 }
155
156 /**
157 * Returns true if both arguments are NaN or they are
158 * equal as defined by {@link #equals(float,float) equals(x, y, 1)}.
159 *
160 * @param x first value
161 * @param y second value
162 * @return {@code true} if the values are equal or both are NaN.
163 */
164 public static boolean equalsIncludingNaN(float x, float y) {
165 final boolean xIsNan = Float.isNaN(x);
166 final boolean yIsNan = Float.isNaN(y);
167 // Combine the booleans with bitwise OR
168 return (xIsNan | yIsNan) ?
169 xIsNan == yIsNan :
170 equals(x, y, 1);
171 }
172
173 /**
174 * Returns {@code true} if there is no float value strictly between the
175 * arguments or the difference between them is within the range of allowed
176 * error (inclusive). Returns {@code false} if either of the arguments
177 * is NaN.
178 *
179 * @param x first value
180 * @param y second value
181 * @param eps the amount of absolute error to allow.
182 * @return {@code true} if the values are equal or within range of each other.
183 */
184 public static boolean equals(float x, float y, float eps) {
185 return equals(x, y, 1) || Math.abs(y - x) <= eps;
186 }
187
188 /**
189 * Returns true if the arguments are both NaN, there are no float value strictly
190 * between the arguments or the difference between them is within the range of allowed
191 * error (inclusive).
192 *
193 * @param x first value
194 * @param y second value
195 * @param eps the amount of absolute error to allow.
196 * @return {@code true} if the values are equal or within range of each other,
197 * or both are NaN.
198 */
199 public static boolean equalsIncludingNaN(float x, float y, float eps) {
200 return equalsIncludingNaN(x, y, 1) || Math.abs(y - x) <= eps;
201 }
202
203 /**
204 * Returns true if the arguments are equal or within the range of allowed
205 * error (inclusive). Returns {@code false} if either of the arguments is NaN.
206 * <p>
207 * Two double numbers are considered equal if there are {@code (maxUlps - 1)}
208 * (or fewer) floating point numbers between them, i.e. two adjacent
209 * floating point numbers are considered equal.
210 * </p>
211 * <p>
212 * Adapted from <a
213 * href="http://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/">
214 * Bruce Dawson</a>.
215 * </p>
216 *
217 * @param x first value
218 * @param y second value
219 * @param maxUlps {@code (maxUlps - 1)} is the number of floating point
220 * values between {@code x} and {@code y}.
221 * @return {@code true} if there are fewer than {@code maxUlps} floating
222 * point values between {@code x} and {@code y}.
223 */
224 public static boolean equals(final float x, final float y, final int maxUlps) {
225 final int xInt = Float.floatToRawIntBits(x);
226 final int yInt = Float.floatToRawIntBits(y);
227
228 final boolean isEqual;
229 if ((xInt ^ yInt) < 0) {
230 // Numbers have opposite signs, take care of overflow.
231 // Remove the sign bit to obtain the absolute ULP above zero.
232 final int deltaPlus = xInt & Integer.MAX_VALUE;
233 final int deltaMinus = yInt & Integer.MAX_VALUE;
234
235 // Avoid possible overflow from adding the deltas by using a long.
236 isEqual = (long) deltaPlus + deltaMinus <= maxUlps;
237 } else {
238 // Numbers have same sign, there is no risk of overflow.
239 isEqual = Math.abs(xInt - yInt) <= maxUlps;
240 }
241
242 return isEqual && !Float.isNaN(x) && !Float.isNaN(y);
243 }
244
245 /**
246 * Returns true if both arguments are NaN or if they are equal as defined
247 * by {@link #equals(float,float,int) equals(x, y, maxUlps)}.
248 *
249 * @param x first value
250 * @param y second value
251 * @param maxUlps {@code (maxUlps - 1)} is the number of floating point
252 * values between {@code x} and {@code y}.
253 * @return {@code true} if both arguments are NaN or if there are less than
254 * {@code maxUlps} floating point values between {@code x} and {@code y}.
255 */
256 public static boolean equalsIncludingNaN(float x, float y, int maxUlps) {
257 final boolean xIsNan = Float.isNaN(x);
258 final boolean yIsNan = Float.isNaN(y);
259 // Combine the booleans with bitwise OR
260 return (xIsNan | yIsNan) ?
261 xIsNan == yIsNan :
262 equals(x, y, maxUlps);
263 }
264
265 /**
266 * Returns true iff they are equal as defined by
267 * {@link #equals(double,double,int) equals(x, y, 1)}.
268 *
269 * @param x first value
270 * @param y second value
271 * @return {@code true} if the values are equal.
272 */
273 public static boolean equals(double x, double y) {
274 return equals(x, y, 1);
275 }
276
277 /**
278 * Returns true if the arguments are both NaN or they are
279 * equal as defined by {@link #equals(double,double) equals(x, y, 1)}.
280 *
281 * @param x first value
282 * @param y second value
283 * @return {@code true} if the values are equal or both are NaN.
284 */
285 public static boolean equalsIncludingNaN(double x, double y) {
286 final boolean xIsNan = Double.isNaN(x);
287 final boolean yIsNan = Double.isNaN(y);
288 // Combine the booleans with bitwise OR
289 return (xIsNan | yIsNan) ?
290 xIsNan == yIsNan :
291 equals(x, y, 1);
292 }
293
294 /**
295 * Returns {@code true} if there is no double value strictly between the
296 * arguments or the difference between them is within the range of allowed
297 * error (inclusive). Returns {@code false} if either of the arguments
298 * is NaN.
299 *
300 * @param x First value.
301 * @param y Second value.
302 * @param eps Amount of allowed absolute error.
303 * @return {@code true} if the values are equal or within range of each other.
304 */
305 public static boolean equals(double x, double y, double eps) {
306 return equals(x, y, 1) || Math.abs(y - x) <= eps;
307 }
308
309 /**
310 * Returns {@code true} if there is no double value strictly between the
311 * arguments or the relative difference between them is less than or equal
312 * to the given tolerance. Returns {@code false} if either of the arguments
313 * is NaN.
314 *
315 * @param x First value.
316 * @param y Second value.
317 * @param eps Amount of allowed relative error.
318 * @return {@code true} if the values are two adjacent floating point
319 * numbers or they are within range of each other.
320 */
321 public static boolean equalsWithRelativeTolerance(double x, double y, double eps) {
322 if (equals(x, y, 1)) {
323 return true;
324 }
325
326 final double absoluteMax = Math.max(Math.abs(x), Math.abs(y));
327 final double relativeDifference = Math.abs((x - y) / absoluteMax);
328
329 return relativeDifference <= eps;
330 }
331
332 /**
333 * Returns true if the arguments are both NaN, there are no double value strictly
334 * between the arguments or the difference between them is within the range of allowed
335 * error (inclusive).
336 *
337 * @param x first value
338 * @param y second value
339 * @param eps the amount of absolute error to allow.
340 * @return {@code true} if the values are equal or within range of each other,
341 * or both are NaN.
342 */
343 public static boolean equalsIncludingNaN(double x, double y, double eps) {
344 return equalsIncludingNaN(x, y) || Math.abs(y - x) <= eps;
345 }
346
347 /**
348 * Returns true if the arguments are equal or within the range of allowed
349 * error (inclusive). Returns {@code false} if either of the arguments is NaN.
350 * <p>
351 * Two float numbers are considered equal if there are {@code (maxUlps - 1)}
352 * (or fewer) floating point numbers between them, i.e. two adjacent
353 * floating point numbers are considered equal.
354 * </p>
355 * <p>
356 * Adapted from <a
357 * href="http://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/">
358 * Bruce Dawson</a>.
359 * </p>
360 *
361 * @param x first value
362 * @param y second value
363 * @param maxUlps {@code (maxUlps - 1)} is the number of floating point
364 * values between {@code x} and {@code y}.
365 * @return {@code true} if there are fewer than {@code maxUlps} floating
366 * point values between {@code x} and {@code y}.
367 */
368 public static boolean equals(final double x, final double y, final int maxUlps) {
369 final long xInt = Double.doubleToRawLongBits(x);
370 final long yInt = Double.doubleToRawLongBits(y);
371
372 if ((xInt ^ yInt) < 0) {
373 // Numbers have opposite signs, take care of overflow.
374 // Remove the sign bit to obtain the absolute ULP above zero.
375 final long deltaPlus = xInt & Long.MAX_VALUE;
376 final long deltaMinus = yInt & Long.MAX_VALUE;
377
378 // Note:
379 // If either value is NaN, the exponent bits are set to (2047L << 52) and the
380 // distance above 0.0 is always above an integer ULP error. So omit the test
381 // for NaN and return directly.
382
383 // Avoid possible overflow from adding the deltas by splitting the comparison
384 return deltaPlus <= maxUlps && deltaMinus <= (maxUlps - deltaPlus);
385 }
386
387 // Numbers have same sign, there is no risk of overflow.
388 return Math.abs(xInt - yInt) <= maxUlps && !Double.isNaN(x) && !Double.isNaN(y);
389 }
390
391 /**
392 * Returns true if both arguments are NaN or if they are equal as defined
393 * by {@link #equals(double,double,int) equals(x, y, maxUlps)}.
394 *
395 * @param x first value
396 * @param y second value
397 * @param maxUlps {@code (maxUlps - 1)} is the number of floating point
398 * values between {@code x} and {@code y}.
399 * @return {@code true} if both arguments are NaN or if there are less than
400 * {@code maxUlps} floating point values between {@code x} and {@code y}.
401 */
402 public static boolean equalsIncludingNaN(double x, double y, int maxUlps) {
403 final boolean xIsNan = Double.isNaN(x);
404 final boolean yIsNan = Double.isNaN(y);
405 // Combine the booleans with bitwise OR
406 return (xIsNan | yIsNan) ?
407 xIsNan == yIsNan :
408 equals(x, y, maxUlps);
409 }
410
411 /**
412 * Rounds the given value to the specified number of decimal places.
413 * The value is rounded using the {@link RoundingMode#HALF_UP} method.
414 *
415 * @param x Value to round.
416 * @param scale Number of digits to the right of the decimal point.
417 * @return the rounded value.
418 */
419 public static double round(double x, int scale) {
420 return round(x, scale, RoundingMode.HALF_UP);
421 }
422
423 /**
424 * Rounds the given value to the specified number of decimal places.
425 * The value is rounded using the given method which is any method defined
426 * in {@link BigDecimal}.
427 * If {@code x} is infinite or {@code NaN}, then the value of {@code x} is
428 * returned unchanged, regardless of the other parameters.
429 *
430 * @param x Value to round.
431 * @param scale Number of digits to the right of the decimal point.
432 * @param roundingMethod Rounding method as defined in {@link BigDecimal}.
433 * @return the rounded value.
434 * @throws ArithmeticException if {@code roundingMethod} is
435 * {@link RoundingMode#UNNECESSARY} and the specified scaling operation
436 * would require rounding.
437 */
438 public static double round(double x,
439 int scale,
440 RoundingMode roundingMethod) {
441 try {
442 final double rounded = (new BigDecimal(Double.toString(x))
443 .setScale(scale, roundingMethod))
444 .doubleValue();
445 // MATH-1089: negative values rounded to zero should result in negative zero
446 return rounded == POSITIVE_ZERO ? POSITIVE_ZERO * x : rounded;
447 } catch (NumberFormatException ex) {
448 if (Double.isInfinite(x)) {
449 return x;
450 }
451 return Double.NaN;
452 }
453 }
454
455 /**
456 * Computes a number close to {@code delta} with the property that
457 * {@code (x + delta - x)} is exactly machine-representable.
458 * This is useful when computing numerical derivatives, in order to
459 * reduce roundoff errors.
460 *
461 * @param x Value.
462 * @param delta Offset value.
463 * @return the machine-representable floating number closest to the
464 * difference between {@code x + delta} and {@code x}.
465 */
466 public static double representableDelta(double x,
467 double delta) {
468 return x + delta - x;
469 }
470
471 /**
472 * Creates a {@link DoubleEquivalence} instance that uses the given epsilon
473 * value for determining equality.
474 *
475 * @param eps Value to use for determining equality.
476 * @return a new instance.
477 */
478 public static DoubleEquivalence doubleEquivalenceOfEpsilon(final double eps) {
479 if (!Double.isFinite(eps) ||
480 eps < 0d) {
481 throw new IllegalArgumentException("Invalid epsilon value: " + eps);
482 }
483
484 return new DoubleEquivalence() {
485 /** Epsilon value. */
486 private final double epsilon = eps;
487
488 /** {@inheritDoc} */
489 @Override
490 public int compare(double a,
491 double b) {
492 return Precision.compareTo(a, b, epsilon);
493 }
494 };
495 }
496
497 /**
498 * Interface containing comparison operations for doubles that allow values
499 * to be <em>considered</em> equal even if they are not exactly equal.
500 * It is intended for comparing outputs of a computation where floating
501 * point errors may have occurred.
502 */
503 public interface DoubleEquivalence {
504 /**
505 * Indicates whether given values are considered equal to each other.
506 *
507 * @param a Value.
508 * @param b Value.
509 * @return true if the given values are considered equal.
510 */
511 default boolean eq(double a, double b) {
512 return compare(a, b) == 0;
513 }
514
515 /**
516 * Indicates whether the given value is considered equal to zero.
517 * It is a shortcut for {@code eq(a, 0.0)}.
518 *
519 * @param a Value.
520 * @return true if the argument is considered equal to zero.
521 */
522 default boolean eqZero(double a) {
523 return eq(a, 0d);
524 }
525
526 /**
527 * Indicates whether the first argument is strictly smaller than the second.
528 *
529 * @param a Value.
530 * @param b Value.
531 * @return true if {@code a < b}
532 */
533 default boolean lt(double a, double b) {
534 return compare(a, b) < 0;
535 }
536
537 /**
538 * Indicates whether the first argument is smaller or considered equal to the second.
539 *
540 * @param a Value.
541 * @param b Value.
542 * @return true if {@code a <= b}
543 */
544 default boolean lte(double a, double b) {
545 return compare(a, b) <= 0;
546 }
547
548 /**
549 * Indicates whether the first argument is strictly greater than the second.
550 *
551 * @param a Value.
552 * @param b Value.
553 * @return true if {@code a > b}
554 */
555 default boolean gt(double a, double b) {
556 return compare(a, b) > 0;
557 }
558
559 /**
560 * Indicates whether the first argument is greater than or considered equal to the second.
561 *
562 * @param a Value.
563 * @param b Value.
564 * @return true if {@code a >= b}
565 */
566 default boolean gte(double a, double b) {
567 return compare(a, b) >= 0;
568 }
569
570 /**
571 * Returns the {@link Math#signum(double) sign} of the argument.
572 * The returned value is
573 * <ul>
574 * <li>{@code -0.0} if {@code a} is considered equal to zero and negatively signed,</li>
575 * <li>{@code +0.0} if {@code a} is considered equal to zero and positively signed,</li>
576 * <li>{@code -1.0} if {@code a} is considered less than zero,</li>
577 * <li>{@code +1.0} if {@code a} is considered greater than zero.</li>
578 * </ul>
579 *
580 * <p>The equality with zero uses the {@link #eqZero(double) eqZero} method.
581 *
582 * @param a Value.
583 * @return the sign (or {@code a} if {@code a == 0} or
584 * {@code a} is NaN).
585 * @see #eqZero(double)
586 */
587 default double signum(double a) {
588 if (a == 0d || Double.isNaN(a)) {
589 return a;
590 }
591 return eqZero(a) ?
592 Math.copySign(0d, a) :
593 Math.copySign(1d, a);
594 }
595
596 /**
597 * Compares two values.
598 * The returned value is
599 * <ul>
600 * <li>{@code 0} if the arguments are considered equal,</li>
601 * <li>{@code -1} if {@code a < b},</li>
602 * <li>{@code +1} if {@code a > b} or if either value is NaN.</li>
603 * </ul>
604 *
605 * @param a Value.
606 * @param b Value.
607 * @return {@code 0} if the values are considered equal, {@code -1}
608 * if the first is smaller than the second, {@code 1} is the first
609 * is larger than the second or either value is NaN.
610 */
611 int compare(double a, double b);
612 }
613 }