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 public static final double SAFE_MIN; 048 049 /** Exponent offset in IEEE754 representation. */ 050 private static final long EXPONENT_OFFSET = 1023L; 051 052 /** Offset to order signed double numbers lexicographically. */ 053 private static final long SGN_MASK = 0x8000000000000000L; 054 /** Offset to order signed double numbers lexicographically. */ 055 private static final int SGN_MASK_FLOAT = 0x80000000; 056 /** Positive zero. */ 057 private static final double POSITIVE_ZERO = 0d; 058 /** Positive zero bits. */ 059 private static final long POSITIVE_ZERO_DOUBLE_BITS = Double.doubleToRawLongBits(+0.0); 060 /** Negative zero bits. */ 061 private static final long NEGATIVE_ZERO_DOUBLE_BITS = Double.doubleToRawLongBits(-0.0); 062 /** Positive zero bits. */ 063 private static final int POSITIVE_ZERO_FLOAT_BITS = Float.floatToRawIntBits(+0.0f); 064 /** Negative zero bits. */ 065 private static final int NEGATIVE_ZERO_FLOAT_BITS = Float.floatToRawIntBits(-0.0f); 066 067 static { 068 /* 069 * This was previously expressed as = 0x1.0p-53 070 * However, OpenJDK (Sparc Solaris) cannot handle such small 071 * constants: MATH-721 072 */ 073 EPSILON = Double.longBitsToDouble((EXPONENT_OFFSET - 53L) << 52); 074 075 /* 076 * This was previously expressed as = 0x1.0p-1022 077 * However, OpenJDK (Sparc Solaris) cannot handle such small 078 * constants: MATH-721 079 */ 080 SAFE_MIN = Double.longBitsToDouble((EXPONENT_OFFSET - 1022L) << 52); 081 } 082 083 /** 084 * Private constructor. 085 */ 086 private Precision() {} 087 088 /** 089 * Compares two numbers given some amount of allowed error. 090 * The returned value is 091 * <ul> 092 * <li> 093 * 0 if {@link #equals(double,double,double) equals(x, y, eps)}, 094 * </li> 095 * <li> 096 * negative if !{@link #equals(double,double,double) equals(x, y, eps)} and {@code x < y}, 097 * </li> 098 * <li> 099 * positive if !{@link #equals(double,double,double) equals(x, y, eps)} and {@code x > y} or 100 * either argument is {@code NaN}. 101 * </li> 102 * </ul> 103 * 104 * @param x First value. 105 * @param y Second value. 106 * @param eps Allowed error when checking for equality. 107 * @return 0 if the value are considered equal, -1 if the first is smaller than 108 * the second, 1 is the first is larger than the second. 109 */ 110 public static int compareTo(double x, double y, double eps) { 111 if (equals(x, y, eps)) { 112 return 0; 113 } else if (x < y) { 114 return -1; 115 } else if (x > y) { 116 return 1; 117 } 118 // NaN input. 119 return Double.compare(x, y); 120 } 121 122 /** 123 * Compares two numbers given some amount of allowed error. 124 * Two float numbers are considered equal if there are {@code (maxUlps - 1)} 125 * (or fewer) floating point numbers between them, i.e. two adjacent floating 126 * point numbers are considered equal. 127 * Adapted from 128 * <a href="http://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/"> 129 * Bruce Dawson</a>. Returns {@code false} if either of the arguments is NaN. 130 * The returned value is 131 * <ul> 132 * <li> 133 * zero if {@link #equals(double,double,int) equals(x, y, maxUlps)}, 134 * </li> 135 * <li> 136 * negative if !{@link #equals(double,double,int) equals(x, y, maxUlps)} and {@code x < y}, 137 * </li> 138 * <li> 139 * positive if !{@link #equals(double,double,int) equals(x, y, maxUlps)} and {@code x > y} 140 * or either argument is {@code NaN}. 141 * </li> 142 * </ul> 143 * 144 * @param x First value. 145 * @param y Second value. 146 * @param maxUlps {@code (maxUlps - 1)} is the number of floating point 147 * values between {@code x} and {@code y}. 148 * @return 0 if the value are considered equal, -1 if the first is smaller than 149 * the second, 1 is the first is larger than the second. 150 */ 151 public static int compareTo(final double x, final double y, final int maxUlps) { 152 if (equals(x, y, maxUlps)) { 153 return 0; 154 } else if (x < y) { 155 return -1; 156 } 157 return 1; 158 } 159 160 /** 161 * Returns true iff they are equal as defined by 162 * {@link #equals(float,float,int) equals(x, y, 1)}. 163 * 164 * @param x first value 165 * @param y second value 166 * @return {@code true} if the values are equal. 167 */ 168 public static boolean equals(float x, float y) { 169 return equals(x, y, 1); 170 } 171 172 /** 173 * Returns true if both arguments are NaN or they are 174 * equal as defined by {@link #equals(float,float) equals(x, y, 1)}. 175 * 176 * @param x first value 177 * @param y second value 178 * @return {@code true} if the values are equal or both are NaN. 179 */ 180 public static boolean equalsIncludingNaN(float x, float y) { 181 final boolean xIsNan = Float.isNaN(x); 182 final boolean yIsNan = Float.isNaN(y); 183 // Combine the booleans with bitwise OR 184 return (xIsNan | yIsNan) ? 185 !(xIsNan ^ yIsNan) : 186 equals(x, y, 1); 187 } 188 189 /** 190 * Returns {@code true} if there is no float value strictly between the 191 * arguments or the difference between them is within the range of allowed 192 * error (inclusive). Returns {@code false} if either of the arguments 193 * is NaN. 194 * 195 * @param x first value 196 * @param y second value 197 * @param eps the amount of absolute error to allow. 198 * @return {@code true} if the values are equal or within range of each other. 199 */ 200 public static boolean equals(float x, float y, float eps) { 201 return equals(x, y, 1) || Math.abs(y - x) <= eps; 202 } 203 204 /** 205 * Returns true if the arguments are both NaN, there are no float value strictly 206 * between the arguments or the difference between them is within the range of allowed 207 * error (inclusive). 208 * 209 * @param x first value 210 * @param y second value 211 * @param eps the amount of absolute error to allow. 212 * @return {@code true} if the values are equal or within range of each other, 213 * or both are NaN. 214 */ 215 public static boolean equalsIncludingNaN(float x, float y, float eps) { 216 return equalsIncludingNaN(x, y, 1) || (Math.abs(y - x) <= eps); 217 } 218 219 /** 220 * Returns true if the arguments are equal or within the range of allowed 221 * error (inclusive). Returns {@code false} if either of the arguments is NaN. 222 * <p> 223 * Two double numbers are considered equal if there are {@code (maxUlps - 1)} 224 * (or fewer) floating point numbers between them, i.e. two adjacent 225 * floating point numbers are considered equal. 226 * </p> 227 * <p> 228 * Adapted from <a 229 * href="http://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/"> 230 * Bruce Dawson</a>. 231 * </p> 232 * 233 * @param x first value 234 * @param y second value 235 * @param maxUlps {@code (maxUlps - 1)} is the number of floating point 236 * values between {@code x} and {@code y}. 237 * @return {@code true} if there are fewer than {@code maxUlps} floating 238 * point values between {@code x} and {@code y}. 239 */ 240 public static boolean equals(final float x, final float y, final int maxUlps) { 241 242 final int xInt = Float.floatToRawIntBits(x); 243 final int yInt = Float.floatToRawIntBits(y); 244 245 final boolean isEqual; 246 if (((xInt ^ yInt) & SGN_MASK_FLOAT) == 0) { 247 // number have same sign, there is no risk of overflow 248 isEqual = Math.abs(xInt - yInt) <= maxUlps; 249 } else { 250 // number have opposite signs, take care of overflow 251 final int deltaPlus; 252 final int deltaMinus; 253 if (xInt < yInt) { 254 deltaPlus = yInt - POSITIVE_ZERO_FLOAT_BITS; 255 deltaMinus = xInt - NEGATIVE_ZERO_FLOAT_BITS; 256 } else { 257 deltaPlus = xInt - POSITIVE_ZERO_FLOAT_BITS; 258 deltaMinus = yInt - NEGATIVE_ZERO_FLOAT_BITS; 259 } 260 261 if (deltaPlus > maxUlps) { 262 isEqual = false; 263 } else { 264 isEqual = deltaMinus <= (maxUlps - deltaPlus); 265 } 266 267 } 268 269 return isEqual && !Float.isNaN(x) && !Float.isNaN(y); 270 } 271 272 /** 273 * Returns true if both arguments are NaN or if they are equal as defined 274 * by {@link #equals(float,float,int) equals(x, y, maxUlps)}. 275 * 276 * @param x first value 277 * @param y second value 278 * @param maxUlps {@code (maxUlps - 1)} is the number of floating point 279 * values between {@code x} and {@code y}. 280 * @return {@code true} if both arguments are NaN or if there are less than 281 * {@code maxUlps} floating point values between {@code x} and {@code y}. 282 */ 283 public static boolean equalsIncludingNaN(float x, float y, int maxUlps) { 284 final boolean xIsNan = Float.isNaN(x); 285 final boolean yIsNan = Float.isNaN(y); 286 // Combine the booleans with bitwise OR 287 return (xIsNan | yIsNan) ? 288 !(xIsNan ^ yIsNan) : 289 equals(x, y, maxUlps); 290 } 291 292 /** 293 * Returns true iff they are equal as defined by 294 * {@link #equals(double,double,int) equals(x, y, 1)}. 295 * 296 * @param x first value 297 * @param y second value 298 * @return {@code true} if the values are equal. 299 */ 300 public static boolean equals(double x, double y) { 301 return equals(x, y, 1); 302 } 303 304 /** 305 * Returns true if the arguments are both NaN or they are 306 * equal as defined by {@link #equals(double,double) equals(x, y, 1)}. 307 * 308 * @param x first value 309 * @param y second value 310 * @return {@code true} if the values are equal or both are NaN. 311 */ 312 public static boolean equalsIncludingNaN(double x, double y) { 313 final boolean xIsNan = Double.isNaN(x); 314 final boolean yIsNan = Double.isNaN(y); 315 // Combine the booleans with bitwise OR 316 return (xIsNan | yIsNan) ? 317 !(xIsNan ^ yIsNan) : 318 equals(x, y, 1); 319 } 320 321 /** 322 * Returns {@code true} if there is no double value strictly between the 323 * arguments or the difference between them is within the range of allowed 324 * error (inclusive). Returns {@code false} if either of the arguments 325 * is NaN. 326 * 327 * @param x First value. 328 * @param y Second value. 329 * @param eps Amount of allowed absolute error. 330 * @return {@code true} if the values are equal or within range of each other. 331 */ 332 public static boolean equals(double x, double y, double eps) { 333 return equals(x, y, 1) || Math.abs(y - x) <= eps; 334 } 335 336 /** 337 * Returns {@code true} if there is no double value strictly between the 338 * arguments or the relative difference between them is less than or equal 339 * to the given tolerance. Returns {@code false} if either of the arguments 340 * is NaN. 341 * 342 * @param x First value. 343 * @param y Second value. 344 * @param eps Amount of allowed relative error. 345 * @return {@code true} if the values are two adjacent floating point 346 * numbers or they are within range of each other. 347 */ 348 public static boolean equalsWithRelativeTolerance(double x, double y, double eps) { 349 if (equals(x, y, 1)) { 350 return true; 351 } 352 353 final double absoluteMax = Math.max(Math.abs(x), Math.abs(y)); 354 final double relativeDifference = Math.abs((x - y) / absoluteMax); 355 356 return relativeDifference <= eps; 357 } 358 359 /** 360 * Returns true if the arguments are both NaN, there are no double value strictly 361 * between the arguments or the difference between them is within the range of allowed 362 * error (inclusive). 363 * 364 * @param x first value 365 * @param y second value 366 * @param eps the amount of absolute error to allow. 367 * @return {@code true} if the values are equal or within range of each other, 368 * or both are NaN. 369 */ 370 public static boolean equalsIncludingNaN(double x, double y, double eps) { 371 return equalsIncludingNaN(x, y) || (Math.abs(y - x) <= eps); 372 } 373 374 /** 375 * Returns true if the arguments are equal or within the range of allowed 376 * error (inclusive). Returns {@code false} if either of the arguments is NaN. 377 * <p> 378 * Two float numbers are considered equal if there are {@code (maxUlps - 1)} 379 * (or fewer) floating point numbers between them, i.e. two adjacent 380 * floating point numbers are considered equal. 381 * </p> 382 * <p> 383 * Adapted from <a 384 * href="http://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/"> 385 * Bruce Dawson</a>. 386 * </p> 387 * 388 * @param x first value 389 * @param y second value 390 * @param maxUlps {@code (maxUlps - 1)} is the number of floating point 391 * values between {@code x} and {@code y}. 392 * @return {@code true} if there are fewer than {@code maxUlps} floating 393 * point values between {@code x} and {@code y}. 394 */ 395 public static boolean equals(final double x, final double y, final int maxUlps) { 396 397 final long xInt = Double.doubleToRawLongBits(x); 398 final long yInt = Double.doubleToRawLongBits(y); 399 400 final boolean isEqual; 401 if (((xInt ^ yInt) & SGN_MASK) == 0L) { 402 // number have same sign, there is no risk of overflow 403 isEqual = Math.abs(xInt - yInt) <= maxUlps; 404 } else { 405 // number have opposite signs, take care of overflow 406 final long deltaPlus; 407 final long deltaMinus; 408 if (xInt < yInt) { 409 deltaPlus = yInt - POSITIVE_ZERO_DOUBLE_BITS; 410 deltaMinus = xInt - NEGATIVE_ZERO_DOUBLE_BITS; 411 } else { 412 deltaPlus = xInt - POSITIVE_ZERO_DOUBLE_BITS; 413 deltaMinus = yInt - NEGATIVE_ZERO_DOUBLE_BITS; 414 } 415 416 if (deltaPlus > maxUlps) { 417 isEqual = false; 418 } else { 419 isEqual = deltaMinus <= (maxUlps - deltaPlus); 420 } 421 422 } 423 424 return isEqual && !Double.isNaN(x) && !Double.isNaN(y); 425 } 426 427 /** 428 * Returns true if both arguments are NaN or if they are equal as defined 429 * by {@link #equals(double,double,int) equals(x, y, maxUlps)}. 430 * 431 * @param x first value 432 * @param y second value 433 * @param maxUlps {@code (maxUlps - 1)} is the number of floating point 434 * values between {@code x} and {@code y}. 435 * @return {@code true} if both arguments are NaN or if there are less than 436 * {@code maxUlps} floating point values between {@code x} and {@code y}. 437 */ 438 public static boolean equalsIncludingNaN(double x, double y, int maxUlps) { 439 final boolean xIsNan = Double.isNaN(x); 440 final boolean yIsNan = Double.isNaN(y); 441 // Combine the booleans with bitwise OR 442 return (xIsNan | yIsNan) ? 443 !(xIsNan ^ yIsNan) : 444 equals(x, y, maxUlps); 445 } 446 447 /** 448 * Rounds the given value to the specified number of decimal places. 449 * The value is rounded using the {@link RoundingMode#HALF_UP} method. 450 * 451 * @param x Value to round. 452 * @param scale Number of digits to the right of the decimal point. 453 * @return the rounded value. 454 */ 455 public static double round(double x, int scale) { 456 return round(x, scale, RoundingMode.HALF_UP); 457 } 458 459 /** 460 * Rounds the given value to the specified number of decimal places. 461 * The value is rounded using the given method which is any method defined 462 * in {@link BigDecimal}. 463 * If {@code x} is infinite or {@code NaN}, then the value of {@code x} is 464 * returned unchanged, regardless of the other parameters. 465 * 466 * @param x Value to round. 467 * @param scale Number of digits to the right of the decimal point. 468 * @param roundingMethod Rounding method as defined in {@link BigDecimal}. 469 * @return the rounded value. 470 * @throws ArithmeticException if {@code roundingMethod} is 471 * {@link RoundingMode#UNNECESSARY} and the specified scaling operation 472 * would require rounding. 473 */ 474 public static double round(double x, 475 int scale, 476 RoundingMode roundingMethod) { 477 try { 478 final double rounded = (new BigDecimal(Double.toString(x)) 479 .setScale(scale, roundingMethod)) 480 .doubleValue(); 481 // MATH-1089: negative values rounded to zero should result in negative zero 482 return rounded == POSITIVE_ZERO ? POSITIVE_ZERO * x : rounded; 483 } catch (NumberFormatException ex) { 484 if (Double.isInfinite(x)) { 485 return x; 486 } 487 return Double.NaN; 488 } 489 } 490 491 /** 492 * Computes a number close to {@code delta} with the property that 493 * {@code (x + delta - x)} is exactly machine-representable. 494 * This is useful when computing numerical derivatives, in order to 495 * reduce roundoff errors. 496 * 497 * @param x Value. 498 * @param delta Offset value. 499 * @return the machine-representable floating number closest to the 500 * difference between {@code x + delta} and {@code x}. 501 */ 502 public static double representableDelta(double x, 503 double delta) { 504 return x + delta - x; 505 } 506 507 /** 508 * Creates a {@link DoubleEquivalence} instance that uses the given epsilon 509 * value for determining equality. 510 * 511 * @param eps Value to use for determining equality. 512 * @return a new instance. 513 */ 514 public static DoubleEquivalence doubleEquivalenceOfEpsilon(final double eps) { 515 if (!Double.isFinite(eps) || 516 eps < 0d) { 517 throw new IllegalArgumentException("Invalid epsilon value: " + eps); 518 } 519 520 return new DoubleEquivalence() { 521 /** Epsilon value. */ 522 private final double epsilon = eps; 523 524 /** {@inheritDoc} */ 525 @Override 526 public int compare(double a, 527 double b) { 528 return Precision.compareTo(a, b, epsilon); 529 } 530 }; 531 } 532 533 /** 534 * Interface containing comparison operations for doubles that allow values 535 * to be <em>considered</em> equal even if they are not exactly equal. 536 * It is intended for comparing outputs of a computation where floating 537 * point errors may have occurred. 538 */ 539 public interface DoubleEquivalence { 540 /** 541 * Indicates whether given values are considered equal to each other. 542 * 543 * @param a Value. 544 * @param b Value. 545 * @return true if the given values are considered equal. 546 */ 547 default boolean eq(double a, double b) { 548 return compare(a, b) == 0; 549 } 550 551 /** 552 * Indicates whether the given value is considered equal to zero. 553 * It is a shortcut for {@code eq(a, 0.0)}. 554 * 555 * @param a Value. 556 * @return true if the argument is considered equal to zero. 557 */ 558 default boolean eqZero(double a) { 559 return eq(a, 0d); 560 } 561 562 /** 563 * Indicates whether the first argument is strictly smaller than the second. 564 * 565 * @param a Value. 566 * @param b Value. 567 * @return true if {@code a < b} 568 */ 569 default boolean lt(double a, double b) { 570 return compare(a, b) < 0; 571 } 572 573 /** 574 * Indicates whether the first argument is smaller or considered equal to the second. 575 * 576 * @param a Value. 577 * @param b Value. 578 * @return true if {@code a <= b} 579 */ 580 default boolean lte(double a, double b) { 581 return compare(a, b) <= 0; 582 } 583 584 /** 585 * Indicates whether the first argument is strictly greater than the second. 586 * 587 * @param a Value. 588 * @param b Value. 589 * @return true if {@code a > b} 590 */ 591 default boolean gt(double a, double b) { 592 return compare(a, b) > 0; 593 } 594 595 /** 596 * Indicates whether the first argument is greater than or considered equal to the second. 597 * 598 * @param a Value. 599 * @param b Value. 600 * @return true if {@code a >= b} 601 */ 602 default boolean gte(double a, double b) { 603 return compare(a, b) >= 0; 604 } 605 606 /** 607 * Returns the {@link Math#signum(double) sign} of the argument. 608 * The returned value is 609 * <ul> 610 * <li>{@code -0.0} if {@code a} is considered equal to zero and negatively signed,</li> 611 * <li>{@code +0.0} if {@code a} is considered equal to zero and positively signed,</li> 612 * <li>{@code -1.0} if {@code a} is considered less than zero,</li> 613 * <li>{@code +1.0} if {@code a} is considered greater than zero.</li> 614 * </ul> 615 * 616 * <p>The equality with zero uses the {@link #eqZero(double) eqZero} method. 617 * 618 * @param a Value. 619 * @return the sign (or {@code a} if {@code a == 0} or 620 * {@code a} is NaN). 621 * @see #eqZero(double) 622 */ 623 default double signum(double a) { 624 if (a == 0d || Double.isNaN(a)) { 625 return a; 626 } 627 return eqZero(a) ? 628 Math.copySign(0d, a) : 629 Math.copySign(1d, a); 630 } 631 632 /** 633 * Compares two values. 634 * The returned value is 635 * <ul> 636 * <li>{@code 0} if the arguments are considered equal,</li> 637 * <li>{@code -1} if {@code a < b},</li> 638 * <li>{@code +1} if {@code a > b} or if either value is NaN.</li> 639 * </ul> 640 * 641 * @param a Value. 642 * @param b Value. 643 * @return {@code 0} if the values are considered equal, {@code -1} 644 * if the first is smaller than the second, {@code 1} is the first 645 * is larger than the second or either value is NaN. 646 */ 647 int compare(double a, double b); 648 } 649}