View Javadoc

1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  
18  package org.apache.commons.math.util;
19  
20  import java.math.BigDecimal;
21  import java.math.BigInteger;
22  import java.util.Arrays;
23  
24  import org.apache.commons.math.MathRuntimeException;
25  
26  /**
27   * Some useful additions to the built-in functions in {@link Math}.
28   * @version $Revision: 772119 $ $Date: 2009-05-06 05:43:28 -0400 (Wed, 06 May 2009) $
29   */
30  public final class MathUtils {
31  
32      /** Smallest positive number such that 1 - EPSILON is not numerically equal to 1. */
33      public static final double EPSILON = 0x1.0p-53;
34  
35      /** Safe minimum, such that 1 / SAFE_MIN does not overflow.
36       * <p>In IEEE 754 arithmetic, this is also the smallest normalized
37       * number 2<sup>-1022</sup>.</p>
38       */
39      public static final double SAFE_MIN = 0x1.0p-1022;
40  
41      /** -1.0 cast as a byte. */
42      private static final byte  NB = (byte)-1;
43  
44      /** -1.0 cast as a short. */
45      private static final short NS = (short)-1;
46  
47      /** 1.0 cast as a byte. */
48      private static final byte  PB = (byte)1;
49  
50      /** 1.0 cast as a short. */
51      private static final short PS = (short)1;
52  
53      /** 0.0 cast as a byte. */
54      private static final byte  ZB = (byte)0;
55  
56      /** 0.0 cast as a short. */
57      private static final short ZS = (short)0;
58  
59      /** 2 &pi;. */
60      private static final double TWO_PI = 2 * Math.PI;
61  
62      /** Gap between NaN and regular numbers. */
63      private static final int NAN_GAP = 4 * 1024 * 1024;
64  
65      /** Offset to order signed double numbers lexicographically. */
66      private static final long SGN_MASK = 0x8000000000000000L;
67  
68      /**
69       * Private Constructor
70       */
71      private MathUtils() {
72          super();
73      }
74  
75      /**
76       * Add two integers, checking for overflow.
77       * 
78       * @param x an addend
79       * @param y an addend
80       * @return the sum <code>x+y</code>
81       * @throws ArithmeticException if the result can not be represented as an
82       *         int
83       * @since 1.1
84       */
85      public static int addAndCheck(int x, int y) {
86          long s = (long)x + (long)y;
87          if (s < Integer.MIN_VALUE || s > Integer.MAX_VALUE) {
88              throw new ArithmeticException("overflow: add");
89          }
90          return (int)s;
91      }
92  
93      /**
94       * Add two long integers, checking for overflow.
95       * 
96       * @param a an addend
97       * @param b an addend
98       * @return the sum <code>a+b</code>
99       * @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 }