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 }