001/* 002 * Licensed to the Apache Software Foundation (ASF) under one or more 003 * contributor license agreements. See the NOTICE file distributed with 004 * this work for additional information regarding copyright ownership. 005 * The ASF licenses this file to You under the Apache License, Version 2.0 006 * (the "License"); you may not use this file except in compliance with 007 * the License. You may obtain a copy of the License at 008 * 009 * http://www.apache.org/licenses/LICENSE-2.0 010 * 011 * Unless required by applicable law or agreed to in writing, software 012 * distributed under the License is distributed on an "AS IS" BASIS, 013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 014 * See the License for the specific language governing permissions and 015 * limitations under the License. 016 */ 017 018package org.apache.commons.numbers.core; 019 020import java.math.BigDecimal; 021import java.math.RoundingMode; 022 023/** 024 * Utilities for comparing numbers. 025 */ 026public final class Precision { 027 /** 028 * <p> 029 * Largest double-precision floating-point number such that 030 * {@code 1 + EPSILON} is numerically equal to 1. This value is an upper 031 * bound on the relative error due to rounding real numbers to double 032 * precision floating-point numbers. 033 * </p> 034 * <p> 035 * In IEEE 754 arithmetic, this is 2<sup>-53</sup>. 036 * </p> 037 * 038 * @see <a href="http://en.wikipedia.org/wiki/Machine_epsilon">Machine epsilon</a> 039 */ 040 public static final double EPSILON; 041 042 /** 043 * Safe minimum, such that {@code 1 / SAFE_MIN} does not overflow. 044 * In IEEE 754 arithmetic, this is also the smallest normalized 045 * number 2<sup>-1022</sup>. 046 * 047 * @see Double#MIN_NORMAL 048 */ 049 public static final double SAFE_MIN = Double.MIN_NORMAL; 050 051 /** Exponent offset in IEEE754 representation. */ 052 private static final long EXPONENT_OFFSET = 1023L; 053 054 /** Positive zero. */ 055 private static final double POSITIVE_ZERO = 0d; 056 057 static { 058 /* 059 * This was previously expressed as = 0x1.0p-53 060 * However, OpenJDK (Sparc Solaris) cannot handle such small 061 * constants: MATH-721 062 */ 063 EPSILON = Double.longBitsToDouble((EXPONENT_OFFSET - 53L) << 52); 064 } 065 066 /** 067 * Private constructor. 068 */ 069 private Precision() {} 070 071 /** 072 * Compares two numbers given some amount of allowed error. 073 * The returned value is: 074 * <ul> 075 * <li>zero if considered equal using {@link #equals(double,double,double) equals(x, y, eps)} 076 * <li>negative if not equal and {@code x < y} 077 * <li>positive if not equal and {@code x > y} 078 * </ul> 079 * 080 * <p>NaN values are handled as if using {@link Double#compare(double, double)} where the 081 * returned value is: 082 * <ul> 083 * <li>zero if {@code NaN, NaN} 084 * <li>negative if {@code !NaN, NaN} 085 * <li>positive if {@code NaN, !NaN} 086 * </ul> 087 * 088 * @param x First value. 089 * @param y Second value. 090 * @param eps Allowed error when checking for equality. 091 * @return 0 if the value are considered equal, -1 if the first is smaller than 092 * the second, 1 if the first is larger than the second. 093 * @see #equals(double, double, double) 094 */ 095 public static int compareTo(double x, double y, double eps) { 096 if (equals(x, y, eps)) { 097 return 0; 098 } else if (x < y) { 099 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 {@link RoundingMode#HALF_UP}. 414 * 415 * <p>Note: This method is intended to act on the String representation 416 * of the {@code double} argument. See {@link #round(double, int, RoundingMode)} 417 * for details. 418 * 419 * @param x Value to round. 420 * @param scale Number of digits to the right of the decimal point. 421 * @return the rounded value. 422 * @see #round(double, int, RoundingMode) 423 */ 424 public static double round(double x, int scale) { 425 return round(x, scale, RoundingMode.HALF_UP); 426 } 427 428 /** 429 * Rounds the given value to the specified number of decimal places. 430 * The value is rounded using the given {@link RoundingMode rounding mode}. 431 * 432 * <p>If {@code x} is infinite or {@code NaN}, then the value of {@code x} is 433 * returned unchanged, regardless of the other parameters. 434 * 435 * <p><b>Note</b> 436 * 437 * <p>Rounding of a 64-bit base-2 format {@code double} using a decimal representation 438 * may result in rounding during conversion to and/or from a base-10 representation. 439 * 440 * <p>This method is intended to act on the String representation of the {@code double} 441 * argument, i.e. the closest representable decimal value. The argument is converted to 442 * a String (with possible rounding), rounding is performed on the decimal representation, 443 * and the resulting String is returned as the closest representable {@code double}. 444 * 445 * <p>Conversion from base-2 to base-10 format uses the {@link BigDecimal#valueOf(double)} 446 * method. The alternative would be to use 447 * {@link BigDecimal#BigDecimal(double) new BigDecimal(x)} to construct an exact decimal 448 * representation of the value. 449 * 450 * <p>Conversion from base-10 to base-2 format uses the equivalent of 451 * {@link Double#parseDouble(String)} to create the closest representable {@code double} 452 * to the decimal value. 453 * 454 * <p>The following code demonstrates how to eliminate rounding during conversion to a 455 * decimal format. The return conversion to a {@code double} may be inexact: 456 * <pre> 457 * double rounded = new BigDecimal(x).setScale(scale, roundingMode).doubleValue(); 458 * </pre> 459 * 460 * <p>Acting on the String representation of the {@code double} allows this method to 461 * return expected output when rounding {@code double} representations of decimal text. 462 * <pre> 463 * Precision.round(39.245, 2) == 39.25 464 * Precision.round(30.095, 2) == 30.1 465 * Precision.round(30.645, 2) == 30.65 466 * </pre> 467 * 468 * @param x Value to round. 469 * @param scale Number of digits to the right of the decimal point. 470 * @param roundingMode Rounding mode as defined in {@link BigDecimal}. 471 * @return the rounded value. 472 * @see BigDecimal#doubleValue() 473 * @throws ArithmeticException if {@code roundingMode} is 474 * {@link RoundingMode#UNNECESSARY} and the specified scaling operation 475 * would require rounding. 476 */ 477 public static double round(double x, 478 int scale, 479 RoundingMode roundingMode) { 480 try { 481 final double rounded = (new BigDecimal(Double.toString(x)) 482 .setScale(scale, roundingMode)) 483 .doubleValue(); 484 // MATH-1089: negative values rounded to zero should result in negative zero 485 return rounded == POSITIVE_ZERO ? POSITIVE_ZERO * x : rounded; 486 } catch (NumberFormatException ex) { 487 if (Double.isInfinite(x)) { 488 return x; 489 } 490 return Double.NaN; 491 } 492 } 493 494 /** 495 * Computes a number close to {@code delta} with the property that 496 * {@code (x + delta - x)} is exactly machine-representable. 497 * This is useful when computing numerical derivatives, in order to 498 * reduce roundoff errors. 499 * 500 * @param x Value. 501 * @param delta Offset value. 502 * @return the machine-representable floating number closest to the 503 * difference between {@code x + delta} and {@code x}. 504 */ 505 public static double representableDelta(double x, 506 double delta) { 507 return x + delta - x; 508 } 509 510 /** 511 * Creates a {@link DoubleEquivalence} instance that uses the given epsilon 512 * value for determining equality. 513 * 514 * @param eps Value to use for determining equality. 515 * @return a new instance. 516 */ 517 public static DoubleEquivalence doubleEquivalenceOfEpsilon(final double eps) { 518 if (!Double.isFinite(eps) || 519 eps < 0d) { 520 throw new IllegalArgumentException("Invalid epsilon value: " + eps); 521 } 522 523 return new DoubleEquivalence() { 524 /** Epsilon value. */ 525 private final double epsilon = eps; 526 527 /** {@inheritDoc} */ 528 @Override 529 public int compare(double a, 530 double b) { 531 return Precision.compareTo(a, b, epsilon); 532 } 533 }; 534 } 535 536 /** 537 * Interface containing comparison operations for doubles that allow values 538 * to be <em>considered</em> equal even if they are not exactly equal. 539 * It is intended for comparing outputs of a computation where floating 540 * point errors may have occurred. 541 */ 542 public interface DoubleEquivalence { 543 /** 544 * Indicates whether given values are considered equal to each other. 545 * 546 * @param a Value. 547 * @param b Value. 548 * @return true if the given values are considered equal. 549 */ 550 default boolean eq(double a, double b) { 551 return compare(a, b) == 0; 552 } 553 554 /** 555 * Indicates whether the given value is considered equal to zero. 556 * It is a shortcut for {@code eq(a, 0.0)}. 557 * 558 * @param a Value. 559 * @return true if the argument is considered equal to zero. 560 */ 561 default boolean eqZero(double a) { 562 return eq(a, 0d); 563 } 564 565 /** 566 * Indicates whether the first argument is strictly smaller than the second. 567 * 568 * @param a Value. 569 * @param b Value. 570 * @return true if {@code a < b} 571 */ 572 default boolean lt(double a, double b) { 573 return compare(a, b) < 0; 574 } 575 576 /** 577 * Indicates whether the first argument is smaller or considered equal to the second. 578 * 579 * @param a Value. 580 * @param b Value. 581 * @return true if {@code a <= b} 582 */ 583 default boolean lte(double a, double b) { 584 return compare(a, b) <= 0; 585 } 586 587 /** 588 * Indicates whether the first argument is strictly greater than the second. 589 * 590 * @param a Value. 591 * @param b Value. 592 * @return true if {@code a > b} 593 */ 594 default boolean gt(double a, double b) { 595 return compare(a, b) > 0; 596 } 597 598 /** 599 * Indicates whether the first argument is greater than or considered equal to the second. 600 * 601 * @param a Value. 602 * @param b Value. 603 * @return true if {@code a >= b} 604 */ 605 default boolean gte(double a, double b) { 606 return compare(a, b) >= 0; 607 } 608 609 /** 610 * Returns the {@link Math#signum(double) sign} of the argument. 611 * The returned value is 612 * <ul> 613 * <li>{@code -0.0} if {@code a} is considered equal to zero and negatively signed,</li> 614 * <li>{@code +0.0} if {@code a} is considered equal to zero and positively signed,</li> 615 * <li>{@code -1.0} if {@code a} is considered less than zero,</li> 616 * <li>{@code +1.0} if {@code a} is considered greater than zero.</li> 617 * </ul> 618 * 619 * <p>The equality with zero uses the {@link #eqZero(double) eqZero} method. 620 * 621 * @param a Value. 622 * @return the sign (or {@code a} if {@code a == 0} or 623 * {@code a} is NaN). 624 * @see #eqZero(double) 625 */ 626 default double signum(double a) { 627 if (a == 0d || Double.isNaN(a)) { 628 return a; 629 } 630 return eqZero(a) ? 631 Math.copySign(0d, a) : 632 Math.copySign(1d, a); 633 } 634 635 /** 636 * Compares two values. 637 * The returned value is 638 * <ul> 639 * <li>{@code 0} if the arguments are considered equal,</li> 640 * <li>{@code -1} if {@code a < b},</li> 641 * <li>{@code +1} if {@code a > b} or if either value is NaN.</li> 642 * </ul> 643 * 644 * @param a Value. 645 * @param b Value. 646 * @return {@code 0} if the values are considered equal, {@code -1} 647 * if the first is smaller than the second, {@code 1} is the first 648 * is larger than the second or either value is NaN. 649 */ 650 int compare(double a, double b); 651 } 652}