    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    }