001    /*
002     * Licensed to the Apache Software Foundation (ASF) under one or more
003     * contributor license agreements.  See the NOTICE file distributed with
004     * this work for additional information regarding copyright ownership.
005     * The ASF licenses this file to You under the Apache License, Version 2.0
006     * (the "License"); you may not use this file except in compliance with
007     * the License.  You may obtain a copy of the License at
008     *
009     *      http://www.apache.org/licenses/LICENSE-2.0
010     *
011     * Unless required by applicable law or agreed to in writing, software
012     * distributed under the License is distributed on an "AS IS" BASIS,
013     * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014     * See the License for the specific language governing permissions and
015     * limitations under the License.
016     */
017    
018    package org.apache.commons.math.util;
019    
020    import java.math.BigDecimal;
021    import java.math.BigInteger;
022    import java.util.Arrays;
023    
024    import org.apache.commons.math.MathRuntimeException;
025    
026    /**
027     * Some useful additions to the built-in functions in {@link Math}.
028     * @version $Revision: 772119 $ $Date: 2009-05-06 05:43:28 -0400 (Wed, 06 May 2009) $
029     */
030    public final class MathUtils {
031    
032        /** Smallest positive number such that 1 - EPSILON is not numerically equal to 1. */
033        public static final double EPSILON = 0x1.0p-53;
034    
035        /** Safe minimum, such that 1 / SAFE_MIN does not overflow.
036         * <p>In IEEE 754 arithmetic, this is also the smallest normalized
037         * number 2<sup>-1022</sup>.</p>
038         */
039        public static final double SAFE_MIN = 0x1.0p-1022;
040    
041        /** -1.0 cast as a byte. */
042        private static final byte  NB = (byte)-1;
043    
044        /** -1.0 cast as a short. */
045        private static final short NS = (short)-1;
046    
047        /** 1.0 cast as a byte. */
048        private static final byte  PB = (byte)1;
049    
050        /** 1.0 cast as a short. */
051        private static final short PS = (short)1;
052    
053        /** 0.0 cast as a byte. */
054        private static final byte  ZB = (byte)0;
055    
056        /** 0.0 cast as a short. */
057        private static final short ZS = (short)0;
058    
059        /** 2 &pi;. */
060        private static final double TWO_PI = 2 * Math.PI;
061    
062        /** Gap between NaN and regular numbers. */
063        private static final int NAN_GAP = 4 * 1024 * 1024;
064    
065        /** Offset to order signed double numbers lexicographically. */
066        private static final long SGN_MASK = 0x8000000000000000L;
067    
068        /**
069         * Private Constructor
070         */
071        private MathUtils() {
072            super();
073        }
074    
075        /**
076         * Add two integers, checking for overflow.
077         * 
078         * @param x an addend
079         * @param y an addend
080         * @return the sum <code>x+y</code>
081         * @throws ArithmeticException if the result can not be represented as an
082         *         int
083         * @since 1.1
084         */
085        public static int addAndCheck(int x, int y) {
086            long s = (long)x + (long)y;
087            if (s < Integer.MIN_VALUE || s > Integer.MAX_VALUE) {
088                throw new ArithmeticException("overflow: add");
089            }
090            return (int)s;
091        }
092    
093        /**
094         * Add two long integers, checking for overflow.
095         * 
096         * @param a an addend
097         * @param b an addend
098         * @return the sum <code>a+b</code>
099         * @throws ArithmeticException if the result can not be represented as an
100         *         long
101         * @since 1.2
102         */
103        public static long addAndCheck(long a, long b) {
104            return addAndCheck(a, b, "overflow: add");
105        }
106        
107        /**
108         * Add two long integers, checking for overflow.
109         * 
110         * @param a an addend
111         * @param b an addend
112         * @param msg the message to use for any thrown exception.
113         * @return the sum <code>a+b</code>
114         * @throws ArithmeticException if the result can not be represented as an
115         *         long
116         * @since 1.2
117         */
118        private static long addAndCheck(long a, long b, String msg) {
119            long ret;
120            if (a > b) {
121                // use symmetry to reduce boundary cases
122                ret = addAndCheck(b, a, msg);
123            } else {
124                // assert a <= b
125                
126                if (a < 0) {
127                    if (b < 0) {
128                        // check for negative overflow
129                        if (Long.MIN_VALUE - b <= a) {
130                            ret = a + b;
131                        } else {
132                            throw new ArithmeticException(msg);
133                        }
134                    } else {
135                        // opposite sign addition is always safe
136                        ret = a + b;
137                    }
138                } else {
139                    // assert a >= 0
140                    // assert b >= 0
141    
142                    // check for positive overflow
143                    if (a <= Long.MAX_VALUE - b) {
144                        ret = a + b;
145                    } else {
146                        throw new ArithmeticException(msg);
147                    }
148                }
149            }
150            return ret;
151        }
152        
153        /**
154         * Returns an exact representation of the <a
155         * href="http://mathworld.wolfram.com/BinomialCoefficient.html"> Binomial
156         * Coefficient</a>, "<code>n choose k</code>", the number of
157         * <code>k</code>-element subsets that can be selected from an
158         * <code>n</code>-element set.
159         * <p>
160         * <Strong>Preconditions</strong>:
161         * <ul>
162         * <li> <code>0 <= k <= n </code> (otherwise
163         * <code>IllegalArgumentException</code> is thrown)</li>
164         * <li> The result is small enough to fit into a <code>long</code>. The
165         * largest value of <code>n</code> for which all coefficients are
166         * <code> < Long.MAX_VALUE</code> is 66. If the computed value exceeds
167         * <code>Long.MAX_VALUE</code> an <code>ArithMeticException</code> is
168         * thrown.</li>
169         * </ul></p>
170         * 
171         * @param n the size of the set
172         * @param k the size of the subsets to be counted
173         * @return <code>n choose k</code>
174         * @throws IllegalArgumentException if preconditions are not met.
175         * @throws ArithmeticException if the result is too large to be represented
176         *         by a long integer.
177         */
178        public static long binomialCoefficient(final int n, final int k) {
179            checkBinomial(n, k);
180            if ((n == k) || (k == 0)) {
181                return 1;
182            }
183            if ((k == 1) || (k == n - 1)) {
184                return n;
185            }
186            // Use symmetry for large k
187            if (k > n / 2)
188                return binomialCoefficient(n, n - k);
189            
190            // We use the formula
191            // (n choose k) = n! / (n-k)! / k!
192            // (n choose k) == ((n-k+1)*...*n) / (1*...*k)
193            // which could be written
194            // (n choose k) == (n-1 choose k-1) * n / k
195            long result = 1;
196            if (n <= 61) {
197                // For n <= 61, the naive implementation cannot overflow.
198                for (int j = 1, i = n - k + 1; j <= k; i++, j++) {
199                    result = result * i / j;
200                }
201            } else if (n <= 66) {
202                // For n > 61 but n <= 66, the result cannot overflow,
203                // but we must take care not to overflow intermediate values.
204                for (int j = 1, i = n - k + 1; j <= k; i++, j++) {
205                    // We know that (result * i) is divisible by j,
206                    // but (result * i) may overflow, so we split j:
207                    // Filter out the gcd, d, so j/d and i/d are integer.
208                    // result is divisible by (j/d) because (j/d)
209                    // is relative prime to (i/d) and is a divisor of
210                    // result * (i/d).
211                    long d = gcd(i, j);
212                    result = (result / (j / d)) * (i / d);
213                }
214            } else {
215                // For n > 66, a result overflow might occur, so we check
216                // the multiplication, taking care to not overflow
217                // unnecessary.
218                for (int j = 1, i = n - k + 1; j <= k; i++, j++) {
219                    long d = gcd(i, j);
220                    result = mulAndCheck((result / (j / d)), (i / d));
221                }
222            }
223            return result;
224        }
225    
226        /**
227         * Returns a <code>double</code> representation of the <a
228         * href="http://mathworld.wolfram.com/BinomialCoefficient.html"> Binomial
229         * Coefficient</a>, "<code>n choose k</code>", the number of
230         * <code>k</code>-element subsets that can be selected from an
231         * <code>n</code>-element set.
232         * <p>
233         * <Strong>Preconditions</strong>:
234         * <ul>
235         * <li> <code>0 <= k <= n </code> (otherwise
236         * <code>IllegalArgumentException</code> is thrown)</li>
237         * <li> The result is small enough to fit into a <code>double</code>. The
238         * largest value of <code>n</code> for which all coefficients are <
239         * Double.MAX_VALUE is 1029. If the computed value exceeds Double.MAX_VALUE,
240         * Double.POSITIVE_INFINITY is returned</li>
241         * </ul></p>
242         * 
243         * @param n the size of the set
244         * @param k the size of the subsets to be counted
245         * @return <code>n choose k</code>
246         * @throws IllegalArgumentException if preconditions are not met.
247         */
248        public static double binomialCoefficientDouble(final int n, final int k) {
249            checkBinomial(n, k);
250            if ((n == k) || (k == 0)) {
251                return 1d;
252            }
253            if ((k == 1) || (k == n - 1)) {
254                return n;
255            }
256            if (k > n/2) {
257                return binomialCoefficientDouble(n, n - k);
258            }
259            if (n < 67) {
260                return binomialCoefficient(n,k);
261            }
262            
263            double result = 1d;
264            for (int i = 1; i <= k; i++) {
265                 result *= (double)(n - k + i) / (double)i;
266            }
267      
268            return Math.floor(result + 0.5);
269        }
270        
271        /**
272         * Returns the natural <code>log</code> of the <a
273         * href="http://mathworld.wolfram.com/BinomialCoefficient.html"> Binomial
274         * Coefficient</a>, "<code>n choose k</code>", the number of
275         * <code>k</code>-element subsets that can be selected from an
276         * <code>n</code>-element set.
277         * <p>
278         * <Strong>Preconditions</strong>:
279         * <ul>
280         * <li> <code>0 <= k <= n </code> (otherwise
281         * <code>IllegalArgumentException</code> is thrown)</li>
282         * </ul></p>
283         * 
284         * @param n the size of the set
285         * @param k the size of the subsets to be counted
286         * @return <code>n choose k</code>
287         * @throws IllegalArgumentException if preconditions are not met.
288         */
289        public static double binomialCoefficientLog(final int n, final int k) {
290            checkBinomial(n, k);
291            if ((n == k) || (k == 0)) {
292                return 0;
293            }
294            if ((k == 1) || (k == n - 1)) {
295                return Math.log(n);
296            }
297            
298            /*
299             * For values small enough to do exact integer computation,
300             * return the log of the exact value 
301             */
302            if (n < 67) {  
303                return Math.log(binomialCoefficient(n,k));
304            }
305            
306            /*
307             * Return the log of binomialCoefficientDouble for values that will not
308             * overflow binomialCoefficientDouble
309             */
310            if (n < 1030) { 
311                return Math.log(binomialCoefficientDouble(n, k));
312            } 
313    
314            if (k > n / 2) {
315                return binomialCoefficientLog(n, n - k);
316            }
317    
318            /*
319             * Sum logs for values that could overflow
320             */
321            double logSum = 0;
322    
323            // n!/(n-k)!
324            for (int i = n - k + 1; i <= n; i++) {
325                logSum += Math.log(i);
326            }
327    
328            // divide by k!
329            for (int i = 2; i <= k; i++) {
330                logSum -= Math.log(i);
331            }
332    
333            return logSum;      
334        }
335    
336        /**
337         * Check binomial preconditions.
338         * @param n the size of the set
339         * @param k the size of the subsets to be counted
340         * @exception IllegalArgumentException if preconditions are not met.
341         */
342        private static void checkBinomial(final int n, final int k)
343            throws IllegalArgumentException {
344            if (n < k) {
345                throw MathRuntimeException.createIllegalArgumentException(
346                    "must have n >= k for binomial coefficient (n,k), got n = {0}, k = {1}",
347                    n, k);
348            }
349            if (n < 0) {
350                throw MathRuntimeException.createIllegalArgumentException(
351                      "must have n >= 0 for binomial coefficient (n,k), got n = {0}",
352                      n);
353            }
354        }
355        
356        /**
357         * Compares two numbers given some amount of allowed error.
358         * 
359         * @param x the first number
360         * @param y the second number
361         * @param eps the amount of error to allow when checking for equality
362         * @return <ul><li>0 if  {@link #equals(double, double, double) equals(x, y, eps)}</li>
363         *       <li>&lt; 0 if !{@link #equals(double, double, double) equals(x, y, eps)} &amp;&amp; x &lt; y</li>
364         *       <li>> 0 if !{@link #equals(double, double, double) equals(x, y, eps)} &amp;&amp; x > y</li></ul>
365         */
366        public static int compareTo(double x, double y, double eps) {
367            if (equals(x, y, eps)) {
368                return 0;
369            } else if (x < y) {
370              return -1;
371            }
372            return 1;
373        }
374        
375        /**
376         * Returns the <a href="http://mathworld.wolfram.com/HyperbolicCosine.html">
377         * hyperbolic cosine</a> of x.
378         * 
379         * @param x double value for which to find the hyperbolic cosine
380         * @return hyperbolic cosine of x
381         */
382        public static double cosh(double x) {
383            return (Math.exp(x) + Math.exp(-x)) / 2.0;
384        }
385        
386        /**
387         * Returns true iff both arguments are NaN or neither is NaN and they are
388         * equal
389         * 
390         * @param x first value
391         * @param y second value
392         * @return true if the values are equal or both are NaN
393         */
394        public static boolean equals(double x, double y) {
395            return ((Double.isNaN(x) && Double.isNaN(y)) || x == y);
396        }
397    
398        /**
399         * Returns true iff both arguments are equal or within the range of allowed
400         * error (inclusive).
401         * <p>
402         * Two NaNs are considered equals, as are two infinities with same sign.
403         * </p>
404         * 
405         * @param x first value
406         * @param y second value
407         * @param eps the amount of absolute error to allow
408         * @return true if the values are equal or within range of each other
409         */
410        public static boolean equals(double x, double y, double eps) {
411          return equals(x, y) || (Math.abs(y - x) <= eps);
412        }
413        
414        /**
415         * Returns true iff both arguments are equal or within the range of allowed
416         * error (inclusive).
417         * Adapted from <a
418         * href="http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm">
419         * Bruce Dawson</a>
420         *
421         * @param x first value
422         * @param y second value
423         * @param maxUlps {@code (maxUlps - 1)} is the number of floating point
424         * values between {@code x} and {@code y}.
425         * @return {@code true} if there are less than {@code maxUlps} floating
426         * point values between {@code x} and {@code y}
427         */
428        public static boolean equals(double x, double y, int maxUlps) {
429            // Check that "maxUlps" is non-negative and small enough so that the
430            // default NAN won't compare as equal to anything.
431            assert maxUlps > 0 && maxUlps < NAN_GAP;
432    
433            long xInt = Double.doubleToLongBits(x);
434            long yInt = Double.doubleToLongBits(y);
435    
436            // Make lexicographically ordered as a two's-complement integer.
437            if (xInt < 0) {
438                xInt = SGN_MASK - xInt;
439            }
440            if (yInt < 0) {
441                yInt = SGN_MASK - yInt;
442            }
443    
444            return Math.abs(xInt - yInt) <= maxUlps;
445        }
446    
447        /**
448         * Returns true iff both arguments are null or have same dimensions
449         * and all their elements are {@link #equals(double,double) equals}
450         * 
451         * @param x first array
452         * @param y second array
453         * @return true if the values are both null or have same dimension
454         * and equal elements
455         * @since 1.2
456         */
457        public static boolean equals(double[] x, double[] y) {
458            if ((x == null) || (y == null)) {
459                return !((x == null) ^ (y == null));
460            }
461            if (x.length != y.length) {
462                return false;
463            }
464            for (int i = 0; i < x.length; ++i) {
465                if (!equals(x[i], y[i])) {
466                    return false;
467                }
468            }
469            return true;
470        }
471        
472        /** All long-representable factorials */
473        private static final long[] factorials = new long[] 
474           {1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800, 39916800,
475            479001600, 6227020800l, 87178291200l, 1307674368000l, 20922789888000l,
476            355687428096000l, 6402373705728000l, 121645100408832000l,
477            2432902008176640000l};
478    
479        /**
480         * Returns n!. Shorthand for <code>n</code> <a
481         * href="http://mathworld.wolfram.com/Factorial.html"> Factorial</a>, the
482         * product of the numbers <code>1,...,n</code>.
483         * <p>
484         * <Strong>Preconditions</strong>:
485         * <ul>
486         * <li> <code>n >= 0</code> (otherwise
487         * <code>IllegalArgumentException</code> is thrown)</li>
488         * <li> The result is small enough to fit into a <code>long</code>. The
489         * largest value of <code>n</code> for which <code>n!</code> <
490         * Long.MAX_VALUE</code> is 20. If the computed value exceeds <code>Long.MAX_VALUE</code>
491         * an <code>ArithMeticException </code> is thrown.</li>
492         * </ul>
493         * </p>
494         * 
495         * @param n argument
496         * @return <code>n!</code>
497         * @throws ArithmeticException if the result is too large to be represented
498         *         by a long integer.
499         * @throws IllegalArgumentException if n < 0
500         */
501        public static long factorial(final int n) {
502            if (n < 0) {
503                throw MathRuntimeException.createIllegalArgumentException(
504                      "must have n >= 0 for n!, got n = {0}",
505                      n);
506            }
507            if (n > 20) {
508                throw new ArithmeticException(
509                        "factorial value is too large to fit in a long");
510            }
511            return factorials[n];
512        }
513    
514        /**
515         * Returns n!. Shorthand for <code>n</code> <a
516         * href="http://mathworld.wolfram.com/Factorial.html"> Factorial</a>, the
517         * product of the numbers <code>1,...,n</code> as a <code>double</code>.
518         * <p>
519         * <Strong>Preconditions</strong>:
520         * <ul>
521         * <li> <code>n >= 0</code> (otherwise
522         * <code>IllegalArgumentException</code> is thrown)</li>
523         * <li> The result is small enough to fit into a <code>double</code>. The
524         * largest value of <code>n</code> for which <code>n!</code> <
525         * Double.MAX_VALUE</code> is 170. If the computed value exceeds
526         * Double.MAX_VALUE, Double.POSITIVE_INFINITY is returned</li>
527         * </ul>
528         * </p>
529         * 
530         * @param n argument
531         * @return <code>n!</code>
532         * @throws IllegalArgumentException if n < 0
533         */
534        public static double factorialDouble(final int n) {
535            if (n < 0) {
536                throw MathRuntimeException.createIllegalArgumentException(
537                      "must have n >= 0 for n!, got n = {0}",
538                      n);
539            }
540            if (n < 21) {
541                return factorial(n);
542            }
543            return Math.floor(Math.exp(factorialLog(n)) + 0.5);
544        }
545    
546        /**
547         * Returns the natural logarithm of n!.
548         * <p>
549         * <Strong>Preconditions</strong>:
550         * <ul>
551         * <li> <code>n >= 0</code> (otherwise
552         * <code>IllegalArgumentException</code> is thrown)</li>
553         * </ul></p>
554         * 
555         * @param n argument
556         * @return <code>n!</code>
557         * @throws IllegalArgumentException if preconditions are not met.
558         */
559        public static double factorialLog(final int n) {
560            if (n < 0) {
561                throw MathRuntimeException.createIllegalArgumentException(
562                      "must have n >= 0 for n!, got n = {0}",
563                      n);
564            }
565            if (n < 21) {
566                return Math.log(factorial(n));
567            }
568            double logSum = 0;
569            for (int i = 2; i <= n; i++) {
570                logSum += Math.log(i);
571            }
572            return logSum;
573        }
574    
575        /**
576         * <p>
577         * Gets the greatest common divisor of the absolute value of two numbers,
578         * using the "binary gcd" method which avoids division and modulo
579         * operations. See Knuth 4.5.2 algorithm B. This algorithm is due to Josef
580         * Stein (1961).
581         * </p>
582         * Special cases:
583         * <ul>
584         * <li>The invocations
585         * <code>gcd(Integer.MIN_VALUE, Integer.MIN_VALUE)</code>,
586         * <code>gcd(Integer.MIN_VALUE, 0)</code> and
587         * <code>gcd(0, Integer.MIN_VALUE)</code> throw an
588         * <code>ArithmeticException</code>, because the result would be 2^31, which
589         * is too large for an int value.</li>
590         * <li>The result of <code>gcd(x, x)</code>, <code>gcd(0, x)</code> and
591         * <code>gcd(x, 0)</code> is the absolute value of <code>x</code>, except
592         * for the special cases above.
593         * <li>The invocation <code>gcd(0, 0)</code> is the only one which returns
594         * <code>0</code>.</li>
595         * </ul>
596         * 
597         * @param p any number
598         * @param q any number
599         * @return the greatest common divisor, never negative
600         * @throws ArithmeticException
601         *             if the result cannot be represented as a nonnegative int
602         *             value
603         * @since 1.1
604         */
605        public static int gcd(final int p, final int q) {
606            int u = p;
607            int v = q;
608            if ((u == 0) || (v == 0)) {
609                if ((u == Integer.MIN_VALUE) || (v == Integer.MIN_VALUE)) {
610                    throw MathRuntimeException.createArithmeticException(
611                            "overflow: gcd({0}, {1}) is 2^31",
612                            p, q);
613                }
614                return (Math.abs(u) + Math.abs(v));
615            }
616            // keep u and v negative, as negative integers range down to
617            // -2^31, while positive numbers can only be as large as 2^31-1
618            // (i.e. we can't necessarily negate a negative number without
619            // overflow)
620            /* assert u!=0 && v!=0; */
621            if (u > 0) {
622                u = -u;
623            } // make u negative
624            if (v > 0) {
625                v = -v;
626            } // make v negative
627            // B1. [Find power of 2]
628            int k = 0;
629            while ((u & 1) == 0 && (v & 1) == 0 && k < 31) { // while u and v are
630                                                                // both even...
631                u /= 2;
632                v /= 2;
633                k++; // cast out twos.
634            }
635            if (k == 31) {
636                throw MathRuntimeException.createArithmeticException(
637                        "overflow: gcd({0}, {1}) is 2^31",
638                        p, q);
639            }
640            // B2. Initialize: u and v have been divided by 2^k and at least
641            // one is odd.
642            int t = ((u & 1) == 1) ? v : -(u / 2)/* B3 */;
643            // t negative: u was odd, v may be even (t replaces v)
644            // t positive: u was even, v is odd (t replaces u)
645            do {
646                /* assert u<0 && v<0; */
647                // B4/B3: cast out twos from t.
648                while ((t & 1) == 0) { // while t is even..
649                    t /= 2; // cast out twos
650                }
651                // B5 [reset max(u,v)]
652                if (t > 0) {
653                    u = -t;
654                } else {
655                    v = t;
656                }
657                // B6/B3. at this point both u and v should be odd.
658                t = (v - u) / 2;
659                // |u| larger: t positive (replace u)
660                // |v| larger: t negative (replace v)
661            } while (t != 0);
662            return -u * (1 << k); // gcd is u*2^k
663        }
664    
665        /**
666         * Returns an integer hash code representing the given double value.
667         * 
668         * @param value the value to be hashed
669         * @return the hash code
670         */
671        public static int hash(double value) {
672            return new Double(value).hashCode();
673        }
674    
675        /**
676         * Returns an integer hash code representing the given double array.
677         * 
678         * @param value the value to be hashed (may be null)
679         * @return the hash code
680         * @since 1.2
681         */
682        public static int hash(double[] value) {
683            return Arrays.hashCode(value);
684        }
685    
686        /**
687         * For a byte value x, this method returns (byte)(+1) if x >= 0 and
688         * (byte)(-1) if x < 0.
689         * 
690         * @param x the value, a byte
691         * @return (byte)(+1) or (byte)(-1), depending on the sign of x
692         */
693        public static byte indicator(final byte x) {
694            return (x >= ZB) ? PB : NB;
695        }
696    
697        /**
698         * For a double precision value x, this method returns +1.0 if x >= 0 and
699         * -1.0 if x < 0. Returns <code>NaN</code> if <code>x</code> is
700         * <code>NaN</code>.
701         * 
702         * @param x the value, a double
703         * @return +1.0 or -1.0, depending on the sign of x
704         */
705        public static double indicator(final double x) {
706            if (Double.isNaN(x)) {
707                return Double.NaN;
708            }
709            return (x >= 0.0) ? 1.0 : -1.0;
710        }
711    
712        /**
713         * For a float value x, this method returns +1.0F if x >= 0 and -1.0F if x <
714         * 0. Returns <code>NaN</code> if <code>x</code> is <code>NaN</code>.
715         * 
716         * @param x the value, a float
717         * @return +1.0F or -1.0F, depending on the sign of x
718         */
719        public static float indicator(final float x) {
720            if (Float.isNaN(x)) {
721                return Float.NaN;
722            }
723            return (x >= 0.0F) ? 1.0F : -1.0F;
724        }
725    
726        /**
727         * For an int value x, this method returns +1 if x >= 0 and -1 if x < 0.
728         * 
729         * @param x the value, an int
730         * @return +1 or -1, depending on the sign of x
731         */
732        public static int indicator(final int x) {
733            return (x >= 0) ? 1 : -1;
734        }
735    
736        /**
737         * For a long value x, this method returns +1L if x >= 0 and -1L if x < 0.
738         * 
739         * @param x the value, a long
740         * @return +1L or -1L, depending on the sign of x
741         */
742        public static long indicator(final long x) {
743            return (x >= 0L) ? 1L : -1L;
744        }
745    
746        /**
747         * For a short value x, this method returns (short)(+1) if x >= 0 and
748         * (short)(-1) if x < 0.
749         * 
750         * @param x the value, a short
751         * @return (short)(+1) or (short)(-1), depending on the sign of x
752         */
753        public static short indicator(final short x) {
754            return (x >= ZS) ? PS : NS;
755        }
756    
757        /**
758         * <p>
759         * Returns the least common multiple of the absolute value of two numbers,
760         * using the formula <code>lcm(a,b) = (a / gcd(a,b)) * b</code>.
761         * </p>
762         * Special cases:
763         * <ul>
764         * <li>The invocations <code>lcm(Integer.MIN_VALUE, n)</code> and
765         * <code>lcm(n, Integer.MIN_VALUE)</code>, where <code>abs(n)</code> is a
766         * power of 2, throw an <code>ArithmeticException</code>, because the result
767         * would be 2^31, which is too large for an int value.</li>
768         * <li>The result of <code>lcm(0, x)</code> and <code>lcm(x, 0)</code> is
769         * <code>0</code> for any <code>x</code>.
770         * </ul>
771         * 
772         * @param a any number
773         * @param b any number
774         * @return the least common multiple, never negative
775         * @throws ArithmeticException
776         *             if the result cannot be represented as a nonnegative int
777         *             value
778         * @since 1.1
779         */
780        public static int lcm(int a, int b) {
781            if (a==0 || b==0){
782                return 0;
783            }
784            int lcm = Math.abs(mulAndCheck(a / gcd(a, b), b));
785            if (lcm == Integer.MIN_VALUE){
786                throw new ArithmeticException("overflow: lcm is 2^31");
787            }
788            return lcm;
789        }
790    
791        /** 
792         * <p>Returns the 
793         * <a href="http://mathworld.wolfram.com/Logarithm.html">logarithm</a>
794         * for base <code>b</code> of <code>x</code>.
795         * </p>
796         * <p>Returns <code>NaN<code> if either argument is negative.  If 
797         * <code>base</code> is 0 and <code>x</code> is positive, 0 is returned.
798         * If <code>base</code> is positive and <code>x</code> is 0, 
799         * <code>Double.NEGATIVE_INFINITY</code> is returned.  If both arguments
800         * are 0, the result is <code>NaN</code>.</p>
801         * 
802         * @param base the base of the logarithm, must be greater than 0
803         * @param x argument, must be greater than 0
804         * @return the value of the logarithm - the number y such that base^y = x.
805         * @since 1.2
806         */ 
807        public static double log(double base, double x) {
808            return Math.log(x)/Math.log(base);
809        }
810    
811        /**
812         * Multiply two integers, checking for overflow.
813         * 
814         * @param x a factor
815         * @param y a factor
816         * @return the product <code>x*y</code>
817         * @throws ArithmeticException if the result can not be represented as an
818         *         int
819         * @since 1.1
820         */
821        public static int mulAndCheck(int x, int y) {
822            long m = ((long)x) * ((long)y);
823            if (m < Integer.MIN_VALUE || m > Integer.MAX_VALUE) {
824                throw new ArithmeticException("overflow: mul");
825            }
826            return (int)m;
827        }
828    
829        /**
830         * Multiply two long integers, checking for overflow.
831         * 
832         * @param a first value
833         * @param b second value
834         * @return the product <code>a * b</code>
835         * @throws ArithmeticException if the result can not be represented as an
836         *         long
837         * @since 1.2
838         */
839        public static long mulAndCheck(long a, long b) {
840            long ret;
841            String msg = "overflow: multiply";
842            if (a > b) {
843                // use symmetry to reduce boundary cases
844                ret = mulAndCheck(b, a);
845            } else {
846                if (a < 0) {
847                    if (b < 0) {
848                        // check for positive overflow with negative a, negative b
849                        if (a >= Long.MAX_VALUE / b) {
850                            ret = a * b;
851                        } else {
852                            throw new ArithmeticException(msg);
853                        }
854                    } else if (b > 0) {
855                        // check for negative overflow with negative a, positive b
856                        if (Long.MIN_VALUE / b <= a) {
857                            ret = a * b;
858                        } else {
859                            throw new ArithmeticException(msg);
860                            
861                        }
862                    } else {
863                        // assert b == 0
864                        ret = 0;
865                    }
866                } else if (a > 0) {
867                    // assert a > 0
868                    // assert b > 0
869                    
870                    // check for positive overflow with positive a, positive b
871                    if (a <= Long.MAX_VALUE / b) {
872                        ret = a * b;
873                    } else {
874                        throw new ArithmeticException(msg);
875                    }
876                } else {
877                    // assert a == 0
878                    ret = 0;
879                }
880            }
881            return ret;
882        }
883    
884        /**
885         * Get the next machine representable number after a number, moving
886         * in the direction of another number.
887         * <p>
888         * If <code>direction</code> is greater than or equal to<code>d</code>,
889         * the smallest machine representable number strictly greater than
890         * <code>d</code> is returned; otherwise the largest representable number
891         * strictly less than <code>d</code> is returned.</p>
892         * <p>
893         * If <code>d</code> is NaN or Infinite, it is returned unchanged.</p>
894         * 
895         * @param d base number
896         * @param direction (the only important thing is whether
897         * direction is greater or smaller than d)
898         * @return the next machine representable number in the specified direction
899         * @since 1.2
900         */
901        public static double nextAfter(double d, double direction) {
902    
903            // handling of some important special cases
904            if (Double.isNaN(d) || Double.isInfinite(d)) {
905                    return d;
906            } else if (d == 0) {
907                    return (direction < 0) ? -Double.MIN_VALUE : Double.MIN_VALUE;
908            }
909            // special cases MAX_VALUE to infinity and  MIN_VALUE to 0
910            // are handled just as normal numbers
911    
912            // split the double in raw components
913            long bits     = Double.doubleToLongBits(d);
914            long sign     = bits & 0x8000000000000000L;
915            long exponent = bits & 0x7ff0000000000000L;
916            long mantissa = bits & 0x000fffffffffffffL;
917    
918            if (d * (direction - d) >= 0) {
919                    // we should increase the mantissa
920                    if (mantissa == 0x000fffffffffffffL) {
921                            return Double.longBitsToDouble(sign |
922                                            (exponent + 0x0010000000000000L));
923                    } else {
924                            return Double.longBitsToDouble(sign |
925                                            exponent | (mantissa + 1));
926                    }
927            } else {
928                    // we should decrease the mantissa
929                    if (mantissa == 0L) {
930                            return Double.longBitsToDouble(sign |
931                                            (exponent - 0x0010000000000000L) |
932                                            0x000fffffffffffffL);
933                    } else {
934                            return Double.longBitsToDouble(sign |
935                                            exponent | (mantissa - 1));
936                    }
937            }
938    
939        }
940    
941        /**
942         * Scale a number by 2<sup>scaleFactor</sup>.
943         * <p>If <code>d</code> is 0 or NaN or Infinite, it is returned unchanged.</p>
944         * 
945         * @param d base number
946         * @param scaleFactor power of two by which d sould be multiplied
947         * @return d &times; 2<sup>scaleFactor</sup>
948         * @since 2.0
949         */
950        public static double scalb(final double d, final int scaleFactor) {
951    
952            // handling of some important special cases
953            if ((d == 0) || Double.isNaN(d) || Double.isInfinite(d)) {
954                return d;
955            }
956    
957            // split the double in raw components
958            final long bits     = Double.doubleToLongBits(d);
959            final long exponent = bits & 0x7ff0000000000000L;
960            final long rest     = bits & 0x800fffffffffffffL;
961    
962            // shift the exponent
963            final long newBits = rest | (exponent + (((long) scaleFactor) << 52));
964            return Double.longBitsToDouble(newBits);
965    
966        }
967    
968        /**
969         * Normalize an angle in a 2&pi wide interval around a center value.
970         * <p>This method has three main uses:</p>
971         * <ul>
972         *   <li>normalize an angle between 0 and 2&pi;:<br/>
973         *       <code>a = MathUtils.normalizeAngle(a, Math.PI);</code></li>
974         *   <li>normalize an angle between -&pi; and +&pi;<br/>
975         *       <code>a = MathUtils.normalizeAngle(a, 0.0);</code></li>
976         *   <li>compute the angle between two defining angular positions:<br>
977         *       <code>angle = MathUtils.normalizeAngle(end, start) - start;</code></li>
978         * </ul>
979         * <p>Note that due to numerical accuracy and since &pi; cannot be represented
980         * exactly, the result interval is <em>closed</em>, it cannot be half-closed
981         * as would be more satisfactory in a purely mathematical view.</p>
982         * @param a angle to normalize
983         * @param center center of the desired 2&pi; interval for the result
984         * @return a-2k&pi; with integer k and center-&pi; &lt;= a-2k&pi; &lt;= center+&pi;
985         * @since 1.2
986         */
987         public static double normalizeAngle(double a, double center) {
988             return a - TWO_PI * Math.floor((a + Math.PI - center) / TWO_PI);
989         }
990    
991        /**
992         * Round the given value to the specified number of decimal places. The
993         * value is rounded using the {@link BigDecimal#ROUND_HALF_UP} method.
994         * 
995         * @param x the value to round.
996         * @param scale the number of digits to the right of the decimal point.
997         * @return the rounded value.
998         * @since 1.1
999         */
1000        public static double round(double x, int scale) {
1001            return round(x, scale, BigDecimal.ROUND_HALF_UP);
1002        }
1003    
1004        /**
1005         * Round the given value to the specified number of decimal places. The
1006         * value is rounded using the given method which is any method defined in
1007         * {@link BigDecimal}.
1008         * 
1009         * @param x the value to round.
1010         * @param scale the number of digits to the right of the decimal point.
1011         * @param roundingMethod the rounding method as defined in
1012         *        {@link BigDecimal}.
1013         * @return the rounded value.
1014         * @since 1.1
1015         */
1016        public static double round(double x, int scale, int roundingMethod) {
1017            try {
1018                return (new BigDecimal
1019                       (Double.toString(x))
1020                       .setScale(scale, roundingMethod))
1021                       .doubleValue();
1022            } catch (NumberFormatException ex) {
1023                if (Double.isInfinite(x)) {
1024                    return x;          
1025                } else {
1026                    return Double.NaN;
1027                }
1028            }
1029        }
1030    
1031        /**
1032         * Round the given value to the specified number of decimal places. The
1033         * value is rounding using the {@link BigDecimal#ROUND_HALF_UP} method.
1034         * 
1035         * @param x the value to round.
1036         * @param scale the number of digits to the right of the decimal point.
1037         * @return the rounded value.
1038         * @since 1.1
1039         */
1040        public static float round(float x, int scale) {
1041            return round(x, scale, BigDecimal.ROUND_HALF_UP);
1042        }
1043    
1044        /**
1045         * Round the given value to the specified number of decimal places. The
1046         * value is rounded using the given method which is any method defined in
1047         * {@link BigDecimal}.
1048         * 
1049         * @param x the value to round.
1050         * @param scale the number of digits to the right of the decimal point.
1051         * @param roundingMethod the rounding method as defined in
1052         *        {@link BigDecimal}.
1053         * @return the rounded value.
1054         * @since 1.1
1055         */
1056        public static float round(float x, int scale, int roundingMethod) {
1057            float sign = indicator(x);
1058            float factor = (float)Math.pow(10.0f, scale) * sign;
1059            return (float)roundUnscaled(x * factor, sign, roundingMethod) / factor;
1060        }
1061    
1062        /**
1063         * Round the given non-negative, value to the "nearest" integer. Nearest is
1064         * determined by the rounding method specified. Rounding methods are defined
1065         * in {@link BigDecimal}.
1066         * 
1067         * @param unscaled the value to round.
1068         * @param sign the sign of the original, scaled value.
1069         * @param roundingMethod the rounding method as defined in
1070         *        {@link BigDecimal}.
1071         * @return the rounded value.
1072         * @since 1.1
1073         */
1074        private static double roundUnscaled(double unscaled, double sign,
1075            int roundingMethod) {
1076            switch (roundingMethod) {
1077            case BigDecimal.ROUND_CEILING :
1078                if (sign == -1) {
1079                    unscaled = Math.floor(nextAfter(unscaled, Double.NEGATIVE_INFINITY));
1080                } else {
1081                    unscaled = Math.ceil(nextAfter(unscaled, Double.POSITIVE_INFINITY));
1082                }
1083                break;
1084            case BigDecimal.ROUND_DOWN :
1085                unscaled = Math.floor(nextAfter(unscaled, Double.NEGATIVE_INFINITY));
1086                break;
1087            case BigDecimal.ROUND_FLOOR :
1088                if (sign == -1) {
1089                    unscaled = Math.ceil(nextAfter(unscaled, Double.POSITIVE_INFINITY));
1090                } else {
1091                    unscaled = Math.floor(nextAfter(unscaled, Double.NEGATIVE_INFINITY));
1092                }
1093                break;
1094            case BigDecimal.ROUND_HALF_DOWN : {
1095                unscaled = nextAfter(unscaled, Double.NEGATIVE_INFINITY);
1096                double fraction = unscaled - Math.floor(unscaled);
1097                if (fraction > 0.5) {
1098                    unscaled = Math.ceil(unscaled);
1099                } else {
1100                    unscaled = Math.floor(unscaled);
1101                }
1102                break;
1103            }
1104            case BigDecimal.ROUND_HALF_EVEN : {
1105                double fraction = unscaled - Math.floor(unscaled);
1106                if (fraction > 0.5) {
1107                    unscaled = Math.ceil(unscaled);
1108                } else if (fraction < 0.5) {
1109                    unscaled = Math.floor(unscaled);
1110                } else {
1111                    // The following equality test is intentional and needed for rounding purposes
1112                    if (Math.floor(unscaled) / 2.0 == Math.floor(Math
1113                        .floor(unscaled) / 2.0)) { // even
1114                        unscaled = Math.floor(unscaled);
1115                    } else { // odd
1116                        unscaled = Math.ceil(unscaled);
1117                    }
1118                }
1119                break;
1120            }
1121            case BigDecimal.ROUND_HALF_UP : {
1122                unscaled = nextAfter(unscaled, Double.POSITIVE_INFINITY);
1123                double fraction = unscaled - Math.floor(unscaled);
1124                if (fraction >= 0.5) {
1125                    unscaled = Math.ceil(unscaled);
1126                } else {
1127                    unscaled = Math.floor(unscaled);
1128                }
1129                break;
1130            }
1131            case BigDecimal.ROUND_UNNECESSARY :
1132                if (unscaled != Math.floor(unscaled)) {
1133                    throw new ArithmeticException("Inexact result from rounding");
1134                }
1135                break;
1136            case BigDecimal.ROUND_UP :
1137                unscaled = Math.ceil(nextAfter(unscaled,  Double.POSITIVE_INFINITY));
1138                break;
1139            default :
1140                throw MathRuntimeException.createIllegalArgumentException(
1141                      "invalid rounding method {0}, valid methods: {1} ({2}), {3} ({4})," +
1142                      " {5} ({6}), {7} ({8}), {9} ({10}), {11} ({12}), {13} ({14}), {15} ({16})",
1143                      roundingMethod,
1144                      "ROUND_CEILING",     BigDecimal.ROUND_CEILING,
1145                      "ROUND_DOWN",        BigDecimal.ROUND_DOWN,
1146                      "ROUND_FLOOR",       BigDecimal.ROUND_FLOOR,
1147                      "ROUND_HALF_DOWN",   BigDecimal.ROUND_HALF_DOWN,
1148                      "ROUND_HALF_EVEN",   BigDecimal.ROUND_HALF_EVEN,
1149                      "ROUND_HALF_UP",     BigDecimal.ROUND_HALF_UP,
1150                      "ROUND_UNNECESSARY", BigDecimal.ROUND_UNNECESSARY,
1151                      "ROUND_UP",          BigDecimal.ROUND_UP);
1152            }
1153            return unscaled;
1154        }
1155    
1156        /**
1157         * Returns the <a href="http://mathworld.wolfram.com/Sign.html"> sign</a>
1158         * for byte value <code>x</code>.
1159         * <p>
1160         * For a byte value x, this method returns (byte)(+1) if x > 0, (byte)(0) if
1161         * x = 0, and (byte)(-1) if x < 0.</p>
1162         * 
1163         * @param x the value, a byte
1164         * @return (byte)(+1), (byte)(0), or (byte)(-1), depending on the sign of x
1165         */
1166        public static byte sign(final byte x) {
1167            return (x == ZB) ? ZB : (x > ZB) ? PB : NB;
1168        }
1169    
1170        /**
1171         * Returns the <a href="http://mathworld.wolfram.com/Sign.html"> sign</a>
1172         * for double precision <code>x</code>.
1173         * <p>
1174         * For a double value <code>x</code>, this method returns
1175         * <code>+1.0</code> if <code>x > 0</code>, <code>0.0</code> if
1176         * <code>x = 0.0</code>, and <code>-1.0</code> if <code>x < 0</code>.
1177         * Returns <code>NaN</code> if <code>x</code> is <code>NaN</code>.</p>
1178         * 
1179         * @param x the value, a double
1180         * @return +1.0, 0.0, or -1.0, depending on the sign of x
1181         */
1182        public static double sign(final double x) {
1183            if (Double.isNaN(x)) {
1184                return Double.NaN;
1185            }
1186            return (x == 0.0) ? 0.0 : (x > 0.0) ? 1.0 : -1.0;
1187        }
1188    
1189        /**
1190         * Returns the <a href="http://mathworld.wolfram.com/Sign.html"> sign</a>
1191         * for float value <code>x</code>.
1192         * <p>
1193         * For a float value x, this method returns +1.0F if x > 0, 0.0F if x =
1194         * 0.0F, and -1.0F if x < 0. Returns <code>NaN</code> if <code>x</code>
1195         * is <code>NaN</code>.</p>
1196         * 
1197         * @param x the value, a float
1198         * @return +1.0F, 0.0F, or -1.0F, depending on the sign of x
1199         */
1200        public static float sign(final float x) {
1201            if (Float.isNaN(x)) {
1202                return Float.NaN;
1203            }
1204            return (x == 0.0F) ? 0.0F : (x > 0.0F) ? 1.0F : -1.0F;
1205        }
1206    
1207        /**
1208         * Returns the <a href="http://mathworld.wolfram.com/Sign.html"> sign</a>
1209         * for int value <code>x</code>.
1210         * <p>
1211         * For an int value x, this method returns +1 if x > 0, 0 if x = 0, and -1
1212         * if x < 0.</p>
1213         * 
1214         * @param x the value, an int
1215         * @return +1, 0, or -1, depending on the sign of x
1216         */
1217        public static int sign(final int x) {
1218            return (x == 0) ? 0 : (x > 0) ? 1 : -1;
1219        }
1220    
1221        /**
1222         * Returns the <a href="http://mathworld.wolfram.com/Sign.html"> sign</a>
1223         * for long value <code>x</code>.
1224         * <p>
1225         * For a long value x, this method returns +1L if x > 0, 0L if x = 0, and
1226         * -1L if x < 0.</p>
1227         * 
1228         * @param x the value, a long
1229         * @return +1L, 0L, or -1L, depending on the sign of x
1230         */
1231        public static long sign(final long x) {
1232            return (x == 0L) ? 0L : (x > 0L) ? 1L : -1L;
1233        }
1234    
1235        /**
1236         * Returns the <a href="http://mathworld.wolfram.com/Sign.html"> sign</a>
1237         * for short value <code>x</code>.
1238         * <p>
1239         * For a short value x, this method returns (short)(+1) if x > 0, (short)(0)
1240         * if x = 0, and (short)(-1) if x < 0.</p>
1241         * 
1242         * @param x the value, a short
1243         * @return (short)(+1), (short)(0), or (short)(-1), depending on the sign of
1244         *         x
1245         */
1246        public static short sign(final short x) {
1247            return (x == ZS) ? ZS : (x > ZS) ? PS : NS;
1248        }
1249    
1250        /**
1251         * Returns the <a href="http://mathworld.wolfram.com/HyperbolicSine.html">
1252         * hyperbolic sine</a> of x.
1253         * 
1254         * @param x double value for which to find the hyperbolic sine
1255         * @return hyperbolic sine of x
1256         */
1257        public static double sinh(double x) {
1258            return (Math.exp(x) - Math.exp(-x)) / 2.0;
1259        }
1260    
1261        /**
1262         * Subtract two integers, checking for overflow.
1263         * 
1264         * @param x the minuend
1265         * @param y the subtrahend
1266         * @return the difference <code>x-y</code>
1267         * @throws ArithmeticException if the result can not be represented as an
1268         *         int
1269         * @since 1.1
1270         */
1271        public static int subAndCheck(int x, int y) {
1272            long s = (long)x - (long)y;
1273            if (s < Integer.MIN_VALUE || s > Integer.MAX_VALUE) {
1274                throw new ArithmeticException("overflow: subtract");
1275            }
1276            return (int)s;
1277        }
1278    
1279        /**
1280         * Subtract two long integers, checking for overflow.
1281         * 
1282         * @param a first value
1283         * @param b second value
1284         * @return the difference <code>a-b</code>
1285         * @throws ArithmeticException if the result can not be represented as an
1286         *         long
1287         * @since 1.2
1288         */
1289        public static long subAndCheck(long a, long b) {
1290            long ret;
1291            String msg = "overflow: subtract";
1292            if (b == Long.MIN_VALUE) {
1293                if (a < 0) {
1294                    ret = a - b;
1295                } else {
1296                    throw new ArithmeticException(msg);
1297                }
1298            } else {
1299                // use additive inverse
1300                ret = addAndCheck(a, -b, msg);
1301            }
1302            return ret;
1303        }
1304    
1305        /**
1306         * Raise an int to an int power.
1307         * @param k number to raise
1308         * @param e exponent (must be positive or null)
1309         * @return k<sup>e</sup>
1310         * @exception IllegalArgumentException if e is negative
1311         */
1312        public static int pow(final int k, int e)
1313            throws IllegalArgumentException {
1314    
1315            if (e < 0) {
1316                throw MathRuntimeException.createIllegalArgumentException(
1317                    "cannot raise an integral value to a negative power ({0}^{1})",
1318                    k, e);
1319            }
1320    
1321            int result = 1;
1322            int k2p    = k;
1323            while (e != 0) {
1324                if ((e & 0x1) != 0) {
1325                    result *= k2p;
1326                }
1327                k2p *= k2p;
1328                e = e >> 1;
1329            }
1330    
1331            return result;
1332    
1333        }
1334    
1335        /**
1336         * Raise an int to a long power.
1337         * @param k number to raise
1338         * @param e exponent (must be positive or null)
1339         * @return k<sup>e</sup>
1340         * @exception IllegalArgumentException if e is negative
1341         */
1342        public static int pow(final int k, long e)
1343            throws IllegalArgumentException {
1344    
1345            if (e < 0) {
1346                throw MathRuntimeException.createIllegalArgumentException(
1347                    "cannot raise an integral value to a negative power ({0}^{1})",
1348                    k, e);
1349            }
1350    
1351            int result = 1;
1352            int k2p    = k;
1353            while (e != 0) {
1354                if ((e & 0x1) != 0) {
1355                    result *= k2p;
1356                }
1357                k2p *= k2p;
1358                e = e >> 1;
1359            }
1360    
1361            return result;
1362    
1363        }
1364    
1365        /**
1366         * Raise a long to an int power.
1367         * @param k number to raise
1368         * @param e exponent (must be positive or null)
1369         * @return k<sup>e</sup>
1370         * @exception IllegalArgumentException if e is negative
1371         */
1372        public static long pow(final long k, int e)
1373            throws IllegalArgumentException {
1374    
1375            if (e < 0) {
1376                throw MathRuntimeException.createIllegalArgumentException(
1377                    "cannot raise an integral value to a negative power ({0}^{1})",
1378                    k, e);
1379            }
1380    
1381            long result = 1l;
1382            long k2p    = k;
1383            while (e != 0) {
1384                if ((e & 0x1) != 0) {
1385                    result *= k2p;
1386                }
1387                k2p *= k2p;
1388                e = e >> 1;
1389            }
1390    
1391            return result;
1392    
1393        }
1394    
1395        /**
1396         * Raise a long to a long power.
1397         * @param k number to raise
1398         * @param e exponent (must be positive or null)
1399         * @return k<sup>e</sup>
1400         * @exception IllegalArgumentException if e is negative
1401         */
1402        public static long pow(final long k, long e)
1403            throws IllegalArgumentException {
1404    
1405            if (e < 0) {
1406                throw MathRuntimeException.createIllegalArgumentException(
1407                    "cannot raise an integral value to a negative power ({0}^{1})",
1408                    k, e);
1409            }
1410    
1411            long result = 1l;
1412            long k2p    = k;
1413            while (e != 0) {
1414                if ((e & 0x1) != 0) {
1415                    result *= k2p;
1416                }
1417                k2p *= k2p;
1418                e = e >> 1;
1419            }
1420    
1421            return result;
1422    
1423        }
1424    
1425        /**
1426         * Raise a BigInteger to an int power.
1427         * @param k number to raise
1428         * @param e exponent (must be positive or null)
1429         * @return k<sup>e</sup>
1430         * @exception IllegalArgumentException if e is negative
1431         */
1432        public static BigInteger pow(final BigInteger k, int e)
1433            throws IllegalArgumentException {
1434    
1435            if (e < 0) {
1436                throw MathRuntimeException.createIllegalArgumentException(
1437                    "cannot raise an integral value to a negative power ({0}^{1})",
1438                    k, e);
1439            }
1440    
1441            return k.pow(e);
1442    
1443        }
1444    
1445        /**
1446         * Raise a BigInteger to a long power.
1447         * @param k number to raise
1448         * @param e exponent (must be positive or null)
1449         * @return k<sup>e</sup>
1450         * @exception IllegalArgumentException if e is negative
1451         */
1452        public static BigInteger pow(final BigInteger k, long e)
1453            throws IllegalArgumentException {
1454    
1455            if (e < 0) {
1456                throw MathRuntimeException.createIllegalArgumentException(
1457                    "cannot raise an integral value to a negative power ({0}^{1})",
1458                    k, e);
1459            }
1460    
1461            BigInteger result = BigInteger.ONE;
1462            BigInteger k2p    = k;
1463            while (e != 0) {
1464                if ((e & 0x1) != 0) {
1465                    result = result.multiply(k2p);
1466                }
1467                k2p = k2p.multiply(k2p);
1468                e = e >> 1;
1469            }
1470    
1471            return result;
1472    
1473        }
1474    
1475        /**
1476         * Raise a BigInteger to a BigInteger power.
1477         * @param k number to raise
1478         * @param e exponent (must be positive or null)
1479         * @return k<sup>e</sup>
1480         * @exception IllegalArgumentException if e is negative
1481         */
1482        public static BigInteger pow(final BigInteger k, BigInteger e)
1483            throws IllegalArgumentException {
1484    
1485            if (e.compareTo(BigInteger.ZERO) < 0) {
1486                throw MathRuntimeException.createIllegalArgumentException(
1487                    "cannot raise an integral value to a negative power ({0}^{1})",
1488                    k, e);
1489            }
1490    
1491            BigInteger result = BigInteger.ONE;
1492            BigInteger k2p    = k;
1493            while (!BigInteger.ZERO.equals(e)) {
1494                if (e.testBit(0)) {
1495                    result = result.multiply(k2p);
1496                }
1497                k2p = k2p.multiply(k2p);
1498                e = e.shiftRight(1);
1499            }
1500    
1501            return result;
1502    
1503        }
1504    
1505        /**
1506         * Calculates the L<sub>1</sub> (sum of abs) distance between two points.
1507         *
1508         * @param p1 the first point
1509         * @param p2 the second point
1510         * @return the L<sub>1</sub> distance between the two points
1511         */
1512        public static final double distance1(double[] p1, double[] p2) {
1513            double sum = 0;
1514            for (int i = 0; i < p1.length; i++) {
1515                sum += Math.abs(p1[i] - p2[i]);
1516            }
1517            return sum;
1518        }
1519        
1520        /**
1521         * Calculates the L<sub>1</sub> (sum of abs) distance between two points.
1522         *
1523         * @param p1 the first point
1524         * @param p2 the second point
1525         * @return the L<sub>1</sub> distance between the two points
1526         */
1527        public static final int distance1(int[] p1, int[] p2) {
1528          int sum = 0;
1529          for (int i = 0; i < p1.length; i++) {
1530              sum += Math.abs(p1[i] - p2[i]);
1531          }
1532          return sum;
1533        }
1534    
1535        /**
1536         * Calculates the L<sub>2</sub> (Euclidean) distance between two points.
1537         *
1538         * @param p1 the first point
1539         * @param p2 the second point
1540         * @return the L<sub>2</sub> distance between the two points
1541         */
1542        public static final double distance(double[] p1, double[] p2) {
1543            double sum = 0;
1544            for (int i = 0; i < p1.length; i++) {
1545                final double dp = p1[i] - p2[i];
1546                sum += dp * dp;
1547            }
1548            return Math.sqrt(sum);
1549        }
1550        
1551        /**
1552         * Calculates the L<sub>2</sub> (Euclidean) distance between two points.
1553         *
1554         * @param p1 the first point
1555         * @param p2 the second point
1556         * @return the L<sub>2</sub> distance between the two points
1557         */
1558        public static final double distance(int[] p1, int[] p2) {
1559          int sum = 0;
1560          for (int i = 0; i < p1.length; i++) {
1561              final int dp = p1[i] - p2[i];
1562              sum += dp * dp;
1563          }
1564          return Math.sqrt(sum);
1565        }
1566        
1567        /**
1568         * Calculates the L<sub>&infin;</sub> (max of abs) distance between two points.
1569         *
1570         * @param p1 the first point
1571         * @param p2 the second point
1572         * @return the L<sub>&infin;</sub> distance between the two points
1573         */
1574        public static final double distanceInf(double[] p1, double[] p2) {
1575            double max = 0;
1576            for (int i = 0; i < p1.length; i++) {
1577                max = Math.max(max, Math.abs(p1[i] - p2[i]));
1578            }
1579            return max;
1580        }
1581        
1582        /**
1583         * Calculates the L<sub>&infin;</sub> (max of abs) distance between two points.
1584         *
1585         * @param p1 the first point
1586         * @param p2 the second point
1587         * @return the L<sub>&infin;</sub> distance between the two points
1588         */
1589        public static final int distanceInf(int[] p1, int[] p2) {
1590            int max = 0;
1591            for (int i = 0; i < p1.length; i++) {
1592                max = Math.max(max, Math.abs(p1[i] - p2[i]));
1593            }
1594            return max;
1595        }
1596    
1597        
1598    }