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 */
017package org.apache.commons.math3.util;
018
019import java.math.BigInteger;
020
021import org.apache.commons.math3.exception.MathArithmeticException;
022import org.apache.commons.math3.exception.NotPositiveException;
023import org.apache.commons.math3.exception.NumberIsTooLargeException;
024import org.apache.commons.math3.exception.util.Localizable;
025import org.apache.commons.math3.exception.util.LocalizedFormats;
026
027/**
028 * Some useful, arithmetics related, additions to the built-in functions in
029 * {@link Math}.
030 *
031 */
032public final class ArithmeticUtils {
033
034    /** Private constructor. */
035    private ArithmeticUtils() {
036        super();
037    }
038
039    /**
040     * Add two integers, checking for overflow.
041     *
042     * @param x an addend
043     * @param y an addend
044     * @return the sum {@code x+y}
045     * @throws MathArithmeticException if the result can not be represented
046     * as an {@code int}.
047     * @since 1.1
048     */
049    public static int addAndCheck(int x, int y)
050            throws MathArithmeticException {
051        long s = (long)x + (long)y;
052        if (s < Integer.MIN_VALUE || s > Integer.MAX_VALUE) {
053            throw new MathArithmeticException(LocalizedFormats.OVERFLOW_IN_ADDITION, x, y);
054        }
055        return (int)s;
056    }
057
058    /**
059     * Add two long integers, checking for overflow.
060     *
061     * @param a an addend
062     * @param b an addend
063     * @return the sum {@code a+b}
064     * @throws MathArithmeticException if the result can not be represented as an long
065     * @since 1.2
066     */
067    public static long addAndCheck(long a, long b) throws MathArithmeticException {
068        return addAndCheck(a, b, LocalizedFormats.OVERFLOW_IN_ADDITION);
069    }
070
071    /**
072     * Returns an exact representation of the <a
073     * href="http://mathworld.wolfram.com/BinomialCoefficient.html"> Binomial
074     * Coefficient</a>, "{@code n choose k}", the number of
075     * {@code k}-element subsets that can be selected from an
076     * {@code n}-element set.
077     * <p>
078     * <Strong>Preconditions</strong>:
079     * <ul>
080     * <li> {@code 0 <= k <= n } (otherwise
081     * {@code IllegalArgumentException} is thrown)</li>
082     * <li> The result is small enough to fit into a {@code long}. The
083     * largest value of {@code n} for which all coefficients are
084     * {@code  < Long.MAX_VALUE} is 66. If the computed value exceeds
085     * {@code Long.MAX_VALUE} an {@code ArithMeticException} is
086     * thrown.</li>
087     * </ul></p>
088     *
089     * @param n the size of the set
090     * @param k the size of the subsets to be counted
091     * @return {@code n choose k}
092     * @throws NotPositiveException if {@code n < 0}.
093     * @throws NumberIsTooLargeException if {@code k > n}.
094     * @throws MathArithmeticException if the result is too large to be
095     * represented by a long integer.
096     * @deprecated use {@link CombinatoricsUtils#binomialCoefficient(int, int)}
097     */
098    @Deprecated
099    public static long binomialCoefficient(final int n, final int k)
100        throws NotPositiveException, NumberIsTooLargeException, MathArithmeticException {
101       return CombinatoricsUtils.binomialCoefficient(n, k);
102    }
103
104    /**
105     * Returns a {@code double} representation of the <a
106     * href="http://mathworld.wolfram.com/BinomialCoefficient.html"> Binomial
107     * Coefficient</a>, "{@code n choose k}", the number of
108     * {@code k}-element subsets that can be selected from an
109     * {@code n}-element set.
110     * <p>
111     * <Strong>Preconditions</strong>:
112     * <ul>
113     * <li> {@code 0 <= k <= n } (otherwise
114     * {@code IllegalArgumentException} is thrown)</li>
115     * <li> The result is small enough to fit into a {@code double}. The
116     * largest value of {@code n} for which all coefficients are <
117     * Double.MAX_VALUE is 1029. If the computed value exceeds Double.MAX_VALUE,
118     * Double.POSITIVE_INFINITY is returned</li>
119     * </ul></p>
120     *
121     * @param n the size of the set
122     * @param k the size of the subsets to be counted
123     * @return {@code n choose k}
124     * @throws NotPositiveException if {@code n < 0}.
125     * @throws NumberIsTooLargeException if {@code k > n}.
126     * @throws MathArithmeticException if the result is too large to be
127     * represented by a long integer.
128     * @deprecated use {@link CombinatoricsUtils#binomialCoefficientDouble(int, int)}
129     */
130    @Deprecated
131    public static double binomialCoefficientDouble(final int n, final int k)
132        throws NotPositiveException, NumberIsTooLargeException, MathArithmeticException {
133        return CombinatoricsUtils.binomialCoefficientDouble(n, k);
134    }
135
136    /**
137     * Returns the natural {@code log} of the <a
138     * href="http://mathworld.wolfram.com/BinomialCoefficient.html"> Binomial
139     * Coefficient</a>, "{@code n choose k}", the number of
140     * {@code k}-element subsets that can be selected from an
141     * {@code n}-element set.
142     * <p>
143     * <Strong>Preconditions</strong>:
144     * <ul>
145     * <li> {@code 0 <= k <= n } (otherwise
146     * {@code IllegalArgumentException} is thrown)</li>
147     * </ul></p>
148     *
149     * @param n the size of the set
150     * @param k the size of the subsets to be counted
151     * @return {@code n choose k}
152     * @throws NotPositiveException if {@code n < 0}.
153     * @throws NumberIsTooLargeException if {@code k > n}.
154     * @throws MathArithmeticException if the result is too large to be
155     * represented by a long integer.
156     * @deprecated use {@link CombinatoricsUtils#binomialCoefficientLog(int, int)}
157     */
158    @Deprecated
159    public static double binomialCoefficientLog(final int n, final int k)
160        throws NotPositiveException, NumberIsTooLargeException, MathArithmeticException {
161        return CombinatoricsUtils.binomialCoefficientLog(n, k);
162    }
163
164    /**
165     * Returns n!. Shorthand for {@code n} <a
166     * href="http://mathworld.wolfram.com/Factorial.html"> Factorial</a>, the
167     * product of the numbers {@code 1,...,n}.
168     * <p>
169     * <Strong>Preconditions</strong>:
170     * <ul>
171     * <li> {@code n >= 0} (otherwise
172     * {@code IllegalArgumentException} is thrown)</li>
173     * <li> The result is small enough to fit into a {@code long}. The
174     * largest value of {@code n} for which {@code n!} <
175     * Long.MAX_VALUE} is 20. If the computed value exceeds {@code Long.MAX_VALUE}
176     * an {@code ArithMeticException } is thrown.</li>
177     * </ul>
178     * </p>
179     *
180     * @param n argument
181     * @return {@code n!}
182     * @throws MathArithmeticException if the result is too large to be represented
183     * by a {@code long}.
184     * @throws NotPositiveException if {@code n < 0}.
185     * @throws MathArithmeticException if {@code n > 20}: The factorial value is too
186     * large to fit in a {@code long}.
187     * @deprecated use {@link CombinatoricsUtils#factorial(int)}
188     */
189    @Deprecated
190    public static long factorial(final int n) throws NotPositiveException, MathArithmeticException {
191        return CombinatoricsUtils.factorial(n);
192    }
193
194    /**
195     * Compute n!, the<a href="http://mathworld.wolfram.com/Factorial.html">
196     * factorial</a> of {@code n} (the product of the numbers 1 to n), as a
197     * {@code double}.
198     * The result should be small enough to fit into a {@code double}: The
199     * largest {@code n} for which {@code n! < Double.MAX_VALUE} is 170.
200     * If the computed value exceeds {@code Double.MAX_VALUE},
201     * {@code Double.POSITIVE_INFINITY} is returned.
202     *
203     * @param n Argument.
204     * @return {@code n!}
205     * @throws NotPositiveException if {@code n < 0}.
206     * @deprecated use {@link CombinatoricsUtils#factorialDouble(int)}
207     */
208    @Deprecated
209    public static double factorialDouble(final int n) throws NotPositiveException {
210         return CombinatoricsUtils.factorialDouble(n);
211    }
212
213    /**
214     * Compute the natural logarithm of the factorial of {@code n}.
215     *
216     * @param n Argument.
217     * @return {@code n!}
218     * @throws NotPositiveException if {@code n < 0}.
219     * @deprecated use {@link CombinatoricsUtils#factorialLog(int)}
220     */
221    @Deprecated
222    public static double factorialLog(final int n) throws NotPositiveException {
223        return CombinatoricsUtils.factorialLog(n);
224    }
225
226    /**
227     * Computes the greatest common divisor of the absolute value of two
228     * numbers, using a modified version of the "binary gcd" method.
229     * See Knuth 4.5.2 algorithm B.
230     * The algorithm is due to Josef Stein (1961).
231     * <br/>
232     * Special cases:
233     * <ul>
234     *  <li>The invocations
235     *   {@code gcd(Integer.MIN_VALUE, Integer.MIN_VALUE)},
236     *   {@code gcd(Integer.MIN_VALUE, 0)} and
237     *   {@code gcd(0, Integer.MIN_VALUE)} throw an
238     *   {@code ArithmeticException}, because the result would be 2^31, which
239     *   is too large for an int value.</li>
240     *  <li>The result of {@code gcd(x, x)}, {@code gcd(0, x)} and
241     *   {@code gcd(x, 0)} is the absolute value of {@code x}, except
242     *   for the special cases above.</li>
243     *  <li>The invocation {@code gcd(0, 0)} is the only one which returns
244     *   {@code 0}.</li>
245     * </ul>
246     *
247     * @param p Number.
248     * @param q Number.
249     * @return the greatest common divisor (never negative).
250     * @throws MathArithmeticException if the result cannot be represented as
251     * a non-negative {@code int} value.
252     * @since 1.1
253     */
254    public static int gcd(int p, int q) throws MathArithmeticException {
255        int a = p;
256        int b = q;
257        if (a == 0 ||
258            b == 0) {
259            if (a == Integer.MIN_VALUE ||
260                b == Integer.MIN_VALUE) {
261                throw new MathArithmeticException(LocalizedFormats.GCD_OVERFLOW_32_BITS,
262                                                  p, q);
263            }
264            return FastMath.abs(a + b);
265        }
266
267        long al = a;
268        long bl = b;
269        boolean useLong = false;
270        if (a < 0) {
271            if(Integer.MIN_VALUE == a) {
272                useLong = true;
273            } else {
274                a = -a;
275            }
276            al = -al;
277        }
278        if (b < 0) {
279            if (Integer.MIN_VALUE == b) {
280                useLong = true;
281            } else {
282                b = -b;
283            }
284            bl = -bl;
285        }
286        if (useLong) {
287            if(al == bl) {
288                throw new MathArithmeticException(LocalizedFormats.GCD_OVERFLOW_32_BITS,
289                                                  p, q);
290            }
291            long blbu = bl;
292            bl = al;
293            al = blbu % al;
294            if (al == 0) {
295                if (bl > Integer.MAX_VALUE) {
296                    throw new MathArithmeticException(LocalizedFormats.GCD_OVERFLOW_32_BITS,
297                                                      p, q);
298                }
299                return (int) bl;
300            }
301            blbu = bl;
302
303            // Now "al" and "bl" fit in an "int".
304            b = (int) al;
305            a = (int) (blbu % al);
306        }
307
308        return gcdPositive(a, b);
309    }
310
311    /**
312     * Computes the greatest common divisor of two <em>positive</em> numbers
313     * (this precondition is <em>not</em> checked and the result is undefined
314     * if not fulfilled) using the "binary gcd" method which avoids division
315     * and modulo operations.
316     * See Knuth 4.5.2 algorithm B.
317     * The algorithm is due to Josef Stein (1961).
318     * <br/>
319     * Special cases:
320     * <ul>
321     *  <li>The result of {@code gcd(x, x)}, {@code gcd(0, x)} and
322     *   {@code gcd(x, 0)} is the value of {@code x}.</li>
323     *  <li>The invocation {@code gcd(0, 0)} is the only one which returns
324     *   {@code 0}.</li>
325     * </ul>
326     *
327     * @param a Positive number.
328     * @param b Positive number.
329     * @return the greatest common divisor.
330     */
331    private static int gcdPositive(int a, int b) {
332        if (a == 0) {
333            return b;
334        }
335        else if (b == 0) {
336            return a;
337        }
338
339        // Make "a" and "b" odd, keeping track of common power of 2.
340        final int aTwos = Integer.numberOfTrailingZeros(a);
341        a >>= aTwos;
342        final int bTwos = Integer.numberOfTrailingZeros(b);
343        b >>= bTwos;
344        final int shift = FastMath.min(aTwos, bTwos);
345
346        // "a" and "b" are positive.
347        // If a > b then "gdc(a, b)" is equal to "gcd(a - b, b)".
348        // If a < b then "gcd(a, b)" is equal to "gcd(b - a, a)".
349        // Hence, in the successive iterations:
350        //  "a" becomes the absolute difference of the current values,
351        //  "b" becomes the minimum of the current values.
352        while (a != b) {
353            final int delta = a - b;
354            b = Math.min(a, b);
355            a = Math.abs(delta);
356
357            // Remove any power of 2 in "a" ("b" is guaranteed to be odd).
358            a >>= Integer.numberOfTrailingZeros(a);
359        }
360
361        // Recover the common power of 2.
362        return a << shift;
363    }
364
365    /**
366     * <p>
367     * Gets the greatest common divisor of the absolute value of two numbers,
368     * using the "binary gcd" method which avoids division and modulo
369     * operations. See Knuth 4.5.2 algorithm B. This algorithm is due to Josef
370     * Stein (1961).
371     * </p>
372     * Special cases:
373     * <ul>
374     * <li>The invocations
375     * {@code gcd(Long.MIN_VALUE, Long.MIN_VALUE)},
376     * {@code gcd(Long.MIN_VALUE, 0L)} and
377     * {@code gcd(0L, Long.MIN_VALUE)} throw an
378     * {@code ArithmeticException}, because the result would be 2^63, which
379     * is too large for a long value.</li>
380     * <li>The result of {@code gcd(x, x)}, {@code gcd(0L, x)} and
381     * {@code gcd(x, 0L)} is the absolute value of {@code x}, except
382     * for the special cases above.
383     * <li>The invocation {@code gcd(0L, 0L)} is the only one which returns
384     * {@code 0L}.</li>
385     * </ul>
386     *
387     * @param p Number.
388     * @param q Number.
389     * @return the greatest common divisor, never negative.
390     * @throws MathArithmeticException if the result cannot be represented as
391     * a non-negative {@code long} value.
392     * @since 2.1
393     */
394    public static long gcd(final long p, final long q) throws MathArithmeticException {
395        long u = p;
396        long v = q;
397        if ((u == 0) || (v == 0)) {
398            if ((u == Long.MIN_VALUE) || (v == Long.MIN_VALUE)){
399                throw new MathArithmeticException(LocalizedFormats.GCD_OVERFLOW_64_BITS,
400                                                  p, q);
401            }
402            return FastMath.abs(u) + FastMath.abs(v);
403        }
404        // keep u and v negative, as negative integers range down to
405        // -2^63, while positive numbers can only be as large as 2^63-1
406        // (i.e. we can't necessarily negate a negative number without
407        // overflow)
408        /* assert u!=0 && v!=0; */
409        if (u > 0) {
410            u = -u;
411        } // make u negative
412        if (v > 0) {
413            v = -v;
414        } // make v negative
415        // B1. [Find power of 2]
416        int k = 0;
417        while ((u & 1) == 0 && (v & 1) == 0 && k < 63) { // while u and v are
418                                                            // both even...
419            u /= 2;
420            v /= 2;
421            k++; // cast out twos.
422        }
423        if (k == 63) {
424            throw new MathArithmeticException(LocalizedFormats.GCD_OVERFLOW_64_BITS,
425                                              p, q);
426        }
427        // B2. Initialize: u and v have been divided by 2^k and at least
428        // one is odd.
429        long t = ((u & 1) == 1) ? v : -(u / 2)/* B3 */;
430        // t negative: u was odd, v may be even (t replaces v)
431        // t positive: u was even, v is odd (t replaces u)
432        do {
433            /* assert u<0 && v<0; */
434            // B4/B3: cast out twos from t.
435            while ((t & 1) == 0) { // while t is even..
436                t /= 2; // cast out twos
437            }
438            // B5 [reset max(u,v)]
439            if (t > 0) {
440                u = -t;
441            } else {
442                v = t;
443            }
444            // B6/B3. at this point both u and v should be odd.
445            t = (v - u) / 2;
446            // |u| larger: t positive (replace u)
447            // |v| larger: t negative (replace v)
448        } while (t != 0);
449        return -u * (1L << k); // gcd is u*2^k
450    }
451
452    /**
453     * <p>
454     * Returns the least common multiple of the absolute value of two numbers,
455     * using the formula {@code lcm(a,b) = (a / gcd(a,b)) * b}.
456     * </p>
457     * Special cases:
458     * <ul>
459     * <li>The invocations {@code lcm(Integer.MIN_VALUE, n)} and
460     * {@code lcm(n, Integer.MIN_VALUE)}, where {@code abs(n)} is a
461     * power of 2, throw an {@code ArithmeticException}, because the result
462     * would be 2^31, which is too large for an int value.</li>
463     * <li>The result of {@code lcm(0, x)} and {@code lcm(x, 0)} is
464     * {@code 0} for any {@code x}.
465     * </ul>
466     *
467     * @param a Number.
468     * @param b Number.
469     * @return the least common multiple, never negative.
470     * @throws MathArithmeticException if the result cannot be represented as
471     * a non-negative {@code int} value.
472     * @since 1.1
473     */
474    public static int lcm(int a, int b) throws MathArithmeticException {
475        if (a == 0 || b == 0){
476            return 0;
477        }
478        int lcm = FastMath.abs(ArithmeticUtils.mulAndCheck(a / gcd(a, b), b));
479        if (lcm == Integer.MIN_VALUE) {
480            throw new MathArithmeticException(LocalizedFormats.LCM_OVERFLOW_32_BITS,
481                                              a, b);
482        }
483        return lcm;
484    }
485
486    /**
487     * <p>
488     * Returns the least common multiple of the absolute value of two numbers,
489     * using the formula {@code lcm(a,b) = (a / gcd(a,b)) * b}.
490     * </p>
491     * Special cases:
492     * <ul>
493     * <li>The invocations {@code lcm(Long.MIN_VALUE, n)} and
494     * {@code lcm(n, Long.MIN_VALUE)}, where {@code abs(n)} is a
495     * power of 2, throw an {@code ArithmeticException}, because the result
496     * would be 2^63, which is too large for an int value.</li>
497     * <li>The result of {@code lcm(0L, x)} and {@code lcm(x, 0L)} is
498     * {@code 0L} for any {@code x}.
499     * </ul>
500     *
501     * @param a Number.
502     * @param b Number.
503     * @return the least common multiple, never negative.
504     * @throws MathArithmeticException if the result cannot be represented
505     * as a non-negative {@code long} value.
506     * @since 2.1
507     */
508    public static long lcm(long a, long b) throws MathArithmeticException {
509        if (a == 0 || b == 0){
510            return 0;
511        }
512        long lcm = FastMath.abs(ArithmeticUtils.mulAndCheck(a / gcd(a, b), b));
513        if (lcm == Long.MIN_VALUE){
514            throw new MathArithmeticException(LocalizedFormats.LCM_OVERFLOW_64_BITS,
515                                              a, b);
516        }
517        return lcm;
518    }
519
520    /**
521     * Multiply two integers, checking for overflow.
522     *
523     * @param x Factor.
524     * @param y Factor.
525     * @return the product {@code x * y}.
526     * @throws MathArithmeticException if the result can not be
527     * represented as an {@code int}.
528     * @since 1.1
529     */
530    public static int mulAndCheck(int x, int y) throws MathArithmeticException {
531        long m = ((long)x) * ((long)y);
532        if (m < Integer.MIN_VALUE || m > Integer.MAX_VALUE) {
533            throw new MathArithmeticException();
534        }
535        return (int)m;
536    }
537
538    /**
539     * Multiply two long integers, checking for overflow.
540     *
541     * @param a Factor.
542     * @param b Factor.
543     * @return the product {@code a * b}.
544     * @throws MathArithmeticException if the result can not be represented
545     * as a {@code long}.
546     * @since 1.2
547     */
548    public static long mulAndCheck(long a, long b) throws MathArithmeticException {
549        long ret;
550        if (a > b) {
551            // use symmetry to reduce boundary cases
552            ret = mulAndCheck(b, a);
553        } else {
554            if (a < 0) {
555                if (b < 0) {
556                    // check for positive overflow with negative a, negative b
557                    if (a >= Long.MAX_VALUE / b) {
558                        ret = a * b;
559                    } else {
560                        throw new MathArithmeticException();
561                    }
562                } else if (b > 0) {
563                    // check for negative overflow with negative a, positive b
564                    if (Long.MIN_VALUE / b <= a) {
565                        ret = a * b;
566                    } else {
567                        throw new MathArithmeticException();
568
569                    }
570                } else {
571                    // assert b == 0
572                    ret = 0;
573                }
574            } else if (a > 0) {
575                // assert a > 0
576                // assert b > 0
577
578                // check for positive overflow with positive a, positive b
579                if (a <= Long.MAX_VALUE / b) {
580                    ret = a * b;
581                } else {
582                    throw new MathArithmeticException();
583                }
584            } else {
585                // assert a == 0
586                ret = 0;
587            }
588        }
589        return ret;
590    }
591
592    /**
593     * Subtract two integers, checking for overflow.
594     *
595     * @param x Minuend.
596     * @param y Subtrahend.
597     * @return the difference {@code x - y}.
598     * @throws MathArithmeticException if the result can not be represented
599     * as an {@code int}.
600     * @since 1.1
601     */
602    public static int subAndCheck(int x, int y) throws MathArithmeticException {
603        long s = (long)x - (long)y;
604        if (s < Integer.MIN_VALUE || s > Integer.MAX_VALUE) {
605            throw new MathArithmeticException(LocalizedFormats.OVERFLOW_IN_SUBTRACTION, x, y);
606        }
607        return (int)s;
608    }
609
610    /**
611     * Subtract two long integers, checking for overflow.
612     *
613     * @param a Value.
614     * @param b Value.
615     * @return the difference {@code a - b}.
616     * @throws MathArithmeticException if the result can not be represented as a
617     * {@code long}.
618     * @since 1.2
619     */
620    public static long subAndCheck(long a, long b) throws MathArithmeticException {
621        long ret;
622        if (b == Long.MIN_VALUE) {
623            if (a < 0) {
624                ret = a - b;
625            } else {
626                throw new MathArithmeticException(LocalizedFormats.OVERFLOW_IN_ADDITION, a, -b);
627            }
628        } else {
629            // use additive inverse
630            ret = addAndCheck(a, -b, LocalizedFormats.OVERFLOW_IN_ADDITION);
631        }
632        return ret;
633    }
634
635    /**
636     * Raise an int to an int power.
637     *
638     * @param k Number to raise.
639     * @param e Exponent (must be positive or zero).
640     * @return \( k^e \)
641     * @throws NotPositiveException if {@code e < 0}.
642     * @throws MathArithmeticException if the result would overflow.
643     */
644    public static int pow(final int k,
645                          final int e)
646        throws NotPositiveException,
647               MathArithmeticException {
648        if (e < 0) {
649            throw new NotPositiveException(LocalizedFormats.EXPONENT, e);
650        }
651
652        try {
653            int exp = e;
654            int result = 1;
655            int k2p    = k;
656            while (true) {
657                if ((exp & 0x1) != 0) {
658                    result = mulAndCheck(result, k2p);
659                }
660
661                exp >>= 1;
662                if (exp == 0) {
663                    break;
664                }
665
666                k2p = mulAndCheck(k2p, k2p);
667            }
668
669            return result;
670        } catch (MathArithmeticException mae) {
671            // Add context information.
672            mae.getContext().addMessage(LocalizedFormats.OVERFLOW);
673            mae.getContext().addMessage(LocalizedFormats.BASE, k);
674            mae.getContext().addMessage(LocalizedFormats.EXPONENT, e);
675
676            // Rethrow.
677            throw mae;
678        }
679    }
680
681    /**
682     * Raise an int to a long power.
683     *
684     * @param k Number to raise.
685     * @param e Exponent (must be positive or zero).
686     * @return k<sup>e</sup>
687     * @throws NotPositiveException if {@code e < 0}.
688     * @deprecated As of 3.3. Please use {@link #pow(int,int)} instead.
689     */
690    @Deprecated
691    public static int pow(final int k, long e) throws NotPositiveException {
692        if (e < 0) {
693            throw new NotPositiveException(LocalizedFormats.EXPONENT, e);
694        }
695
696        int result = 1;
697        int k2p    = k;
698        while (e != 0) {
699            if ((e & 0x1) != 0) {
700                result *= k2p;
701            }
702            k2p *= k2p;
703            e >>= 1;
704        }
705
706        return result;
707    }
708
709    /**
710     * Raise a long to an int power.
711     *
712     * @param k Number to raise.
713     * @param e Exponent (must be positive or zero).
714     * @return \( k^e \)
715     * @throws NotPositiveException if {@code e < 0}.
716     * @throws MathArithmeticException if the result would overflow.
717     */
718    public static long pow(final long k,
719                           final int e)
720        throws NotPositiveException,
721               MathArithmeticException {
722        if (e < 0) {
723            throw new NotPositiveException(LocalizedFormats.EXPONENT, e);
724        }
725
726        try {
727            int exp = e;
728            long result = 1;
729            long k2p    = k;
730            while (true) {
731                if ((exp & 0x1) != 0) {
732                    result = mulAndCheck(result, k2p);
733                }
734
735                exp >>= 1;
736                if (exp == 0) {
737                    break;
738                }
739
740                k2p = mulAndCheck(k2p, k2p);
741            }
742
743            return result;
744        } catch (MathArithmeticException mae) {
745            // Add context information.
746            mae.getContext().addMessage(LocalizedFormats.OVERFLOW);
747            mae.getContext().addMessage(LocalizedFormats.BASE, k);
748            mae.getContext().addMessage(LocalizedFormats.EXPONENT, e);
749
750            // Rethrow.
751            throw mae;
752        }
753    }
754
755    /**
756     * Raise a long to a long power.
757     *
758     * @param k Number to raise.
759     * @param e Exponent (must be positive or zero).
760     * @return k<sup>e</sup>
761     * @throws NotPositiveException if {@code e < 0}.
762     * @deprecated As of 3.3. Please use {@link #pow(long,int)} instead.
763     */
764    @Deprecated
765    public static long pow(final long k, long e) throws NotPositiveException {
766        if (e < 0) {
767            throw new NotPositiveException(LocalizedFormats.EXPONENT, e);
768        }
769
770        long result = 1l;
771        long k2p    = k;
772        while (e != 0) {
773            if ((e & 0x1) != 0) {
774                result *= k2p;
775            }
776            k2p *= k2p;
777            e >>= 1;
778        }
779
780        return result;
781    }
782
783    /**
784     * Raise a BigInteger to an int power.
785     *
786     * @param k Number to raise.
787     * @param e Exponent (must be positive or zero).
788     * @return k<sup>e</sup>
789     * @throws NotPositiveException if {@code e < 0}.
790     */
791    public static BigInteger pow(final BigInteger k, int e) throws NotPositiveException {
792        if (e < 0) {
793            throw new NotPositiveException(LocalizedFormats.EXPONENT, e);
794        }
795
796        return k.pow(e);
797    }
798
799    /**
800     * Raise a BigInteger to a long power.
801     *
802     * @param k Number to raise.
803     * @param e Exponent (must be positive or zero).
804     * @return k<sup>e</sup>
805     * @throws NotPositiveException if {@code e < 0}.
806     */
807    public static BigInteger pow(final BigInteger k, long e) throws NotPositiveException {
808        if (e < 0) {
809            throw new NotPositiveException(LocalizedFormats.EXPONENT, e);
810        }
811
812        BigInteger result = BigInteger.ONE;
813        BigInteger k2p    = k;
814        while (e != 0) {
815            if ((e & 0x1) != 0) {
816                result = result.multiply(k2p);
817            }
818            k2p = k2p.multiply(k2p);
819            e >>= 1;
820        }
821
822        return result;
823
824    }
825
826    /**
827     * Raise a BigInteger to a BigInteger power.
828     *
829     * @param k Number to raise.
830     * @param e Exponent (must be positive or zero).
831     * @return k<sup>e</sup>
832     * @throws NotPositiveException if {@code e < 0}.
833     */
834    public static BigInteger pow(final BigInteger k, BigInteger e) throws NotPositiveException {
835        if (e.compareTo(BigInteger.ZERO) < 0) {
836            throw new NotPositiveException(LocalizedFormats.EXPONENT, e);
837        }
838
839        BigInteger result = BigInteger.ONE;
840        BigInteger k2p    = k;
841        while (!BigInteger.ZERO.equals(e)) {
842            if (e.testBit(0)) {
843                result = result.multiply(k2p);
844            }
845            k2p = k2p.multiply(k2p);
846            e = e.shiftRight(1);
847        }
848
849        return result;
850    }
851
852    /**
853     * Returns the <a
854     * href="http://mathworld.wolfram.com/StirlingNumberoftheSecondKind.html">
855     * Stirling number of the second kind</a>, "{@code S(n,k)}", the number of
856     * ways of partitioning an {@code n}-element set into {@code k} non-empty
857     * subsets.
858     * <p>
859     * The preconditions are {@code 0 <= k <= n } (otherwise
860     * {@code NotPositiveException} is thrown)
861     * </p>
862     * @param n the size of the set
863     * @param k the number of non-empty subsets
864     * @return {@code S(n,k)}
865     * @throws NotPositiveException if {@code k < 0}.
866     * @throws NumberIsTooLargeException if {@code k > n}.
867     * @throws MathArithmeticException if some overflow happens, typically for n exceeding 25 and
868     * k between 20 and n-2 (S(n,n-1) is handled specifically and does not overflow)
869     * @since 3.1
870     * @deprecated use {@link CombinatoricsUtils#stirlingS2(int, int)}
871     */
872    @Deprecated
873    public static long stirlingS2(final int n, final int k)
874        throws NotPositiveException, NumberIsTooLargeException, MathArithmeticException {
875        return CombinatoricsUtils.stirlingS2(n, k);
876
877    }
878
879    /**
880     * Add two long integers, checking for overflow.
881     *
882     * @param a Addend.
883     * @param b Addend.
884     * @param pattern Pattern to use for any thrown exception.
885     * @return the sum {@code a + b}.
886     * @throws MathArithmeticException if the result cannot be represented
887     * as a {@code long}.
888     * @since 1.2
889     */
890     private static long addAndCheck(long a, long b, Localizable pattern) throws MathArithmeticException {
891         final long result = a + b;
892         if (!((a ^ b) < 0 | (a ^ result) >= 0)) {
893             throw new MathArithmeticException(pattern, a, b);
894         }
895         return result;
896    }
897
898    /**
899     * Returns true if the argument is a power of two.
900     *
901     * @param n the number to test
902     * @return true if the argument is a power of two
903     */
904    public static boolean isPowerOfTwo(long n) {
905        return (n > 0) && ((n & (n - 1)) == 0);
906    }
907}