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