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 package org.apache.commons.math3.util; 018 019 import java.math.BigInteger; 020 import java.util.concurrent.atomic.AtomicReference; 021 022 import org.apache.commons.math3.exception.MathArithmeticException; 023 import org.apache.commons.math3.exception.NotPositiveException; 024 import org.apache.commons.math3.exception.NumberIsTooLargeException; 025 import org.apache.commons.math3.exception.util.Localizable; 026 import org.apache.commons.math3.exception.util.LocalizedFormats; 027 028 /** 029 * Some useful, arithmetics related, additions to the built-in functions in 030 * {@link Math}. 031 * 032 * @version $Id: ArithmeticUtils.java 1422313 2012-12-15 18:53:41Z psteitz $ 033 */ 034 public final class ArithmeticUtils { 035 036 /** All long-representable factorials */ 037 static final long[] FACTORIALS = new long[] { 038 1l, 1l, 2l, 039 6l, 24l, 120l, 040 720l, 5040l, 40320l, 041 362880l, 3628800l, 39916800l, 042 479001600l, 6227020800l, 87178291200l, 043 1307674368000l, 20922789888000l, 355687428096000l, 044 6402373705728000l, 121645100408832000l, 2432902008176640000l }; 045 046 /** Stirling numbers of the second kind. */ 047 static final AtomicReference<long[][]> STIRLING_S2 = new AtomicReference<long[][]> (null); 048 049 /** Private constructor. */ 050 private ArithmeticUtils() { 051 super(); 052 } 053 054 /** 055 * Add two integers, checking for overflow. 056 * 057 * @param x an addend 058 * @param y an addend 059 * @return the sum {@code x+y} 060 * @throws MathArithmeticException if the result can not be represented 061 * as an {@code int}. 062 * @since 1.1 063 */ 064 public static int addAndCheck(int x, int y) 065 throws MathArithmeticException { 066 long s = (long)x + (long)y; 067 if (s < Integer.MIN_VALUE || s > Integer.MAX_VALUE) { 068 throw new MathArithmeticException(LocalizedFormats.OVERFLOW_IN_ADDITION, x, y); 069 } 070 return (int)s; 071 } 072 073 /** 074 * Add two long integers, checking for overflow. 075 * 076 * @param a an addend 077 * @param b an addend 078 * @return the sum {@code a+b} 079 * @throws MathArithmeticException if the result can not be represented as an 080 * long 081 * @since 1.2 082 */ 083 public static long addAndCheck(long a, long b) throws MathArithmeticException { 084 return ArithmeticUtils.addAndCheck(a, b, LocalizedFormats.OVERFLOW_IN_ADDITION); 085 } 086 087 /** 088 * Returns an exact representation of the <a 089 * href="http://mathworld.wolfram.com/BinomialCoefficient.html"> Binomial 090 * Coefficient</a>, "{@code n choose k}", the number of 091 * {@code k}-element subsets that can be selected from an 092 * {@code n}-element set. 093 * <p> 094 * <Strong>Preconditions</strong>: 095 * <ul> 096 * <li> {@code 0 <= k <= n } (otherwise 097 * {@code IllegalArgumentException} is thrown)</li> 098 * <li> The result is small enough to fit into a {@code long}. The 099 * largest value of {@code n} for which all coefficients are 100 * {@code < Long.MAX_VALUE} is 66. If the computed value exceeds 101 * {@code Long.MAX_VALUE} an {@code ArithMeticException} is 102 * thrown.</li> 103 * </ul></p> 104 * 105 * @param n the size of the set 106 * @param k the size of the subsets to be counted 107 * @return {@code n choose k} 108 * @throws NotPositiveException if {@code n < 0}. 109 * @throws NumberIsTooLargeException if {@code k > n}. 110 * @throws MathArithmeticException if the result is too large to be 111 * represented by a long integer. 112 */ 113 public static long binomialCoefficient(final int n, final int k) 114 throws NotPositiveException, NumberIsTooLargeException, MathArithmeticException { 115 ArithmeticUtils.checkBinomial(n, k); 116 if ((n == k) || (k == 0)) { 117 return 1; 118 } 119 if ((k == 1) || (k == n - 1)) { 120 return n; 121 } 122 // Use symmetry for large k 123 if (k > n / 2) { 124 return binomialCoefficient(n, n - k); 125 } 126 127 // We use the formula 128 // (n choose k) = n! / (n-k)! / k! 129 // (n choose k) == ((n-k+1)*...*n) / (1*...*k) 130 // which could be written 131 // (n choose k) == (n-1 choose k-1) * n / k 132 long result = 1; 133 if (n <= 61) { 134 // For n <= 61, the naive implementation cannot overflow. 135 int i = n - k + 1; 136 for (int j = 1; j <= k; j++) { 137 result = result * i / j; 138 i++; 139 } 140 } else if (n <= 66) { 141 // For n > 61 but n <= 66, the result cannot overflow, 142 // but we must take care not to overflow intermediate values. 143 int i = n - k + 1; 144 for (int j = 1; j <= k; j++) { 145 // We know that (result * i) is divisible by j, 146 // but (result * i) may overflow, so we split j: 147 // Filter out the gcd, d, so j/d and i/d are integer. 148 // result is divisible by (j/d) because (j/d) 149 // is relative prime to (i/d) and is a divisor of 150 // result * (i/d). 151 final long d = gcd(i, j); 152 result = (result / (j / d)) * (i / d); 153 i++; 154 } 155 } else { 156 // For n > 66, a result overflow might occur, so we check 157 // the multiplication, taking care to not overflow 158 // unnecessary. 159 int i = n - k + 1; 160 for (int j = 1; j <= k; j++) { 161 final long d = gcd(i, j); 162 result = mulAndCheck(result / (j / d), i / d); 163 i++; 164 } 165 } 166 return result; 167 } 168 169 /** 170 * Returns a {@code double} representation of the <a 171 * href="http://mathworld.wolfram.com/BinomialCoefficient.html"> Binomial 172 * Coefficient</a>, "{@code n choose k}", the number of 173 * {@code k}-element subsets that can be selected from an 174 * {@code n}-element set. 175 * <p> 176 * <Strong>Preconditions</strong>: 177 * <ul> 178 * <li> {@code 0 <= k <= n } (otherwise 179 * {@code IllegalArgumentException} is thrown)</li> 180 * <li> The result is small enough to fit into a {@code double}. The 181 * largest value of {@code n} for which all coefficients are < 182 * Double.MAX_VALUE is 1029. If the computed value exceeds Double.MAX_VALUE, 183 * Double.POSITIVE_INFINITY is returned</li> 184 * </ul></p> 185 * 186 * @param n the size of the set 187 * @param k the size of the subsets to be counted 188 * @return {@code n choose k} 189 * @throws NotPositiveException if {@code n < 0}. 190 * @throws NumberIsTooLargeException if {@code k > n}. 191 * @throws MathArithmeticException if the result is too large to be 192 * represented by a long integer. 193 */ 194 public static double binomialCoefficientDouble(final int n, final int k) 195 throws NotPositiveException, NumberIsTooLargeException, MathArithmeticException { 196 ArithmeticUtils.checkBinomial(n, k); 197 if ((n == k) || (k == 0)) { 198 return 1d; 199 } 200 if ((k == 1) || (k == n - 1)) { 201 return n; 202 } 203 if (k > n/2) { 204 return binomialCoefficientDouble(n, n - k); 205 } 206 if (n < 67) { 207 return binomialCoefficient(n,k); 208 } 209 210 double result = 1d; 211 for (int i = 1; i <= k; i++) { 212 result *= (double)(n - k + i) / (double)i; 213 } 214 215 return FastMath.floor(result + 0.5); 216 } 217 218 /** 219 * Returns the natural {@code log} of the <a 220 * href="http://mathworld.wolfram.com/BinomialCoefficient.html"> Binomial 221 * Coefficient</a>, "{@code n choose k}", the number of 222 * {@code k}-element subsets that can be selected from an 223 * {@code n}-element set. 224 * <p> 225 * <Strong>Preconditions</strong>: 226 * <ul> 227 * <li> {@code 0 <= k <= n } (otherwise 228 * {@code IllegalArgumentException} is thrown)</li> 229 * </ul></p> 230 * 231 * @param n the size of the set 232 * @param k the size of the subsets to be counted 233 * @return {@code n choose k} 234 * @throws NotPositiveException if {@code n < 0}. 235 * @throws NumberIsTooLargeException if {@code k > n}. 236 * @throws MathArithmeticException if the result is too large to be 237 * represented by a long integer. 238 */ 239 public static double binomialCoefficientLog(final int n, final int k) 240 throws NotPositiveException, NumberIsTooLargeException, MathArithmeticException { 241 ArithmeticUtils.checkBinomial(n, k); 242 if ((n == k) || (k == 0)) { 243 return 0; 244 } 245 if ((k == 1) || (k == n - 1)) { 246 return FastMath.log(n); 247 } 248 249 /* 250 * For values small enough to do exact integer computation, 251 * return the log of the exact value 252 */ 253 if (n < 67) { 254 return FastMath.log(binomialCoefficient(n,k)); 255 } 256 257 /* 258 * Return the log of binomialCoefficientDouble for values that will not 259 * overflow binomialCoefficientDouble 260 */ 261 if (n < 1030) { 262 return FastMath.log(binomialCoefficientDouble(n, k)); 263 } 264 265 if (k > n / 2) { 266 return binomialCoefficientLog(n, n - k); 267 } 268 269 /* 270 * Sum logs for values that could overflow 271 */ 272 double logSum = 0; 273 274 // n!/(n-k)! 275 for (int i = n - k + 1; i <= n; i++) { 276 logSum += FastMath.log(i); 277 } 278 279 // divide by k! 280 for (int i = 2; i <= k; i++) { 281 logSum -= FastMath.log(i); 282 } 283 284 return logSum; 285 } 286 287 /** 288 * Returns n!. Shorthand for {@code n} <a 289 * href="http://mathworld.wolfram.com/Factorial.html"> Factorial</a>, the 290 * product of the numbers {@code 1,...,n}. 291 * <p> 292 * <Strong>Preconditions</strong>: 293 * <ul> 294 * <li> {@code n >= 0} (otherwise 295 * {@code IllegalArgumentException} is thrown)</li> 296 * <li> The result is small enough to fit into a {@code long}. The 297 * largest value of {@code n} for which {@code n!} < 298 * Long.MAX_VALUE} is 20. If the computed value exceeds {@code Long.MAX_VALUE} 299 * an {@code ArithMeticException } is thrown.</li> 300 * </ul> 301 * </p> 302 * 303 * @param n argument 304 * @return {@code n!} 305 * @throws MathArithmeticException if the result is too large to be represented 306 * by a {@code long}. 307 * @throws NotPositiveException if {@code n < 0}. 308 * @throws MathArithmeticException if {@code n > 20}: The factorial value is too 309 * large to fit in a {@code long}. 310 */ 311 public static long factorial(final int n) throws NotPositiveException, MathArithmeticException { 312 if (n < 0) { 313 throw new NotPositiveException(LocalizedFormats.FACTORIAL_NEGATIVE_PARAMETER, 314 n); 315 } 316 if (n > 20) { 317 throw new MathArithmeticException(); 318 } 319 return FACTORIALS[n]; 320 } 321 322 /** 323 * Compute n!, the<a href="http://mathworld.wolfram.com/Factorial.html"> 324 * factorial</a> of {@code n} (the product of the numbers 1 to n), as a 325 * {@code double}. 326 * The result should be small enough to fit into a {@code double}: The 327 * largest {@code n} for which {@code n! < Double.MAX_VALUE} is 170. 328 * If the computed value exceeds {@code Double.MAX_VALUE}, 329 * {@code Double.POSITIVE_INFINITY} is returned. 330 * 331 * @param n Argument. 332 * @return {@code n!} 333 * @throws NotPositiveException if {@code n < 0}. 334 */ 335 public static double factorialDouble(final int n) throws NotPositiveException { 336 if (n < 0) { 337 throw new NotPositiveException(LocalizedFormats.FACTORIAL_NEGATIVE_PARAMETER, 338 n); 339 } 340 if (n < 21) { 341 return FACTORIALS[n]; 342 } 343 return FastMath.floor(FastMath.exp(ArithmeticUtils.factorialLog(n)) + 0.5); 344 } 345 346 /** 347 * Compute the natural logarithm of the factorial of {@code n}. 348 * 349 * @param n Argument. 350 * @return {@code n!} 351 * @throws NotPositiveException if {@code n < 0}. 352 */ 353 public static double factorialLog(final int n) throws NotPositiveException { 354 if (n < 0) { 355 throw new NotPositiveException(LocalizedFormats.FACTORIAL_NEGATIVE_PARAMETER, 356 n); 357 } 358 if (n < 21) { 359 return FastMath.log(FACTORIALS[n]); 360 } 361 double logSum = 0; 362 for (int i = 2; i <= n; i++) { 363 logSum += FastMath.log(i); 364 } 365 return logSum; 366 } 367 368 /** 369 * Computes the greatest common divisor of the absolute value of two 370 * numbers, using a modified version of the "binary gcd" method. 371 * See Knuth 4.5.2 algorithm B. 372 * The algorithm is due to Josef Stein (1961). 373 * <br/> 374 * Special cases: 375 * <ul> 376 * <li>The invocations 377 * {@code gcd(Integer.MIN_VALUE, Integer.MIN_VALUE)}, 378 * {@code gcd(Integer.MIN_VALUE, 0)} and 379 * {@code gcd(0, Integer.MIN_VALUE)} throw an 380 * {@code ArithmeticException}, because the result would be 2^31, which 381 * is too large for an int value.</li> 382 * <li>The result of {@code gcd(x, x)}, {@code gcd(0, x)} and 383 * {@code gcd(x, 0)} is the absolute value of {@code x}, except 384 * for the special cases above.</li> 385 * <li>The invocation {@code gcd(0, 0)} is the only one which returns 386 * {@code 0}.</li> 387 * </ul> 388 * 389 * @param p Number. 390 * @param q Number. 391 * @return the greatest common divisor (never negative). 392 * @throws MathArithmeticException if the result cannot be represented as 393 * a non-negative {@code int} value. 394 * @since 1.1 395 */ 396 public static int gcd(int p, 397 int q) 398 throws MathArithmeticException { 399 int a = p; 400 int b = q; 401 if (a == 0 || 402 b == 0) { 403 if (a == Integer.MIN_VALUE || 404 b == Integer.MIN_VALUE) { 405 throw new MathArithmeticException(LocalizedFormats.GCD_OVERFLOW_32_BITS, 406 p, q); 407 } 408 return FastMath.abs(a + b); 409 } 410 411 long al = a; 412 long bl = b; 413 boolean useLong = false; 414 if (a < 0) { 415 if(Integer.MIN_VALUE == a) { 416 useLong = true; 417 } else { 418 a = -a; 419 } 420 al = -al; 421 } 422 if (b < 0) { 423 if (Integer.MIN_VALUE == b) { 424 useLong = true; 425 } else { 426 b = -b; 427 } 428 bl = -bl; 429 } 430 if (useLong) { 431 if(al == bl) { 432 throw new MathArithmeticException(LocalizedFormats.GCD_OVERFLOW_32_BITS, 433 p, q); 434 } 435 long blbu = bl; 436 bl = al; 437 al = blbu % al; 438 if (al == 0) { 439 if (bl > Integer.MAX_VALUE) { 440 throw new MathArithmeticException(LocalizedFormats.GCD_OVERFLOW_32_BITS, 441 p, q); 442 } 443 return (int) bl; 444 } 445 blbu = bl; 446 447 // Now "al" and "bl" fit in an "int". 448 b = (int) al; 449 a = (int) (blbu % al); 450 } 451 452 return gcdPositive(a, b); 453 } 454 455 /** 456 * Computes the greatest common divisor of two <em>positive</em> numbers 457 * (this precondition is <em>not</em> checked and the result is undefined 458 * if not fulfilled) using the "binary gcd" method which avoids division 459 * and modulo operations. 460 * See Knuth 4.5.2 algorithm B. 461 * The algorithm is due to Josef Stein (1961). 462 * <br/> 463 * Special cases: 464 * <ul> 465 * <li>The result of {@code gcd(x, x)}, {@code gcd(0, x)} and 466 * {@code gcd(x, 0)} is the value of {@code x}.</li> 467 * <li>The invocation {@code gcd(0, 0)} is the only one which returns 468 * {@code 0}.</li> 469 * </ul> 470 * 471 * @param a Positive number. 472 * @param b Positive number. 473 * @return the greatest common divisor. 474 */ 475 private static int gcdPositive(int a, 476 int b) { 477 if (a == 0) { 478 return b; 479 } 480 else if (b == 0) { 481 return a; 482 } 483 484 // Make "a" and "b" odd, keeping track of common power of 2. 485 final int aTwos = Integer.numberOfTrailingZeros(a); 486 a >>= aTwos; 487 final int bTwos = Integer.numberOfTrailingZeros(b); 488 b >>= bTwos; 489 final int shift = Math.min(aTwos, bTwos); 490 491 // "a" and "b" are positive. 492 // If a > b then "gdc(a, b)" is equal to "gcd(a - b, b)". 493 // If a < b then "gcd(a, b)" is equal to "gcd(b - a, a)". 494 // Hence, in the successive iterations: 495 // "a" becomes the absolute difference of the current values, 496 // "b" becomes the minimum of the current values. 497 while (a != b) { 498 final int delta = a - b; 499 b = Math.min(a, b); 500 a = Math.abs(delta); 501 502 // Remove any power of 2 in "a" ("b" is guaranteed to be odd). 503 a >>= Integer.numberOfTrailingZeros(a); 504 } 505 506 // Recover the common power of 2. 507 return a << shift; 508 } 509 510 /** 511 * <p> 512 * Gets the greatest common divisor of the absolute value of two numbers, 513 * using the "binary gcd" method which avoids division and modulo 514 * operations. See Knuth 4.5.2 algorithm B. This algorithm is due to Josef 515 * Stein (1961). 516 * </p> 517 * Special cases: 518 * <ul> 519 * <li>The invocations 520 * {@code gcd(Long.MIN_VALUE, Long.MIN_VALUE)}, 521 * {@code gcd(Long.MIN_VALUE, 0L)} and 522 * {@code gcd(0L, Long.MIN_VALUE)} throw an 523 * {@code ArithmeticException}, because the result would be 2^63, which 524 * is too large for a long value.</li> 525 * <li>The result of {@code gcd(x, x)}, {@code gcd(0L, x)} and 526 * {@code gcd(x, 0L)} is the absolute value of {@code x}, except 527 * for the special cases above. 528 * <li>The invocation {@code gcd(0L, 0L)} is the only one which returns 529 * {@code 0L}.</li> 530 * </ul> 531 * 532 * @param p Number. 533 * @param q Number. 534 * @return the greatest common divisor, never negative. 535 * @throws MathArithmeticException if the result cannot be represented as 536 * a non-negative {@code long} value. 537 * @since 2.1 538 */ 539 public static long gcd(final long p, final long q) throws MathArithmeticException { 540 long u = p; 541 long v = q; 542 if ((u == 0) || (v == 0)) { 543 if ((u == Long.MIN_VALUE) || (v == Long.MIN_VALUE)){ 544 throw new MathArithmeticException(LocalizedFormats.GCD_OVERFLOW_64_BITS, 545 p, q); 546 } 547 return FastMath.abs(u) + FastMath.abs(v); 548 } 549 // keep u and v negative, as negative integers range down to 550 // -2^63, while positive numbers can only be as large as 2^63-1 551 // (i.e. we can't necessarily negate a negative number without 552 // overflow) 553 /* assert u!=0 && v!=0; */ 554 if (u > 0) { 555 u = -u; 556 } // make u negative 557 if (v > 0) { 558 v = -v; 559 } // make v negative 560 // B1. [Find power of 2] 561 int k = 0; 562 while ((u & 1) == 0 && (v & 1) == 0 && k < 63) { // while u and v are 563 // both even... 564 u /= 2; 565 v /= 2; 566 k++; // cast out twos. 567 } 568 if (k == 63) { 569 throw new MathArithmeticException(LocalizedFormats.GCD_OVERFLOW_64_BITS, 570 p, q); 571 } 572 // B2. Initialize: u and v have been divided by 2^k and at least 573 // one is odd. 574 long t = ((u & 1) == 1) ? v : -(u / 2)/* B3 */; 575 // t negative: u was odd, v may be even (t replaces v) 576 // t positive: u was even, v is odd (t replaces u) 577 do { 578 /* assert u<0 && v<0; */ 579 // B4/B3: cast out twos from t. 580 while ((t & 1) == 0) { // while t is even.. 581 t /= 2; // cast out twos 582 } 583 // B5 [reset max(u,v)] 584 if (t > 0) { 585 u = -t; 586 } else { 587 v = t; 588 } 589 // B6/B3. at this point both u and v should be odd. 590 t = (v - u) / 2; 591 // |u| larger: t positive (replace u) 592 // |v| larger: t negative (replace v) 593 } while (t != 0); 594 return -u * (1L << k); // gcd is u*2^k 595 } 596 597 /** 598 * <p> 599 * Returns the least common multiple of the absolute value of two numbers, 600 * using the formula {@code lcm(a,b) = (a / gcd(a,b)) * b}. 601 * </p> 602 * Special cases: 603 * <ul> 604 * <li>The invocations {@code lcm(Integer.MIN_VALUE, n)} and 605 * {@code lcm(n, Integer.MIN_VALUE)}, where {@code abs(n)} is a 606 * power of 2, throw an {@code ArithmeticException}, because the result 607 * would be 2^31, which is too large for an int value.</li> 608 * <li>The result of {@code lcm(0, x)} and {@code lcm(x, 0)} is 609 * {@code 0} for any {@code x}. 610 * </ul> 611 * 612 * @param a Number. 613 * @param b Number. 614 * @return the least common multiple, never negative. 615 * @throws MathArithmeticException if the result cannot be represented as 616 * a non-negative {@code int} value. 617 * @since 1.1 618 */ 619 public static int lcm(int a, int b) throws MathArithmeticException { 620 if (a == 0 || b == 0){ 621 return 0; 622 } 623 int lcm = FastMath.abs(ArithmeticUtils.mulAndCheck(a / gcd(a, b), b)); 624 if (lcm == Integer.MIN_VALUE) { 625 throw new MathArithmeticException(LocalizedFormats.LCM_OVERFLOW_32_BITS, 626 a, b); 627 } 628 return lcm; 629 } 630 631 /** 632 * <p> 633 * Returns the least common multiple of the absolute value of two numbers, 634 * using the formula {@code lcm(a,b) = (a / gcd(a,b)) * b}. 635 * </p> 636 * Special cases: 637 * <ul> 638 * <li>The invocations {@code lcm(Long.MIN_VALUE, n)} and 639 * {@code lcm(n, Long.MIN_VALUE)}, where {@code abs(n)} is a 640 * power of 2, throw an {@code ArithmeticException}, because the result 641 * would be 2^63, which is too large for an int value.</li> 642 * <li>The result of {@code lcm(0L, x)} and {@code lcm(x, 0L)} is 643 * {@code 0L} for any {@code x}. 644 * </ul> 645 * 646 * @param a Number. 647 * @param b Number. 648 * @return the least common multiple, never negative. 649 * @throws MathArithmeticException if the result cannot be represented 650 * as a non-negative {@code long} value. 651 * @since 2.1 652 */ 653 public static long lcm(long a, long b) throws MathArithmeticException { 654 if (a == 0 || b == 0){ 655 return 0; 656 } 657 long lcm = FastMath.abs(ArithmeticUtils.mulAndCheck(a / gcd(a, b), b)); 658 if (lcm == Long.MIN_VALUE){ 659 throw new MathArithmeticException(LocalizedFormats.LCM_OVERFLOW_64_BITS, 660 a, b); 661 } 662 return lcm; 663 } 664 665 /** 666 * Multiply two integers, checking for overflow. 667 * 668 * @param x Factor. 669 * @param y Factor. 670 * @return the product {@code x * y}. 671 * @throws MathArithmeticException if the result can not be 672 * represented as an {@code int}. 673 * @since 1.1 674 */ 675 public static int mulAndCheck(int x, int y) throws MathArithmeticException { 676 long m = ((long)x) * ((long)y); 677 if (m < Integer.MIN_VALUE || m > Integer.MAX_VALUE) { 678 throw new MathArithmeticException(); 679 } 680 return (int)m; 681 } 682 683 /** 684 * Multiply two long integers, checking for overflow. 685 * 686 * @param a Factor. 687 * @param b Factor. 688 * @return the product {@code a * b}. 689 * @throws MathArithmeticException if the result can not be represented 690 * as a {@code long}. 691 * @since 1.2 692 */ 693 public static long mulAndCheck(long a, long b) throws MathArithmeticException { 694 long ret; 695 if (a > b) { 696 // use symmetry to reduce boundary cases 697 ret = mulAndCheck(b, a); 698 } else { 699 if (a < 0) { 700 if (b < 0) { 701 // check for positive overflow with negative a, negative b 702 if (a >= Long.MAX_VALUE / b) { 703 ret = a * b; 704 } else { 705 throw new MathArithmeticException(); 706 } 707 } else if (b > 0) { 708 // check for negative overflow with negative a, positive b 709 if (Long.MIN_VALUE / b <= a) { 710 ret = a * b; 711 } else { 712 throw new MathArithmeticException(); 713 714 } 715 } else { 716 // assert b == 0 717 ret = 0; 718 } 719 } else if (a > 0) { 720 // assert a > 0 721 // assert b > 0 722 723 // check for positive overflow with positive a, positive b 724 if (a <= Long.MAX_VALUE / b) { 725 ret = a * b; 726 } else { 727 throw new MathArithmeticException(); 728 } 729 } else { 730 // assert a == 0 731 ret = 0; 732 } 733 } 734 return ret; 735 } 736 737 /** 738 * Subtract two integers, checking for overflow. 739 * 740 * @param x Minuend. 741 * @param y Subtrahend. 742 * @return the difference {@code x - y}. 743 * @throws MathArithmeticException if the result can not be represented 744 * as an {@code int}. 745 * @since 1.1 746 */ 747 public static int subAndCheck(int x, int y) throws MathArithmeticException { 748 long s = (long)x - (long)y; 749 if (s < Integer.MIN_VALUE || s > Integer.MAX_VALUE) { 750 throw new MathArithmeticException(LocalizedFormats.OVERFLOW_IN_SUBTRACTION, x, y); 751 } 752 return (int)s; 753 } 754 755 /** 756 * Subtract two long integers, checking for overflow. 757 * 758 * @param a Value. 759 * @param b Value. 760 * @return the difference {@code a - b}. 761 * @throws MathArithmeticException if the result can not be represented as a 762 * {@code long}. 763 * @since 1.2 764 */ 765 public static long subAndCheck(long a, long b) throws MathArithmeticException { 766 long ret; 767 if (b == Long.MIN_VALUE) { 768 if (a < 0) { 769 ret = a - b; 770 } else { 771 throw new MathArithmeticException(LocalizedFormats.OVERFLOW_IN_ADDITION, a, -b); 772 } 773 } else { 774 // use additive inverse 775 ret = addAndCheck(a, -b, LocalizedFormats.OVERFLOW_IN_ADDITION); 776 } 777 return ret; 778 } 779 780 /** 781 * Raise an int to an int power. 782 * 783 * @param k Number to raise. 784 * @param e Exponent (must be positive or zero). 785 * @return k<sup>e</sup> 786 * @throws NotPositiveException if {@code e < 0}. 787 */ 788 public static int pow(final int k, int e) throws NotPositiveException { 789 if (e < 0) { 790 throw new NotPositiveException(LocalizedFormats.EXPONENT, e); 791 } 792 793 int result = 1; 794 int k2p = k; 795 while (e != 0) { 796 if ((e & 0x1) != 0) { 797 result *= k2p; 798 } 799 k2p *= k2p; 800 e = e >> 1; 801 } 802 803 return result; 804 } 805 806 /** 807 * Raise an int to a long power. 808 * 809 * @param k Number to raise. 810 * @param e Exponent (must be positive or zero). 811 * @return k<sup>e</sup> 812 * @throws NotPositiveException if {@code e < 0}. 813 */ 814 public static int pow(final int k, long e) throws NotPositiveException { 815 if (e < 0) { 816 throw new NotPositiveException(LocalizedFormats.EXPONENT, e); 817 } 818 819 int result = 1; 820 int k2p = k; 821 while (e != 0) { 822 if ((e & 0x1) != 0) { 823 result *= k2p; 824 } 825 k2p *= k2p; 826 e = e >> 1; 827 } 828 829 return result; 830 } 831 832 /** 833 * Raise a long to an int power. 834 * 835 * @param k Number to raise. 836 * @param e Exponent (must be positive or zero). 837 * @return k<sup>e</sup> 838 * @throws NotPositiveException if {@code e < 0}. 839 */ 840 public static long pow(final long k, int e) throws NotPositiveException { 841 if (e < 0) { 842 throw new NotPositiveException(LocalizedFormats.EXPONENT, e); 843 } 844 845 long result = 1l; 846 long k2p = k; 847 while (e != 0) { 848 if ((e & 0x1) != 0) { 849 result *= k2p; 850 } 851 k2p *= k2p; 852 e = e >> 1; 853 } 854 855 return result; 856 } 857 858 /** 859 * Raise a long to a long power. 860 * 861 * @param k Number to raise. 862 * @param e Exponent (must be positive or zero). 863 * @return k<sup>e</sup> 864 * @throws NotPositiveException if {@code e < 0}. 865 */ 866 public static long pow(final long k, long e) throws NotPositiveException { 867 if (e < 0) { 868 throw new NotPositiveException(LocalizedFormats.EXPONENT, e); 869 } 870 871 long result = 1l; 872 long k2p = k; 873 while (e != 0) { 874 if ((e & 0x1) != 0) { 875 result *= k2p; 876 } 877 k2p *= k2p; 878 e = e >> 1; 879 } 880 881 return result; 882 } 883 884 /** 885 * Raise a BigInteger to an int power. 886 * 887 * @param k Number to raise. 888 * @param e Exponent (must be positive or zero). 889 * @return k<sup>e</sup> 890 * @throws NotPositiveException if {@code e < 0}. 891 */ 892 public static BigInteger pow(final BigInteger k, int e) throws NotPositiveException { 893 if (e < 0) { 894 throw new NotPositiveException(LocalizedFormats.EXPONENT, e); 895 } 896 897 return k.pow(e); 898 } 899 900 /** 901 * Raise a BigInteger to a long power. 902 * 903 * @param k Number to raise. 904 * @param e Exponent (must be positive or zero). 905 * @return k<sup>e</sup> 906 * @throws NotPositiveException if {@code e < 0}. 907 */ 908 public static BigInteger pow(final BigInteger k, long e) throws NotPositiveException { 909 if (e < 0) { 910 throw new NotPositiveException(LocalizedFormats.EXPONENT, e); 911 } 912 913 BigInteger result = BigInteger.ONE; 914 BigInteger k2p = k; 915 while (e != 0) { 916 if ((e & 0x1) != 0) { 917 result = result.multiply(k2p); 918 } 919 k2p = k2p.multiply(k2p); 920 e = e >> 1; 921 } 922 923 return result; 924 925 } 926 927 /** 928 * Raise a BigInteger to a BigInteger power. 929 * 930 * @param k Number to raise. 931 * @param e Exponent (must be positive or zero). 932 * @return k<sup>e</sup> 933 * @throws NotPositiveException if {@code e < 0}. 934 */ 935 public static BigInteger pow(final BigInteger k, BigInteger e) throws NotPositiveException { 936 if (e.compareTo(BigInteger.ZERO) < 0) { 937 throw new NotPositiveException(LocalizedFormats.EXPONENT, e); 938 } 939 940 BigInteger result = BigInteger.ONE; 941 BigInteger k2p = k; 942 while (!BigInteger.ZERO.equals(e)) { 943 if (e.testBit(0)) { 944 result = result.multiply(k2p); 945 } 946 k2p = k2p.multiply(k2p); 947 e = e.shiftRight(1); 948 } 949 950 return result; 951 } 952 953 /** 954 * Returns the <a 955 * href="http://mathworld.wolfram.com/StirlingNumberoftheSecondKind.html"> 956 * Stirling number of the second kind</a>, "{@code S(n,k)}", the number of 957 * ways of partitioning an {@code n}-element set into {@code k} non-empty 958 * subsets. 959 * <p> 960 * The preconditions are {@code 0 <= k <= n } (otherwise 961 * {@code NotPositiveException} is thrown) 962 * </p> 963 * @param n the size of the set 964 * @param k the number of non-empty subsets 965 * @return {@code S(n,k)} 966 * @throws NotPositiveException if {@code k < 0}. 967 * @throws NumberIsTooLargeException if {@code k > n}. 968 * @throws MathArithmeticException if some overflow happens, typically for n exceeding 25 and 969 * k between 20 and n-2 (S(n,n-1) is handled specifically and does not overflow) 970 * @since 3.1 971 */ 972 public static long stirlingS2(final int n, final int k) 973 throws NotPositiveException, NumberIsTooLargeException, MathArithmeticException { 974 if (k < 0) { 975 throw new NotPositiveException(k); 976 } 977 if (k > n) { 978 throw new NumberIsTooLargeException(k, n, true); 979 } 980 981 long[][] stirlingS2 = STIRLING_S2.get(); 982 983 if (stirlingS2 == null) { 984 // the cache has never been initialized, compute the first numbers 985 // by direct recurrence relation 986 987 // as S(26,9) = 11201516780955125625 is larger than Long.MAX_VALUE 988 // we must stop computation at row 26 989 final int maxIndex = 26; 990 stirlingS2 = new long[maxIndex][]; 991 stirlingS2[0] = new long[] { 1l }; 992 for (int i = 1; i < stirlingS2.length; ++i) { 993 stirlingS2[i] = new long[i + 1]; 994 stirlingS2[i][0] = 0; 995 stirlingS2[i][1] = 1; 996 stirlingS2[i][i] = 1; 997 for (int j = 2; j < i; ++j) { 998 stirlingS2[i][j] = j * stirlingS2[i - 1][j] + stirlingS2[i - 1][j - 1]; 999 } 1000 } 1001 1002 // atomically save the cache 1003 STIRLING_S2.compareAndSet(null, stirlingS2); 1004 1005 } 1006 1007 if (n < stirlingS2.length) { 1008 // the number is in the small cache 1009 return stirlingS2[n][k]; 1010 } else { 1011 // use explicit formula to compute the number without caching it 1012 if (k == 0) { 1013 return 0; 1014 } else if (k == 1 || k == n) { 1015 return 1; 1016 } else if (k == 2) { 1017 return (1l << (n - 1)) - 1l; 1018 } else if (k == n - 1) { 1019 return binomialCoefficient(n, 2); 1020 } else { 1021 // definition formula: note that this may trigger some overflow 1022 long sum = 0; 1023 long sign = ((k & 0x1) == 0) ? 1 : -1; 1024 for (int j = 1; j <= k; ++j) { 1025 sign = -sign; 1026 sum += sign * binomialCoefficient(k, j) * pow(j, n); 1027 if (sum < 0) { 1028 // there was an overflow somewhere 1029 throw new MathArithmeticException(LocalizedFormats.ARGUMENT_OUTSIDE_DOMAIN, 1030 n, 0, stirlingS2.length - 1); 1031 } 1032 } 1033 return sum / factorial(k); 1034 } 1035 } 1036 1037 } 1038 1039 /** 1040 * Add two long integers, checking for overflow. 1041 * 1042 * @param a Addend. 1043 * @param b Addend. 1044 * @param pattern Pattern to use for any thrown exception. 1045 * @return the sum {@code a + b}. 1046 * @throws MathArithmeticException if the result cannot be represented 1047 * as a {@code long}. 1048 * @since 1.2 1049 */ 1050 private static long addAndCheck(long a, long b, Localizable pattern) throws MathArithmeticException { 1051 long ret; 1052 if (a > b) { 1053 // use symmetry to reduce boundary cases 1054 ret = addAndCheck(b, a, pattern); 1055 } else { 1056 // assert a <= b 1057 1058 if (a < 0) { 1059 if (b < 0) { 1060 // check for negative overflow 1061 if (Long.MIN_VALUE - b <= a) { 1062 ret = a + b; 1063 } else { 1064 throw new MathArithmeticException(pattern, a, b); 1065 } 1066 } else { 1067 // opposite sign addition is always safe 1068 ret = a + b; 1069 } 1070 } else { 1071 // assert a >= 0 1072 // assert b >= 0 1073 1074 // check for positive overflow 1075 if (a <= Long.MAX_VALUE - b) { 1076 ret = a + b; 1077 } else { 1078 throw new MathArithmeticException(pattern, a, b); 1079 } 1080 } 1081 } 1082 return ret; 1083 } 1084 1085 /** 1086 * Check binomial preconditions. 1087 * 1088 * @param n Size of the set. 1089 * @param k Size of the subsets to be counted. 1090 * @throws NotPositiveException if {@code n < 0}. 1091 * @throws NumberIsTooLargeException if {@code k > n}. 1092 */ 1093 private static void checkBinomial(final int n, final int k) throws NumberIsTooLargeException, NotPositiveException { 1094 if (n < k) { 1095 throw new NumberIsTooLargeException(LocalizedFormats.BINOMIAL_INVALID_PARAMETERS_ORDER, 1096 k, n, true); 1097 } 1098 if (n < 0) { 1099 throw new NotPositiveException(LocalizedFormats.BINOMIAL_NEGATIVE_PARAMETER, n); 1100 } 1101 } 1102 1103 /** 1104 * Returns true if the argument is a power of two. 1105 * 1106 * @param n the number to test 1107 * @return true if the argument is a power of two 1108 */ 1109 public static boolean isPowerOfTwo(long n) { 1110 return (n > 0) && ((n & (n - 1)) == 0); 1111 } 1112 }