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    }