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.math3.util; 019 020import java.math.BigDecimal; 021 022import org.apache.commons.math3.exception.MathArithmeticException; 023import org.apache.commons.math3.exception.MathIllegalArgumentException; 024import org.apache.commons.math3.exception.util.LocalizedFormats; 025 026/** 027 * Utilities for comparing numbers. 028 * 029 * @since 3.0 030 */ 031public class Precision { 032 /** 033 * <p> 034 * Largest double-precision floating-point number such that 035 * {@code 1 + EPSILON} is numerically equal to 1. This value is an upper 036 * bound on the relative error due to rounding real numbers to double 037 * precision floating-point numbers. 038 * </p> 039 * <p> 040 * In IEEE 754 arithmetic, this is 2<sup>-53</sup>. 041 * </p> 042 * 043 * @see <a href="http://en.wikipedia.org/wiki/Machine_epsilon">Machine epsilon</a> 044 */ 045 public static final double EPSILON; 046 047 /** 048 * Safe minimum, such that {@code 1 / SAFE_MIN} does not overflow. 049 * <br/> 050 * In IEEE 754 arithmetic, this is also the smallest normalized 051 * number 2<sup>-1022</sup>. 052 */ 053 public static final double SAFE_MIN; 054 055 /** Exponent offset in IEEE754 representation. */ 056 private static final long EXPONENT_OFFSET = 1023l; 057 058 /** Offset to order signed double numbers lexicographically. */ 059 private static final long SGN_MASK = 0x8000000000000000L; 060 /** Offset to order signed double numbers lexicographically. */ 061 private static final int SGN_MASK_FLOAT = 0x80000000; 062 /** Positive zero. */ 063 private static final double POSITIVE_ZERO = 0d; 064 /** Positive zero bits. */ 065 private static final long POSITIVE_ZERO_DOUBLE_BITS = Double.doubleToRawLongBits(+0.0); 066 /** Negative zero bits. */ 067 private static final long NEGATIVE_ZERO_DOUBLE_BITS = Double.doubleToRawLongBits(-0.0); 068 /** Positive zero bits. */ 069 private static final int POSITIVE_ZERO_FLOAT_BITS = Float.floatToRawIntBits(+0.0f); 070 /** Negative zero bits. */ 071 private static final int NEGATIVE_ZERO_FLOAT_BITS = Float.floatToRawIntBits(-0.0f); 072 073 static { 074 /* 075 * This was previously expressed as = 0x1.0p-53; 076 * However, OpenJDK (Sparc Solaris) cannot handle such small 077 * constants: MATH-721 078 */ 079 EPSILON = Double.longBitsToDouble((EXPONENT_OFFSET - 53l) << 52); 080 081 /* 082 * This was previously expressed as = 0x1.0p-1022; 083 * However, OpenJDK (Sparc Solaris) cannot handle such small 084 * constants: MATH-721 085 */ 086 SAFE_MIN = Double.longBitsToDouble((EXPONENT_OFFSET - 1022l) << 52); 087 } 088 089 /** 090 * Private constructor. 091 */ 092 private Precision() {} 093 094 /** 095 * Compares two numbers given some amount of allowed error. 096 * 097 * @param x the first number 098 * @param y the second number 099 * @param eps the amount of error to allow when checking for equality 100 * @return <ul><li>0 if {@link #equals(double, double, double) equals(x, y, eps)}</li> 101 * <li>< 0 if !{@link #equals(double, double, double) equals(x, y, eps)} && x < y</li> 102 * <li>> 0 if !{@link #equals(double, double, double) equals(x, y, eps)} && x > y or 103 * either argument is NaN</li></ul> 104 */ 105 public static int compareTo(double x, double y, double eps) { 106 if (equals(x, y, eps)) { 107 return 0; 108 } else if (x < y) { 109 return -1; 110 } 111 return 1; 112 } 113 114 /** 115 * Compares two numbers given some amount of allowed error. 116 * Two float numbers are considered equal if there are {@code (maxUlps - 1)} 117 * (or fewer) floating point numbers between them, i.e. two adjacent floating 118 * point numbers are considered equal. 119 * Adapted from <a 120 * href="http://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/"> 121 * Bruce Dawson</a>. Returns {@code false} if either of the arguments is NaN. 122 * 123 * @param x first value 124 * @param y second value 125 * @param maxUlps {@code (maxUlps - 1)} is the number of floating point 126 * values between {@code x} and {@code y}. 127 * @return <ul><li>0 if {@link #equals(double, double, int) equals(x, y, maxUlps)}</li> 128 * <li>< 0 if !{@link #equals(double, double, int) equals(x, y, maxUlps)} && x < y</li> 129 * <li>> 0 if !{@link #equals(double, double, int) equals(x, y, maxUlps)} && x > y 130 * or either argument is NaN</li></ul> 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 } 138 return 1; 139 } 140 141 /** 142 * Returns true iff they are equal as defined by 143 * {@link #equals(float,float,int) equals(x, y, 1)}. 144 * 145 * @param x first value 146 * @param y second value 147 * @return {@code true} if the values are equal. 148 */ 149 public static boolean equals(float x, float y) { 150 return equals(x, y, 1); 151 } 152 153 /** 154 * Returns true if both arguments are NaN or they are 155 * equal as defined by {@link #equals(float,float) equals(x, y, 1)}. 156 * 157 * @param x first value 158 * @param y second value 159 * @return {@code true} if the values are equal or both are NaN. 160 * @since 2.2 161 */ 162 public static boolean equalsIncludingNaN(float x, float y) { 163 return (x != x || y != y) ? !(x != x ^ y != y) : equals(x, y, 1); 164 } 165 166 /** 167 * Returns true if the arguments are equal or within the range of allowed 168 * error (inclusive). Returns {@code false} if either of the arguments 169 * is NaN. 170 * 171 * @param x first value 172 * @param y second value 173 * @param eps the amount of absolute error to allow. 174 * @return {@code true} if the values are equal or within range of each other. 175 * @since 2.2 176 */ 177 public static boolean equals(float x, float y, float eps) { 178 return equals(x, y, 1) || FastMath.abs(y - x) <= eps; 179 } 180 181 /** 182 * Returns true if the arguments are both NaN, are equal, or are within the range 183 * of allowed error (inclusive). 184 * 185 * @param x first value 186 * @param y second value 187 * @param eps the amount of absolute error to allow. 188 * @return {@code true} if the values are equal or within range of each other, 189 * or both are NaN. 190 * @since 2.2 191 */ 192 public static boolean equalsIncludingNaN(float x, float y, float eps) { 193 return equalsIncludingNaN(x, y) || (FastMath.abs(y - x) <= eps); 194 } 195 196 /** 197 * Returns true if the arguments are equal or within the range of allowed 198 * error (inclusive). 199 * Two float numbers are considered equal if there are {@code (maxUlps - 1)} 200 * (or fewer) floating point numbers between them, i.e. two adjacent floating 201 * point numbers are considered equal. 202 * Adapted from <a 203 * href="http://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/"> 204 * Bruce Dawson</a>. Returns {@code false} if either of the arguments is NaN. 205 * 206 * @param x first value 207 * @param y second value 208 * @param maxUlps {@code (maxUlps - 1)} is the number of floating point 209 * values between {@code x} and {@code y}. 210 * @return {@code true} if there are fewer than {@code maxUlps} floating 211 * point values between {@code x} and {@code y}. 212 * @since 2.2 213 */ 214 public static boolean equals(final float x, final float y, final int maxUlps) { 215 216 final int xInt = Float.floatToRawIntBits(x); 217 final int yInt = Float.floatToRawIntBits(y); 218 219 final boolean isEqual; 220 if (((xInt ^ yInt) & SGN_MASK_FLOAT) == 0) { 221 // number have same sign, there is no risk of overflow 222 isEqual = FastMath.abs(xInt - yInt) <= maxUlps; 223 } else { 224 // number have opposite signs, take care of overflow 225 final int deltaPlus; 226 final int deltaMinus; 227 if (xInt < yInt) { 228 deltaPlus = yInt - POSITIVE_ZERO_FLOAT_BITS; 229 deltaMinus = xInt - NEGATIVE_ZERO_FLOAT_BITS; 230 } else { 231 deltaPlus = xInt - POSITIVE_ZERO_FLOAT_BITS; 232 deltaMinus = yInt - NEGATIVE_ZERO_FLOAT_BITS; 233 } 234 235 if (deltaPlus > maxUlps) { 236 isEqual = false; 237 } else { 238 isEqual = deltaMinus <= (maxUlps - deltaPlus); 239 } 240 241 } 242 243 return isEqual && !Float.isNaN(x) && !Float.isNaN(y); 244 245 } 246 247 /** 248 * Returns true if the arguments are both NaN or if they are equal as defined 249 * by {@link #equals(float,float,int) equals(x, y, maxUlps)}. 250 * 251 * @param x first value 252 * @param y second value 253 * @param maxUlps {@code (maxUlps - 1)} is the number of floating point 254 * values between {@code x} and {@code y}. 255 * @return {@code true} if both arguments are NaN or if there are less than 256 * {@code maxUlps} floating point values between {@code x} and {@code y}. 257 * @since 2.2 258 */ 259 public static boolean equalsIncludingNaN(float x, float y, int maxUlps) { 260 return (x != x || y != y) ? !(x != x ^ y != y) : equals(x, y, maxUlps); 261 } 262 263 /** 264 * Returns true iff they are equal as defined by 265 * {@link #equals(double,double,int) equals(x, y, 1)}. 266 * 267 * @param x first value 268 * @param y second value 269 * @return {@code true} if the values are equal. 270 */ 271 public static boolean equals(double x, double y) { 272 return equals(x, y, 1); 273 } 274 275 /** 276 * Returns true if the arguments are both NaN or they are 277 * equal as defined by {@link #equals(double,double) equals(x, y, 1)}. 278 * 279 * @param x first value 280 * @param y second value 281 * @return {@code true} if the values are equal or both are NaN. 282 * @since 2.2 283 */ 284 public static boolean equalsIncludingNaN(double x, double y) { 285 return (x != x || y != y) ? !(x != x ^ y != y) : equals(x, y, 1); 286 } 287 288 /** 289 * Returns {@code true} if there is no double value strictly between the 290 * arguments or the difference between them is within the range of allowed 291 * error (inclusive). Returns {@code false} if either of the arguments 292 * is NaN. 293 * 294 * @param x First value. 295 * @param y Second value. 296 * @param eps Amount of allowed absolute error. 297 * @return {@code true} if the values are two adjacent floating point 298 * numbers or they are within range of each other. 299 */ 300 public static boolean equals(double x, double y, double eps) { 301 return equals(x, y, 1) || FastMath.abs(y - x) <= eps; 302 } 303 304 /** 305 * Returns {@code true} if there is no double value strictly between the 306 * arguments or the relative difference between them is less than or equal 307 * to the given tolerance. Returns {@code false} if either of the arguments 308 * is NaN. 309 * 310 * @param x First value. 311 * @param y Second value. 312 * @param eps Amount of allowed relative error. 313 * @return {@code true} if the values are two adjacent floating point 314 * numbers or they are within range of each other. 315 * @since 3.1 316 */ 317 public static boolean equalsWithRelativeTolerance(double x, double y, double eps) { 318 if (equals(x, y, 1)) { 319 return true; 320 } 321 322 final double absoluteMax = FastMath.max(FastMath.abs(x), FastMath.abs(y)); 323 final double relativeDifference = FastMath.abs((x - y) / absoluteMax); 324 325 return relativeDifference <= eps; 326 } 327 328 /** 329 * Returns true if the arguments are both NaN, are equal or are within the range 330 * of allowed error (inclusive). 331 * 332 * @param x first value 333 * @param y second value 334 * @param eps the amount of absolute error to allow. 335 * @return {@code true} if the values are equal or within range of each other, 336 * or both are NaN. 337 * @since 2.2 338 */ 339 public static boolean equalsIncludingNaN(double x, double y, double eps) { 340 return equalsIncludingNaN(x, y) || (FastMath.abs(y - x) <= eps); 341 } 342 343 /** 344 * Returns true if the arguments are equal or within the range of allowed 345 * error (inclusive). 346 * <p> 347 * Two float numbers are considered equal if there are {@code (maxUlps - 1)} 348 * (or fewer) floating point numbers between them, i.e. two adjacent 349 * floating point numbers are considered equal. 350 * </p> 351 * <p> 352 * Adapted from <a 353 * href="http://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/"> 354 * Bruce Dawson</a>. Returns {@code false} if either of the arguments is NaN. 355 * </p> 356 * 357 * @param x first value 358 * @param y second value 359 * @param maxUlps {@code (maxUlps - 1)} is the number of floating point 360 * values between {@code x} and {@code y}. 361 * @return {@code true} if there are fewer than {@code maxUlps} floating 362 * point values between {@code x} and {@code y}. 363 */ 364 public static boolean equals(final double x, final double y, final int maxUlps) { 365 366 final long xInt = Double.doubleToRawLongBits(x); 367 final long yInt = Double.doubleToRawLongBits(y); 368 369 final boolean isEqual; 370 if (((xInt ^ yInt) & SGN_MASK) == 0l) { 371 // number have same sign, there is no risk of overflow 372 isEqual = FastMath.abs(xInt - yInt) <= maxUlps; 373 } else { 374 // number have opposite signs, take care of overflow 375 final long deltaPlus; 376 final long deltaMinus; 377 if (xInt < yInt) { 378 deltaPlus = yInt - POSITIVE_ZERO_DOUBLE_BITS; 379 deltaMinus = xInt - NEGATIVE_ZERO_DOUBLE_BITS; 380 } else { 381 deltaPlus = xInt - POSITIVE_ZERO_DOUBLE_BITS; 382 deltaMinus = yInt - NEGATIVE_ZERO_DOUBLE_BITS; 383 } 384 385 if (deltaPlus > maxUlps) { 386 isEqual = false; 387 } else { 388 isEqual = deltaMinus <= (maxUlps - deltaPlus); 389 } 390 391 } 392 393 return isEqual && !Double.isNaN(x) && !Double.isNaN(y); 394 395 } 396 397 /** 398 * Returns true if both arguments are NaN or if they are equal as defined 399 * by {@link #equals(double,double,int) equals(x, y, maxUlps)}. 400 * 401 * @param x first value 402 * @param y second value 403 * @param maxUlps {@code (maxUlps - 1)} is the number of floating point 404 * values between {@code x} and {@code y}. 405 * @return {@code true} if both arguments are NaN or if there are less than 406 * {@code maxUlps} floating point values between {@code x} and {@code y}. 407 * @since 2.2 408 */ 409 public static boolean equalsIncludingNaN(double x, double y, int maxUlps) { 410 return (x != x || y != y) ? !(x != x ^ y != y) : equals(x, y, maxUlps); 411 } 412 413 /** 414 * Rounds the given value to the specified number of decimal places. 415 * The value is rounded using the {@link BigDecimal#ROUND_HALF_UP} method. 416 * 417 * @param x Value to round. 418 * @param scale Number of digits to the right of the decimal point. 419 * @return the rounded value. 420 * @since 1.1 (previously in {@code MathUtils}, moved as of version 3.0) 421 */ 422 public static double round(double x, int scale) { 423 return round(x, scale, BigDecimal.ROUND_HALF_UP); 424 } 425 426 /** 427 * Rounds the given value to the specified number of decimal places. 428 * The value is rounded using the given method which is any method defined 429 * in {@link BigDecimal}. 430 * If {@code x} is infinite or {@code NaN}, then the value of {@code x} is 431 * returned unchanged, regardless of the other parameters. 432 * 433 * @param x Value to round. 434 * @param scale Number of digits to the right of the decimal point. 435 * @param roundingMethod Rounding method as defined in {@link BigDecimal}. 436 * @return the rounded value. 437 * @throws ArithmeticException if {@code roundingMethod == ROUND_UNNECESSARY} 438 * and the specified scaling operation would require rounding. 439 * @throws IllegalArgumentException if {@code roundingMethod} does not 440 * represent a valid rounding mode. 441 * @since 1.1 (previously in {@code MathUtils}, moved as of version 3.0) 442 */ 443 public static double round(double x, int scale, int roundingMethod) { 444 try { 445 final double rounded = (new BigDecimal(Double.toString(x)) 446 .setScale(scale, roundingMethod)) 447 .doubleValue(); 448 // MATH-1089: negative values rounded to zero should result in negative zero 449 return rounded == POSITIVE_ZERO ? POSITIVE_ZERO * x : rounded; 450 } catch (NumberFormatException ex) { 451 if (Double.isInfinite(x)) { 452 return x; 453 } else { 454 return Double.NaN; 455 } 456 } 457 } 458 459 /** 460 * Rounds the given value to the specified number of decimal places. 461 * The value is rounded using the {@link BigDecimal#ROUND_HALF_UP} method. 462 * 463 * @param x Value to round. 464 * @param scale Number of digits to the right of the decimal point. 465 * @return the rounded value. 466 * @since 1.1 (previously in {@code MathUtils}, moved as of version 3.0) 467 */ 468 public static float round(float x, int scale) { 469 return round(x, scale, BigDecimal.ROUND_HALF_UP); 470 } 471 472 /** 473 * Rounds the given value to the specified number of decimal places. 474 * The value is rounded using the given method which is any method defined 475 * in {@link BigDecimal}. 476 * 477 * @param x Value to round. 478 * @param scale Number of digits to the right of the decimal point. 479 * @param roundingMethod Rounding method as defined in {@link BigDecimal}. 480 * @return the rounded value. 481 * @since 1.1 (previously in {@code MathUtils}, moved as of version 3.0) 482 * @throws MathArithmeticException if an exact operation is required but result is not exact 483 * @throws MathIllegalArgumentException if {@code roundingMethod} is not a valid rounding method. 484 */ 485 public static float round(float x, int scale, int roundingMethod) 486 throws MathArithmeticException, MathIllegalArgumentException { 487 final float sign = FastMath.copySign(1f, x); 488 final float factor = (float) FastMath.pow(10.0f, scale) * sign; 489 return (float) roundUnscaled(x * factor, sign, roundingMethod) / factor; 490 } 491 492 /** 493 * Rounds the given non-negative value to the "nearest" integer. Nearest is 494 * determined by the rounding method specified. Rounding methods are defined 495 * in {@link BigDecimal}. 496 * 497 * @param unscaled Value to round. 498 * @param sign Sign of the original, scaled value. 499 * @param roundingMethod Rounding method, as defined in {@link BigDecimal}. 500 * @return the rounded value. 501 * @throws MathArithmeticException if an exact operation is required but result is not exact 502 * @throws MathIllegalArgumentException if {@code roundingMethod} is not a valid rounding method. 503 * @since 1.1 (previously in {@code MathUtils}, moved as of version 3.0) 504 */ 505 private static double roundUnscaled(double unscaled, 506 double sign, 507 int roundingMethod) 508 throws MathArithmeticException, MathIllegalArgumentException { 509 switch (roundingMethod) { 510 case BigDecimal.ROUND_CEILING : 511 if (sign == -1) { 512 unscaled = FastMath.floor(FastMath.nextAfter(unscaled, Double.NEGATIVE_INFINITY)); 513 } else { 514 unscaled = FastMath.ceil(FastMath.nextAfter(unscaled, Double.POSITIVE_INFINITY)); 515 } 516 break; 517 case BigDecimal.ROUND_DOWN : 518 unscaled = FastMath.floor(FastMath.nextAfter(unscaled, Double.NEGATIVE_INFINITY)); 519 break; 520 case BigDecimal.ROUND_FLOOR : 521 if (sign == -1) { 522 unscaled = FastMath.ceil(FastMath.nextAfter(unscaled, Double.POSITIVE_INFINITY)); 523 } else { 524 unscaled = FastMath.floor(FastMath.nextAfter(unscaled, Double.NEGATIVE_INFINITY)); 525 } 526 break; 527 case BigDecimal.ROUND_HALF_DOWN : { 528 unscaled = FastMath.nextAfter(unscaled, Double.NEGATIVE_INFINITY); 529 double fraction = unscaled - FastMath.floor(unscaled); 530 if (fraction > 0.5) { 531 unscaled = FastMath.ceil(unscaled); 532 } else { 533 unscaled = FastMath.floor(unscaled); 534 } 535 break; 536 } 537 case BigDecimal.ROUND_HALF_EVEN : { 538 double fraction = unscaled - FastMath.floor(unscaled); 539 if (fraction > 0.5) { 540 unscaled = FastMath.ceil(unscaled); 541 } else if (fraction < 0.5) { 542 unscaled = FastMath.floor(unscaled); 543 } else { 544 // The following equality test is intentional and needed for rounding purposes 545 if (FastMath.floor(unscaled) / 2.0 == FastMath.floor(FastMath.floor(unscaled) / 2.0)) { // even 546 unscaled = FastMath.floor(unscaled); 547 } else { // odd 548 unscaled = FastMath.ceil(unscaled); 549 } 550 } 551 break; 552 } 553 case BigDecimal.ROUND_HALF_UP : { 554 unscaled = FastMath.nextAfter(unscaled, Double.POSITIVE_INFINITY); 555 double fraction = unscaled - FastMath.floor(unscaled); 556 if (fraction >= 0.5) { 557 unscaled = FastMath.ceil(unscaled); 558 } else { 559 unscaled = FastMath.floor(unscaled); 560 } 561 break; 562 } 563 case BigDecimal.ROUND_UNNECESSARY : 564 if (unscaled != FastMath.floor(unscaled)) { 565 throw new MathArithmeticException(); 566 } 567 break; 568 case BigDecimal.ROUND_UP : 569 // do not round if the discarded fraction is equal to zero 570 if (unscaled != FastMath.floor(unscaled)) { 571 unscaled = FastMath.ceil(FastMath.nextAfter(unscaled, Double.POSITIVE_INFINITY)); 572 } 573 break; 574 default : 575 throw new MathIllegalArgumentException(LocalizedFormats.INVALID_ROUNDING_METHOD, 576 roundingMethod, 577 "ROUND_CEILING", BigDecimal.ROUND_CEILING, 578 "ROUND_DOWN", BigDecimal.ROUND_DOWN, 579 "ROUND_FLOOR", BigDecimal.ROUND_FLOOR, 580 "ROUND_HALF_DOWN", BigDecimal.ROUND_HALF_DOWN, 581 "ROUND_HALF_EVEN", BigDecimal.ROUND_HALF_EVEN, 582 "ROUND_HALF_UP", BigDecimal.ROUND_HALF_UP, 583 "ROUND_UNNECESSARY", BigDecimal.ROUND_UNNECESSARY, 584 "ROUND_UP", BigDecimal.ROUND_UP); 585 } 586 return unscaled; 587 } 588 589 590 /** 591 * Computes a number {@code delta} close to {@code originalDelta} with 592 * the property that <pre><code> 593 * x + delta - x 594 * </code></pre> 595 * is exactly machine-representable. 596 * This is useful when computing numerical derivatives, in order to reduce 597 * roundoff errors. 598 * 599 * @param x Value. 600 * @param originalDelta Offset value. 601 * @return a number {@code delta} so that {@code x + delta} and {@code x} 602 * differ by a representable floating number. 603 */ 604 public static double representableDelta(double x, 605 double originalDelta) { 606 return x + originalDelta - x; 607 } 608}