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