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