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