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 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}