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.math.util; 019 020 import java.math.BigDecimal; 021 import java.math.BigInteger; 022 import java.util.Arrays; 023 024 import org.apache.commons.math.MathRuntimeException; 025 026 /** 027 * Some useful additions to the built-in functions in {@link Math}. 028 * @version $Revision: 927249 $ $Date: 2010-03-24 21:06:51 -0400 (Wed, 24 Mar 2010) $ 029 */ 030 public final class MathUtils { 031 032 /** Smallest positive number such that 1 - EPSILON is not numerically equal to 1. */ 033 public static final double EPSILON = 0x1.0p-53; 034 035 /** Safe minimum, such that 1 / SAFE_MIN does not overflow. 036 * <p>In IEEE 754 arithmetic, this is also the smallest normalized 037 * number 2<sup>-1022</sup>.</p> 038 */ 039 public static final double SAFE_MIN = 0x1.0p-1022; 040 041 /** 042 * 2 π. 043 * @since 2.1 044 */ 045 public static final double TWO_PI = 2 * Math.PI; 046 047 /** -1.0 cast as a byte. */ 048 private static final byte NB = (byte)-1; 049 050 /** -1.0 cast as a short. */ 051 private static final short NS = (short)-1; 052 053 /** 1.0 cast as a byte. */ 054 private static final byte PB = (byte)1; 055 056 /** 1.0 cast as a short. */ 057 private static final short PS = (short)1; 058 059 /** 0.0 cast as a byte. */ 060 private static final byte ZB = (byte)0; 061 062 /** 0.0 cast as a short. */ 063 private static final short ZS = (short)0; 064 065 /** Gap between NaN and regular numbers. */ 066 private static final int NAN_GAP = 4 * 1024 * 1024; 067 068 /** Offset to order signed double numbers lexicographically. */ 069 private static final long SGN_MASK = 0x8000000000000000L; 070 071 /** All long-representable factorials */ 072 private static final long[] FACTORIALS = new long[] { 073 1l, 1l, 2l, 074 6l, 24l, 120l, 075 720l, 5040l, 40320l, 076 362880l, 3628800l, 39916800l, 077 479001600l, 6227020800l, 87178291200l, 078 1307674368000l, 20922789888000l, 355687428096000l, 079 6402373705728000l, 121645100408832000l, 2432902008176640000l }; 080 081 /** 082 * Private Constructor 083 */ 084 private MathUtils() { 085 super(); 086 } 087 088 /** 089 * Add two integers, checking for overflow. 090 * 091 * @param x an addend 092 * @param y an addend 093 * @return the sum <code>x+y</code> 094 * @throws ArithmeticException if the result can not be represented as an 095 * int 096 * @since 1.1 097 */ 098 public static int addAndCheck(int x, int y) { 099 long s = (long)x + (long)y; 100 if (s < Integer.MIN_VALUE || s > Integer.MAX_VALUE) { 101 throw new ArithmeticException("overflow: add"); 102 } 103 return (int)s; 104 } 105 106 /** 107 * Add two long integers, checking for overflow. 108 * 109 * @param a an addend 110 * @param b an addend 111 * @return the sum <code>a+b</code> 112 * @throws ArithmeticException if the result can not be represented as an 113 * long 114 * @since 1.2 115 */ 116 public static long addAndCheck(long a, long b) { 117 return addAndCheck(a, b, "overflow: add"); 118 } 119 120 /** 121 * Add two long integers, checking for overflow. 122 * 123 * @param a an addend 124 * @param b an addend 125 * @param msg the message to use for any thrown exception. 126 * @return the sum <code>a+b</code> 127 * @throws ArithmeticException if the result can not be represented as an 128 * long 129 * @since 1.2 130 */ 131 private static long addAndCheck(long a, long b, String msg) { 132 long ret; 133 if (a > b) { 134 // use symmetry to reduce boundary cases 135 ret = addAndCheck(b, a, msg); 136 } else { 137 // assert a <= b 138 139 if (a < 0) { 140 if (b < 0) { 141 // check for negative overflow 142 if (Long.MIN_VALUE - b <= a) { 143 ret = a + b; 144 } else { 145 throw new ArithmeticException(msg); 146 } 147 } else { 148 // opposite sign addition is always safe 149 ret = a + b; 150 } 151 } else { 152 // assert a >= 0 153 // assert b >= 0 154 155 // check for positive overflow 156 if (a <= Long.MAX_VALUE - b) { 157 ret = a + b; 158 } else { 159 throw new ArithmeticException(msg); 160 } 161 } 162 } 163 return ret; 164 } 165 166 /** 167 * Returns an exact representation of the <a 168 * href="http://mathworld.wolfram.com/BinomialCoefficient.html"> Binomial 169 * Coefficient</a>, "<code>n choose k</code>", the number of 170 * <code>k</code>-element subsets that can be selected from an 171 * <code>n</code>-element set. 172 * <p> 173 * <Strong>Preconditions</strong>: 174 * <ul> 175 * <li> <code>0 <= k <= n </code> (otherwise 176 * <code>IllegalArgumentException</code> is thrown)</li> 177 * <li> The result is small enough to fit into a <code>long</code>. The 178 * largest value of <code>n</code> for which all coefficients are 179 * <code> < Long.MAX_VALUE</code> is 66. If the computed value exceeds 180 * <code>Long.MAX_VALUE</code> an <code>ArithMeticException</code> is 181 * thrown.</li> 182 * </ul></p> 183 * 184 * @param n the size of the set 185 * @param k the size of the subsets to be counted 186 * @return <code>n choose k</code> 187 * @throws IllegalArgumentException if preconditions are not met. 188 * @throws ArithmeticException if the result is too large to be represented 189 * by a long integer. 190 */ 191 public static long binomialCoefficient(final int n, final int k) { 192 checkBinomial(n, k); 193 if ((n == k) || (k == 0)) { 194 return 1; 195 } 196 if ((k == 1) || (k == n - 1)) { 197 return n; 198 } 199 // Use symmetry for large k 200 if (k > n / 2) 201 return binomialCoefficient(n, n - k); 202 203 // We use the formula 204 // (n choose k) = n! / (n-k)! / k! 205 // (n choose k) == ((n-k+1)*...*n) / (1*...*k) 206 // which could be written 207 // (n choose k) == (n-1 choose k-1) * n / k 208 long result = 1; 209 if (n <= 61) { 210 // For n <= 61, the naive implementation cannot overflow. 211 int i = n - k + 1; 212 for (int j = 1; j <= k; j++) { 213 result = result * i / j; 214 i++; 215 } 216 } else if (n <= 66) { 217 // For n > 61 but n <= 66, the result cannot overflow, 218 // but we must take care not to overflow intermediate values. 219 int i = n - k + 1; 220 for (int j = 1; j <= k; j++) { 221 // We know that (result * i) is divisible by j, 222 // but (result * i) may overflow, so we split j: 223 // Filter out the gcd, d, so j/d and i/d are integer. 224 // result is divisible by (j/d) because (j/d) 225 // is relative prime to (i/d) and is a divisor of 226 // result * (i/d). 227 final long d = gcd(i, j); 228 result = (result / (j / d)) * (i / d); 229 i++; 230 } 231 } else { 232 // For n > 66, a result overflow might occur, so we check 233 // the multiplication, taking care to not overflow 234 // unnecessary. 235 int i = n - k + 1; 236 for (int j = 1; j <= k; j++) { 237 final long d = gcd(i, j); 238 result = mulAndCheck(result / (j / d), i / d); 239 i++; 240 } 241 } 242 return result; 243 } 244 245 /** 246 * Returns a <code>double</code> representation of the <a 247 * href="http://mathworld.wolfram.com/BinomialCoefficient.html"> Binomial 248 * Coefficient</a>, "<code>n choose k</code>", the number of 249 * <code>k</code>-element subsets that can be selected from an 250 * <code>n</code>-element set. 251 * <p> 252 * <Strong>Preconditions</strong>: 253 * <ul> 254 * <li> <code>0 <= k <= n </code> (otherwise 255 * <code>IllegalArgumentException</code> is thrown)</li> 256 * <li> The result is small enough to fit into a <code>double</code>. The 257 * largest value of <code>n</code> for which all coefficients are < 258 * Double.MAX_VALUE is 1029. If the computed value exceeds Double.MAX_VALUE, 259 * Double.POSITIVE_INFINITY is returned</li> 260 * </ul></p> 261 * 262 * @param n the size of the set 263 * @param k the size of the subsets to be counted 264 * @return <code>n choose k</code> 265 * @throws IllegalArgumentException if preconditions are not met. 266 */ 267 public static double binomialCoefficientDouble(final int n, final int k) { 268 checkBinomial(n, k); 269 if ((n == k) || (k == 0)) { 270 return 1d; 271 } 272 if ((k == 1) || (k == n - 1)) { 273 return n; 274 } 275 if (k > n/2) { 276 return binomialCoefficientDouble(n, n - k); 277 } 278 if (n < 67) { 279 return binomialCoefficient(n,k); 280 } 281 282 double result = 1d; 283 for (int i = 1; i <= k; i++) { 284 result *= (double)(n - k + i) / (double)i; 285 } 286 287 return Math.floor(result + 0.5); 288 } 289 290 /** 291 * Returns the natural <code>log</code> of the <a 292 * href="http://mathworld.wolfram.com/BinomialCoefficient.html"> Binomial 293 * Coefficient</a>, "<code>n choose k</code>", the number of 294 * <code>k</code>-element subsets that can be selected from an 295 * <code>n</code>-element set. 296 * <p> 297 * <Strong>Preconditions</strong>: 298 * <ul> 299 * <li> <code>0 <= k <= n </code> (otherwise 300 * <code>IllegalArgumentException</code> is thrown)</li> 301 * </ul></p> 302 * 303 * @param n the size of the set 304 * @param k the size of the subsets to be counted 305 * @return <code>n choose k</code> 306 * @throws IllegalArgumentException if preconditions are not met. 307 */ 308 public static double binomialCoefficientLog(final int n, final int k) { 309 checkBinomial(n, k); 310 if ((n == k) || (k == 0)) { 311 return 0; 312 } 313 if ((k == 1) || (k == n - 1)) { 314 return Math.log(n); 315 } 316 317 /* 318 * For values small enough to do exact integer computation, 319 * return the log of the exact value 320 */ 321 if (n < 67) { 322 return Math.log(binomialCoefficient(n,k)); 323 } 324 325 /* 326 * Return the log of binomialCoefficientDouble for values that will not 327 * overflow binomialCoefficientDouble 328 */ 329 if (n < 1030) { 330 return Math.log(binomialCoefficientDouble(n, k)); 331 } 332 333 if (k > n / 2) { 334 return binomialCoefficientLog(n, n - k); 335 } 336 337 /* 338 * Sum logs for values that could overflow 339 */ 340 double logSum = 0; 341 342 // n!/(n-k)! 343 for (int i = n - k + 1; i <= n; i++) { 344 logSum += Math.log(i); 345 } 346 347 // divide by k! 348 for (int i = 2; i <= k; i++) { 349 logSum -= Math.log(i); 350 } 351 352 return logSum; 353 } 354 355 /** 356 * Check binomial preconditions. 357 * @param n the size of the set 358 * @param k the size of the subsets to be counted 359 * @exception IllegalArgumentException if preconditions are not met. 360 */ 361 private static void checkBinomial(final int n, final int k) 362 throws IllegalArgumentException { 363 if (n < k) { 364 throw MathRuntimeException.createIllegalArgumentException( 365 "must have n >= k for binomial coefficient (n,k), got n = {0}, k = {1}", 366 n, k); 367 } 368 if (n < 0) { 369 throw MathRuntimeException.createIllegalArgumentException( 370 "must have n >= 0 for binomial coefficient (n,k), got n = {0}", 371 n); 372 } 373 } 374 375 /** 376 * Compares two numbers given some amount of allowed error. 377 * 378 * @param x the first number 379 * @param y the second number 380 * @param eps the amount of error to allow when checking for equality 381 * @return <ul><li>0 if {@link #equals(double, double, double) equals(x, y, eps)}</li> 382 * <li>< 0 if !{@link #equals(double, double, double) equals(x, y, eps)} && x < y</li> 383 * <li>> 0 if !{@link #equals(double, double, double) equals(x, y, eps)} && x > y</li></ul> 384 */ 385 public static int compareTo(double x, double y, double eps) { 386 if (equals(x, y, eps)) { 387 return 0; 388 } else if (x < y) { 389 return -1; 390 } 391 return 1; 392 } 393 394 /** 395 * Returns the <a href="http://mathworld.wolfram.com/HyperbolicCosine.html"> 396 * hyperbolic cosine</a> of x. 397 * 398 * @param x double value for which to find the hyperbolic cosine 399 * @return hyperbolic cosine of x 400 */ 401 public static double cosh(double x) { 402 return (Math.exp(x) + Math.exp(-x)) / 2.0; 403 } 404 405 /** 406 * Returns true iff both arguments are NaN or neither is NaN and they are 407 * equal 408 * 409 * @param x first value 410 * @param y second value 411 * @return true if the values are equal or both are NaN 412 */ 413 public static boolean equals(double x, double y) { 414 return (Double.isNaN(x) && Double.isNaN(y)) || x == y; 415 } 416 417 /** 418 * Returns true iff both arguments are equal or within the range of allowed 419 * error (inclusive). 420 * <p> 421 * Two NaNs are considered equals, as are two infinities with same sign. 422 * </p> 423 * 424 * @param x first value 425 * @param y second value 426 * @param eps the amount of absolute error to allow 427 * @return true if the values are equal or within range of each other 428 */ 429 public static boolean equals(double x, double y, double eps) { 430 return equals(x, y) || (Math.abs(y - x) <= eps); 431 } 432 433 /** 434 * Returns true iff both arguments are equal or within the range of allowed 435 * error (inclusive). 436 * Adapted from <a 437 * href="http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm"> 438 * Bruce Dawson</a> 439 * 440 * @param x first value 441 * @param y second value 442 * @param maxUlps {@code (maxUlps - 1)} is the number of floating point 443 * values between {@code x} and {@code y}. 444 * @return {@code true} if there are less than {@code maxUlps} floating 445 * point values between {@code x} and {@code y} 446 */ 447 public static boolean equals(double x, double y, int maxUlps) { 448 // Check that "maxUlps" is non-negative and small enough so that the 449 // default NAN won't compare as equal to anything. 450 assert maxUlps > 0 && maxUlps < NAN_GAP; 451 452 long xInt = Double.doubleToLongBits(x); 453 long yInt = Double.doubleToLongBits(y); 454 455 // Make lexicographically ordered as a two's-complement integer. 456 if (xInt < 0) { 457 xInt = SGN_MASK - xInt; 458 } 459 if (yInt < 0) { 460 yInt = SGN_MASK - yInt; 461 } 462 463 return Math.abs(xInt - yInt) <= maxUlps; 464 } 465 466 /** 467 * Returns true iff both arguments are null or have same dimensions 468 * and all their elements are {@link #equals(double,double) equals} 469 * 470 * @param x first array 471 * @param y second array 472 * @return true if the values are both null or have same dimension 473 * and equal elements 474 * @since 1.2 475 */ 476 public static boolean equals(double[] x, double[] y) { 477 if ((x == null) || (y == null)) { 478 return !((x == null) ^ (y == null)); 479 } 480 if (x.length != y.length) { 481 return false; 482 } 483 for (int i = 0; i < x.length; ++i) { 484 if (!equals(x[i], y[i])) { 485 return false; 486 } 487 } 488 return true; 489 } 490 491 /** 492 * Returns n!. Shorthand for <code>n</code> <a 493 * href="http://mathworld.wolfram.com/Factorial.html"> Factorial</a>, the 494 * product of the numbers <code>1,...,n</code>. 495 * <p> 496 * <Strong>Preconditions</strong>: 497 * <ul> 498 * <li> <code>n >= 0</code> (otherwise 499 * <code>IllegalArgumentException</code> is thrown)</li> 500 * <li> The result is small enough to fit into a <code>long</code>. The 501 * largest value of <code>n</code> for which <code>n!</code> < 502 * Long.MAX_VALUE</code> is 20. If the computed value exceeds <code>Long.MAX_VALUE</code> 503 * an <code>ArithMeticException </code> is thrown.</li> 504 * </ul> 505 * </p> 506 * 507 * @param n argument 508 * @return <code>n!</code> 509 * @throws ArithmeticException if the result is too large to be represented 510 * by a long integer. 511 * @throws IllegalArgumentException if n < 0 512 */ 513 public static long factorial(final int n) { 514 if (n < 0) { 515 throw MathRuntimeException.createIllegalArgumentException( 516 "must have n >= 0 for n!, got n = {0}", 517 n); 518 } 519 if (n > 20) { 520 throw new ArithmeticException( 521 "factorial value is too large to fit in a long"); 522 } 523 return FACTORIALS[n]; 524 } 525 526 /** 527 * Returns n!. Shorthand for <code>n</code> <a 528 * href="http://mathworld.wolfram.com/Factorial.html"> Factorial</a>, the 529 * product of the numbers <code>1,...,n</code> as a <code>double</code>. 530 * <p> 531 * <Strong>Preconditions</strong>: 532 * <ul> 533 * <li> <code>n >= 0</code> (otherwise 534 * <code>IllegalArgumentException</code> is thrown)</li> 535 * <li> The result is small enough to fit into a <code>double</code>. The 536 * largest value of <code>n</code> for which <code>n!</code> < 537 * Double.MAX_VALUE</code> is 170. If the computed value exceeds 538 * Double.MAX_VALUE, Double.POSITIVE_INFINITY is returned</li> 539 * </ul> 540 * </p> 541 * 542 * @param n argument 543 * @return <code>n!</code> 544 * @throws IllegalArgumentException if n < 0 545 */ 546 public static double factorialDouble(final int n) { 547 if (n < 0) { 548 throw MathRuntimeException.createIllegalArgumentException( 549 "must have n >= 0 for n!, got n = {0}", 550 n); 551 } 552 if (n < 21) { 553 return factorial(n); 554 } 555 return Math.floor(Math.exp(factorialLog(n)) + 0.5); 556 } 557 558 /** 559 * Returns the natural logarithm of n!. 560 * <p> 561 * <Strong>Preconditions</strong>: 562 * <ul> 563 * <li> <code>n >= 0</code> (otherwise 564 * <code>IllegalArgumentException</code> is thrown)</li> 565 * </ul></p> 566 * 567 * @param n argument 568 * @return <code>n!</code> 569 * @throws IllegalArgumentException if preconditions are not met. 570 */ 571 public static double factorialLog(final int n) { 572 if (n < 0) { 573 throw MathRuntimeException.createIllegalArgumentException( 574 "must have n >= 0 for n!, got n = {0}", 575 n); 576 } 577 if (n < 21) { 578 return Math.log(factorial(n)); 579 } 580 double logSum = 0; 581 for (int i = 2; i <= n; i++) { 582 logSum += Math.log(i); 583 } 584 return logSum; 585 } 586 587 /** 588 * <p> 589 * Gets the greatest common divisor of the absolute value of two numbers, 590 * using the "binary gcd" method which avoids division and modulo 591 * operations. See Knuth 4.5.2 algorithm B. This algorithm is due to Josef 592 * Stein (1961). 593 * </p> 594 * Special cases: 595 * <ul> 596 * <li>The invocations 597 * <code>gcd(Integer.MIN_VALUE, Integer.MIN_VALUE)</code>, 598 * <code>gcd(Integer.MIN_VALUE, 0)</code> and 599 * <code>gcd(0, Integer.MIN_VALUE)</code> throw an 600 * <code>ArithmeticException</code>, because the result would be 2^31, which 601 * is too large for an int value.</li> 602 * <li>The result of <code>gcd(x, x)</code>, <code>gcd(0, x)</code> and 603 * <code>gcd(x, 0)</code> is the absolute value of <code>x</code>, except 604 * for the special cases above. 605 * <li>The invocation <code>gcd(0, 0)</code> is the only one which returns 606 * <code>0</code>.</li> 607 * </ul> 608 * 609 * @param p any number 610 * @param q any number 611 * @return the greatest common divisor, never negative 612 * @throws ArithmeticException if the result cannot be represented as a 613 * nonnegative int value 614 * @since 1.1 615 */ 616 public static int gcd(final int p, final int q) { 617 int u = p; 618 int v = q; 619 if ((u == 0) || (v == 0)) { 620 if ((u == Integer.MIN_VALUE) || (v == Integer.MIN_VALUE)) { 621 throw MathRuntimeException.createArithmeticException( 622 "overflow: gcd({0}, {1}) is 2^31", 623 p, q); 624 } 625 return Math.abs(u) + Math.abs(v); 626 } 627 // keep u and v negative, as negative integers range down to 628 // -2^31, while positive numbers can only be as large as 2^31-1 629 // (i.e. we can't necessarily negate a negative number without 630 // overflow) 631 /* assert u!=0 && v!=0; */ 632 if (u > 0) { 633 u = -u; 634 } // make u negative 635 if (v > 0) { 636 v = -v; 637 } // make v negative 638 // B1. [Find power of 2] 639 int k = 0; 640 while ((u & 1) == 0 && (v & 1) == 0 && k < 31) { // while u and v are 641 // both even... 642 u /= 2; 643 v /= 2; 644 k++; // cast out twos. 645 } 646 if (k == 31) { 647 throw MathRuntimeException.createArithmeticException( 648 "overflow: gcd({0}, {1}) is 2^31", 649 p, q); 650 } 651 // B2. Initialize: u and v have been divided by 2^k and at least 652 // one is odd. 653 int t = ((u & 1) == 1) ? v : -(u / 2)/* B3 */; 654 // t negative: u was odd, v may be even (t replaces v) 655 // t positive: u was even, v is odd (t replaces u) 656 do { 657 /* assert u<0 && v<0; */ 658 // B4/B3: cast out twos from t. 659 while ((t & 1) == 0) { // while t is even.. 660 t /= 2; // cast out twos 661 } 662 // B5 [reset max(u,v)] 663 if (t > 0) { 664 u = -t; 665 } else { 666 v = t; 667 } 668 // B6/B3. at this point both u and v should be odd. 669 t = (v - u) / 2; 670 // |u| larger: t positive (replace u) 671 // |v| larger: t negative (replace v) 672 } while (t != 0); 673 return -u * (1 << k); // gcd is u*2^k 674 } 675 676 /** 677 * <p> 678 * Gets the greatest common divisor of the absolute value of two numbers, 679 * using the "binary gcd" method which avoids division and modulo 680 * operations. See Knuth 4.5.2 algorithm B. This algorithm is due to Josef 681 * Stein (1961). 682 * </p> 683 * Special cases: 684 * <ul> 685 * <li>The invocations 686 * <code>gcd(Long.MIN_VALUE, Long.MIN_VALUE)</code>, 687 * <code>gcd(Long.MIN_VALUE, 0L)</code> and 688 * <code>gcd(0L, Long.MIN_VALUE)</code> throw an 689 * <code>ArithmeticException</code>, because the result would be 2^63, which 690 * is too large for a long value.</li> 691 * <li>The result of <code>gcd(x, x)</code>, <code>gcd(0L, x)</code> and 692 * <code>gcd(x, 0L)</code> is the absolute value of <code>x</code>, except 693 * for the special cases above. 694 * <li>The invocation <code>gcd(0L, 0L)</code> is the only one which returns 695 * <code>0L</code>.</li> 696 * </ul> 697 * 698 * @param p any number 699 * @param q any number 700 * @return the greatest common divisor, never negative 701 * @throws ArithmeticException if the result cannot be represented as a nonnegative long 702 * value 703 * @since 2.1 704 */ 705 public static long gcd(final long p, final long q) { 706 long u = p; 707 long v = q; 708 if ((u == 0) || (v == 0)) { 709 if ((u == Long.MIN_VALUE) || (v == Long.MIN_VALUE)){ 710 throw MathRuntimeException.createArithmeticException( 711 "overflow: gcd({0}, {1}) is 2^63", 712 p, q); 713 } 714 return Math.abs(u) + Math.abs(v); 715 } 716 // keep u and v negative, as negative integers range down to 717 // -2^63, while positive numbers can only be as large as 2^63-1 718 // (i.e. we can't necessarily negate a negative number without 719 // overflow) 720 /* assert u!=0 && v!=0; */ 721 if (u > 0) { 722 u = -u; 723 } // make u negative 724 if (v > 0) { 725 v = -v; 726 } // make v negative 727 // B1. [Find power of 2] 728 int k = 0; 729 while ((u & 1) == 0 && (v & 1) == 0 && k < 63) { // while u and v are 730 // both even... 731 u /= 2; 732 v /= 2; 733 k++; // cast out twos. 734 } 735 if (k == 63) { 736 throw MathRuntimeException.createArithmeticException( 737 "overflow: gcd({0}, {1}) is 2^63", 738 p, q); 739 } 740 // B2. Initialize: u and v have been divided by 2^k and at least 741 // one is odd. 742 long t = ((u & 1) == 1) ? v : -(u / 2)/* B3 */; 743 // t negative: u was odd, v may be even (t replaces v) 744 // t positive: u was even, v is odd (t replaces u) 745 do { 746 /* assert u<0 && v<0; */ 747 // B4/B3: cast out twos from t. 748 while ((t & 1) == 0) { // while t is even.. 749 t /= 2; // cast out twos 750 } 751 // B5 [reset max(u,v)] 752 if (t > 0) { 753 u = -t; 754 } else { 755 v = t; 756 } 757 // B6/B3. at this point both u and v should be odd. 758 t = (v - u) / 2; 759 // |u| larger: t positive (replace u) 760 // |v| larger: t negative (replace v) 761 } while (t != 0); 762 return -u * (1L << k); // gcd is u*2^k 763 } 764 765 /** 766 * Returns an integer hash code representing the given double value. 767 * 768 * @param value the value to be hashed 769 * @return the hash code 770 */ 771 public static int hash(double value) { 772 return new Double(value).hashCode(); 773 } 774 775 /** 776 * Returns an integer hash code representing the given double array. 777 * 778 * @param value the value to be hashed (may be null) 779 * @return the hash code 780 * @since 1.2 781 */ 782 public static int hash(double[] value) { 783 return Arrays.hashCode(value); 784 } 785 786 /** 787 * For a byte value x, this method returns (byte)(+1) if x >= 0 and 788 * (byte)(-1) if x < 0. 789 * 790 * @param x the value, a byte 791 * @return (byte)(+1) or (byte)(-1), depending on the sign of x 792 */ 793 public static byte indicator(final byte x) { 794 return (x >= ZB) ? PB : NB; 795 } 796 797 /** 798 * For a double precision value x, this method returns +1.0 if x >= 0 and 799 * -1.0 if x < 0. Returns <code>NaN</code> if <code>x</code> is 800 * <code>NaN</code>. 801 * 802 * @param x the value, a double 803 * @return +1.0 or -1.0, depending on the sign of x 804 */ 805 public static double indicator(final double x) { 806 if (Double.isNaN(x)) { 807 return Double.NaN; 808 } 809 return (x >= 0.0) ? 1.0 : -1.0; 810 } 811 812 /** 813 * For a float value x, this method returns +1.0F if x >= 0 and -1.0F if x < 814 * 0. Returns <code>NaN</code> if <code>x</code> is <code>NaN</code>. 815 * 816 * @param x the value, a float 817 * @return +1.0F or -1.0F, depending on the sign of x 818 */ 819 public static float indicator(final float x) { 820 if (Float.isNaN(x)) { 821 return Float.NaN; 822 } 823 return (x >= 0.0F) ? 1.0F : -1.0F; 824 } 825 826 /** 827 * For an int value x, this method returns +1 if x >= 0 and -1 if x < 0. 828 * 829 * @param x the value, an int 830 * @return +1 or -1, depending on the sign of x 831 */ 832 public static int indicator(final int x) { 833 return (x >= 0) ? 1 : -1; 834 } 835 836 /** 837 * For a long value x, this method returns +1L if x >= 0 and -1L if x < 0. 838 * 839 * @param x the value, a long 840 * @return +1L or -1L, depending on the sign of x 841 */ 842 public static long indicator(final long x) { 843 return (x >= 0L) ? 1L : -1L; 844 } 845 846 /** 847 * For a short value x, this method returns (short)(+1) if x >= 0 and 848 * (short)(-1) if x < 0. 849 * 850 * @param x the value, a short 851 * @return (short)(+1) or (short)(-1), depending on the sign of x 852 */ 853 public static short indicator(final short x) { 854 return (x >= ZS) ? PS : NS; 855 } 856 857 /** 858 * <p> 859 * Returns the least common multiple of the absolute value of two numbers, 860 * using the formula <code>lcm(a,b) = (a / gcd(a,b)) * b</code>. 861 * </p> 862 * Special cases: 863 * <ul> 864 * <li>The invocations <code>lcm(Integer.MIN_VALUE, n)</code> and 865 * <code>lcm(n, Integer.MIN_VALUE)</code>, where <code>abs(n)</code> is a 866 * power of 2, throw an <code>ArithmeticException</code>, because the result 867 * would be 2^31, which is too large for an int value.</li> 868 * <li>The result of <code>lcm(0, x)</code> and <code>lcm(x, 0)</code> is 869 * <code>0</code> for any <code>x</code>. 870 * </ul> 871 * 872 * @param a any number 873 * @param b any number 874 * @return the least common multiple, never negative 875 * @throws ArithmeticException 876 * if the result cannot be represented as a nonnegative int 877 * value 878 * @since 1.1 879 */ 880 public static int lcm(int a, int b) { 881 if (a==0 || b==0){ 882 return 0; 883 } 884 int lcm = Math.abs(mulAndCheck(a / gcd(a, b), b)); 885 if (lcm == Integer.MIN_VALUE) { 886 throw MathRuntimeException.createArithmeticException( 887 "overflow: lcm({0}, {1}) is 2^31", 888 a, b); 889 } 890 return lcm; 891 } 892 893 /** 894 * <p> 895 * Returns the least common multiple of the absolute value of two numbers, 896 * using the formula <code>lcm(a,b) = (a / gcd(a,b)) * b</code>. 897 * </p> 898 * Special cases: 899 * <ul> 900 * <li>The invocations <code>lcm(Long.MIN_VALUE, n)</code> and 901 * <code>lcm(n, Long.MIN_VALUE)</code>, where <code>abs(n)</code> is a 902 * power of 2, throw an <code>ArithmeticException</code>, because the result 903 * would be 2^63, which is too large for an int value.</li> 904 * <li>The result of <code>lcm(0L, x)</code> and <code>lcm(x, 0L)</code> is 905 * <code>0L</code> for any <code>x</code>. 906 * </ul> 907 * 908 * @param a any number 909 * @param b any number 910 * @return the least common multiple, never negative 911 * @throws ArithmeticException if the result cannot be represented as a nonnegative long 912 * value 913 * @since 2.1 914 */ 915 public static long lcm(long a, long b) { 916 if (a==0 || b==0){ 917 return 0; 918 } 919 long lcm = Math.abs(mulAndCheck(a / gcd(a, b), b)); 920 if (lcm == Long.MIN_VALUE){ 921 throw MathRuntimeException.createArithmeticException( 922 "overflow: lcm({0}, {1}) is 2^63", 923 a, b); 924 } 925 return lcm; 926 } 927 928 /** 929 * <p>Returns the 930 * <a href="http://mathworld.wolfram.com/Logarithm.html">logarithm</a> 931 * for base <code>b</code> of <code>x</code>. 932 * </p> 933 * <p>Returns <code>NaN<code> if either argument is negative. If 934 * <code>base</code> is 0 and <code>x</code> is positive, 0 is returned. 935 * If <code>base</code> is positive and <code>x</code> is 0, 936 * <code>Double.NEGATIVE_INFINITY</code> is returned. If both arguments 937 * are 0, the result is <code>NaN</code>.</p> 938 * 939 * @param base the base of the logarithm, must be greater than 0 940 * @param x argument, must be greater than 0 941 * @return the value of the logarithm - the number y such that base^y = x. 942 * @since 1.2 943 */ 944 public static double log(double base, double x) { 945 return Math.log(x)/Math.log(base); 946 } 947 948 /** 949 * Multiply two integers, checking for overflow. 950 * 951 * @param x a factor 952 * @param y a factor 953 * @return the product <code>x*y</code> 954 * @throws ArithmeticException if the result can not be represented as an 955 * int 956 * @since 1.1 957 */ 958 public static int mulAndCheck(int x, int y) { 959 long m = ((long)x) * ((long)y); 960 if (m < Integer.MIN_VALUE || m > Integer.MAX_VALUE) { 961 throw new ArithmeticException("overflow: mul"); 962 } 963 return (int)m; 964 } 965 966 /** 967 * Multiply two long integers, checking for overflow. 968 * 969 * @param a first value 970 * @param b second value 971 * @return the product <code>a * b</code> 972 * @throws ArithmeticException if the result can not be represented as an 973 * long 974 * @since 1.2 975 */ 976 public static long mulAndCheck(long a, long b) { 977 long ret; 978 String msg = "overflow: multiply"; 979 if (a > b) { 980 // use symmetry to reduce boundary cases 981 ret = mulAndCheck(b, a); 982 } else { 983 if (a < 0) { 984 if (b < 0) { 985 // check for positive overflow with negative a, negative b 986 if (a >= Long.MAX_VALUE / b) { 987 ret = a * b; 988 } else { 989 throw new ArithmeticException(msg); 990 } 991 } else if (b > 0) { 992 // check for negative overflow with negative a, positive b 993 if (Long.MIN_VALUE / b <= a) { 994 ret = a * b; 995 } else { 996 throw new ArithmeticException(msg); 997 998 } 999 } else { 1000 // assert b == 0 1001 ret = 0; 1002 } 1003 } else if (a > 0) { 1004 // assert a > 0 1005 // assert b > 0 1006 1007 // check for positive overflow with positive a, positive b 1008 if (a <= Long.MAX_VALUE / b) { 1009 ret = a * b; 1010 } else { 1011 throw new ArithmeticException(msg); 1012 } 1013 } else { 1014 // assert a == 0 1015 ret = 0; 1016 } 1017 } 1018 return ret; 1019 } 1020 1021 /** 1022 * Get the next machine representable number after a number, moving 1023 * in the direction of another number. 1024 * <p> 1025 * If <code>direction</code> is greater than or equal to<code>d</code>, 1026 * the smallest machine representable number strictly greater than 1027 * <code>d</code> is returned; otherwise the largest representable number 1028 * strictly less than <code>d</code> is returned.</p> 1029 * <p> 1030 * If <code>d</code> is NaN or Infinite, it is returned unchanged.</p> 1031 * 1032 * @param d base number 1033 * @param direction (the only important thing is whether 1034 * direction is greater or smaller than d) 1035 * @return the next machine representable number in the specified direction 1036 * @since 1.2 1037 */ 1038 public static double nextAfter(double d, double direction) { 1039 1040 // handling of some important special cases 1041 if (Double.isNaN(d) || Double.isInfinite(d)) { 1042 return d; 1043 } else if (d == 0) { 1044 return (direction < 0) ? -Double.MIN_VALUE : Double.MIN_VALUE; 1045 } 1046 // special cases MAX_VALUE to infinity and MIN_VALUE to 0 1047 // are handled just as normal numbers 1048 1049 // split the double in raw components 1050 long bits = Double.doubleToLongBits(d); 1051 long sign = bits & 0x8000000000000000L; 1052 long exponent = bits & 0x7ff0000000000000L; 1053 long mantissa = bits & 0x000fffffffffffffL; 1054 1055 if (d * (direction - d) >= 0) { 1056 // we should increase the mantissa 1057 if (mantissa == 0x000fffffffffffffL) { 1058 return Double.longBitsToDouble(sign | 1059 (exponent + 0x0010000000000000L)); 1060 } else { 1061 return Double.longBitsToDouble(sign | 1062 exponent | (mantissa + 1)); 1063 } 1064 } else { 1065 // we should decrease the mantissa 1066 if (mantissa == 0L) { 1067 return Double.longBitsToDouble(sign | 1068 (exponent - 0x0010000000000000L) | 1069 0x000fffffffffffffL); 1070 } else { 1071 return Double.longBitsToDouble(sign | 1072 exponent | (mantissa - 1)); 1073 } 1074 } 1075 1076 } 1077 1078 /** 1079 * Scale a number by 2<sup>scaleFactor</sup>. 1080 * <p>If <code>d</code> is 0 or NaN or Infinite, it is returned unchanged.</p> 1081 * 1082 * @param d base number 1083 * @param scaleFactor power of two by which d sould be multiplied 1084 * @return d × 2<sup>scaleFactor</sup> 1085 * @since 2.0 1086 */ 1087 public static double scalb(final double d, final int scaleFactor) { 1088 1089 // handling of some important special cases 1090 if ((d == 0) || Double.isNaN(d) || Double.isInfinite(d)) { 1091 return d; 1092 } 1093 1094 // split the double in raw components 1095 final long bits = Double.doubleToLongBits(d); 1096 final long exponent = bits & 0x7ff0000000000000L; 1097 final long rest = bits & 0x800fffffffffffffL; 1098 1099 // shift the exponent 1100 final long newBits = rest | (exponent + (((long) scaleFactor) << 52)); 1101 return Double.longBitsToDouble(newBits); 1102 1103 } 1104 1105 /** 1106 * Normalize an angle in a 2&pi wide interval around a center value. 1107 * <p>This method has three main uses:</p> 1108 * <ul> 1109 * <li>normalize an angle between 0 and 2π:<br/> 1110 * <code>a = MathUtils.normalizeAngle(a, Math.PI);</code></li> 1111 * <li>normalize an angle between -π and +π<br/> 1112 * <code>a = MathUtils.normalizeAngle(a, 0.0);</code></li> 1113 * <li>compute the angle between two defining angular positions:<br> 1114 * <code>angle = MathUtils.normalizeAngle(end, start) - start;</code></li> 1115 * </ul> 1116 * <p>Note that due to numerical accuracy and since π cannot be represented 1117 * exactly, the result interval is <em>closed</em>, it cannot be half-closed 1118 * as would be more satisfactory in a purely mathematical view.</p> 1119 * @param a angle to normalize 1120 * @param center center of the desired 2π interval for the result 1121 * @return a-2kπ with integer k and center-π <= a-2kπ <= center+π 1122 * @since 1.2 1123 */ 1124 public static double normalizeAngle(double a, double center) { 1125 return a - TWO_PI * Math.floor((a + Math.PI - center) / TWO_PI); 1126 } 1127 1128 /** 1129 * <p>Normalizes an array to make it sum to a specified value. 1130 * Returns the result of the transformation <pre> 1131 * x |-> x * normalizedSum / sum 1132 * </pre> 1133 * applied to each non-NaN element x of the input array, where sum is the 1134 * sum of the non-NaN entries in the input array.</p> 1135 * 1136 * <p>Throws IllegalArgumentException if <code>normalizedSum</code> is infinite 1137 * or NaN and ArithmeticException if the input array contains any infinite elements 1138 * or sums to 0</p> 1139 * 1140 * <p>Ignores (i.e., copies unchanged to the output array) NaNs in the input array.</p> 1141 * 1142 * @param values input array to be normalized 1143 * @param normalizedSum target sum for the normalized array 1144 * @return normalized array 1145 * @throws ArithmeticException if the input array contains infinite elements or sums to zero 1146 * @throws IllegalArgumentException if the target sum is infinite or NaN 1147 * @since 2.1 1148 */ 1149 public static double[] normalizeArray(double[] values, double normalizedSum) 1150 throws ArithmeticException, IllegalArgumentException { 1151 if (Double.isInfinite(normalizedSum)) { 1152 throw MathRuntimeException.createIllegalArgumentException( 1153 "Cannot normalize to an infinite value"); 1154 } 1155 if (Double.isNaN(normalizedSum)) { 1156 throw MathRuntimeException.createIllegalArgumentException( 1157 "Cannot normalize to NaN"); 1158 } 1159 double sum = 0d; 1160 final int len = values.length; 1161 double[] out = new double[len]; 1162 for (int i = 0; i < len; i++) { 1163 if (Double.isInfinite(values[i])) { 1164 throw MathRuntimeException.createArithmeticException( 1165 "Array contains an infinite element, {0} at index {1}", values[i], i); 1166 } 1167 if (!Double.isNaN(values[i])) { 1168 sum += values[i]; 1169 } 1170 } 1171 if (sum == 0) { 1172 throw MathRuntimeException.createArithmeticException( 1173 "Array sums to zero"); 1174 } 1175 for (int i = 0; i < len; i++) { 1176 if (Double.isNaN(values[i])) { 1177 out[i] = Double.NaN; 1178 } else { 1179 out[i] = values[i] * normalizedSum / sum; 1180 } 1181 } 1182 return out; 1183 } 1184 1185 /** 1186 * Round the given value to the specified number of decimal places. The 1187 * value is rounded using the {@link BigDecimal#ROUND_HALF_UP} method. 1188 * 1189 * @param x the value to round. 1190 * @param scale the number of digits to the right of the decimal point. 1191 * @return the rounded value. 1192 * @since 1.1 1193 */ 1194 public static double round(double x, int scale) { 1195 return round(x, scale, BigDecimal.ROUND_HALF_UP); 1196 } 1197 1198 /** 1199 * Round the given value to the specified number of decimal places. The 1200 * value is rounded using the given method which is any method defined in 1201 * {@link BigDecimal}. 1202 * 1203 * @param x the value to round. 1204 * @param scale the number of digits to the right of the decimal point. 1205 * @param roundingMethod the rounding method as defined in 1206 * {@link BigDecimal}. 1207 * @return the rounded value. 1208 * @since 1.1 1209 */ 1210 public static double round(double x, int scale, int roundingMethod) { 1211 try { 1212 return (new BigDecimal 1213 (Double.toString(x)) 1214 .setScale(scale, roundingMethod)) 1215 .doubleValue(); 1216 } catch (NumberFormatException ex) { 1217 if (Double.isInfinite(x)) { 1218 return x; 1219 } else { 1220 return Double.NaN; 1221 } 1222 } 1223 } 1224 1225 /** 1226 * Round the given value to the specified number of decimal places. The 1227 * value is rounding using the {@link BigDecimal#ROUND_HALF_UP} method. 1228 * 1229 * @param x the value to round. 1230 * @param scale the number of digits to the right of the decimal point. 1231 * @return the rounded value. 1232 * @since 1.1 1233 */ 1234 public static float round(float x, int scale) { 1235 return round(x, scale, BigDecimal.ROUND_HALF_UP); 1236 } 1237 1238 /** 1239 * Round the given value to the specified number of decimal places. The 1240 * value is rounded using the given method which is any method defined in 1241 * {@link BigDecimal}. 1242 * 1243 * @param x the value to round. 1244 * @param scale the number of digits to the right of the decimal point. 1245 * @param roundingMethod the rounding method as defined in 1246 * {@link BigDecimal}. 1247 * @return the rounded value. 1248 * @since 1.1 1249 */ 1250 public static float round(float x, int scale, int roundingMethod) { 1251 float sign = indicator(x); 1252 float factor = (float)Math.pow(10.0f, scale) * sign; 1253 return (float)roundUnscaled(x * factor, sign, roundingMethod) / factor; 1254 } 1255 1256 /** 1257 * Round the given non-negative, value to the "nearest" integer. Nearest is 1258 * determined by the rounding method specified. Rounding methods are defined 1259 * in {@link BigDecimal}. 1260 * 1261 * @param unscaled the value to round. 1262 * @param sign the sign of the original, scaled value. 1263 * @param roundingMethod the rounding method as defined in 1264 * {@link BigDecimal}. 1265 * @return the rounded value. 1266 * @since 1.1 1267 */ 1268 private static double roundUnscaled(double unscaled, double sign, 1269 int roundingMethod) { 1270 switch (roundingMethod) { 1271 case BigDecimal.ROUND_CEILING : 1272 if (sign == -1) { 1273 unscaled = Math.floor(nextAfter(unscaled, Double.NEGATIVE_INFINITY)); 1274 } else { 1275 unscaled = Math.ceil(nextAfter(unscaled, Double.POSITIVE_INFINITY)); 1276 } 1277 break; 1278 case BigDecimal.ROUND_DOWN : 1279 unscaled = Math.floor(nextAfter(unscaled, Double.NEGATIVE_INFINITY)); 1280 break; 1281 case BigDecimal.ROUND_FLOOR : 1282 if (sign == -1) { 1283 unscaled = Math.ceil(nextAfter(unscaled, Double.POSITIVE_INFINITY)); 1284 } else { 1285 unscaled = Math.floor(nextAfter(unscaled, Double.NEGATIVE_INFINITY)); 1286 } 1287 break; 1288 case BigDecimal.ROUND_HALF_DOWN : { 1289 unscaled = nextAfter(unscaled, Double.NEGATIVE_INFINITY); 1290 double fraction = unscaled - Math.floor(unscaled); 1291 if (fraction > 0.5) { 1292 unscaled = Math.ceil(unscaled); 1293 } else { 1294 unscaled = Math.floor(unscaled); 1295 } 1296 break; 1297 } 1298 case BigDecimal.ROUND_HALF_EVEN : { 1299 double fraction = unscaled - Math.floor(unscaled); 1300 if (fraction > 0.5) { 1301 unscaled = Math.ceil(unscaled); 1302 } else if (fraction < 0.5) { 1303 unscaled = Math.floor(unscaled); 1304 } else { 1305 // The following equality test is intentional and needed for rounding purposes 1306 if (Math.floor(unscaled) / 2.0 == Math.floor(Math 1307 .floor(unscaled) / 2.0)) { // even 1308 unscaled = Math.floor(unscaled); 1309 } else { // odd 1310 unscaled = Math.ceil(unscaled); 1311 } 1312 } 1313 break; 1314 } 1315 case BigDecimal.ROUND_HALF_UP : { 1316 unscaled = nextAfter(unscaled, Double.POSITIVE_INFINITY); 1317 double fraction = unscaled - Math.floor(unscaled); 1318 if (fraction >= 0.5) { 1319 unscaled = Math.ceil(unscaled); 1320 } else { 1321 unscaled = Math.floor(unscaled); 1322 } 1323 break; 1324 } 1325 case BigDecimal.ROUND_UNNECESSARY : 1326 if (unscaled != Math.floor(unscaled)) { 1327 throw new ArithmeticException("Inexact result from rounding"); 1328 } 1329 break; 1330 case BigDecimal.ROUND_UP : 1331 unscaled = Math.ceil(nextAfter(unscaled, Double.POSITIVE_INFINITY)); 1332 break; 1333 default : 1334 throw MathRuntimeException.createIllegalArgumentException( 1335 "invalid rounding method {0}, valid methods: {1} ({2}), {3} ({4})," + 1336 " {5} ({6}), {7} ({8}), {9} ({10}), {11} ({12}), {13} ({14}), {15} ({16})", 1337 roundingMethod, 1338 "ROUND_CEILING", BigDecimal.ROUND_CEILING, 1339 "ROUND_DOWN", BigDecimal.ROUND_DOWN, 1340 "ROUND_FLOOR", BigDecimal.ROUND_FLOOR, 1341 "ROUND_HALF_DOWN", BigDecimal.ROUND_HALF_DOWN, 1342 "ROUND_HALF_EVEN", BigDecimal.ROUND_HALF_EVEN, 1343 "ROUND_HALF_UP", BigDecimal.ROUND_HALF_UP, 1344 "ROUND_UNNECESSARY", BigDecimal.ROUND_UNNECESSARY, 1345 "ROUND_UP", BigDecimal.ROUND_UP); 1346 } 1347 return unscaled; 1348 } 1349 1350 /** 1351 * Returns the <a href="http://mathworld.wolfram.com/Sign.html"> sign</a> 1352 * for byte value <code>x</code>. 1353 * <p> 1354 * For a byte value x, this method returns (byte)(+1) if x > 0, (byte)(0) if 1355 * x = 0, and (byte)(-1) if x < 0.</p> 1356 * 1357 * @param x the value, a byte 1358 * @return (byte)(+1), (byte)(0), or (byte)(-1), depending on the sign of x 1359 */ 1360 public static byte sign(final byte x) { 1361 return (x == ZB) ? ZB : (x > ZB) ? PB : NB; 1362 } 1363 1364 /** 1365 * Returns the <a href="http://mathworld.wolfram.com/Sign.html"> sign</a> 1366 * for double precision <code>x</code>. 1367 * <p> 1368 * For a double value <code>x</code>, this method returns 1369 * <code>+1.0</code> if <code>x > 0</code>, <code>0.0</code> if 1370 * <code>x = 0.0</code>, and <code>-1.0</code> if <code>x < 0</code>. 1371 * Returns <code>NaN</code> if <code>x</code> is <code>NaN</code>.</p> 1372 * 1373 * @param x the value, a double 1374 * @return +1.0, 0.0, or -1.0, depending on the sign of x 1375 */ 1376 public static double sign(final double x) { 1377 if (Double.isNaN(x)) { 1378 return Double.NaN; 1379 } 1380 return (x == 0.0) ? 0.0 : (x > 0.0) ? 1.0 : -1.0; 1381 } 1382 1383 /** 1384 * Returns the <a href="http://mathworld.wolfram.com/Sign.html"> sign</a> 1385 * for float value <code>x</code>. 1386 * <p> 1387 * For a float value x, this method returns +1.0F if x > 0, 0.0F if x = 1388 * 0.0F, and -1.0F if x < 0. Returns <code>NaN</code> if <code>x</code> 1389 * is <code>NaN</code>.</p> 1390 * 1391 * @param x the value, a float 1392 * @return +1.0F, 0.0F, or -1.0F, depending on the sign of x 1393 */ 1394 public static float sign(final float x) { 1395 if (Float.isNaN(x)) { 1396 return Float.NaN; 1397 } 1398 return (x == 0.0F) ? 0.0F : (x > 0.0F) ? 1.0F : -1.0F; 1399 } 1400 1401 /** 1402 * Returns the <a href="http://mathworld.wolfram.com/Sign.html"> sign</a> 1403 * for int value <code>x</code>. 1404 * <p> 1405 * For an int value x, this method returns +1 if x > 0, 0 if x = 0, and -1 1406 * if x < 0.</p> 1407 * 1408 * @param x the value, an int 1409 * @return +1, 0, or -1, depending on the sign of x 1410 */ 1411 public static int sign(final int x) { 1412 return (x == 0) ? 0 : (x > 0) ? 1 : -1; 1413 } 1414 1415 /** 1416 * Returns the <a href="http://mathworld.wolfram.com/Sign.html"> sign</a> 1417 * for long value <code>x</code>. 1418 * <p> 1419 * For a long value x, this method returns +1L if x > 0, 0L if x = 0, and 1420 * -1L if x < 0.</p> 1421 * 1422 * @param x the value, a long 1423 * @return +1L, 0L, or -1L, depending on the sign of x 1424 */ 1425 public static long sign(final long x) { 1426 return (x == 0L) ? 0L : (x > 0L) ? 1L : -1L; 1427 } 1428 1429 /** 1430 * Returns the <a href="http://mathworld.wolfram.com/Sign.html"> sign</a> 1431 * for short value <code>x</code>. 1432 * <p> 1433 * For a short value x, this method returns (short)(+1) if x > 0, (short)(0) 1434 * if x = 0, and (short)(-1) if x < 0.</p> 1435 * 1436 * @param x the value, a short 1437 * @return (short)(+1), (short)(0), or (short)(-1), depending on the sign of 1438 * x 1439 */ 1440 public static short sign(final short x) { 1441 return (x == ZS) ? ZS : (x > ZS) ? PS : NS; 1442 } 1443 1444 /** 1445 * Returns the <a href="http://mathworld.wolfram.com/HyperbolicSine.html"> 1446 * hyperbolic sine</a> of x. 1447 * 1448 * @param x double value for which to find the hyperbolic sine 1449 * @return hyperbolic sine of x 1450 */ 1451 public static double sinh(double x) { 1452 return (Math.exp(x) - Math.exp(-x)) / 2.0; 1453 } 1454 1455 /** 1456 * Subtract two integers, checking for overflow. 1457 * 1458 * @param x the minuend 1459 * @param y the subtrahend 1460 * @return the difference <code>x-y</code> 1461 * @throws ArithmeticException if the result can not be represented as an 1462 * int 1463 * @since 1.1 1464 */ 1465 public static int subAndCheck(int x, int y) { 1466 long s = (long)x - (long)y; 1467 if (s < Integer.MIN_VALUE || s > Integer.MAX_VALUE) { 1468 throw new ArithmeticException("overflow: subtract"); 1469 } 1470 return (int)s; 1471 } 1472 1473 /** 1474 * Subtract two long integers, checking for overflow. 1475 * 1476 * @param a first value 1477 * @param b second value 1478 * @return the difference <code>a-b</code> 1479 * @throws ArithmeticException if the result can not be represented as an 1480 * long 1481 * @since 1.2 1482 */ 1483 public static long subAndCheck(long a, long b) { 1484 long ret; 1485 String msg = "overflow: subtract"; 1486 if (b == Long.MIN_VALUE) { 1487 if (a < 0) { 1488 ret = a - b; 1489 } else { 1490 throw new ArithmeticException(msg); 1491 } 1492 } else { 1493 // use additive inverse 1494 ret = addAndCheck(a, -b, msg); 1495 } 1496 return ret; 1497 } 1498 1499 /** 1500 * Raise an int to an int power. 1501 * @param k number to raise 1502 * @param e exponent (must be positive or null) 1503 * @return k<sup>e</sup> 1504 * @exception IllegalArgumentException if e is negative 1505 */ 1506 public static int pow(final int k, int e) 1507 throws IllegalArgumentException { 1508 1509 if (e < 0) { 1510 throw MathRuntimeException.createIllegalArgumentException( 1511 "cannot raise an integral value to a negative power ({0}^{1})", 1512 k, e); 1513 } 1514 1515 int result = 1; 1516 int k2p = k; 1517 while (e != 0) { 1518 if ((e & 0x1) != 0) { 1519 result *= k2p; 1520 } 1521 k2p *= k2p; 1522 e = e >> 1; 1523 } 1524 1525 return result; 1526 1527 } 1528 1529 /** 1530 * Raise an int to a long power. 1531 * @param k number to raise 1532 * @param e exponent (must be positive or null) 1533 * @return k<sup>e</sup> 1534 * @exception IllegalArgumentException if e is negative 1535 */ 1536 public static int pow(final int k, long e) 1537 throws IllegalArgumentException { 1538 1539 if (e < 0) { 1540 throw MathRuntimeException.createIllegalArgumentException( 1541 "cannot raise an integral value to a negative power ({0}^{1})", 1542 k, e); 1543 } 1544 1545 int result = 1; 1546 int k2p = k; 1547 while (e != 0) { 1548 if ((e & 0x1) != 0) { 1549 result *= k2p; 1550 } 1551 k2p *= k2p; 1552 e = e >> 1; 1553 } 1554 1555 return result; 1556 1557 } 1558 1559 /** 1560 * Raise a long to an int power. 1561 * @param k number to raise 1562 * @param e exponent (must be positive or null) 1563 * @return k<sup>e</sup> 1564 * @exception IllegalArgumentException if e is negative 1565 */ 1566 public static long pow(final long k, int e) 1567 throws IllegalArgumentException { 1568 1569 if (e < 0) { 1570 throw MathRuntimeException.createIllegalArgumentException( 1571 "cannot raise an integral value to a negative power ({0}^{1})", 1572 k, e); 1573 } 1574 1575 long result = 1l; 1576 long k2p = k; 1577 while (e != 0) { 1578 if ((e & 0x1) != 0) { 1579 result *= k2p; 1580 } 1581 k2p *= k2p; 1582 e = e >> 1; 1583 } 1584 1585 return result; 1586 1587 } 1588 1589 /** 1590 * Raise a long to a long power. 1591 * @param k number to raise 1592 * @param e exponent (must be positive or null) 1593 * @return k<sup>e</sup> 1594 * @exception IllegalArgumentException if e is negative 1595 */ 1596 public static long pow(final long k, long e) 1597 throws IllegalArgumentException { 1598 1599 if (e < 0) { 1600 throw MathRuntimeException.createIllegalArgumentException( 1601 "cannot raise an integral value to a negative power ({0}^{1})", 1602 k, e); 1603 } 1604 1605 long result = 1l; 1606 long k2p = k; 1607 while (e != 0) { 1608 if ((e & 0x1) != 0) { 1609 result *= k2p; 1610 } 1611 k2p *= k2p; 1612 e = e >> 1; 1613 } 1614 1615 return result; 1616 1617 } 1618 1619 /** 1620 * Raise a BigInteger to an int power. 1621 * @param k number to raise 1622 * @param e exponent (must be positive or null) 1623 * @return k<sup>e</sup> 1624 * @exception IllegalArgumentException if e is negative 1625 */ 1626 public static BigInteger pow(final BigInteger k, int e) 1627 throws IllegalArgumentException { 1628 1629 if (e < 0) { 1630 throw MathRuntimeException.createIllegalArgumentException( 1631 "cannot raise an integral value to a negative power ({0}^{1})", 1632 k, e); 1633 } 1634 1635 return k.pow(e); 1636 1637 } 1638 1639 /** 1640 * Raise a BigInteger to a long power. 1641 * @param k number to raise 1642 * @param e exponent (must be positive or null) 1643 * @return k<sup>e</sup> 1644 * @exception IllegalArgumentException if e is negative 1645 */ 1646 public static BigInteger pow(final BigInteger k, long e) 1647 throws IllegalArgumentException { 1648 1649 if (e < 0) { 1650 throw MathRuntimeException.createIllegalArgumentException( 1651 "cannot raise an integral value to a negative power ({0}^{1})", 1652 k, e); 1653 } 1654 1655 BigInteger result = BigInteger.ONE; 1656 BigInteger k2p = k; 1657 while (e != 0) { 1658 if ((e & 0x1) != 0) { 1659 result = result.multiply(k2p); 1660 } 1661 k2p = k2p.multiply(k2p); 1662 e = e >> 1; 1663 } 1664 1665 return result; 1666 1667 } 1668 1669 /** 1670 * Raise a BigInteger to a BigInteger power. 1671 * @param k number to raise 1672 * @param e exponent (must be positive or null) 1673 * @return k<sup>e</sup> 1674 * @exception IllegalArgumentException if e is negative 1675 */ 1676 public static BigInteger pow(final BigInteger k, BigInteger e) 1677 throws IllegalArgumentException { 1678 1679 if (e.compareTo(BigInteger.ZERO) < 0) { 1680 throw MathRuntimeException.createIllegalArgumentException( 1681 "cannot raise an integral value to a negative power ({0}^{1})", 1682 k, e); 1683 } 1684 1685 BigInteger result = BigInteger.ONE; 1686 BigInteger k2p = k; 1687 while (!BigInteger.ZERO.equals(e)) { 1688 if (e.testBit(0)) { 1689 result = result.multiply(k2p); 1690 } 1691 k2p = k2p.multiply(k2p); 1692 e = e.shiftRight(1); 1693 } 1694 1695 return result; 1696 1697 } 1698 1699 /** 1700 * Calculates the L<sub>1</sub> (sum of abs) distance between two points. 1701 * 1702 * @param p1 the first point 1703 * @param p2 the second point 1704 * @return the L<sub>1</sub> distance between the two points 1705 */ 1706 public static double distance1(double[] p1, double[] p2) { 1707 double sum = 0; 1708 for (int i = 0; i < p1.length; i++) { 1709 sum += Math.abs(p1[i] - p2[i]); 1710 } 1711 return sum; 1712 } 1713 1714 /** 1715 * Calculates the L<sub>1</sub> (sum of abs) distance between two points. 1716 * 1717 * @param p1 the first point 1718 * @param p2 the second point 1719 * @return the L<sub>1</sub> distance between the two points 1720 */ 1721 public static int distance1(int[] p1, int[] p2) { 1722 int sum = 0; 1723 for (int i = 0; i < p1.length; i++) { 1724 sum += Math.abs(p1[i] - p2[i]); 1725 } 1726 return sum; 1727 } 1728 1729 /** 1730 * Calculates the L<sub>2</sub> (Euclidean) distance between two points. 1731 * 1732 * @param p1 the first point 1733 * @param p2 the second point 1734 * @return the L<sub>2</sub> distance between the two points 1735 */ 1736 public static double distance(double[] p1, double[] p2) { 1737 double sum = 0; 1738 for (int i = 0; i < p1.length; i++) { 1739 final double dp = p1[i] - p2[i]; 1740 sum += dp * dp; 1741 } 1742 return Math.sqrt(sum); 1743 } 1744 1745 /** 1746 * Calculates the L<sub>2</sub> (Euclidean) distance between two points. 1747 * 1748 * @param p1 the first point 1749 * @param p2 the second point 1750 * @return the L<sub>2</sub> distance between the two points 1751 */ 1752 public static double distance(int[] p1, int[] p2) { 1753 double sum = 0; 1754 for (int i = 0; i < p1.length; i++) { 1755 final double dp = p1[i] - p2[i]; 1756 sum += dp * dp; 1757 } 1758 return Math.sqrt(sum); 1759 } 1760 1761 /** 1762 * Calculates the L<sub>∞</sub> (max of abs) distance between two points. 1763 * 1764 * @param p1 the first point 1765 * @param p2 the second point 1766 * @return the L<sub>∞</sub> distance between the two points 1767 */ 1768 public static double distanceInf(double[] p1, double[] p2) { 1769 double max = 0; 1770 for (int i = 0; i < p1.length; i++) { 1771 max = Math.max(max, Math.abs(p1[i] - p2[i])); 1772 } 1773 return max; 1774 } 1775 1776 /** 1777 * Calculates the L<sub>∞</sub> (max of abs) distance between two points. 1778 * 1779 * @param p1 the first point 1780 * @param p2 the second point 1781 * @return the L<sub>∞</sub> distance between the two points 1782 */ 1783 public static int distanceInf(int[] p1, int[] p2) { 1784 int max = 0; 1785 for (int i = 0; i < p1.length; i++) { 1786 max = Math.max(max, Math.abs(p1[i] - p2[i])); 1787 } 1788 return max; 1789 } 1790 1791 /** 1792 * Checks that the given array is sorted. 1793 * 1794 * @param val Values 1795 * @param dir Order direction (-1 for decreasing, 1 for increasing) 1796 * @param strict Whether the order should be strict 1797 * @throws IllegalArgumentException if the array is not sorted. 1798 */ 1799 public static void checkOrder(double[] val, int dir, boolean strict) { 1800 double previous = val[0]; 1801 1802 int max = val.length; 1803 for (int i = 1; i < max; i++) { 1804 if (dir > 0) { 1805 if (strict) { 1806 if (val[i] <= previous) { 1807 throw MathRuntimeException.createIllegalArgumentException("points {0} and {1} are not strictly increasing ({2} >= {3})", 1808 i - 1, i, previous, val[i]); 1809 } 1810 } else { 1811 if (val[i] < previous) { 1812 throw MathRuntimeException.createIllegalArgumentException("points {0} and {1} are not increasing ({2} > {3})", 1813 i - 1, i, previous, val[i]); 1814 } 1815 } 1816 } else { 1817 if (strict) { 1818 if (val[i] >= previous) { 1819 throw MathRuntimeException.createIllegalArgumentException("points {0} and {1} are not strictly decreasing ({2} <= {3})", 1820 i - 1, i, previous, val[i]); 1821 } 1822 } else { 1823 if (val[i] > previous) { 1824 throw MathRuntimeException.createIllegalArgumentException("points {0} and {1} are not decreasing ({2} < {3})", 1825 i - 1, i, previous, val[i]); 1826 } 1827 } 1828 } 1829 1830 previous = val[i]; 1831 } 1832 } 1833 }