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