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