1 /* 2 * Licensed to the Apache Software Foundation (ASF) under one or more 3 * contributor license agreements. See the NOTICE file distributed with 4 * this work for additional information regarding copyright ownership. 5 * The ASF licenses this file to You under the Apache License, Version 2.0 6 * (the "License"); you may not use this file except in compliance with 7 * the License. You may obtain a copy of the License at 8 * 9 * http://www.apache.org/licenses/LICENSE-2.0 10 * 11 * Unless required by applicable law or agreed to in writing, software 12 * distributed under the License is distributed on an "AS IS" BASIS, 13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 14 * See the License for the specific language governing permissions and 15 * limitations under the License. 16 */ 17 18 package org.apache.commons.math.util; 19 20 import java.math.BigDecimal; 21 import java.math.BigInteger; 22 import java.util.Arrays; 23 24 import org.apache.commons.math.MathRuntimeException; 25 import org.apache.commons.math.exception.util.Localizable; 26 import org.apache.commons.math.exception.util.LocalizedFormats; 27 import org.apache.commons.math.exception.NonMonotonousSequenceException; 28 29 /** 30 * Some useful additions to the built-in functions in {@link Math}. 31 * @version $Revision: 1073472 $ $Date: 2011-02-22 20:49:07 +0100 (mar. 22 fvr. 2011) $ 32 */ 33 public final class MathUtils { 34 35 /** Smallest positive number such that 1 - EPSILON is not numerically equal to 1. */ 36 public static final double EPSILON = 0x1.0p-53; 37 38 /** Safe minimum, such that 1 / SAFE_MIN does not overflow. 39 * <p>In IEEE 754 arithmetic, this is also the smallest normalized 40 * number 2<sup>-1022</sup>.</p> 41 */ 42 public static final double SAFE_MIN = 0x1.0p-1022; 43 44 /** 45 * 2 π. 46 * @since 2.1 47 */ 48 public static final double TWO_PI = 2 * FastMath.PI; 49 50 /** -1.0 cast as a byte. */ 51 private static final byte NB = (byte)-1; 52 53 /** -1.0 cast as a short. */ 54 private static final short NS = (short)-1; 55 56 /** 1.0 cast as a byte. */ 57 private static final byte PB = (byte)1; 58 59 /** 1.0 cast as a short. */ 60 private static final short PS = (short)1; 61 62 /** 0.0 cast as a byte. */ 63 private static final byte ZB = (byte)0; 64 65 /** 0.0 cast as a short. */ 66 private static final short ZS = (short)0; 67 68 /** Gap between NaN and regular numbers. */ 69 private static final int NAN_GAP = 4 * 1024 * 1024; 70 71 /** Offset to order signed double numbers lexicographically. */ 72 private static final long SGN_MASK = 0x8000000000000000L; 73 74 /** Offset to order signed double numbers lexicographically. */ 75 private static final int SGN_MASK_FLOAT = 0x80000000; 76 77 /** All long-representable factorials */ 78 private static final long[] FACTORIALS = new long[] { 79 1l, 1l, 2l, 80 6l, 24l, 120l, 81 720l, 5040l, 40320l, 82 362880l, 3628800l, 39916800l, 83 479001600l, 6227020800l, 87178291200l, 84 1307674368000l, 20922789888000l, 355687428096000l, 85 6402373705728000l, 121645100408832000l, 2432902008176640000l }; 86 87 /** 88 * Private Constructor 89 */ 90 private MathUtils() { 91 super(); 92 } 93 94 /** 95 * Add two integers, checking for overflow. 96 * 97 * @param x an addend 98 * @param y an addend 99 * @return the sum <code>x+y</code> 100 * @throws ArithmeticException if the result can not be represented as an 101 * int 102 * @since 1.1 103 */ 104 public static int addAndCheck(int x, int y) { 105 long s = (long)x + (long)y; 106 if (s < Integer.MIN_VALUE || s > Integer.MAX_VALUE) { 107 throw MathRuntimeException.createArithmeticException(LocalizedFormats.OVERFLOW_IN_ADDITION, x, y); 108 } 109 return (int)s; 110 } 111 112 /** 113 * Add two long integers, checking for overflow. 114 * 115 * @param a an addend 116 * @param b an addend 117 * @return the sum <code>a+b</code> 118 * @throws ArithmeticException if the result can not be represented as an 119 * long 120 * @since 1.2 121 */ 122 public static long addAndCheck(long a, long b) { 123 return addAndCheck(a, b, LocalizedFormats.OVERFLOW_IN_ADDITION); 124 } 125 126 /** 127 * Add two long integers, checking for overflow. 128 * 129 * @param a an addend 130 * @param b an addend 131 * @param pattern the pattern to use for any thrown exception. 132 * @return the sum <code>a+b</code> 133 * @throws ArithmeticException if the result can not be represented as an 134 * long 135 * @since 1.2 136 */ 137 private static long addAndCheck(long a, long b, Localizable pattern) { 138 long ret; 139 if (a > b) { 140 // use symmetry to reduce boundary cases 141 ret = addAndCheck(b, a, pattern); 142 } else { 143 // assert a <= b 144 145 if (a < 0) { 146 if (b < 0) { 147 // check for negative overflow 148 if (Long.MIN_VALUE - b <= a) { 149 ret = a + b; 150 } else { 151 throw MathRuntimeException.createArithmeticException(pattern, a, b); 152 } 153 } else { 154 // opposite sign addition is always safe 155 ret = a + b; 156 } 157 } else { 158 // assert a >= 0 159 // assert b >= 0 160 161 // check for positive overflow 162 if (a <= Long.MAX_VALUE - b) { 163 ret = a + b; 164 } else { 165 throw MathRuntimeException.createArithmeticException(pattern, a, b); 166 } 167 } 168 } 169 return ret; 170 } 171 172 /** 173 * Returns an exact representation of the <a 174 * href="http://mathworld.wolfram.com/BinomialCoefficient.html"> Binomial 175 * Coefficient</a>, "<code>n choose k</code>", the number of 176 * <code>k</code>-element subsets that can be selected from an 177 * <code>n</code>-element set. 178 * <p> 179 * <Strong>Preconditions</strong>: 180 * <ul> 181 * <li> <code>0 <= k <= n </code> (otherwise 182 * <code>IllegalArgumentException</code> is thrown)</li> 183 * <li> The result is small enough to fit into a <code>long</code>. The 184 * largest value of <code>n</code> for which all coefficients are 185 * <code> < Long.MAX_VALUE</code> is 66. If the computed value exceeds 186 * <code>Long.MAX_VALUE</code> an <code>ArithMeticException</code> is 187 * thrown.</li> 188 * </ul></p> 189 * 190 * @param n the size of the set 191 * @param k the size of the subsets to be counted 192 * @return <code>n choose k</code> 193 * @throws IllegalArgumentException if preconditions are not met. 194 * @throws ArithmeticException if the result is too large to be represented 195 * by a long integer. 196 */ 197 public static long binomialCoefficient(final int n, final int k) { 198 checkBinomial(n, k); 199 if ((n == k) || (k == 0)) { 200 return 1; 201 } 202 if ((k == 1) || (k == n - 1)) { 203 return n; 204 } 205 // Use symmetry for large k 206 if (k > n / 2) 207 return binomialCoefficient(n, n - k); 208 209 // We use the formula 210 // (n choose k) = n! / (n-k)! / k! 211 // (n choose k) == ((n-k+1)*...*n) / (1*...*k) 212 // which could be written 213 // (n choose k) == (n-1 choose k-1) * n / k 214 long result = 1; 215 if (n <= 61) { 216 // For n <= 61, the naive implementation cannot overflow. 217 int i = n - k + 1; 218 for (int j = 1; j <= k; j++) { 219 result = result * i / j; 220 i++; 221 } 222 } else if (n <= 66) { 223 // For n > 61 but n <= 66, the result cannot overflow, 224 // but we must take care not to overflow intermediate values. 225 int i = n - k + 1; 226 for (int j = 1; j <= k; j++) { 227 // We know that (result * i) is divisible by j, 228 // but (result * i) may overflow, so we split j: 229 // Filter out the gcd, d, so j/d and i/d are integer. 230 // result is divisible by (j/d) because (j/d) 231 // is relative prime to (i/d) and is a divisor of 232 // result * (i/d). 233 final long d = gcd(i, j); 234 result = (result / (j / d)) * (i / d); 235 i++; 236 } 237 } else { 238 // For n > 66, a result overflow might occur, so we check 239 // the multiplication, taking care to not overflow 240 // unnecessary. 241 int i = n - k + 1; 242 for (int j = 1; j <= k; j++) { 243 final long d = gcd(i, j); 244 result = mulAndCheck(result / (j / d), i / d); 245 i++; 246 } 247 } 248 return result; 249 } 250 251 /** 252 * Returns a <code>double</code> representation of the <a 253 * href="http://mathworld.wolfram.com/BinomialCoefficient.html"> Binomial 254 * Coefficient</a>, "<code>n choose k</code>", the number of 255 * <code>k</code>-element subsets that can be selected from an 256 * <code>n</code>-element set. 257 * <p> 258 * <Strong>Preconditions</strong>: 259 * <ul> 260 * <li> <code>0 <= k <= n </code> (otherwise 261 * <code>IllegalArgumentException</code> is thrown)</li> 262 * <li> The result is small enough to fit into a <code>double</code>. The 263 * largest value of <code>n</code> for which all coefficients are < 264 * Double.MAX_VALUE is 1029. If the computed value exceeds Double.MAX_VALUE, 265 * Double.POSITIVE_INFINITY is returned</li> 266 * </ul></p> 267 * 268 * @param n the size of the set 269 * @param k the size of the subsets to be counted 270 * @return <code>n choose k</code> 271 * @throws IllegalArgumentException if preconditions are not met. 272 */ 273 public static double binomialCoefficientDouble(final int n, final int k) { 274 checkBinomial(n, k); 275 if ((n == k) || (k == 0)) { 276 return 1d; 277 } 278 if ((k == 1) || (k == n - 1)) { 279 return n; 280 } 281 if (k > n/2) { 282 return binomialCoefficientDouble(n, n - k); 283 } 284 if (n < 67) { 285 return binomialCoefficient(n,k); 286 } 287 288 double result = 1d; 289 for (int i = 1; i <= k; i++) { 290 result *= (double)(n - k + i) / (double)i; 291 } 292 293 return FastMath.floor(result + 0.5); 294 } 295 296 /** 297 * Returns the natural <code>log</code> of the <a 298 * href="http://mathworld.wolfram.com/BinomialCoefficient.html"> Binomial 299 * Coefficient</a>, "<code>n choose k</code>", the number of 300 * <code>k</code>-element subsets that can be selected from an 301 * <code>n</code>-element set. 302 * <p> 303 * <Strong>Preconditions</strong>: 304 * <ul> 305 * <li> <code>0 <= k <= n </code> (otherwise 306 * <code>IllegalArgumentException</code> is thrown)</li> 307 * </ul></p> 308 * 309 * @param n the size of the set 310 * @param k the size of the subsets to be counted 311 * @return <code>n choose k</code> 312 * @throws IllegalArgumentException if preconditions are not met. 313 */ 314 public static double binomialCoefficientLog(final int n, final int k) { 315 checkBinomial(n, k); 316 if ((n == k) || (k == 0)) { 317 return 0; 318 } 319 if ((k == 1) || (k == n - 1)) { 320 return FastMath.log(n); 321 } 322 323 /* 324 * For values small enough to do exact integer computation, 325 * return the log of the exact value 326 */ 327 if (n < 67) { 328 return FastMath.log(binomialCoefficient(n,k)); 329 } 330 331 /* 332 * Return the log of binomialCoefficientDouble for values that will not 333 * overflow binomialCoefficientDouble 334 */ 335 if (n < 1030) { 336 return FastMath.log(binomialCoefficientDouble(n, k)); 337 } 338 339 if (k > n / 2) { 340 return binomialCoefficientLog(n, n - k); 341 } 342 343 /* 344 * Sum logs for values that could overflow 345 */ 346 double logSum = 0; 347 348 // n!/(n-k)! 349 for (int i = n - k + 1; i <= n; i++) { 350 logSum += FastMath.log(i); 351 } 352 353 // divide by k! 354 for (int i = 2; i <= k; i++) { 355 logSum -= FastMath.log(i); 356 } 357 358 return logSum; 359 } 360 361 /** 362 * Check binomial preconditions. 363 * @param n the size of the set 364 * @param k the size of the subsets to be counted 365 * @exception IllegalArgumentException if preconditions are not met. 366 */ 367 private static void checkBinomial(final int n, final int k) 368 throws IllegalArgumentException { 369 if (n < k) { 370 throw MathRuntimeException.createIllegalArgumentException( 371 LocalizedFormats.BINOMIAL_INVALID_PARAMETERS_ORDER, 372 n, k); 373 } 374 if (n < 0) { 375 throw MathRuntimeException.createIllegalArgumentException( 376 LocalizedFormats.BINOMIAL_NEGATIVE_PARAMETER, 377 n); 378 } 379 } 380 381 /** 382 * Compares two numbers given some amount of allowed error. 383 * 384 * @param x the first number 385 * @param y the second number 386 * @param eps the amount of error to allow when checking for equality 387 * @return <ul><li>0 if {@link #equals(double, double, double) equals(x, y, eps)}</li> 388 * <li>< 0 if !{@link #equals(double, double, double) equals(x, y, eps)} && x < y</li> 389 * <li>> 0 if !{@link #equals(double, double, double) equals(x, y, eps)} && x > y</li></ul> 390 */ 391 public static int compareTo(double x, double y, double eps) { 392 if (equals(x, y, eps)) { 393 return 0; 394 } else if (x < y) { 395 return -1; 396 } 397 return 1; 398 } 399 400 /** 401 * Returns the <a href="http://mathworld.wolfram.com/HyperbolicCosine.html"> 402 * hyperbolic cosine</a> of x. 403 * 404 * @param x double value for which to find the hyperbolic cosine 405 * @return hyperbolic cosine of x 406 */ 407 public static double cosh(double x) { 408 return (FastMath.exp(x) + FastMath.exp(-x)) / 2.0; 409 } 410 411 /** 412 * Returns true iff they are strictly equal. 413 * 414 * @param x first value 415 * @param y second value 416 * @return {@code true} if the values are equal. 417 * @deprecated as of 2.2 his method considers that {@code NaN == NaN}. In release 418 * 3.0, the semantics will change in order to comply with IEEE754 where it 419 * is specified that {@code NaN != NaN}. 420 * New methods have been added for those cases wher the old semantics is 421 * useful (see e.g. {@link #equalsIncludingNaN(float,float) 422 * equalsIncludingNaN}. 423 */ 424 @Deprecated 425 public static boolean equals(float x, float y) { 426 return (Float.isNaN(x) && Float.isNaN(y)) || x == y; 427 } 428 429 /** 430 * Returns true if both arguments are NaN or neither is NaN and they are 431 * equal as defined by {@link #equals(float,float,int) equals(x, y, 1)}. 432 * 433 * @param x first value 434 * @param y second value 435 * @return {@code true} if the values are equal or both are NaN. 436 * @since 2.2 437 */ 438 public static boolean equalsIncludingNaN(float x, float y) { 439 return (Float.isNaN(x) && Float.isNaN(y)) || equals(x, y, 1); 440 } 441 442 /** 443 * Returns true if both arguments are equal or within the range of allowed 444 * error (inclusive). 445 * 446 * @param x first value 447 * @param y second value 448 * @param eps the amount of absolute error to allow. 449 * @return {@code true} if the values are equal or within range of each other. 450 * @since 2.2 451 */ 452 public static boolean equals(float x, float y, float eps) { 453 return equals(x, y, 1) || FastMath.abs(y - x) <= eps; 454 } 455 456 /** 457 * Returns true if both arguments are NaN or are equal or within the range 458 * of allowed error (inclusive). 459 * 460 * @param x first value 461 * @param y second value 462 * @param eps the amount of absolute error to allow. 463 * @return {@code true} if the values are equal or within range of each other, 464 * or both are NaN. 465 * @since 2.2 466 */ 467 public static boolean equalsIncludingNaN(float x, float y, float eps) { 468 return equalsIncludingNaN(x, y) || (FastMath.abs(y - x) <= eps); 469 } 470 471 /** 472 * Returns true if both arguments are equal or within the range of allowed 473 * error (inclusive). 474 * Two float numbers are considered equal if there are {@code (maxUlps - 1)} 475 * (or fewer) floating point numbers between them, i.e. two adjacent floating 476 * point numbers are considered equal. 477 * Adapted from <a 478 * href="http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm"> 479 * Bruce Dawson</a> 480 * 481 * @param x first value 482 * @param y second value 483 * @param maxUlps {@code (maxUlps - 1)} is the number of floating point 484 * values between {@code x} and {@code y}. 485 * @return {@code true} if there are fewer than {@code maxUlps} floating 486 * point values between {@code x} and {@code y}. 487 * @since 2.2 488 */ 489 public static boolean equals(float x, float y, int maxUlps) { 490 // Check that "maxUlps" is non-negative and small enough so that 491 // NaN won't compare as equal to anything (except another NaN). 492 assert maxUlps > 0 && maxUlps < NAN_GAP; 493 494 int xInt = Float.floatToIntBits(x); 495 int yInt = Float.floatToIntBits(y); 496 497 // Make lexicographically ordered as a two's-complement integer. 498 if (xInt < 0) { 499 xInt = SGN_MASK_FLOAT - xInt; 500 } 501 if (yInt < 0) { 502 yInt = SGN_MASK_FLOAT - yInt; 503 } 504 505 final boolean isEqual = FastMath.abs(xInt - yInt) <= maxUlps; 506 507 return isEqual && !Float.isNaN(x) && !Float.isNaN(y); 508 } 509 510 /** 511 * Returns true if both arguments are NaN or if they are equal as defined 512 * by {@link #equals(float,float,int) equals(x, y, maxUlps)}. 513 * 514 * @param x first value 515 * @param y second value 516 * @param maxUlps {@code (maxUlps - 1)} is the number of floating point 517 * values between {@code x} and {@code y}. 518 * @return {@code true} if both arguments are NaN or if there are less than 519 * {@code maxUlps} floating point values between {@code x} and {@code y}. 520 * @since 2.2 521 */ 522 public static boolean equalsIncludingNaN(float x, float y, int maxUlps) { 523 return (Float.isNaN(x) && Float.isNaN(y)) || equals(x, y, maxUlps); 524 } 525 526 /** 527 * Returns true iff both arguments are null or have same dimensions and all 528 * their elements are equal as defined by 529 * {@link #equals(float,float)}. 530 * 531 * @param x first array 532 * @param y second array 533 * @return true if the values are both null or have same dimension 534 * and equal elements. 535 * @deprecated as of 2.2 this method considers that {@code NaN == NaN}. In release 536 * 3.0, the semantics will change in order to comply with IEEE754 where it 537 * is specified that {@code NaN != NaN}. 538 * New methods have been added for those cases where the old semantics is 539 * useful (see e.g. {@link #equalsIncludingNaN(float[],float[]) 540 * equalsIncludingNaN}. 541 */ 542 @Deprecated 543 public static boolean equals(float[] x, float[] y) { 544 if ((x == null) || (y == null)) { 545 return !((x == null) ^ (y == null)); 546 } 547 if (x.length != y.length) { 548 return false; 549 } 550 for (int i = 0; i < x.length; ++i) { 551 if (!equals(x[i], y[i])) { 552 return false; 553 } 554 } 555 return true; 556 } 557 558 /** 559 * Returns true iff both arguments are null or have same dimensions and all 560 * their elements are equal as defined by 561 * {@link #equalsIncludingNaN(float,float)}. 562 * 563 * @param x first array 564 * @param y second array 565 * @return true if the values are both null or have same dimension and 566 * equal elements 567 * @since 2.2 568 */ 569 public static boolean equalsIncludingNaN(float[] x, float[] y) { 570 if ((x == null) || (y == null)) { 571 return !((x == null) ^ (y == null)); 572 } 573 if (x.length != y.length) { 574 return false; 575 } 576 for (int i = 0; i < x.length; ++i) { 577 if (!equalsIncludingNaN(x[i], y[i])) { 578 return false; 579 } 580 } 581 return true; 582 } 583 584 /** 585 * Returns true iff both arguments are NaN or neither is NaN and they are 586 * equal 587 * 588 * <p>This method considers that {@code NaN == NaN}. In release 589 * 3.0, the semantics will change in order to comply with IEEE754 where it 590 * is specified that {@code NaN != NaN}. 591 * New methods have been added for those cases where the old semantics 592 * (w.r.t. NaN) is useful (see e.g. 593 * {@link #equalsIncludingNaN(double,double, double) equalsIncludingNaN}. 594 * </p> 595 * 596 * @param x first value 597 * @param y second value 598 * @return {@code true} if the values are equal. 599 */ 600 public static boolean equals(double x, double y) { 601 return (Double.isNaN(x) && Double.isNaN(y)) || x == y; 602 } 603 604 /** 605 * Returns true if both arguments are NaN or neither is NaN and they are 606 * equal as defined by {@link #equals(double,double,int) equals(x, y, 1)}. 607 * 608 * @param x first value 609 * @param y second value 610 * @return {@code true} if the values are equal or both are NaN. 611 * @since 2.2 612 */ 613 public static boolean equalsIncludingNaN(double x, double y) { 614 return (Double.isNaN(x) && Double.isNaN(y)) || equals(x, y, 1); 615 } 616 617 /** 618 * Returns true if both arguments are equal or within the range of allowed 619 * error (inclusive). 620 * <p> 621 * Two NaNs are considered equals, as are two infinities with same sign. 622 * </p> 623 * <p>This method considers that {@code NaN == NaN}. In release 624 * 3.0, the semantics will change in order to comply with IEEE754 where it 625 * is specified that {@code NaN != NaN}. 626 * New methods have been added for those cases where the old semantics 627 * (w.r.t. NaN) is useful (see e.g. 628 * {@link #equalsIncludingNaN(double,double, double) equalsIncludingNaN}. 629 * </p> 630 * @param x first value 631 * @param y second value 632 * @param eps the amount of absolute error to allow. 633 * @return {@code true} if the values are equal or within range of each other. 634 */ 635 public static boolean equals(double x, double y, double eps) { 636 return equals(x, y) || FastMath.abs(y - x) <= eps; 637 } 638 639 /** 640 * Returns true if both arguments are NaN or are equal or within the range 641 * of allowed error (inclusive). 642 * 643 * @param x first value 644 * @param y second value 645 * @param eps the amount of absolute error to allow. 646 * @return {@code true} if the values are equal or within range of each other, 647 * or both are NaN. 648 * @since 2.2 649 */ 650 public static boolean equalsIncludingNaN(double x, double y, double eps) { 651 return equalsIncludingNaN(x, y) || (FastMath.abs(y - x) <= eps); 652 } 653 654 /** 655 * Returns true if both arguments are equal or within the range of allowed 656 * error (inclusive). 657 * Two float numbers are considered equal if there are {@code (maxUlps - 1)} 658 * (or fewer) floating point numbers between them, i.e. two adjacent floating 659 * point numbers are considered equal. 660 * Adapted from <a 661 * href="http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm"> 662 * Bruce Dawson</a> 663 * 664 * <p>This method considers that {@code NaN == NaN}. In release 665 * 3.0, the semantics will change in order to comply with IEEE754 where it 666 * is specified that {@code NaN != NaN}. 667 * New methods have been added for those cases where the old semantics 668 * (w.r.t. NaN) is useful (see e.g. 669 * {@link #equalsIncludingNaN(double,double, int) equalsIncludingNaN}. 670 * </p> 671 * 672 * @param x first value 673 * @param y second value 674 * @param maxUlps {@code (maxUlps - 1)} is the number of floating point 675 * values between {@code x} and {@code y}. 676 * @return {@code true} if there are fewer than {@code maxUlps} floating 677 * point values between {@code x} and {@code y}. 678 */ 679 public static boolean equals(double x, double y, int maxUlps) { 680 // Check that "maxUlps" is non-negative and small enough so that 681 // NaN won't compare as equal to anything (except another NaN). 682 assert maxUlps > 0 && maxUlps < NAN_GAP; 683 684 long xInt = Double.doubleToLongBits(x); 685 long yInt = Double.doubleToLongBits(y); 686 687 // Make lexicographically ordered as a two's-complement integer. 688 if (xInt < 0) { 689 xInt = SGN_MASK - xInt; 690 } 691 if (yInt < 0) { 692 yInt = SGN_MASK - yInt; 693 } 694 695 return FastMath.abs(xInt - yInt) <= maxUlps; 696 } 697 698 /** 699 * Returns true if both arguments are NaN or if they are equal as defined 700 * by {@link #equals(double,double,int) equals(x, y, maxUlps}. 701 * 702 * @param x first value 703 * @param y second value 704 * @param maxUlps {@code (maxUlps - 1)} is the number of floating point 705 * values between {@code x} and {@code y}. 706 * @return {@code true} if both arguments are NaN or if there are less than 707 * {@code maxUlps} floating point values between {@code x} and {@code y}. 708 * @since 2.2 709 */ 710 public static boolean equalsIncludingNaN(double x, double y, int maxUlps) { 711 return (Double.isNaN(x) && Double.isNaN(y)) || equals(x, y, maxUlps); 712 } 713 714 /** 715 * Returns true iff both arguments are null or have same dimensions and all 716 * their elements are equal as defined by 717 * {@link #equals(double,double)}. 718 * 719 * <p>This method considers that {@code NaN == NaN}. In release 720 * 3.0, the semantics will change in order to comply with IEEE754 where it 721 * is specified that {@code NaN != NaN}. 722 * New methods have been added for those cases wher the old semantics is 723 * useful (see e.g. {@link #equalsIncludingNaN(double[],double[]) 724 * equalsIncludingNaN}. 725 * </p> 726 * @param x first array 727 * @param y second array 728 * @return true if the values are both null or have same dimension 729 * and equal elements. 730 */ 731 public static boolean equals(double[] x, double[] y) { 732 if ((x == null) || (y == null)) { 733 return !((x == null) ^ (y == null)); 734 } 735 if (x.length != y.length) { 736 return false; 737 } 738 for (int i = 0; i < x.length; ++i) { 739 if (!equals(x[i], y[i])) { 740 return false; 741 } 742 } 743 return true; 744 } 745 746 /** 747 * Returns true iff both arguments are null or have same dimensions and all 748 * their elements are equal as defined by 749 * {@link #equalsIncludingNaN(double,double)}. 750 * 751 * @param x first array 752 * @param y second array 753 * @return true if the values are both null or have same dimension and 754 * equal elements 755 * @since 2.2 756 */ 757 public static boolean equalsIncludingNaN(double[] x, double[] y) { 758 if ((x == null) || (y == null)) { 759 return !((x == null) ^ (y == null)); 760 } 761 if (x.length != y.length) { 762 return false; 763 } 764 for (int i = 0; i < x.length; ++i) { 765 if (!equalsIncludingNaN(x[i], y[i])) { 766 return false; 767 } 768 } 769 return true; 770 } 771 772 /** 773 * Returns n!. Shorthand for <code>n</code> <a 774 * href="http://mathworld.wolfram.com/Factorial.html"> Factorial</a>, the 775 * product of the numbers <code>1,...,n</code>. 776 * <p> 777 * <Strong>Preconditions</strong>: 778 * <ul> 779 * <li> <code>n >= 0</code> (otherwise 780 * <code>IllegalArgumentException</code> is thrown)</li> 781 * <li> The result is small enough to fit into a <code>long</code>. The 782 * largest value of <code>n</code> for which <code>n!</code> < 783 * Long.MAX_VALUE</code> is 20. If the computed value exceeds <code>Long.MAX_VALUE</code> 784 * an <code>ArithMeticException </code> is thrown.</li> 785 * </ul> 786 * </p> 787 * 788 * @param n argument 789 * @return <code>n!</code> 790 * @throws ArithmeticException if the result is too large to be represented 791 * by a long integer. 792 * @throws IllegalArgumentException if n < 0 793 */ 794 public static long factorial(final int n) { 795 if (n < 0) { 796 throw MathRuntimeException.createIllegalArgumentException( 797 LocalizedFormats.FACTORIAL_NEGATIVE_PARAMETER, 798 n); 799 } 800 if (n > 20) { 801 throw new ArithmeticException( 802 "factorial value is too large to fit in a long"); 803 } 804 return FACTORIALS[n]; 805 } 806 807 /** 808 * Returns n!. Shorthand for <code>n</code> <a 809 * href="http://mathworld.wolfram.com/Factorial.html"> Factorial</a>, the 810 * product of the numbers <code>1,...,n</code> as a <code>double</code>. 811 * <p> 812 * <Strong>Preconditions</strong>: 813 * <ul> 814 * <li> <code>n >= 0</code> (otherwise 815 * <code>IllegalArgumentException</code> is thrown)</li> 816 * <li> The result is small enough to fit into a <code>double</code>. The 817 * largest value of <code>n</code> for which <code>n!</code> < 818 * Double.MAX_VALUE</code> is 170. If the computed value exceeds 819 * Double.MAX_VALUE, Double.POSITIVE_INFINITY is returned</li> 820 * </ul> 821 * </p> 822 * 823 * @param n argument 824 * @return <code>n!</code> 825 * @throws IllegalArgumentException if n < 0 826 */ 827 public static double factorialDouble(final int n) { 828 if (n < 0) { 829 throw MathRuntimeException.createIllegalArgumentException( 830 LocalizedFormats.FACTORIAL_NEGATIVE_PARAMETER, 831 n); 832 } 833 if (n < 21) { 834 return factorial(n); 835 } 836 return FastMath.floor(FastMath.exp(factorialLog(n)) + 0.5); 837 } 838 839 /** 840 * Returns the natural logarithm of n!. 841 * <p> 842 * <Strong>Preconditions</strong>: 843 * <ul> 844 * <li> <code>n >= 0</code> (otherwise 845 * <code>IllegalArgumentException</code> is thrown)</li> 846 * </ul></p> 847 * 848 * @param n argument 849 * @return <code>n!</code> 850 * @throws IllegalArgumentException if preconditions are not met. 851 */ 852 public static double factorialLog(final int n) { 853 if (n < 0) { 854 throw MathRuntimeException.createIllegalArgumentException( 855 LocalizedFormats.FACTORIAL_NEGATIVE_PARAMETER, 856 n); 857 } 858 if (n < 21) { 859 return FastMath.log(factorial(n)); 860 } 861 double logSum = 0; 862 for (int i = 2; i <= n; i++) { 863 logSum += FastMath.log(i); 864 } 865 return logSum; 866 } 867 868 /** 869 * <p> 870 * Gets the greatest common divisor of the absolute value of two numbers, 871 * using the "binary gcd" method which avoids division and modulo 872 * operations. See Knuth 4.5.2 algorithm B. This algorithm is due to Josef 873 * Stein (1961). 874 * </p> 875 * Special cases: 876 * <ul> 877 * <li>The invocations 878 * <code>gcd(Integer.MIN_VALUE, Integer.MIN_VALUE)</code>, 879 * <code>gcd(Integer.MIN_VALUE, 0)</code> and 880 * <code>gcd(0, Integer.MIN_VALUE)</code> throw an 881 * <code>ArithmeticException</code>, because the result would be 2^31, which 882 * is too large for an int value.</li> 883 * <li>The result of <code>gcd(x, x)</code>, <code>gcd(0, x)</code> and 884 * <code>gcd(x, 0)</code> is the absolute value of <code>x</code>, except 885 * for the special cases above. 886 * <li>The invocation <code>gcd(0, 0)</code> is the only one which returns 887 * <code>0</code>.</li> 888 * </ul> 889 * 890 * @param p any number 891 * @param q any number 892 * @return the greatest common divisor, never negative 893 * @throws ArithmeticException if the result cannot be represented as a 894 * nonnegative int value 895 * @since 1.1 896 */ 897 public static int gcd(final int p, final int q) { 898 int u = p; 899 int v = q; 900 if ((u == 0) || (v == 0)) { 901 if ((u == Integer.MIN_VALUE) || (v == Integer.MIN_VALUE)) { 902 throw MathRuntimeException.createArithmeticException( 903 LocalizedFormats.GCD_OVERFLOW_32_BITS, 904 p, q); 905 } 906 return FastMath.abs(u) + FastMath.abs(v); 907 } 908 // keep u and v negative, as negative integers range down to 909 // -2^31, while positive numbers can only be as large as 2^31-1 910 // (i.e. we can't necessarily negate a negative number without 911 // overflow) 912 /* assert u!=0 && v!=0; */ 913 if (u > 0) { 914 u = -u; 915 } // make u negative 916 if (v > 0) { 917 v = -v; 918 } // make v negative 919 // B1. [Find power of 2] 920 int k = 0; 921 while ((u & 1) == 0 && (v & 1) == 0 && k < 31) { // while u and v are 922 // both even... 923 u /= 2; 924 v /= 2; 925 k++; // cast out twos. 926 } 927 if (k == 31) { 928 throw MathRuntimeException.createArithmeticException( 929 LocalizedFormats.GCD_OVERFLOW_32_BITS, 930 p, q); 931 } 932 // B2. Initialize: u and v have been divided by 2^k and at least 933 // one is odd. 934 int t = ((u & 1) == 1) ? v : -(u / 2)/* B3 */; 935 // t negative: u was odd, v may be even (t replaces v) 936 // t positive: u was even, v is odd (t replaces u) 937 do { 938 /* assert u<0 && v<0; */ 939 // B4/B3: cast out twos from t. 940 while ((t & 1) == 0) { // while t is even.. 941 t /= 2; // cast out twos 942 } 943 // B5 [reset max(u,v)] 944 if (t > 0) { 945 u = -t; 946 } else { 947 v = t; 948 } 949 // B6/B3. at this point both u and v should be odd. 950 t = (v - u) / 2; 951 // |u| larger: t positive (replace u) 952 // |v| larger: t negative (replace v) 953 } while (t != 0); 954 return -u * (1 << k); // gcd is u*2^k 955 } 956 957 /** 958 * <p> 959 * Gets the greatest common divisor of the absolute value of two numbers, 960 * using the "binary gcd" method which avoids division and modulo 961 * operations. See Knuth 4.5.2 algorithm B. This algorithm is due to Josef 962 * Stein (1961). 963 * </p> 964 * Special cases: 965 * <ul> 966 * <li>The invocations 967 * <code>gcd(Long.MIN_VALUE, Long.MIN_VALUE)</code>, 968 * <code>gcd(Long.MIN_VALUE, 0L)</code> and 969 * <code>gcd(0L, Long.MIN_VALUE)</code> throw an 970 * <code>ArithmeticException</code>, because the result would be 2^63, which 971 * is too large for a long value.</li> 972 * <li>The result of <code>gcd(x, x)</code>, <code>gcd(0L, x)</code> and 973 * <code>gcd(x, 0L)</code> is the absolute value of <code>x</code>, except 974 * for the special cases above. 975 * <li>The invocation <code>gcd(0L, 0L)</code> is the only one which returns 976 * <code>0L</code>.</li> 977 * </ul> 978 * 979 * @param p any number 980 * @param q any number 981 * @return the greatest common divisor, never negative 982 * @throws ArithmeticException if the result cannot be represented as a nonnegative long 983 * value 984 * @since 2.1 985 */ 986 public static long gcd(final long p, final long q) { 987 long u = p; 988 long v = q; 989 if ((u == 0) || (v == 0)) { 990 if ((u == Long.MIN_VALUE) || (v == Long.MIN_VALUE)){ 991 throw MathRuntimeException.createArithmeticException( 992 LocalizedFormats.GCD_OVERFLOW_64_BITS, 993 p, q); 994 } 995 return FastMath.abs(u) + FastMath.abs(v); 996 } 997 // keep u and v negative, as negative integers range down to 998 // -2^63, while positive numbers can only be as large as 2^63-1 999 // (i.e. we can't necessarily negate a negative number without 1000 // overflow) 1001 /* assert u!=0 && v!=0; */ 1002 if (u > 0) { 1003 u = -u; 1004 } // make u negative 1005 if (v > 0) { 1006 v = -v; 1007 } // make v negative 1008 // B1. [Find power of 2] 1009 int k = 0; 1010 while ((u & 1) == 0 && (v & 1) == 0 && k < 63) { // while u and v are 1011 // both even... 1012 u /= 2; 1013 v /= 2; 1014 k++; // cast out twos. 1015 } 1016 if (k == 63) { 1017 throw MathRuntimeException.createArithmeticException( 1018 LocalizedFormats.GCD_OVERFLOW_64_BITS, 1019 p, q); 1020 } 1021 // B2. Initialize: u and v have been divided by 2^k and at least 1022 // one is odd. 1023 long t = ((u & 1) == 1) ? v : -(u / 2)/* B3 */; 1024 // t negative: u was odd, v may be even (t replaces v) 1025 // t positive: u was even, v is odd (t replaces u) 1026 do { 1027 /* assert u<0 && v<0; */ 1028 // B4/B3: cast out twos from t. 1029 while ((t & 1) == 0) { // while t is even.. 1030 t /= 2; // cast out twos 1031 } 1032 // B5 [reset max(u,v)] 1033 if (t > 0) { 1034 u = -t; 1035 } else { 1036 v = t; 1037 } 1038 // B6/B3. at this point both u and v should be odd. 1039 t = (v - u) / 2; 1040 // |u| larger: t positive (replace u) 1041 // |v| larger: t negative (replace v) 1042 } while (t != 0); 1043 return -u * (1L << k); // gcd is u*2^k 1044 } 1045 1046 /** 1047 * Returns an integer hash code representing the given double value. 1048 * 1049 * @param value the value to be hashed 1050 * @return the hash code 1051 */ 1052 public static int hash(double value) { 1053 return new Double(value).hashCode(); 1054 } 1055 1056 /** 1057 * Returns an integer hash code representing the given double array. 1058 * 1059 * @param value the value to be hashed (may be null) 1060 * @return the hash code 1061 * @since 1.2 1062 */ 1063 public static int hash(double[] value) { 1064 return Arrays.hashCode(value); 1065 } 1066 1067 /** 1068 * For a byte value x, this method returns (byte)(+1) if x >= 0 and 1069 * (byte)(-1) if x < 0. 1070 * 1071 * @param x the value, a byte 1072 * @return (byte)(+1) or (byte)(-1), depending on the sign of x 1073 */ 1074 public static byte indicator(final byte x) { 1075 return (x >= ZB) ? PB : NB; 1076 } 1077 1078 /** 1079 * For a double precision value x, this method returns +1.0 if x >= 0 and 1080 * -1.0 if x < 0. Returns <code>NaN</code> if <code>x</code> is 1081 * <code>NaN</code>. 1082 * 1083 * @param x the value, a double 1084 * @return +1.0 or -1.0, depending on the sign of x 1085 */ 1086 public static double indicator(final double x) { 1087 if (Double.isNaN(x)) { 1088 return Double.NaN; 1089 } 1090 return (x >= 0.0) ? 1.0 : -1.0; 1091 } 1092 1093 /** 1094 * For a float value x, this method returns +1.0F if x >= 0 and -1.0F if x < 1095 * 0. Returns <code>NaN</code> if <code>x</code> is <code>NaN</code>. 1096 * 1097 * @param x the value, a float 1098 * @return +1.0F or -1.0F, depending on the sign of x 1099 */ 1100 public static float indicator(final float x) { 1101 if (Float.isNaN(x)) { 1102 return Float.NaN; 1103 } 1104 return (x >= 0.0F) ? 1.0F : -1.0F; 1105 } 1106 1107 /** 1108 * For an int value x, this method returns +1 if x >= 0 and -1 if x < 0. 1109 * 1110 * @param x the value, an int 1111 * @return +1 or -1, depending on the sign of x 1112 */ 1113 public static int indicator(final int x) { 1114 return (x >= 0) ? 1 : -1; 1115 } 1116 1117 /** 1118 * For a long value x, this method returns +1L if x >= 0 and -1L if x < 0. 1119 * 1120 * @param x the value, a long 1121 * @return +1L or -1L, depending on the sign of x 1122 */ 1123 public static long indicator(final long x) { 1124 return (x >= 0L) ? 1L : -1L; 1125 } 1126 1127 /** 1128 * For a short value x, this method returns (short)(+1) if x >= 0 and 1129 * (short)(-1) if x < 0. 1130 * 1131 * @param x the value, a short 1132 * @return (short)(+1) or (short)(-1), depending on the sign of x 1133 */ 1134 public static short indicator(final short x) { 1135 return (x >= ZS) ? PS : NS; 1136 } 1137 1138 /** 1139 * <p> 1140 * Returns the least common multiple of the absolute value of two numbers, 1141 * using the formula <code>lcm(a,b) = (a / gcd(a,b)) * b</code>. 1142 * </p> 1143 * Special cases: 1144 * <ul> 1145 * <li>The invocations <code>lcm(Integer.MIN_VALUE, n)</code> and 1146 * <code>lcm(n, Integer.MIN_VALUE)</code>, where <code>abs(n)</code> is a 1147 * power of 2, throw an <code>ArithmeticException</code>, because the result 1148 * would be 2^31, which is too large for an int value.</li> 1149 * <li>The result of <code>lcm(0, x)</code> and <code>lcm(x, 0)</code> is 1150 * <code>0</code> for any <code>x</code>. 1151 * </ul> 1152 * 1153 * @param a any number 1154 * @param b any number 1155 * @return the least common multiple, never negative 1156 * @throws ArithmeticException 1157 * if the result cannot be represented as a nonnegative int 1158 * value 1159 * @since 1.1 1160 */ 1161 public static int lcm(int a, int b) { 1162 if (a==0 || b==0){ 1163 return 0; 1164 } 1165 int lcm = FastMath.abs(mulAndCheck(a / gcd(a, b), b)); 1166 if (lcm == Integer.MIN_VALUE) { 1167 throw MathRuntimeException.createArithmeticException( 1168 LocalizedFormats.LCM_OVERFLOW_32_BITS, 1169 a, b); 1170 } 1171 return lcm; 1172 } 1173 1174 /** 1175 * <p> 1176 * Returns the least common multiple of the absolute value of two numbers, 1177 * using the formula <code>lcm(a,b) = (a / gcd(a,b)) * b</code>. 1178 * </p> 1179 * Special cases: 1180 * <ul> 1181 * <li>The invocations <code>lcm(Long.MIN_VALUE, n)</code> and 1182 * <code>lcm(n, Long.MIN_VALUE)</code>, where <code>abs(n)</code> is a 1183 * power of 2, throw an <code>ArithmeticException</code>, because the result 1184 * would be 2^63, which is too large for an int value.</li> 1185 * <li>The result of <code>lcm(0L, x)</code> and <code>lcm(x, 0L)</code> is 1186 * <code>0L</code> for any <code>x</code>. 1187 * </ul> 1188 * 1189 * @param a any number 1190 * @param b any number 1191 * @return the least common multiple, never negative 1192 * @throws ArithmeticException if the result cannot be represented as a nonnegative long 1193 * value 1194 * @since 2.1 1195 */ 1196 public static long lcm(long a, long b) { 1197 if (a==0 || b==0){ 1198 return 0; 1199 } 1200 long lcm = FastMath.abs(mulAndCheck(a / gcd(a, b), b)); 1201 if (lcm == Long.MIN_VALUE){ 1202 throw MathRuntimeException.createArithmeticException( 1203 LocalizedFormats.LCM_OVERFLOW_64_BITS, 1204 a, b); 1205 } 1206 return lcm; 1207 } 1208 1209 /** 1210 * <p>Returns the 1211 * <a href="http://mathworld.wolfram.com/Logarithm.html">logarithm</a> 1212 * for base <code>b</code> of <code>x</code>. 1213 * </p> 1214 * <p>Returns <code>NaN<code> if either argument is negative. If 1215 * <code>base</code> is 0 and <code>x</code> is positive, 0 is returned. 1216 * If <code>base</code> is positive and <code>x</code> is 0, 1217 * <code>Double.NEGATIVE_INFINITY</code> is returned. If both arguments 1218 * are 0, the result is <code>NaN</code>.</p> 1219 * 1220 * @param base the base of the logarithm, must be greater than 0 1221 * @param x argument, must be greater than 0 1222 * @return the value of the logarithm - the number y such that base^y = x. 1223 * @since 1.2 1224 */ 1225 public static double log(double base, double x) { 1226 return FastMath.log(x)/FastMath.log(base); 1227 } 1228 1229 /** 1230 * Multiply two integers, checking for overflow. 1231 * 1232 * @param x a factor 1233 * @param y a factor 1234 * @return the product <code>x*y</code> 1235 * @throws ArithmeticException if the result can not be represented as an 1236 * int 1237 * @since 1.1 1238 */ 1239 public static int mulAndCheck(int x, int y) { 1240 long m = ((long)x) * ((long)y); 1241 if (m < Integer.MIN_VALUE || m > Integer.MAX_VALUE) { 1242 throw new ArithmeticException("overflow: mul"); 1243 } 1244 return (int)m; 1245 } 1246 1247 /** 1248 * Multiply two long integers, checking for overflow. 1249 * 1250 * @param a first value 1251 * @param b second value 1252 * @return the product <code>a * b</code> 1253 * @throws ArithmeticException if the result can not be represented as an 1254 * long 1255 * @since 1.2 1256 */ 1257 public static long mulAndCheck(long a, long b) { 1258 long ret; 1259 String msg = "overflow: multiply"; 1260 if (a > b) { 1261 // use symmetry to reduce boundary cases 1262 ret = mulAndCheck(b, a); 1263 } else { 1264 if (a < 0) { 1265 if (b < 0) { 1266 // check for positive overflow with negative a, negative b 1267 if (a >= Long.MAX_VALUE / b) { 1268 ret = a * b; 1269 } else { 1270 throw new ArithmeticException(msg); 1271 } 1272 } else if (b > 0) { 1273 // check for negative overflow with negative a, positive b 1274 if (Long.MIN_VALUE / b <= a) { 1275 ret = a * b; 1276 } else { 1277 throw new ArithmeticException(msg); 1278 1279 } 1280 } else { 1281 // assert b == 0 1282 ret = 0; 1283 } 1284 } else if (a > 0) { 1285 // assert a > 0 1286 // assert b > 0 1287 1288 // check for positive overflow with positive a, positive b 1289 if (a <= Long.MAX_VALUE / b) { 1290 ret = a * b; 1291 } else { 1292 throw new ArithmeticException(msg); 1293 } 1294 } else { 1295 // assert a == 0 1296 ret = 0; 1297 } 1298 } 1299 return ret; 1300 } 1301 1302 /** 1303 * Get the next machine representable number after a number, moving 1304 * in the direction of another number. 1305 * <p> 1306 * If <code>direction</code> is greater than or equal to<code>d</code>, 1307 * the smallest machine representable number strictly greater than 1308 * <code>d</code> is returned; otherwise the largest representable number 1309 * strictly less than <code>d</code> is returned.</p> 1310 * <p> 1311 * If <code>d</code> is NaN or Infinite, it is returned unchanged.</p> 1312 * 1313 * @param d base number 1314 * @param direction (the only important thing is whether 1315 * direction is greater or smaller than d) 1316 * @return the next machine representable number in the specified direction 1317 * @since 1.2 1318 * @deprecated as of 2.2, replaced by {@link FastMath#nextAfter(double, double)} 1319 * which handles Infinities differently, and returns direction if d and direction compare equal. 1320 */ 1321 @Deprecated 1322 public static double nextAfter(double d, double direction) { 1323 1324 // handling of some important special cases 1325 if (Double.isNaN(d) || Double.isInfinite(d)) { 1326 return d; 1327 } else if (d == 0) { 1328 return (direction < 0) ? -Double.MIN_VALUE : Double.MIN_VALUE; 1329 } 1330 // special cases MAX_VALUE to infinity and MIN_VALUE to 0 1331 // are handled just as normal numbers 1332 1333 // split the double in raw components 1334 long bits = Double.doubleToLongBits(d); 1335 long sign = bits & 0x8000000000000000L; 1336 long exponent = bits & 0x7ff0000000000000L; 1337 long mantissa = bits & 0x000fffffffffffffL; 1338 1339 if (d * (direction - d) >= 0) { 1340 // we should increase the mantissa 1341 if (mantissa == 0x000fffffffffffffL) { 1342 return Double.longBitsToDouble(sign | 1343 (exponent + 0x0010000000000000L)); 1344 } else { 1345 return Double.longBitsToDouble(sign | 1346 exponent | (mantissa + 1)); 1347 } 1348 } else { 1349 // we should decrease the mantissa 1350 if (mantissa == 0L) { 1351 return Double.longBitsToDouble(sign | 1352 (exponent - 0x0010000000000000L) | 1353 0x000fffffffffffffL); 1354 } else { 1355 return Double.longBitsToDouble(sign | 1356 exponent | (mantissa - 1)); 1357 } 1358 } 1359 } 1360 1361 /** 1362 * Scale a number by 2<sup>scaleFactor</sup>. 1363 * <p>If <code>d</code> is 0 or NaN or Infinite, it is returned unchanged.</p> 1364 * 1365 * @param d base number 1366 * @param scaleFactor power of two by which d should be multiplied 1367 * @return d × 2<sup>scaleFactor</sup> 1368 * @since 2.0 1369 * @deprecated as of 2.2, replaced by {@link FastMath#scalb(double, int)} 1370 */ 1371 @Deprecated 1372 public static double scalb(final double d, final int scaleFactor) { 1373 return FastMath.scalb(d, scaleFactor); 1374 } 1375 1376 /** 1377 * Normalize an angle in a 2&pi wide interval around a center value. 1378 * <p>This method has three main uses:</p> 1379 * <ul> 1380 * <li>normalize an angle between 0 and 2π:<br/> 1381 * <code>a = MathUtils.normalizeAngle(a, FastMath.PI);</code></li> 1382 * <li>normalize an angle between -π and +π<br/> 1383 * <code>a = MathUtils.normalizeAngle(a, 0.0);</code></li> 1384 * <li>compute the angle between two defining angular positions:<br> 1385 * <code>angle = MathUtils.normalizeAngle(end, start) - start;</code></li> 1386 * </ul> 1387 * <p>Note that due to numerical accuracy and since π cannot be represented 1388 * exactly, the result interval is <em>closed</em>, it cannot be half-closed 1389 * as would be more satisfactory in a purely mathematical view.</p> 1390 * @param a angle to normalize 1391 * @param center center of the desired 2π interval for the result 1392 * @return a-2kπ with integer k and center-π <= a-2kπ <= center+π 1393 * @since 1.2 1394 */ 1395 public static double normalizeAngle(double a, double center) { 1396 return a - TWO_PI * FastMath.floor((a + FastMath.PI - center) / TWO_PI); 1397 } 1398 1399 /** 1400 * <p>Normalizes an array to make it sum to a specified value. 1401 * Returns the result of the transformation <pre> 1402 * x |-> x * normalizedSum / sum 1403 * </pre> 1404 * applied to each non-NaN element x of the input array, where sum is the 1405 * sum of the non-NaN entries in the input array.</p> 1406 * 1407 * <p>Throws IllegalArgumentException if <code>normalizedSum</code> is infinite 1408 * or NaN and ArithmeticException if the input array contains any infinite elements 1409 * or sums to 0</p> 1410 * 1411 * <p>Ignores (i.e., copies unchanged to the output array) NaNs in the input array.</p> 1412 * 1413 * @param values input array to be normalized 1414 * @param normalizedSum target sum for the normalized array 1415 * @return normalized array 1416 * @throws ArithmeticException if the input array contains infinite elements or sums to zero 1417 * @throws IllegalArgumentException if the target sum is infinite or NaN 1418 * @since 2.1 1419 */ 1420 public static double[] normalizeArray(double[] values, double normalizedSum) 1421 throws ArithmeticException, IllegalArgumentException { 1422 if (Double.isInfinite(normalizedSum)) { 1423 throw MathRuntimeException.createIllegalArgumentException( 1424 LocalizedFormats.NORMALIZE_INFINITE); 1425 } 1426 if (Double.isNaN(normalizedSum)) { 1427 throw MathRuntimeException.createIllegalArgumentException( 1428 LocalizedFormats.NORMALIZE_NAN); 1429 } 1430 double sum = 0d; 1431 final int len = values.length; 1432 double[] out = new double[len]; 1433 for (int i = 0; i < len; i++) { 1434 if (Double.isInfinite(values[i])) { 1435 throw MathRuntimeException.createArithmeticException( 1436 LocalizedFormats.INFINITE_ARRAY_ELEMENT, values[i], i); 1437 } 1438 if (!Double.isNaN(values[i])) { 1439 sum += values[i]; 1440 } 1441 } 1442 if (sum == 0) { 1443 throw MathRuntimeException.createArithmeticException(LocalizedFormats.ARRAY_SUMS_TO_ZERO); 1444 } 1445 for (int i = 0; i < len; i++) { 1446 if (Double.isNaN(values[i])) { 1447 out[i] = Double.NaN; 1448 } else { 1449 out[i] = values[i] * normalizedSum / sum; 1450 } 1451 } 1452 return out; 1453 } 1454 1455 /** 1456 * Round the given value to the specified number of decimal places. The 1457 * value is rounded using the {@link BigDecimal#ROUND_HALF_UP} method. 1458 * 1459 * @param x the value to round. 1460 * @param scale the number of digits to the right of the decimal point. 1461 * @return the rounded value. 1462 * @since 1.1 1463 */ 1464 public static double round(double x, int scale) { 1465 return round(x, scale, BigDecimal.ROUND_HALF_UP); 1466 } 1467 1468 /** 1469 * Round the given value to the specified number of decimal places. The 1470 * value is rounded using the given method which is any method defined in 1471 * {@link BigDecimal}. 1472 * 1473 * @param x the value to round. 1474 * @param scale the number of digits to the right of the decimal point. 1475 * @param roundingMethod the rounding method as defined in 1476 * {@link BigDecimal}. 1477 * @return the rounded value. 1478 * @since 1.1 1479 */ 1480 public static double round(double x, int scale, int roundingMethod) { 1481 try { 1482 return (new BigDecimal 1483 (Double.toString(x)) 1484 .setScale(scale, roundingMethod)) 1485 .doubleValue(); 1486 } catch (NumberFormatException ex) { 1487 if (Double.isInfinite(x)) { 1488 return x; 1489 } else { 1490 return Double.NaN; 1491 } 1492 } 1493 } 1494 1495 /** 1496 * Round the given value to the specified number of decimal places. The 1497 * value is rounding using the {@link BigDecimal#ROUND_HALF_UP} method. 1498 * 1499 * @param x the value to round. 1500 * @param scale the number of digits to the right of the decimal point. 1501 * @return the rounded value. 1502 * @since 1.1 1503 */ 1504 public static float round(float x, int scale) { 1505 return round(x, scale, BigDecimal.ROUND_HALF_UP); 1506 } 1507 1508 /** 1509 * Round the given value to the specified number of decimal places. The 1510 * value is rounded using the given method which is any method defined in 1511 * {@link BigDecimal}. 1512 * 1513 * @param x the value to round. 1514 * @param scale the number of digits to the right of the decimal point. 1515 * @param roundingMethod the rounding method as defined in 1516 * {@link BigDecimal}. 1517 * @return the rounded value. 1518 * @since 1.1 1519 */ 1520 public static float round(float x, int scale, int roundingMethod) { 1521 float sign = indicator(x); 1522 float factor = (float)FastMath.pow(10.0f, scale) * sign; 1523 return (float)roundUnscaled(x * factor, sign, roundingMethod) / factor; 1524 } 1525 1526 /** 1527 * Round the given non-negative, value to the "nearest" integer. Nearest is 1528 * determined by the rounding method specified. Rounding methods are defined 1529 * in {@link BigDecimal}. 1530 * 1531 * @param unscaled the value to round. 1532 * @param sign the sign of the original, scaled value. 1533 * @param roundingMethod the rounding method as defined in 1534 * {@link BigDecimal}. 1535 * @return the rounded value. 1536 * @since 1.1 1537 */ 1538 private static double roundUnscaled(double unscaled, double sign, 1539 int roundingMethod) { 1540 switch (roundingMethod) { 1541 case BigDecimal.ROUND_CEILING : 1542 if (sign == -1) { 1543 unscaled = FastMath.floor(nextAfter(unscaled, Double.NEGATIVE_INFINITY)); 1544 } else { 1545 unscaled = FastMath.ceil(nextAfter(unscaled, Double.POSITIVE_INFINITY)); 1546 } 1547 break; 1548 case BigDecimal.ROUND_DOWN : 1549 unscaled = FastMath.floor(nextAfter(unscaled, Double.NEGATIVE_INFINITY)); 1550 break; 1551 case BigDecimal.ROUND_FLOOR : 1552 if (sign == -1) { 1553 unscaled = FastMath.ceil(nextAfter(unscaled, Double.POSITIVE_INFINITY)); 1554 } else { 1555 unscaled = FastMath.floor(nextAfter(unscaled, Double.NEGATIVE_INFINITY)); 1556 } 1557 break; 1558 case BigDecimal.ROUND_HALF_DOWN : { 1559 unscaled = nextAfter(unscaled, Double.NEGATIVE_INFINITY); 1560 double fraction = unscaled - FastMath.floor(unscaled); 1561 if (fraction > 0.5) { 1562 unscaled = FastMath.ceil(unscaled); 1563 } else { 1564 unscaled = FastMath.floor(unscaled); 1565 } 1566 break; 1567 } 1568 case BigDecimal.ROUND_HALF_EVEN : { 1569 double fraction = unscaled - FastMath.floor(unscaled); 1570 if (fraction > 0.5) { 1571 unscaled = FastMath.ceil(unscaled); 1572 } else if (fraction < 0.5) { 1573 unscaled = FastMath.floor(unscaled); 1574 } else { 1575 // The following equality test is intentional and needed for rounding purposes 1576 if (FastMath.floor(unscaled) / 2.0 == FastMath.floor(Math 1577 .floor(unscaled) / 2.0)) { // even 1578 unscaled = FastMath.floor(unscaled); 1579 } else { // odd 1580 unscaled = FastMath.ceil(unscaled); 1581 } 1582 } 1583 break; 1584 } 1585 case BigDecimal.ROUND_HALF_UP : { 1586 unscaled = nextAfter(unscaled, Double.POSITIVE_INFINITY); 1587 double fraction = unscaled - FastMath.floor(unscaled); 1588 if (fraction >= 0.5) { 1589 unscaled = FastMath.ceil(unscaled); 1590 } else { 1591 unscaled = FastMath.floor(unscaled); 1592 } 1593 break; 1594 } 1595 case BigDecimal.ROUND_UNNECESSARY : 1596 if (unscaled != FastMath.floor(unscaled)) { 1597 throw new ArithmeticException("Inexact result from rounding"); 1598 } 1599 break; 1600 case BigDecimal.ROUND_UP : 1601 unscaled = FastMath.ceil(nextAfter(unscaled, Double.POSITIVE_INFINITY)); 1602 break; 1603 default : 1604 throw MathRuntimeException.createIllegalArgumentException( 1605 LocalizedFormats.INVALID_ROUNDING_METHOD, 1606 roundingMethod, 1607 "ROUND_CEILING", BigDecimal.ROUND_CEILING, 1608 "ROUND_DOWN", BigDecimal.ROUND_DOWN, 1609 "ROUND_FLOOR", BigDecimal.ROUND_FLOOR, 1610 "ROUND_HALF_DOWN", BigDecimal.ROUND_HALF_DOWN, 1611 "ROUND_HALF_EVEN", BigDecimal.ROUND_HALF_EVEN, 1612 "ROUND_HALF_UP", BigDecimal.ROUND_HALF_UP, 1613 "ROUND_UNNECESSARY", BigDecimal.ROUND_UNNECESSARY, 1614 "ROUND_UP", BigDecimal.ROUND_UP); 1615 } 1616 return unscaled; 1617 } 1618 1619 /** 1620 * Returns the <a href="http://mathworld.wolfram.com/Sign.html"> sign</a> 1621 * for byte value <code>x</code>. 1622 * <p> 1623 * For a byte value x, this method returns (byte)(+1) if x > 0, (byte)(0) if 1624 * x = 0, and (byte)(-1) if x < 0.</p> 1625 * 1626 * @param x the value, a byte 1627 * @return (byte)(+1), (byte)(0), or (byte)(-1), depending on the sign of x 1628 */ 1629 public static byte sign(final byte x) { 1630 return (x == ZB) ? ZB : (x > ZB) ? PB : NB; 1631 } 1632 1633 /** 1634 * Returns the <a href="http://mathworld.wolfram.com/Sign.html"> sign</a> 1635 * for double precision <code>x</code>. 1636 * <p> 1637 * For a double value <code>x</code>, this method returns 1638 * <code>+1.0</code> if <code>x > 0</code>, <code>0.0</code> if 1639 * <code>x = 0.0</code>, and <code>-1.0</code> if <code>x < 0</code>. 1640 * Returns <code>NaN</code> if <code>x</code> is <code>NaN</code>.</p> 1641 * 1642 * @param x the value, a double 1643 * @return +1.0, 0.0, or -1.0, depending on the sign of x 1644 */ 1645 public static double sign(final double x) { 1646 if (Double.isNaN(x)) { 1647 return Double.NaN; 1648 } 1649 return (x == 0.0) ? 0.0 : (x > 0.0) ? 1.0 : -1.0; 1650 } 1651 1652 /** 1653 * Returns the <a href="http://mathworld.wolfram.com/Sign.html"> sign</a> 1654 * for float value <code>x</code>. 1655 * <p> 1656 * For a float value x, this method returns +1.0F if x > 0, 0.0F if x = 1657 * 0.0F, and -1.0F if x < 0. Returns <code>NaN</code> if <code>x</code> 1658 * is <code>NaN</code>.</p> 1659 * 1660 * @param x the value, a float 1661 * @return +1.0F, 0.0F, or -1.0F, depending on the sign of x 1662 */ 1663 public static float sign(final float x) { 1664 if (Float.isNaN(x)) { 1665 return Float.NaN; 1666 } 1667 return (x == 0.0F) ? 0.0F : (x > 0.0F) ? 1.0F : -1.0F; 1668 } 1669 1670 /** 1671 * Returns the <a href="http://mathworld.wolfram.com/Sign.html"> sign</a> 1672 * for int value <code>x</code>. 1673 * <p> 1674 * For an int value x, this method returns +1 if x > 0, 0 if x = 0, and -1 1675 * if x < 0.</p> 1676 * 1677 * @param x the value, an int 1678 * @return +1, 0, or -1, depending on the sign of x 1679 */ 1680 public static int sign(final int x) { 1681 return (x == 0) ? 0 : (x > 0) ? 1 : -1; 1682 } 1683 1684 /** 1685 * Returns the <a href="http://mathworld.wolfram.com/Sign.html"> sign</a> 1686 * for long value <code>x</code>. 1687 * <p> 1688 * For a long value x, this method returns +1L if x > 0, 0L if x = 0, and 1689 * -1L if x < 0.</p> 1690 * 1691 * @param x the value, a long 1692 * @return +1L, 0L, or -1L, depending on the sign of x 1693 */ 1694 public static long sign(final long x) { 1695 return (x == 0L) ? 0L : (x > 0L) ? 1L : -1L; 1696 } 1697 1698 /** 1699 * Returns the <a href="http://mathworld.wolfram.com/Sign.html"> sign</a> 1700 * for short value <code>x</code>. 1701 * <p> 1702 * For a short value x, this method returns (short)(+1) if x > 0, (short)(0) 1703 * if x = 0, and (short)(-1) if x < 0.</p> 1704 * 1705 * @param x the value, a short 1706 * @return (short)(+1), (short)(0), or (short)(-1), depending on the sign of 1707 * x 1708 */ 1709 public static short sign(final short x) { 1710 return (x == ZS) ? ZS : (x > ZS) ? PS : NS; 1711 } 1712 1713 /** 1714 * Returns the <a href="http://mathworld.wolfram.com/HyperbolicSine.html"> 1715 * hyperbolic sine</a> of x. 1716 * 1717 * @param x double value for which to find the hyperbolic sine 1718 * @return hyperbolic sine of x 1719 */ 1720 public static double sinh(double x) { 1721 return (FastMath.exp(x) - FastMath.exp(-x)) / 2.0; 1722 } 1723 1724 /** 1725 * Subtract two integers, checking for overflow. 1726 * 1727 * @param x the minuend 1728 * @param y the subtrahend 1729 * @return the difference <code>x-y</code> 1730 * @throws ArithmeticException if the result can not be represented as an 1731 * int 1732 * @since 1.1 1733 */ 1734 public static int subAndCheck(int x, int y) { 1735 long s = (long)x - (long)y; 1736 if (s < Integer.MIN_VALUE || s > Integer.MAX_VALUE) { 1737 throw MathRuntimeException.createArithmeticException(LocalizedFormats.OVERFLOW_IN_SUBTRACTION, x, y); 1738 } 1739 return (int)s; 1740 } 1741 1742 /** 1743 * Subtract two long integers, checking for overflow. 1744 * 1745 * @param a first value 1746 * @param b second value 1747 * @return the difference <code>a-b</code> 1748 * @throws ArithmeticException if the result can not be represented as an 1749 * long 1750 * @since 1.2 1751 */ 1752 public static long subAndCheck(long a, long b) { 1753 long ret; 1754 String msg = "overflow: subtract"; 1755 if (b == Long.MIN_VALUE) { 1756 if (a < 0) { 1757 ret = a - b; 1758 } else { 1759 throw new ArithmeticException(msg); 1760 } 1761 } else { 1762 // use additive inverse 1763 ret = addAndCheck(a, -b, LocalizedFormats.OVERFLOW_IN_ADDITION); 1764 } 1765 return ret; 1766 } 1767 1768 /** 1769 * Raise an int to an int power. 1770 * @param k number to raise 1771 * @param e exponent (must be positive or null) 1772 * @return k<sup>e</sup> 1773 * @exception IllegalArgumentException if e is negative 1774 */ 1775 public static int pow(final int k, int e) 1776 throws IllegalArgumentException { 1777 1778 if (e < 0) { 1779 throw MathRuntimeException.createIllegalArgumentException( 1780 LocalizedFormats.POWER_NEGATIVE_PARAMETERS, 1781 k, e); 1782 } 1783 1784 int result = 1; 1785 int k2p = k; 1786 while (e != 0) { 1787 if ((e & 0x1) != 0) { 1788 result *= k2p; 1789 } 1790 k2p *= k2p; 1791 e = e >> 1; 1792 } 1793 1794 return result; 1795 1796 } 1797 1798 /** 1799 * Raise an int to a long power. 1800 * @param k number to raise 1801 * @param e exponent (must be positive or null) 1802 * @return k<sup>e</sup> 1803 * @exception IllegalArgumentException if e is negative 1804 */ 1805 public static int pow(final int k, long e) 1806 throws IllegalArgumentException { 1807 1808 if (e < 0) { 1809 throw MathRuntimeException.createIllegalArgumentException( 1810 LocalizedFormats.POWER_NEGATIVE_PARAMETERS, 1811 k, e); 1812 } 1813 1814 int result = 1; 1815 int k2p = k; 1816 while (e != 0) { 1817 if ((e & 0x1) != 0) { 1818 result *= k2p; 1819 } 1820 k2p *= k2p; 1821 e = e >> 1; 1822 } 1823 1824 return result; 1825 1826 } 1827 1828 /** 1829 * Raise a long to an int power. 1830 * @param k number to raise 1831 * @param e exponent (must be positive or null) 1832 * @return k<sup>e</sup> 1833 * @exception IllegalArgumentException if e is negative 1834 */ 1835 public static long pow(final long k, int e) 1836 throws IllegalArgumentException { 1837 1838 if (e < 0) { 1839 throw MathRuntimeException.createIllegalArgumentException( 1840 LocalizedFormats.POWER_NEGATIVE_PARAMETERS, 1841 k, e); 1842 } 1843 1844 long result = 1l; 1845 long k2p = k; 1846 while (e != 0) { 1847 if ((e & 0x1) != 0) { 1848 result *= k2p; 1849 } 1850 k2p *= k2p; 1851 e = e >> 1; 1852 } 1853 1854 return result; 1855 1856 } 1857 1858 /** 1859 * Raise a long to a long power. 1860 * @param k number to raise 1861 * @param e exponent (must be positive or null) 1862 * @return k<sup>e</sup> 1863 * @exception IllegalArgumentException if e is negative 1864 */ 1865 public static long pow(final long k, long e) 1866 throws IllegalArgumentException { 1867 1868 if (e < 0) { 1869 throw MathRuntimeException.createIllegalArgumentException( 1870 LocalizedFormats.POWER_NEGATIVE_PARAMETERS, 1871 k, e); 1872 } 1873 1874 long result = 1l; 1875 long k2p = k; 1876 while (e != 0) { 1877 if ((e & 0x1) != 0) { 1878 result *= k2p; 1879 } 1880 k2p *= k2p; 1881 e = e >> 1; 1882 } 1883 1884 return result; 1885 1886 } 1887 1888 /** 1889 * Raise a BigInteger to an int power. 1890 * @param k number to raise 1891 * @param e exponent (must be positive or null) 1892 * @return k<sup>e</sup> 1893 * @exception IllegalArgumentException if e is negative 1894 */ 1895 public static BigInteger pow(final BigInteger k, int e) 1896 throws IllegalArgumentException { 1897 1898 if (e < 0) { 1899 throw MathRuntimeException.createIllegalArgumentException( 1900 LocalizedFormats.POWER_NEGATIVE_PARAMETERS, 1901 k, e); 1902 } 1903 1904 return k.pow(e); 1905 1906 } 1907 1908 /** 1909 * Raise a BigInteger to a long power. 1910 * @param k number to raise 1911 * @param e exponent (must be positive or null) 1912 * @return k<sup>e</sup> 1913 * @exception IllegalArgumentException if e is negative 1914 */ 1915 public static BigInteger pow(final BigInteger k, long e) 1916 throws IllegalArgumentException { 1917 1918 if (e < 0) { 1919 throw MathRuntimeException.createIllegalArgumentException( 1920 LocalizedFormats.POWER_NEGATIVE_PARAMETERS, 1921 k, e); 1922 } 1923 1924 BigInteger result = BigInteger.ONE; 1925 BigInteger k2p = k; 1926 while (e != 0) { 1927 if ((e & 0x1) != 0) { 1928 result = result.multiply(k2p); 1929 } 1930 k2p = k2p.multiply(k2p); 1931 e = e >> 1; 1932 } 1933 1934 return result; 1935 1936 } 1937 1938 /** 1939 * Raise a BigInteger to a BigInteger power. 1940 * @param k number to raise 1941 * @param e exponent (must be positive or null) 1942 * @return k<sup>e</sup> 1943 * @exception IllegalArgumentException if e is negative 1944 */ 1945 public static BigInteger pow(final BigInteger k, BigInteger e) 1946 throws IllegalArgumentException { 1947 1948 if (e.compareTo(BigInteger.ZERO) < 0) { 1949 throw MathRuntimeException.createIllegalArgumentException( 1950 LocalizedFormats.POWER_NEGATIVE_PARAMETERS, 1951 k, e); 1952 } 1953 1954 BigInteger result = BigInteger.ONE; 1955 BigInteger k2p = k; 1956 while (!BigInteger.ZERO.equals(e)) { 1957 if (e.testBit(0)) { 1958 result = result.multiply(k2p); 1959 } 1960 k2p = k2p.multiply(k2p); 1961 e = e.shiftRight(1); 1962 } 1963 1964 return result; 1965 1966 } 1967 1968 /** 1969 * Calculates the L<sub>1</sub> (sum of abs) distance between two points. 1970 * 1971 * @param p1 the first point 1972 * @param p2 the second point 1973 * @return the L<sub>1</sub> distance between the two points 1974 */ 1975 public static double distance1(double[] p1, double[] p2) { 1976 double sum = 0; 1977 for (int i = 0; i < p1.length; i++) { 1978 sum += FastMath.abs(p1[i] - p2[i]); 1979 } 1980 return sum; 1981 } 1982 1983 /** 1984 * Calculates the L<sub>1</sub> (sum of abs) distance between two points. 1985 * 1986 * @param p1 the first point 1987 * @param p2 the second point 1988 * @return the L<sub>1</sub> distance between the two points 1989 */ 1990 public static int distance1(int[] p1, int[] p2) { 1991 int sum = 0; 1992 for (int i = 0; i < p1.length; i++) { 1993 sum += FastMath.abs(p1[i] - p2[i]); 1994 } 1995 return sum; 1996 } 1997 1998 /** 1999 * Calculates the L<sub>2</sub> (Euclidean) distance between two points. 2000 * 2001 * @param p1 the first point 2002 * @param p2 the second point 2003 * @return the L<sub>2</sub> distance between the two points 2004 */ 2005 public static double distance(double[] p1, double[] p2) { 2006 double sum = 0; 2007 for (int i = 0; i < p1.length; i++) { 2008 final double dp = p1[i] - p2[i]; 2009 sum += dp * dp; 2010 } 2011 return FastMath.sqrt(sum); 2012 } 2013 2014 /** 2015 * Calculates the L<sub>2</sub> (Euclidean) distance between two points. 2016 * 2017 * @param p1 the first point 2018 * @param p2 the second point 2019 * @return the L<sub>2</sub> distance between the two points 2020 */ 2021 public static double distance(int[] p1, int[] p2) { 2022 double sum = 0; 2023 for (int i = 0; i < p1.length; i++) { 2024 final double dp = p1[i] - p2[i]; 2025 sum += dp * dp; 2026 } 2027 return FastMath.sqrt(sum); 2028 } 2029 2030 /** 2031 * Calculates the L<sub>∞</sub> (max of abs) distance between two points. 2032 * 2033 * @param p1 the first point 2034 * @param p2 the second point 2035 * @return the L<sub>∞</sub> distance between the two points 2036 */ 2037 public static double distanceInf(double[] p1, double[] p2) { 2038 double max = 0; 2039 for (int i = 0; i < p1.length; i++) { 2040 max = FastMath.max(max, FastMath.abs(p1[i] - p2[i])); 2041 } 2042 return max; 2043 } 2044 2045 /** 2046 * Calculates the L<sub>∞</sub> (max of abs) distance between two points. 2047 * 2048 * @param p1 the first point 2049 * @param p2 the second point 2050 * @return the L<sub>∞</sub> distance between the two points 2051 */ 2052 public static int distanceInf(int[] p1, int[] p2) { 2053 int max = 0; 2054 for (int i = 0; i < p1.length; i++) { 2055 max = FastMath.max(max, FastMath.abs(p1[i] - p2[i])); 2056 } 2057 return max; 2058 } 2059 2060 /** 2061 * Specification of ordering direction. 2062 */ 2063 public static enum OrderDirection { 2064 /** Constant for increasing direction. */ 2065 INCREASING, 2066 /** Constant for decreasing direction. */ 2067 DECREASING 2068 } 2069 2070 /** 2071 * Checks that the given array is sorted. 2072 * 2073 * @param val Values. 2074 * @param dir Ordering direction. 2075 * @param strict Whether the order should be strict. 2076 * @throws NonMonotonousSequenceException if the array is not sorted. 2077 * @since 2.2 2078 */ 2079 public static void checkOrder(double[] val, OrderDirection dir, boolean strict) { 2080 double previous = val[0]; 2081 boolean ok = true; 2082 2083 int max = val.length; 2084 for (int i = 1; i < max; i++) { 2085 switch (dir) { 2086 case INCREASING: 2087 if (strict) { 2088 if (val[i] <= previous) { 2089 ok = false; 2090 } 2091 } else { 2092 if (val[i] < previous) { 2093 ok = false; 2094 } 2095 } 2096 break; 2097 case DECREASING: 2098 if (strict) { 2099 if (val[i] >= previous) { 2100 ok = false; 2101 } 2102 } else { 2103 if (val[i] > previous) { 2104 ok = false; 2105 } 2106 } 2107 break; 2108 default: 2109 // Should never happen. 2110 throw new IllegalArgumentException(); 2111 } 2112 2113 if (!ok) { 2114 throw new NonMonotonousSequenceException(val[i], previous, i, dir, strict); 2115 } 2116 previous = val[i]; 2117 } 2118 } 2119 2120 /** 2121 * Checks that the given array is sorted in strictly increasing order. 2122 * 2123 * @param val Values. 2124 * @throws NonMonotonousSequenceException if the array is not sorted. 2125 * @since 2.2 2126 */ 2127 public static void checkOrder(double[] val) { 2128 checkOrder(val, OrderDirection.INCREASING, true); 2129 } 2130 2131 /** 2132 * Checks that the given array is sorted. 2133 * 2134 * @param val Values 2135 * @param dir Order direction (-1 for decreasing, 1 for increasing) 2136 * @param strict Whether the order should be strict 2137 * @throws NonMonotonousSequenceException if the array is not sorted. 2138 * @deprecated as of 2.2 (please use the new {@link #checkOrder(double[],OrderDirection,boolean) 2139 * checkOrder} method). To be removed in 3.0. 2140 */ 2141 @Deprecated 2142 public static void checkOrder(double[] val, int dir, boolean strict) { 2143 if (dir > 0) { 2144 checkOrder(val, OrderDirection.INCREASING, strict); 2145 } else { 2146 checkOrder(val, OrderDirection.DECREASING, strict); 2147 } 2148 } 2149 2150 /** 2151 * Returns the Cartesian norm (2-norm), handling both overflow and underflow. 2152 * Translation of the minpack enorm subroutine. 2153 * 2154 * The redistribution policy for MINPACK is available <a 2155 * href="http://www.netlib.org/minpack/disclaimer">here</a>, for convenience, it 2156 * is reproduced below.</p> 2157 * 2158 * <table border="0" width="80%" cellpadding="10" align="center" bgcolor="#E0E0E0"> 2159 * <tr><td> 2160 * Minpack Copyright Notice (1999) University of Chicago. 2161 * All rights reserved 2162 * </td></tr> 2163 * <tr><td> 2164 * Redistribution and use in source and binary forms, with or without 2165 * modification, are permitted provided that the following conditions 2166 * are met: 2167 * <ol> 2168 * <li>Redistributions of source code must retain the above copyright 2169 * notice, this list of conditions and the following disclaimer.</li> 2170 * <li>Redistributions in binary form must reproduce the above 2171 * copyright notice, this list of conditions and the following 2172 * disclaimer in the documentation and/or other materials provided 2173 * with the distribution.</li> 2174 * <li>The end-user documentation included with the redistribution, if any, 2175 * must include the following acknowledgment: 2176 * <code>This product includes software developed by the University of 2177 * Chicago, as Operator of Argonne National Laboratory.</code> 2178 * Alternately, this acknowledgment may appear in the software itself, 2179 * if and wherever such third-party acknowledgments normally appear.</li> 2180 * <li><strong>WARRANTY DISCLAIMER. THE SOFTWARE IS SUPPLIED "AS IS" 2181 * WITHOUT WARRANTY OF ANY KIND. THE COPYRIGHT HOLDER, THE 2182 * UNITED STATES, THE UNITED STATES DEPARTMENT OF ENERGY, AND 2183 * THEIR EMPLOYEES: (1) DISCLAIM ANY WARRANTIES, EXPRESS OR 2184 * IMPLIED, INCLUDING BUT NOT LIMITED TO ANY IMPLIED WARRANTIES 2185 * OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, TITLE 2186 * OR NON-INFRINGEMENT, (2) DO NOT ASSUME ANY LEGAL LIABILITY 2187 * OR RESPONSIBILITY FOR THE ACCURACY, COMPLETENESS, OR 2188 * USEFULNESS OF THE SOFTWARE, (3) DO NOT REPRESENT THAT USE OF 2189 * THE SOFTWARE WOULD NOT INFRINGE PRIVATELY OWNED RIGHTS, (4) 2190 * DO NOT WARRANT THAT THE SOFTWARE WILL FUNCTION 2191 * UNINTERRUPTED, THAT IT IS ERROR-FREE OR THAT ANY ERRORS WILL 2192 * BE CORRECTED.</strong></li> 2193 * <li><strong>LIMITATION OF LIABILITY. IN NO EVENT WILL THE COPYRIGHT 2194 * HOLDER, THE UNITED STATES, THE UNITED STATES DEPARTMENT OF 2195 * ENERGY, OR THEIR EMPLOYEES: BE LIABLE FOR ANY INDIRECT, 2196 * INCIDENTAL, CONSEQUENTIAL, SPECIAL OR PUNITIVE DAMAGES OF 2197 * ANY KIND OR NATURE, INCLUDING BUT NOT LIMITED TO LOSS OF 2198 * PROFITS OR LOSS OF DATA, FOR ANY REASON WHATSOEVER, WHETHER 2199 * SUCH LIABILITY IS ASSERTED ON THE BASIS OF CONTRACT, TORT 2200 * (INCLUDING NEGLIGENCE OR STRICT LIABILITY), OR OTHERWISE, 2201 * EVEN IF ANY OF SAID PARTIES HAS BEEN WARNED OF THE 2202 * POSSIBILITY OF SUCH LOSS OR DAMAGES.</strong></li> 2203 * <ol></td></tr> 2204 * </table> 2205 * 2206 * @param v vector of doubles 2207 * @return the 2-norm of the vector 2208 * @since 2.2 2209 */ 2210 public static double safeNorm(double[] v) { 2211 double rdwarf = 3.834e-20; 2212 double rgiant = 1.304e+19; 2213 double s1=0.0; 2214 double s2=0.0; 2215 double s3=0.0; 2216 double x1max = 0.0; 2217 double x3max = 0.0; 2218 double floatn = (double)v.length; 2219 double agiant = rgiant/floatn; 2220 for (int i=0;i<v.length;i++) { 2221 double xabs = Math.abs(v[i]); 2222 if (xabs<rdwarf || xabs>agiant) { 2223 if (xabs>rdwarf) { 2224 if (xabs>x1max) { 2225 double r=x1max/xabs; 2226 s1=1.0+s1*r*r; 2227 x1max=xabs; 2228 } else { 2229 double r=xabs/x1max; 2230 s1+=r*r; 2231 } 2232 } else { 2233 if (xabs>x3max) { 2234 double r=x3max/xabs; 2235 s3=1.0+s3*r*r; 2236 x3max=xabs; 2237 } else { 2238 if (xabs!=0.0) { 2239 double r=xabs/x3max; 2240 s3+=r*r; 2241 } 2242 } 2243 } 2244 } else { 2245 s2+=xabs*xabs; 2246 } 2247 } 2248 double norm; 2249 if (s1!=0.0) { 2250 norm = x1max*Math.sqrt(s1+(s2/x1max)/x1max); 2251 } else { 2252 if (s2==0.0) { 2253 norm = x3max*Math.sqrt(s3); 2254 } else { 2255 if (s2>=x3max) { 2256 norm = Math.sqrt(s2*(1.0+(x3max/s2)*(x3max*s3))); 2257 } else { 2258 norm = Math.sqrt(x3max*((s2/x3max)+(x3max*s3))); 2259 } 2260 } 2261 } 2262 return norm; 2263 } 2264 2265 } 2266
