    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    import org.apache.commons.math.exception.util.Localizable;
    026    import org.apache.commons.math.exception.util.LocalizedFormats;
    027    import org.apache.commons.math.exception.NonMonotonousSequenceException;
    028
    029    /**
    030     * Some useful additions to the built-in functions in {@link Math}.
    031     * @version $Revision: 1073472 $ $Date: 2011-02-22 20:49:07 +0100 (mar. 22 f??vr. 2011) $
    032     */
    033    public final class MathUtils {
    034
    035        /** Smallest positive number such that 1 - EPSILON is not numerically equal to 1. */
    036        public static final double EPSILON = 0x1.0p-53;
    037
    038        /** Safe minimum, such that 1 / SAFE_MIN does not overflow.
    039         * <p>In IEEE 754 arithmetic, this is also the smallest normalized
    040         * number 2<sup>-1022</sup>.</p>
    041         */
    042        public static final double SAFE_MIN = 0x1.0p-1022;
    043
    044        /**
    045         * 2 &pi;.
    046         * @since 2.1
    047         */
    048        public static final double TWO_PI = 2 * FastMath.PI;
    049
    050        /** -1.0 cast as a byte. */
    051        private static final byte  NB = (byte)-1;
    052
    053        /** -1.0 cast as a short. */
    054        private static final short NS = (short)-1;
    055
    056        /** 1.0 cast as a byte. */
    057        private static final byte  PB = (byte)1;
    058
    059        /** 1.0 cast as a short. */
    060        private static final short PS = (short)1;
    061
    062        /** 0.0 cast as a byte. */
    063        private static final byte  ZB = (byte)0;
    064
    065        /** 0.0 cast as a short. */
    066        private static final short ZS = (short)0;
    067
    068        /** Gap between NaN and regular numbers. */
    069        private static final int NAN_GAP = 4 * 1024 * 1024;
    070
    071        /** Offset to order signed double numbers lexicographically. */
    072        private static final long SGN_MASK = 0x8000000000000000L;
    073
    074        /** Offset to order signed double numbers lexicographically. */
    075        private static final int SGN_MASK_FLOAT = 0x80000000;
    076
    077        /** All long-representable factorials */
    078        private static final long[] FACTORIALS = new long[] {
    079                           1l,                  1l,                   2l,
    080                           6l,                 24l,                 120l,
    081                         720l,               5040l,               40320l,
    082                      362880l,            3628800l,            39916800l,
    083                   479001600l,         6227020800l,         87178291200l,
    084               1307674368000l,     20922789888000l,     355687428096000l,
    085            6402373705728000l, 121645100408832000l, 2432902008176640000l };
    086
    087        /**
    088         * Private Constructor
    089         */
    090        private MathUtils() {
    091            super();
    092        }
    093
    094        /**
    095         * Add two integers, checking for overflow.
    096         *
    097         * @param x an addend
    098         * @param y an addend
    099         * @return the sum <code>x+y</code>
    100         * @throws ArithmeticException if the result can not be represented as an
    101         *         int
    102         * @since 1.1
    103         */
    104        public static int addAndCheck(int x, int y) {
    105            long s = (long)x + (long)y;
    106            if (s < Integer.MIN_VALUE || s > Integer.MAX_VALUE) {
    107                throw MathRuntimeException.createArithmeticException(LocalizedFormats.OVERFLOW_IN_ADDITION, x, y);
    108            }
    109            return (int)s;
    110        }
    111
    112        /**
    113         * Add two long integers, checking for overflow.
    114         *
    115         * @param a an addend
    116         * @param b an addend
    117         * @return the sum <code>a+b</code>
    118         * @throws ArithmeticException if the result can not be represented as an
    119         *         long
    120         * @since 1.2
    121         */
    122        public static long addAndCheck(long a, long b) {
    123            return addAndCheck(a, b, LocalizedFormats.OVERFLOW_IN_ADDITION);
    124        }
    125
    126        /**
    127         * Add two long integers, checking for overflow.
    128         *
    129         * @param a an addend
    130         * @param b an addend
    131         * @param pattern the pattern to use for any thrown exception.
    132         * @return the sum <code>a+b</code>
    133         * @throws ArithmeticException if the result can not be represented as an
    134         *         long
    135         * @since 1.2
    136         */
    137        private static long addAndCheck(long a, long b, Localizable pattern) {
    138            long ret;
    139            if (a > b) {
    140                // use symmetry to reduce boundary cases
    141                ret = addAndCheck(b, a, pattern);
    142            } else {
    143                // assert a <= b
    144
    145                if (a < 0) {
    146                    if (b < 0) {
    147                        // check for negative overflow
    148                        if (Long.MIN_VALUE - b <= a) {
    149                            ret = a + b;
    150                        } else {
    151                            throw MathRuntimeException.createArithmeticException(pattern, a, b);
    152                        }
    153                    } else {
    154                        // opposite sign addition is always safe
    155                        ret = a + b;
    156                    }
    157                } else {
    158                    // assert a >= 0
    159                    // assert b >= 0
    160
    161                    // check for positive overflow
    162                    if (a <= Long.MAX_VALUE - b) {
    163                        ret = a + b;
    164                    } else {
    165                        throw MathRuntimeException.createArithmeticException(pattern, a, b);
    166                    }
    167                }
    168            }
    169            return ret;
    170        }
    171
    172        /**
    173         * Returns an exact representation of the <a
    174         * href="http://mathworld.wolfram.com/BinomialCoefficient.html"> Binomial
    175         * Coefficient</a>, "<code>n choose k</code>", the number of
    176         * <code>k</code>-element subsets that can be selected from an
    177         * <code>n</code>-element set.
    178         * <p>
    179         * <Strong>Preconditions</strong>:
    180         * <ul>
    181         * <li> <code>0 <= k <= n </code> (otherwise
    182         * <code>IllegalArgumentException</code> is thrown)</li>
    183         * <li> The result is small enough to fit into a <code>long</code>. The
    184         * largest value of <code>n</code> for which all coefficients are
    185         * <code> < Long.MAX_VALUE</code> is 66. If the computed value exceeds
    186         * <code>Long.MAX_VALUE</code> an <code>ArithMeticException</code> is
    187         * thrown.</li>
    188         * </ul></p>
    189         *
    190         * @param n the size of the set
    191         * @param k the size of the subsets to be counted
    192         * @return <code>n choose k</code>
    193         * @throws IllegalArgumentException if preconditions are not met.
    194         * @throws ArithmeticException if the result is too large to be represented
    195         *         by a long integer.
    196         */
    197        public static long binomialCoefficient(final int n, final int k) {
    198            checkBinomial(n, k);
    199            if ((n == k) || (k == 0)) {
    200                return 1;
    201            }
    202            if ((k == 1) || (k == n - 1)) {
    203                return n;
    204            }
    205            // Use symmetry for large k
    206            if (k > n / 2)
    207                return binomialCoefficient(n, n - k);
    208
    209            // We use the formula
    210            // (n choose k) = n! / (n-k)! / k!
    211            // (n choose k) == ((n-k+1)*...*n) / (1*...*k)
    212            // which could be written
    213            // (n choose k) == (n-1 choose k-1) * n / k
    214            long result = 1;
    215            if (n <= 61) {
    216                // For n <= 61, the naive implementation cannot overflow.
    217                int i = n - k + 1;
    218                for (int j = 1; j <= k; j++) {
    219                    result = result * i / j;
    220                    i++;
    221                }
    222            } else if (n <= 66) {
    223                // For n > 61 but n <= 66, the result cannot overflow,
    224                // but we must take care not to overflow intermediate values.
    225                int i = n - k + 1;
    226                for (int j = 1; j <= k; j++) {
    227                    // We know that (result * i) is divisible by j,
    228                    // but (result * i) may overflow, so we split j:
    229                    // Filter out the gcd, d, so j/d and i/d are integer.
    230                    // result is divisible by (j/d) because (j/d)
    231                    // is relative prime to (i/d) and is a divisor of
    232                    // result * (i/d).
    233                    final long d = gcd(i, j);
    234                    result = (result / (j / d)) * (i / d);
    235                    i++;
    236                }
    237            } else {
    238                // For n > 66, a result overflow might occur, so we check
    239                // the multiplication, taking care to not overflow
    240                // unnecessary.
    241                int i = n - k + 1;
    242                for (int j = 1; j <= k; j++) {
    243                    final long d = gcd(i, j);
    244                    result = mulAndCheck(result / (j / d), i / d);
    245                    i++;
    246                }
    247            }
    248            return result;
    249        }
    250
    251        /**
    252         * Returns a <code>double</code> representation of the <a
    253         * href="http://mathworld.wolfram.com/BinomialCoefficient.html"> Binomial
    254         * Coefficient</a>, "<code>n choose k</code>", the number of
    255         * <code>k</code>-element subsets that can be selected from an
    256         * <code>n</code>-element set.
    257         * <p>
    258         * <Strong>Preconditions</strong>:
    259         * <ul>
    260         * <li> <code>0 <= k <= n </code> (otherwise
    261         * <code>IllegalArgumentException</code> is thrown)</li>
    262         * <li> The result is small enough to fit into a <code>double</code>. The
    263         * largest value of <code>n</code> for which all coefficients are <
    264         * Double.MAX_VALUE is 1029. If the computed value exceeds Double.MAX_VALUE,
    265         * Double.POSITIVE_INFINITY is returned</li>
    266         * </ul></p>
    267         *
    268         * @param n the size of the set
    269         * @param k the size of the subsets to be counted
    270         * @return <code>n choose k</code>
    271         * @throws IllegalArgumentException if preconditions are not met.
    272         */
    273        public static double binomialCoefficientDouble(final int n, final int k) {
    274            checkBinomial(n, k);
    275            if ((n == k) || (k == 0)) {
    276                return 1d;
    277            }
    278            if ((k == 1) || (k == n - 1)) {
    279                return n;
    280            }
    281            if (k > n/2) {
    282                return binomialCoefficientDouble(n, n - k);
    283            }
    284            if (n < 67) {
    285                return binomialCoefficient(n,k);
    286            }
    287
    288            double result = 1d;
    289            for (int i = 1; i <= k; i++) {
    290                 result *= (double)(n - k + i) / (double)i;
    291            }
    292
    293            return FastMath.floor(result + 0.5);
    294        }
    295
    296        /**
    297         * Returns the natural <code>log</code> of the <a
    298         * href="http://mathworld.wolfram.com/BinomialCoefficient.html"> Binomial
    299         * Coefficient</a>, "<code>n choose k</code>", the number of
    300         * <code>k</code>-element subsets that can be selected from an
    301         * <code>n</code>-element set.
    302         * <p>
    303         * <Strong>Preconditions</strong>:
    304         * <ul>
    305         * <li> <code>0 <= k <= n </code> (otherwise
    306         * <code>IllegalArgumentException</code> is thrown)</li>
    307         * </ul></p>
    308         *
    309         * @param n the size of the set
    310         * @param k the size of the subsets to be counted
    311         * @return <code>n choose k</code>
    312         * @throws IllegalArgumentException if preconditions are not met.
    313         */
    314        public static double binomialCoefficientLog(final int n, final int k) {
    315            checkBinomial(n, k);
    316            if ((n == k) || (k == 0)) {
    317                return 0;
    318            }
    319            if ((k == 1) || (k == n - 1)) {
    320                return FastMath.log(n);
    321            }
    322
    323            /*
    324             * For values small enough to do exact integer computation,
    325             * return the log of the exact value
    326             */
    327            if (n < 67) {
    328                return FastMath.log(binomialCoefficient(n,k));
    329            }
    330
    331            /*
    332             * Return the log of binomialCoefficientDouble for values that will not
    333             * overflow binomialCoefficientDouble
    334             */
    335            if (n < 1030) {
    336                return FastMath.log(binomialCoefficientDouble(n, k));
    337            }
    338
    339            if (k > n / 2) {
    340                return binomialCoefficientLog(n, n - k);
    341            }
    342
    343            /*
    344             * Sum logs for values that could overflow
    345             */
    346            double logSum = 0;
    347
    348            // n!/(n-k)!
    349            for (int i = n - k + 1; i <= n; i++) {
    350                logSum += FastMath.log(i);
    351            }
    352
    353            // divide by k!
    354            for (int i = 2; i <= k; i++) {
    355                logSum -= FastMath.log(i);
    356            }
    357
    358            return logSum;
    359        }
    360
    361        /**
    362         * Check binomial preconditions.
    363         * @param n the size of the set
    364         * @param k the size of the subsets to be counted
    365         * @exception IllegalArgumentException if preconditions are not met.
    366         */
    367        private static void checkBinomial(final int n, final int k)
    368            throws IllegalArgumentException {
    369            if (n < k) {
    370                throw MathRuntimeException.createIllegalArgumentException(
    371                    LocalizedFormats.BINOMIAL_INVALID_PARAMETERS_ORDER,
    372                    n, k);
    373            }
    374            if (n < 0) {
    375                throw MathRuntimeException.createIllegalArgumentException(
    376                      LocalizedFormats.BINOMIAL_NEGATIVE_PARAMETER,
    377                      n);
    378            }
    379        }
    380
    381        /**
    382         * Compares two numbers given some amount of allowed error.
    383         *
    384         * @param x the first number
    385         * @param y the second number
    386         * @param eps the amount of error to allow when checking for equality
    387         * @return <ul><li>0 if  {@link #equals(double, double, double) equals(x, y, eps)}</li>
    388         *       <li>&lt; 0 if !{@link #equals(double, double, double) equals(x, y, eps)} &amp;&amp; x &lt; y</li>
    389         *       <li>> 0 if !{@link #equals(double, double, double) equals(x, y, eps)} &amp;&amp; x > y</li></ul>
    390         */
    391        public static int compareTo(double x, double y, double eps) {
    392            if (equals(x, y, eps)) {
    393                return 0;
    394            } else if (x < y) {
    395              return -1;
    396            }
    397            return 1;
    398        }
    399
    400        /**
    401         * Returns the <a href="http://mathworld.wolfram.com/HyperbolicCosine.html">
    402         * hyperbolic cosine</a> of x.
    403         *
    404         * @param x double value for which to find the hyperbolic cosine
    405         * @return hyperbolic cosine of x
    406         */
    407        public static double cosh(double x) {
    408            return (FastMath.exp(x) + FastMath.exp(-x)) / 2.0;
    409        }
    410
    411        /**
    412         * Returns true iff they are strictly equal.
    413         *
    414         * @param x first value
    415         * @param y second value
    416         * @return {@code true} if the values are equal.
    417         * @deprecated as of 2.2 his method considers that {@code NaN == NaN}. In release
    418         * 3.0, the semantics will change in order to comply with IEEE754 where it
    419         * is specified that {@code NaN != NaN}.
    420         * New methods have been added for those cases wher the old semantics is
    421         * useful (see e.g. {@link #equalsIncludingNaN(float,float)
    422         * equalsIncludingNaN}.
    423         */
    424        @Deprecated
    425        public static boolean equals(float x, float y) {
    426            return (Float.isNaN(x) && Float.isNaN(y)) || x == y;
    427        }
    428
    429        /**
    430         * Returns true if both arguments are NaN or neither is NaN and they are
    431         * equal as defined by {@link #equals(float,float,int) equals(x, y, 1)}.
    432         *
    433         * @param x first value
    434         * @param y second value
    435         * @return {@code true} if the values are equal or both are NaN.
    436         * @since 2.2
    437         */
    438        public static boolean equalsIncludingNaN(float x, float y) {
    439            return (Float.isNaN(x) && Float.isNaN(y)) || equals(x, y, 1);
    440        }
    441
    442        /**
    443         * Returns true if both arguments are equal or within the range of allowed
    444         * error (inclusive).
    445         *
    446         * @param x first value
    447         * @param y second value
    448         * @param eps the amount of absolute error to allow.
    449         * @return {@code true} if the values are equal or within range of each other.
    450         * @since 2.2
    451         */
    452        public static boolean equals(float x, float y, float eps) {
    453            return equals(x, y, 1) || FastMath.abs(y - x) <= eps;
    454        }
    455
    456        /**
    457         * Returns true if both arguments are NaN or are equal or within the range
    458         * of allowed error (inclusive).
    459         *
    460         * @param x first value
    461         * @param y second value
    462         * @param eps the amount of absolute error to allow.
    463         * @return {@code true} if the values are equal or within range of each other,
    464         * or both are NaN.
    465         * @since 2.2
    466         */
    467        public static boolean equalsIncludingNaN(float x, float y, float eps) {
    468            return equalsIncludingNaN(x, y) || (FastMath.abs(y - x) <= eps);
    469        }
    470
    471        /**
    472         * Returns true if both arguments are equal or within the range of allowed
    473         * error (inclusive).
    474         * Two float numbers are considered equal if there are {@code (maxUlps - 1)}
    475         * (or fewer) floating point numbers between them, i.e. two adjacent floating
    476         * point numbers are considered equal.
    477         * Adapted from <a
    478         * href="http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm">
    479         * Bruce Dawson</a>
    480         *
    481         * @param x first value
    482         * @param y second value
    483         * @param maxUlps {@code (maxUlps - 1)} is the number of floating point
    484         * values between {@code x} and {@code y}.
    485         * @return {@code true} if there are fewer than {@code maxUlps} floating
    486         * point values between {@code x} and {@code y}.
    487         * @since 2.2
    488         */
    489        public static boolean equals(float x, float y, int maxUlps) {
    490            // Check that "maxUlps" is non-negative and small enough so that
    491            // NaN won't compare as equal to anything (except another NaN).
    492            assert maxUlps > 0 && maxUlps < NAN_GAP;
    493
    494            int xInt = Float.floatToIntBits(x);
    495            int yInt = Float.floatToIntBits(y);
    496
    497            // Make lexicographically ordered as a two's-complement integer.
    498            if (xInt < 0) {
    499                xInt = SGN_MASK_FLOAT - xInt;
    500            }
    501            if (yInt < 0) {
    502                yInt = SGN_MASK_FLOAT - yInt;
    503            }
    504
    505            final boolean isEqual = FastMath.abs(xInt - yInt) <= maxUlps;
    506
    507            return isEqual && !Float.isNaN(x) && !Float.isNaN(y);
    508        }
    509
    510        /**
    511         * Returns true if both arguments are NaN or if they are equal as defined
    512         * by {@link #equals(float,float,int) equals(x, y, maxUlps)}.
    513         *
    514         * @param x first value
    515         * @param y second value
    516         * @param maxUlps {@code (maxUlps - 1)} is the number of floating point
    517         * values between {@code x} and {@code y}.
    518         * @return {@code true} if both arguments are NaN or if there are less than
    519         * {@code maxUlps} floating point values between {@code x} and {@code y}.
    520         * @since 2.2
    521         */
    522        public static boolean equalsIncludingNaN(float x, float y, int maxUlps) {
    523            return (Float.isNaN(x) && Float.isNaN(y)) || equals(x, y, maxUlps);
    524        }
    525
    526        /**
    527         * Returns true iff both arguments are null or have same dimensions and all
    528         * their elements are equal as defined by
    529         * {@link #equals(float,float)}.
    530         *
    531         * @param x first array
    532         * @param y second array
    533         * @return true if the values are both null or have same dimension
    534         * and equal elements.
    535         * @deprecated as of 2.2 this method considers that {@code NaN == NaN}. In release
    536         * 3.0, the semantics will change in order to comply with IEEE754 where it
    537         * is specified that {@code NaN != NaN}.
    538         * New methods have been added for those cases where the old semantics is
    539         * useful (see e.g. {@link #equalsIncludingNaN(float[],float[])
    540         * equalsIncludingNaN}.
    541         */
    542        @Deprecated
    543        public static boolean equals(float[] x, float[] y) {
    544            if ((x == null) || (y == null)) {
    545                return !((x == null) ^ (y == null));
    546            }
    547            if (x.length != y.length) {
    548                return false;
    549            }
    550            for (int i = 0; i < x.length; ++i) {
    551                if (!equals(x[i], y[i])) {
    552                    return false;
    553                }
    554            }
    555            return true;
    556        }
    557
    558        /**
    559         * Returns true iff both arguments are null or have same dimensions and all
    560         * their elements are equal as defined by
    561         * {@link #equalsIncludingNaN(float,float)}.
    562         *
    563         * @param x first array
    564         * @param y second array
    565         * @return true if the values are both null or have same dimension and
    566         * equal elements
    567         * @since 2.2
    568         */
    569        public static boolean equalsIncludingNaN(float[] x, float[] y) {
    570            if ((x == null) || (y == null)) {
    571                return !((x == null) ^ (y == null));
    572            }
    573            if (x.length != y.length) {
    574                return false;
    575            }
    576            for (int i = 0; i < x.length; ++i) {
    577                if (!equalsIncludingNaN(x[i], y[i])) {
    578                    return false;
    579                }
    580            }
    581            return true;
    582        }
    583
    584        /**
    585         * Returns true iff both arguments are NaN or neither is NaN and they are
    586         * equal
    587         *
    588         * <p>This method considers that {@code NaN == NaN}. In release
    589         * 3.0, the semantics will change in order to comply with IEEE754 where it
    590         * is specified that {@code NaN != NaN}.
    591         * New methods have been added for those cases where the old semantics
    592         * (w.r.t. NaN) is useful (see e.g.
    593         * {@link #equalsIncludingNaN(double,double, double) equalsIncludingNaN}.
    594         * </p>
    595         *
    596         * @param x first value
    597         * @param y second value
    598         * @return {@code true} if the values are equal.
    599         */
    600        public static boolean equals(double x, double y) {
    601            return (Double.isNaN(x) && Double.isNaN(y)) || x == y;
    602        }
    603
    604        /**
    605         * Returns true if both arguments are NaN or neither is NaN and they are
    606         * equal as defined by {@link #equals(double,double,int) equals(x, y, 1)}.
    607         *
    608         * @param x first value
    609         * @param y second value
    610         * @return {@code true} if the values are equal or both are NaN.
    611         * @since 2.2
    612         */
    613        public static boolean equalsIncludingNaN(double x, double y) {
    614            return (Double.isNaN(x) && Double.isNaN(y)) || equals(x, y, 1);
    615        }
    616
    617        /**
    618         * Returns true if both arguments are equal or within the range of allowed
    619         * error (inclusive).
    620         * <p>
    621         * Two NaNs are considered equals, as are two infinities with same sign.
    622         * </p>
    623         * <p>This method considers that {@code NaN == NaN}. In release
    624         * 3.0, the semantics will change in order to comply with IEEE754 where it
    625         * is specified that {@code NaN != NaN}.
    626         * New methods have been added for those cases where the old semantics
    627         * (w.r.t. NaN) is useful (see e.g.
    628         * {@link #equalsIncludingNaN(double,double, double) equalsIncludingNaN}.
    629         * </p>
    630         * @param x first value
    631         * @param y second value
    632         * @param eps the amount of absolute error to allow.
    633         * @return {@code true} if the values are equal or within range of each other.
    634         */
    635        public static boolean equals(double x, double y, double eps) {
    636            return equals(x, y) || FastMath.abs(y - x) <= eps;
    637        }
    638
    639        /**
    640         * Returns true if both arguments are NaN or are equal or within the range
    641         * of allowed error (inclusive).
    642         *
    643         * @param x first value
    644         * @param y second value
    645         * @param eps the amount of absolute error to allow.
    646         * @return {@code true} if the values are equal or within range of each other,
    647         * or both are NaN.
    648         * @since 2.2
    649         */
    650        public static boolean equalsIncludingNaN(double x, double y, double eps) {
    651            return equalsIncludingNaN(x, y) || (FastMath.abs(y - x) <= eps);
    652        }
    653
    654        /**
    655         * Returns true if both arguments are equal or within the range of allowed
    656         * error (inclusive).
    657         * Two float numbers are considered equal if there are {@code (maxUlps - 1)}
    658         * (or fewer) floating point numbers between them, i.e. two adjacent floating
    659         * point numbers are considered equal.
    660         * Adapted from <a
    661         * href="http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm">
    662         * Bruce Dawson</a>
    663         *
    664         * <p>This method considers that {@code NaN == NaN}. In release
    665         * 3.0, the semantics will change in order to comply with IEEE754 where it
    666         * is specified that {@code NaN != NaN}.
    667         * New methods have been added for those cases where the old semantics
    668         * (w.r.t. NaN) is useful (see e.g.
    669         * {@link #equalsIncludingNaN(double,double, int) equalsIncludingNaN}.
    670         * </p>
    671         *
    672         * @param x first value
    673         * @param y second value
    674         * @param maxUlps {@code (maxUlps - 1)} is the number of floating point
    675         * values between {@code x} and {@code y}.
    676         * @return {@code true} if there are fewer than {@code maxUlps} floating
    677         * point values between {@code x} and {@code y}.
    678         */
    679        public static boolean equals(double x, double y, int maxUlps) {
    680            // Check that "maxUlps" is non-negative and small enough so that
    681            // NaN won't compare as equal to anything (except another NaN).
    682            assert maxUlps > 0 && maxUlps < NAN_GAP;
    683
    684            long xInt = Double.doubleToLongBits(x);
    685            long yInt = Double.doubleToLongBits(y);
    686
    687            // Make lexicographically ordered as a two's-complement integer.
    688            if (xInt < 0) {
    689                xInt = SGN_MASK - xInt;
    690            }
    691            if (yInt < 0) {
    692                yInt = SGN_MASK - yInt;
    693            }
    694
    695            return FastMath.abs(xInt - yInt) <= maxUlps;
    696        }
    697
    698        /**
    699         * Returns true if both arguments are NaN or if they are equal as defined
    700         * by {@link #equals(double,double,int) equals(x, y, maxUlps}.
    701         *
    702         * @param x first value
    703         * @param y second value
    704         * @param maxUlps {@code (maxUlps - 1)} is the number of floating point
    705         * values between {@code x} and {@code y}.
    706         * @return {@code true} if both arguments are NaN or if there are less than
    707         * {@code maxUlps} floating point values between {@code x} and {@code y}.
    708         * @since 2.2
    709         */
    710        public static boolean equalsIncludingNaN(double x, double y, int maxUlps) {
    711            return (Double.isNaN(x) && Double.isNaN(y)) || equals(x, y, maxUlps);
    712        }
    713
    714        /**
    715         * Returns true iff both arguments are null or have same dimensions and all
    716         * their elements are equal as defined by
    717         * {@link #equals(double,double)}.
    718         *
    719         * <p>This method considers that {@code NaN == NaN}. In release
    720         * 3.0, the semantics will change in order to comply with IEEE754 where it
    721         * is specified that {@code NaN != NaN}.
    722         * New methods have been added for those cases wher the old semantics is
    723         * useful (see e.g. {@link #equalsIncludingNaN(double[],double[])
    724         * equalsIncludingNaN}.
    725         * </p>
    726         * @param x first array
    727         * @param y second array
    728         * @return true if the values are both null or have same dimension
    729         * and equal elements.
    730         */
    731        public static boolean equals(double[] x, double[] y) {
    732            if ((x == null) || (y == null)) {
    733                return !((x == null) ^ (y == null));
    734            }
    735            if (x.length != y.length) {
    736                return false;
    737            }
    738            for (int i = 0; i < x.length; ++i) {
    739                if (!equals(x[i], y[i])) {
    740                    return false;
    741                }
    742            }
    743            return true;
    744        }
    745
    746        /**
    747         * Returns true iff both arguments are null or have same dimensions and all
    748         * their elements are equal as defined by
    749         * {@link #equalsIncludingNaN(double,double)}.
    750         *
    751         * @param x first array
    752         * @param y second array
    753         * @return true if the values are both null or have same dimension and
    754         * equal elements
    755         * @since 2.2
    756         */
    757        public static boolean equalsIncludingNaN(double[] x, double[] y) {
    758            if ((x == null) || (y == null)) {
    759                return !((x == null) ^ (y == null));
    760            }
    761            if (x.length != y.length) {
    762                return false;
    763            }
    764            for (int i = 0; i < x.length; ++i) {
    765                if (!equalsIncludingNaN(x[i], y[i])) {
    766                    return false;
    767                }
    768            }
    769            return true;
    770        }
    771
    772        /**
    773         * Returns n!. Shorthand for <code>n</code> <a
    774         * href="http://mathworld.wolfram.com/Factorial.html"> Factorial</a>, the
    775         * product of the numbers <code>1,...,n</code>.
    776         * <p>
    777         * <Strong>Preconditions</strong>:
    778         * <ul>
    779         * <li> <code>n >= 0</code> (otherwise
    780         * <code>IllegalArgumentException</code> is thrown)</li>
    781         * <li> The result is small enough to fit into a <code>long</code>. The
    782         * largest value of <code>n</code> for which <code>n!</code> <
    783         * Long.MAX_VALUE</code> is 20. If the computed value exceeds <code>Long.MAX_VALUE</code>
    784         * an <code>ArithMeticException </code> is thrown.</li>
    785         * </ul>
    786         * </p>
    787         *
    788         * @param n argument
    789         * @return <code>n!</code>
    790         * @throws ArithmeticException if the result is too large to be represented
    791         *         by a long integer.
    792         * @throws IllegalArgumentException if n < 0
    793         */
    794        public static long factorial(final int n) {
    795            if (n < 0) {
    796                throw MathRuntimeException.createIllegalArgumentException(
    797                      LocalizedFormats.FACTORIAL_NEGATIVE_PARAMETER,
    798                      n);
    799            }
    800            if (n > 20) {
    801                throw new ArithmeticException(
    802                        "factorial value is too large to fit in a long");
    803            }
    804            return FACTORIALS[n];
    805        }
    806
    807        /**
    808         * Returns n!. Shorthand for <code>n</code> <a
    809         * href="http://mathworld.wolfram.com/Factorial.html"> Factorial</a>, the
    810         * product of the numbers <code>1,...,n</code> as a <code>double</code>.
    811         * <p>
    812         * <Strong>Preconditions</strong>:
    813         * <ul>
    814         * <li> <code>n >= 0</code> (otherwise
    815         * <code>IllegalArgumentException</code> is thrown)</li>
    816         * <li> The result is small enough to fit into a <code>double</code>. The
    817         * largest value of <code>n</code> for which <code>n!</code> <
    818         * Double.MAX_VALUE</code> is 170. If the computed value exceeds
    819         * Double.MAX_VALUE, Double.POSITIVE_INFINITY is returned</li>
    820         * </ul>
    821         * </p>
    822         *
    823         * @param n argument
    824         * @return <code>n!</code>
    825         * @throws IllegalArgumentException if n < 0
    826         */
    827        public static double factorialDouble(final int n) {
    828            if (n < 0) {
    829                throw MathRuntimeException.createIllegalArgumentException(
    830                      LocalizedFormats.FACTORIAL_NEGATIVE_PARAMETER,
    831                      n);
    832            }
    833            if (n < 21) {
    834                return factorial(n);
    835            }
    836            return FastMath.floor(FastMath.exp(factorialLog(n)) + 0.5);
    837        }
    838
    839        /**
    840         * Returns the natural logarithm of n!.
    841         * <p>
    842         * <Strong>Preconditions</strong>:
    843         * <ul>
    844         * <li> <code>n >= 0</code> (otherwise
    845         * <code>IllegalArgumentException</code> is thrown)</li>
    846         * </ul></p>
    847         *
    848         * @param n argument
    849         * @return <code>n!</code>
    850         * @throws IllegalArgumentException if preconditions are not met.
    851         */
    852        public static double factorialLog(final int n) {
    853            if (n < 0) {
    854                throw MathRuntimeException.createIllegalArgumentException(
    855                      LocalizedFormats.FACTORIAL_NEGATIVE_PARAMETER,
    856                      n);
    857            }
    858            if (n < 21) {
    859                return FastMath.log(factorial(n));
    860            }
    861            double logSum = 0;
    862            for (int i = 2; i <= n; i++) {
    863                logSum += FastMath.log(i);
    864            }
    865            return logSum;
    866        }
    867
    868        /**
    869         * <p>
    870         * Gets the greatest common divisor of the absolute value of two numbers,
    871         * using the "binary gcd" method which avoids division and modulo
    872         * operations. See Knuth 4.5.2 algorithm B. This algorithm is due to Josef
    873         * Stein (1961).
    874         * </p>
    875         * Special cases:
    876         * <ul>
    877         * <li>The invocations
    878         * <code>gcd(Integer.MIN_VALUE, Integer.MIN_VALUE)</code>,
    879         * <code>gcd(Integer.MIN_VALUE, 0)</code> and
    880         * <code>gcd(0, Integer.MIN_VALUE)</code> throw an
    881         * <code>ArithmeticException</code>, because the result would be 2^31, which
    882         * is too large for an int value.</li>
    883         * <li>The result of <code>gcd(x, x)</code>, <code>gcd(0, x)</code> and
    884         * <code>gcd(x, 0)</code> is the absolute value of <code>x</code>, except
    885         * for the special cases above.
    886         * <li>The invocation <code>gcd(0, 0)</code> is the only one which returns
    887         * <code>0</code>.</li>
    888         * </ul>
    889         *
    890         * @param p any number
    891         * @param q any number
    892         * @return the greatest common divisor, never negative
    893         * @throws ArithmeticException if the result cannot be represented as a
    894         * nonnegative int value
    895         * @since 1.1
    896         */
    897        public static int gcd(final int p, final int q) {
    898            int u = p;
    899            int v = q;
    900            if ((u == 0) || (v == 0)) {
    901                if ((u == Integer.MIN_VALUE) || (v == Integer.MIN_VALUE)) {
    902                    throw MathRuntimeException.createArithmeticException(
    903                            LocalizedFormats.GCD_OVERFLOW_32_BITS,
    904                            p, q);
    905                }
    906                return FastMath.abs(u) + FastMath.abs(v);
    907            }
    908            // keep u and v negative, as negative integers range down to
    909            // -2^31, while positive numbers can only be as large as 2^31-1
    910            // (i.e. we can't necessarily negate a negative number without
    911            // overflow)
    912            /* assert u!=0 && v!=0; */
    913            if (u > 0) {
    914                u = -u;
    915            } // make u negative
    916            if (v > 0) {
    917                v = -v;
    918            } // make v negative
    919            // B1. [Find power of 2]
    920            int k = 0;
    921            while ((u & 1) == 0 && (v & 1) == 0 && k < 31) { // while u and v are
    922                                                                // both even...
    923                u /= 2;
    924                v /= 2;
    925                k++; // cast out twos.
    926            }
    927            if (k == 31) {
    928                throw MathRuntimeException.createArithmeticException(
    929                        LocalizedFormats.GCD_OVERFLOW_32_BITS,
    930                        p, q);
    931            }
    932            // B2. Initialize: u and v have been divided by 2^k and at least
    933            // one is odd.
    934            int t = ((u & 1) == 1) ? v : -(u / 2)/* B3 */;
    935            // t negative: u was odd, v may be even (t replaces v)
    936            // t positive: u was even, v is odd (t replaces u)
    937            do {
    938                /* assert u<0 && v<0; */
    939                // B4/B3: cast out twos from t.
    940                while ((t & 1) == 0) { // while t is even..
    941                    t /= 2; // cast out twos
    942                }
    943                // B5 [reset max(u,v)]
    944                if (t > 0) {
    945                    u = -t;
    946                } else {
    947                    v = t;
    948                }
    949                // B6/B3. at this point both u and v should be odd.
    950                t = (v - u) / 2;
    951                // |u| larger: t positive (replace u)
    952                // |v| larger: t negative (replace v)
    953            } while (t != 0);
    954            return -u * (1 << k); // gcd is u*2^k
    955        }
    956
    957        /**
    958         * <p>
    959         * Gets the greatest common divisor of the absolute value of two numbers,
    960         * using the "binary gcd" method which avoids division and modulo
    961         * operations. See Knuth 4.5.2 algorithm B. This algorithm is due to Josef
    962         * Stein (1961).
    963         * </p>
    964         * Special cases:
    965         * <ul>
    966         * <li>The invocations
    967         * <code>gcd(Long.MIN_VALUE, Long.MIN_VALUE)</code>,
    968         * <code>gcd(Long.MIN_VALUE, 0L)</code> and
    969         * <code>gcd(0L, Long.MIN_VALUE)</code> throw an
    970         * <code>ArithmeticException</code>, because the result would be 2^63, which
    971         * is too large for a long value.</li>
    972         * <li>The result of <code>gcd(x, x)</code>, <code>gcd(0L, x)</code> and
    973         * <code>gcd(x, 0L)</code> is the absolute value of <code>x</code>, except
    974         * for the special cases above.
    975         * <li>The invocation <code>gcd(0L, 0L)</code> is the only one which returns
    976         * <code>0L</code>.</li>
    977         * </ul>
    978         *
    979         * @param p any number
    980         * @param q any number
    981         * @return the greatest common divisor, never negative
    982         * @throws ArithmeticException if the result cannot be represented as a nonnegative long
    983         * value
    984         * @since 2.1
    985         */
    986        public static long gcd(final long p, final long q) {
    987            long u = p;
    988            long v = q;
    989            if ((u == 0) || (v == 0)) {
    990                if ((u == Long.MIN_VALUE) || (v == Long.MIN_VALUE)){
    991                    throw MathRuntimeException.createArithmeticException(
    992                            LocalizedFormats.GCD_OVERFLOW_64_BITS,
    993                            p, q);
    994                }
    995                return FastMath.abs(u) + FastMath.abs(v);
    996            }
    997            // keep u and v negative, as negative integers range down to
    998            // -2^63, while positive numbers can only be as large as 2^63-1
    999            // (i.e. we can't necessarily negate a negative number without
    1000            // overflow)
    1001            /* assert u!=0 && v!=0; */
    1002            if (u > 0) {
    1003                u = -u;
    1004            } // make u negative
    1005            if (v > 0) {
    1006                v = -v;
    1007            } // make v negative
    1008            // B1. [Find power of 2]
    1009            int k = 0;
    1010            while ((u & 1) == 0 && (v & 1) == 0 && k < 63) { // while u and v are
    1011                                                                // both even...
    1012                u /= 2;
    1013                v /= 2;
    1014                k++; // cast out twos.
    1015            }
    1016            if (k == 63) {
    1017                throw MathRuntimeException.createArithmeticException(
    1018                        LocalizedFormats.GCD_OVERFLOW_64_BITS,
    1019                        p, q);
    1020            }
    1021            // B2. Initialize: u and v have been divided by 2^k and at least
    1022            // one is odd.
    1023            long t = ((u & 1) == 1) ? v : -(u / 2)/* B3 */;
    1024            // t negative: u was odd, v may be even (t replaces v)
    1025            // t positive: u was even, v is odd (t replaces u)
    1026            do {
    1027                /* assert u<0 && v<0; */
    1028                // B4/B3: cast out twos from t.
    1029                while ((t & 1) == 0) { // while t is even..
    1030                    t /= 2; // cast out twos
    1031                }
    1032                // B5 [reset max(u,v)]
    1033                if (t > 0) {
    1034                    u = -t;
    1035                } else {
    1036                    v = t;
    1037                }
    1038                // B6/B3. at this point both u and v should be odd.
    1039                t = (v - u) / 2;
    1040                // |u| larger: t positive (replace u)
    1041                // |v| larger: t negative (replace v)
    1042            } while (t != 0);
    1043            return -u * (1L << k); // gcd is u*2^k
    1044        }
    1045
    1046        /**
    1047         * Returns an integer hash code representing the given double value.
    1048         *
    1049         * @param value the value to be hashed
    1050         * @return the hash code
    1051         */
    1052        public static int hash(double value) {
    1053            return new Double(value).hashCode();
    1054        }
    1055
    1056        /**
    1057         * Returns an integer hash code representing the given double array.
    1058         *
    1059         * @param value the value to be hashed (may be null)
    1060         * @return the hash code
    1061         * @since 1.2
    1062         */
    1063        public static int hash(double[] value) {
    1064            return Arrays.hashCode(value);
    1065        }
    1066
    1067        /**
    1068         * For a byte value x, this method returns (byte)(+1) if x >= 0 and
    1069         * (byte)(-1) if x < 0.
    1070         *
    1071         * @param x the value, a byte
    1072         * @return (byte)(+1) or (byte)(-1), depending on the sign of x
    1073         */
    1074        public static byte indicator(final byte x) {
    1075            return (x >= ZB) ? PB : NB;
    1076        }
    1077
    1078        /**
    1079         * For a double precision value x, this method returns +1.0 if x >= 0 and
    1080         * -1.0 if x < 0. Returns <code>NaN</code> if <code>x</code> is
    1081         * <code>NaN</code>.
    1082         *
    1083         * @param x the value, a double
    1084         * @return +1.0 or -1.0, depending on the sign of x
    1085         */
    1086        public static double indicator(final double x) {
    1087            if (Double.isNaN(x)) {
    1088                return Double.NaN;
    1089            }
    1090            return (x >= 0.0) ? 1.0 : -1.0;
    1091        }
    1092
    1093        /**
    1094         * For a float value x, this method returns +1.0F if x >= 0 and -1.0F if x <
    1095         * 0. Returns <code>NaN</code> if <code>x</code> is <code>NaN</code>.
    1096         *
    1097         * @param x the value, a float
    1098         * @return +1.0F or -1.0F, depending on the sign of x
    1099         */
    1100        public static float indicator(final float x) {
    1101            if (Float.isNaN(x)) {
    1102                return Float.NaN;
    1103            }
    1104            return (x >= 0.0F) ? 1.0F : -1.0F;
    1105        }
    1106
    1107        /**
    1108         * For an int value x, this method returns +1 if x >= 0 and -1 if x < 0.
    1109         *
    1110         * @param x the value, an int
    1111         * @return +1 or -1, depending on the sign of x
    1112         */
    1113        public static int indicator(final int x) {
    1114            return (x >= 0) ? 1 : -1;
    1115        }
    1116
    1117        /**
    1118         * For a long value x, this method returns +1L if x >= 0 and -1L if x < 0.
    1119         *
    1120         * @param x the value, a long
    1121         * @return +1L or -1L, depending on the sign of x
    1122         */
    1123        public static long indicator(final long x) {
    1124            return (x >= 0L) ? 1L : -1L;
    1125        }
    1126
    1127        /**
    1128         * For a short value x, this method returns (short)(+1) if x >= 0 and
    1129         * (short)(-1) if x < 0.
    1130         *
    1131         * @param x the value, a short
    1132         * @return (short)(+1) or (short)(-1), depending on the sign of x
    1133         */
    1134        public static short indicator(final short x) {
    1135            return (x >= ZS) ? PS : NS;
    1136        }
    1137
    1138        /**
    1139         * <p>
    1140         * Returns the least common multiple of the absolute value of two numbers,
    1141         * using the formula <code>lcm(a,b) = (a / gcd(a,b)) * b</code>.
    1142         * </p>
    1143         * Special cases:
    1144         * <ul>
    1145         * <li>The invocations <code>lcm(Integer.MIN_VALUE, n)</code> and
    1146         * <code>lcm(n, Integer.MIN_VALUE)</code>, where <code>abs(n)</code> is a
    1147         * power of 2, throw an <code>ArithmeticException</code>, because the result
    1148         * would be 2^31, which is too large for an int value.</li>
    1149         * <li>The result of <code>lcm(0, x)</code> and <code>lcm(x, 0)</code> is
    1150         * <code>0</code> for any <code>x</code>.
    1151         * </ul>
    1152         *
    1153         * @param a any number
    1154         * @param b any number
    1155         * @return the least common multiple, never negative
    1156         * @throws ArithmeticException
    1157         *             if the result cannot be represented as a nonnegative int
    1158         *             value
    1159         * @since 1.1
    1160         */
    1161        public static int lcm(int a, int b) {
    1162            if (a==0 || b==0){
    1163                return 0;
    1164            }
    1165            int lcm = FastMath.abs(mulAndCheck(a / gcd(a, b), b));
    1166            if (lcm == Integer.MIN_VALUE) {
    1167                throw MathRuntimeException.createArithmeticException(
    1168                    LocalizedFormats.LCM_OVERFLOW_32_BITS,
    1169                    a, b);
    1170            }
    1171            return lcm;
    1172        }
    1173
    1174        /**
    1175         * <p>
    1176         * Returns the least common multiple of the absolute value of two numbers,
    1177         * using the formula <code>lcm(a,b) = (a / gcd(a,b)) * b</code>.
    1178         * </p>
    1179         * Special cases:
    1180         * <ul>
    1181         * <li>The invocations <code>lcm(Long.MIN_VALUE, n)</code> and
    1182         * <code>lcm(n, Long.MIN_VALUE)</code>, where <code>abs(n)</code> is a
    1183         * power of 2, throw an <code>ArithmeticException</code>, because the result
    1184         * would be 2^63, which is too large for an int value.</li>
    1185         * <li>The result of <code>lcm(0L, x)</code> and <code>lcm(x, 0L)</code> is
    1186         * <code>0L</code> for any <code>x</code>.
    1187         * </ul>
    1188         *
    1189         * @param a any number
    1190         * @param b any number
    1191         * @return the least common multiple, never negative
    1192         * @throws ArithmeticException if the result cannot be represented as a nonnegative long
    1193         * value
    1194         * @since 2.1
    1195         */
    1196        public static long lcm(long a, long b) {
    1197            if (a==0 || b==0){
    1198                return 0;
    1199            }
    1200            long lcm = FastMath.abs(mulAndCheck(a / gcd(a, b), b));
    1201            if (lcm == Long.MIN_VALUE){
    1202                throw MathRuntimeException.createArithmeticException(
    1203                    LocalizedFormats.LCM_OVERFLOW_64_BITS,
    1204                    a, b);
    1205            }
    1206            return lcm;
    1207        }
    1208
    1209        /**
    1210         * <p>Returns the
    1211         * <a href="http://mathworld.wolfram.com/Logarithm.html">logarithm</a>
    1212         * for base <code>b</code> of <code>x</code>.
    1213         * </p>
    1214         * <p>Returns <code>NaN<code> if either argument is negative.  If
    1215         * <code>base</code> is 0 and <code>x</code> is positive, 0 is returned.
    1216         * If <code>base</code> is positive and <code>x</code> is 0,
    1217         * <code>Double.NEGATIVE_INFINITY</code> is returned.  If both arguments
    1218         * are 0, the result is <code>NaN</code>.</p>
    1219         *
    1220         * @param base the base of the logarithm, must be greater than 0
    1221         * @param x argument, must be greater than 0
    1222         * @return the value of the logarithm - the number y such that base^y = x.
    1223         * @since 1.2
    1224         */
    1225        public static double log(double base, double x) {
    1226            return FastMath.log(x)/FastMath.log(base);
    1227        }
    1228
    1229        /**
    1230         * Multiply two integers, checking for overflow.
    1231         *
    1232         * @param x a factor
    1233         * @param y a factor
    1234         * @return the product <code>x*y</code>
    1235         * @throws ArithmeticException if the result can not be represented as an
    1236         *         int
    1237         * @since 1.1
    1238         */
    1239        public static int mulAndCheck(int x, int y) {
    1240            long m = ((long)x) * ((long)y);
    1241            if (m < Integer.MIN_VALUE || m > Integer.MAX_VALUE) {
    1242                throw new ArithmeticException("overflow: mul");
    1243            }
    1244            return (int)m;
    1245        }
    1246
    1247        /**
    1248         * Multiply two long integers, checking for overflow.
    1249         *
    1250         * @param a first value
    1251         * @param b second value
    1252         * @return the product <code>a * b</code>
    1253         * @throws ArithmeticException if the result can not be represented as an
    1254         *         long
    1255         * @since 1.2
    1256         */
    1257        public static long mulAndCheck(long a, long b) {
    1258            long ret;
    1259            String msg = "overflow: multiply";
    1260            if (a > b) {
    1261                // use symmetry to reduce boundary cases
    1262                ret = mulAndCheck(b, a);
    1263            } else {
    1264                if (a < 0) {
    1265                    if (b < 0) {
    1266                        // check for positive overflow with negative a, negative b
    1267                        if (a >= Long.MAX_VALUE / b) {
    1268                            ret = a * b;
    1269                        } else {
    1270                            throw new ArithmeticException(msg);
    1271                        }
    1272                    } else if (b > 0) {
    1273                        // check for negative overflow with negative a, positive b
    1274                        if (Long.MIN_VALUE / b <= a) {
    1275                            ret = a * b;
    1276                        } else {
    1277                            throw new ArithmeticException(msg);
    1278
    1279                        }
    1280                    } else {
    1281                        // assert b == 0
    1282                        ret = 0;
    1283                    }
    1284                } else if (a > 0) {
    1285                    // assert a > 0
    1286                    // assert b > 0
    1287
    1288                    // check for positive overflow with positive a, positive b
    1289                    if (a <= Long.MAX_VALUE / b) {
    1290                        ret = a * b;
    1291                    } else {
    1292                        throw new ArithmeticException(msg);
    1293                    }
    1294                } else {
    1295                    // assert a == 0
    1296                    ret = 0;
    1297                }
    1298            }
    1299            return ret;
    1300        }
    1301
    1302        /**
    1303         * Get the next machine representable number after a number, moving
    1304         * in the direction of another number.
    1305         * <p>
    1306         * If <code>direction</code> is greater than or equal to<code>d</code>,
    1307         * the smallest machine representable number strictly greater than
    1308         * <code>d</code> is returned; otherwise the largest representable number
    1309         * strictly less than <code>d</code> is returned.</p>
    1310         * <p>
    1311         * If <code>d</code> is NaN or Infinite, it is returned unchanged.</p>
    1312         *
    1313         * @param d base number
    1314         * @param direction (the only important thing is whether
    1315         * direction is greater or smaller than d)
    1316         * @return the next machine representable number in the specified direction
    1317         * @since 1.2
    1318         * @deprecated as of 2.2, replaced by {@link FastMath#nextAfter(double, double)}
    1319         * which handles Infinities differently, and returns direction if d and direction compare equal.
    1320         */
    1321        @Deprecated
    1322        public static double nextAfter(double d, double direction) {
    1323
    1324            // handling of some important special cases
    1325            if (Double.isNaN(d) || Double.isInfinite(d)) {
    1326                    return d;
    1327            } else if (d == 0) {
    1328                    return (direction < 0) ? -Double.MIN_VALUE : Double.MIN_VALUE;
    1329            }
    1330            // special cases MAX_VALUE to infinity and  MIN_VALUE to 0
    1331            // are handled just as normal numbers
    1332
    1333            // split the double in raw components
    1334            long bits     = Double.doubleToLongBits(d);
    1335            long sign     = bits & 0x8000000000000000L;
    1336            long exponent = bits & 0x7ff0000000000000L;
    1337            long mantissa = bits & 0x000fffffffffffffL;
    1338
    1339            if (d * (direction - d) >= 0) {
    1340                    // we should increase the mantissa
    1341                    if (mantissa == 0x000fffffffffffffL) {
    1342                            return Double.longBitsToDouble(sign |
    1343                                            (exponent + 0x0010000000000000L));
    1344                    } else {
    1345                            return Double.longBitsToDouble(sign |
    1346                                            exponent | (mantissa + 1));
    1347                    }
    1348            } else {
    1349                    // we should decrease the mantissa
    1350                    if (mantissa == 0L) {
    1351                            return Double.longBitsToDouble(sign |
    1352                                            (exponent - 0x0010000000000000L) |
    1353                                            0x000fffffffffffffL);
    1354                    } else {
    1355                            return Double.longBitsToDouble(sign |
    1356                                            exponent | (mantissa - 1));
    1357                    }
    1358            }
    1359        }
    1360
    1361        /**
    1362         * Scale a number by 2<sup>scaleFactor</sup>.
    1363         * <p>If <code>d</code> is 0 or NaN or Infinite, it is returned unchanged.</p>
    1364         *
    1365         * @param d base number
    1366         * @param scaleFactor power of two by which d should be multiplied
    1367         * @return d &times; 2<sup>scaleFactor</sup>
    1368         * @since 2.0
    1369         * @deprecated as of 2.2, replaced by {@link FastMath#scalb(double, int)}
    1370         */
    1371        @Deprecated
    1372        public static double scalb(final double d, final int scaleFactor) {
    1373            return FastMath.scalb(d, scaleFactor);
    1374        }
    1375
    1376        /**
    1377         * Normalize an angle in a 2&pi wide interval around a center value.
    1378         * <p>This method has three main uses:</p>
    1379         * <ul>
    1380         *   <li>normalize an angle between 0 and 2&pi;:<br/>
    1381         *       <code>a = MathUtils.normalizeAngle(a, FastMath.PI);</code></li>
    1382         *   <li>normalize an angle between -&pi; and +&pi;<br/>
    1383         *       <code>a = MathUtils.normalizeAngle(a, 0.0);</code></li>
    1384         *   <li>compute the angle between two defining angular positions:<br>
    1385         *       <code>angle = MathUtils.normalizeAngle(end, start) - start;</code></li>
    1386         * </ul>
    1387         * <p>Note that due to numerical accuracy and since &pi; cannot be represented
    1388         * exactly, the result interval is <em>closed</em>, it cannot be half-closed
    1389         * as would be more satisfactory in a purely mathematical view.</p>
    1390         * @param a angle to normalize
    1391         * @param center center of the desired 2&pi; interval for the result
    1392         * @return a-2k&pi; with integer k and center-&pi; &lt;= a-2k&pi; &lt;= center+&pi;
    1393         * @since 1.2
    1394         */
    1395         public static double normalizeAngle(double a, double center) {
    1396             return a - TWO_PI * FastMath.floor((a + FastMath.PI - center) / TWO_PI);
    1397         }
    1398
    1399         /**
    1400          * <p>Normalizes an array to make it sum to a specified value.
    1401          * Returns the result of the transformation <pre>
    1402          *    x |-> x * normalizedSum / sum
    1403          * </pre>
    1404          * applied to each non-NaN element x of the input array, where sum is the
    1405          * sum of the non-NaN entries in the input array.</p>
    1406          *
    1407          * <p>Throws IllegalArgumentException if <code>normalizedSum</code> is infinite
    1408          * or NaN and ArithmeticException if the input array contains any infinite elements
    1409          * or sums to 0</p>
    1410          *
    1411          * <p>Ignores (i.e., copies unchanged to the output array) NaNs in the input array.</p>
    1412          *
    1413          * @param values input array to be normalized
    1414          * @param normalizedSum target sum for the normalized array
    1415          * @return normalized array
    1416          * @throws ArithmeticException if the input array contains infinite elements or sums to zero
    1417          * @throws IllegalArgumentException if the target sum is infinite or NaN
    1418          * @since 2.1
    1419          */
    1420         public static double[] normalizeArray(double[] values, double normalizedSum)
    1421           throws ArithmeticException, IllegalArgumentException {
    1422             if (Double.isInfinite(normalizedSum)) {
    1423                 throw MathRuntimeException.createIllegalArgumentException(
    1424                         LocalizedFormats.NORMALIZE_INFINITE);
    1425             }
    1426             if (Double.isNaN(normalizedSum)) {
    1427                 throw MathRuntimeException.createIllegalArgumentException(
    1428                         LocalizedFormats.NORMALIZE_NAN);
    1429             }
    1430             double sum = 0d;
    1431             final int len = values.length;
    1432             double[] out = new double[len];
    1433             for (int i = 0; i < len; i++) {
    1434                 if (Double.isInfinite(values[i])) {
    1435                     throw MathRuntimeException.createArithmeticException(
    1436                             LocalizedFormats.INFINITE_ARRAY_ELEMENT, values[i], i);
    1437                 }
    1438                 if (!Double.isNaN(values[i])) {
    1439                     sum += values[i];
    1440                 }
    1441             }
    1442             if (sum == 0) {
    1443                 throw MathRuntimeException.createArithmeticException(LocalizedFormats.ARRAY_SUMS_TO_ZERO);
    1444             }
    1445             for (int i = 0; i < len; i++) {
    1446                 if (Double.isNaN(values[i])) {
    1447                     out[i] = Double.NaN;
    1448                 } else {
    1449                     out[i] = values[i] * normalizedSum / sum;
    1450                 }
    1451             }
    1452             return out;
    1453         }
    1454
    1455        /**
    1456         * Round the given value to the specified number of decimal places. The
    1457         * value is rounded using the {@link BigDecimal#ROUND_HALF_UP} method.
    1458         *
    1459         * @param x the value to round.
    1460         * @param scale the number of digits to the right of the decimal point.
    1461         * @return the rounded value.
    1462         * @since 1.1
    1463         */
    1464        public static double round(double x, int scale) {
    1465            return round(x, scale, BigDecimal.ROUND_HALF_UP);
    1466        }
    1467
    1468        /**
    1469         * Round the given value to the specified number of decimal places. The
    1470         * value is rounded using the given method which is any method defined in
    1471         * {@link BigDecimal}.
    1472         *
    1473         * @param x the value to round.
    1474         * @param scale the number of digits to the right of the decimal point.
    1475         * @param roundingMethod the rounding method as defined in
    1476         *        {@link BigDecimal}.
    1477         * @return the rounded value.
    1478         * @since 1.1
    1479         */
    1480        public static double round(double x, int scale, int roundingMethod) {
    1481            try {
    1482                return (new BigDecimal
    1483                       (Double.toString(x))
    1484                       .setScale(scale, roundingMethod))
    1485                       .doubleValue();
    1486            } catch (NumberFormatException ex) {
    1487                if (Double.isInfinite(x)) {
    1488                    return x;
    1489                } else {
    1490                    return Double.NaN;
    1491                }
    1492            }
    1493        }
    1494
    1495        /**
    1496         * Round the given value to the specified number of decimal places. The
    1497         * value is rounding using the {@link BigDecimal#ROUND_HALF_UP} method.
    1498         *
    1499         * @param x the value to round.
    1500         * @param scale the number of digits to the right of the decimal point.
    1501         * @return the rounded value.
    1502         * @since 1.1
    1503         */
    1504        public static float round(float x, int scale) {
    1505            return round(x, scale, BigDecimal.ROUND_HALF_UP);
    1506        }
    1507
    1508        /**
    1509         * Round the given value to the specified number of decimal places. The
    1510         * value is rounded using the given method which is any method defined in
    1511         * {@link BigDecimal}.
    1512         *
    1513         * @param x the value to round.
    1514         * @param scale the number of digits to the right of the decimal point.
    1515         * @param roundingMethod the rounding method as defined in
    1516         *        {@link BigDecimal}.
    1517         * @return the rounded value.
    1518         * @since 1.1
    1519         */
    1520        public static float round(float x, int scale, int roundingMethod) {
    1521            float sign = indicator(x);
    1522            float factor = (float)FastMath.pow(10.0f, scale) * sign;
    1523            return (float)roundUnscaled(x * factor, sign, roundingMethod) / factor;
    1524        }
    1525
    1526        /**
    1527         * Round the given non-negative, value to the "nearest" integer. Nearest is
    1528         * determined by the rounding method specified. Rounding methods are defined
    1529         * in {@link BigDecimal}.
    1530         *
    1531         * @param unscaled the value to round.
    1532         * @param sign the sign of the original, scaled value.
    1533         * @param roundingMethod the rounding method as defined in
    1534         *        {@link BigDecimal}.
    1535         * @return the rounded value.
    1536         * @since 1.1
    1537         */
    1538        private static double roundUnscaled(double unscaled, double sign,
    1539            int roundingMethod) {
    1540            switch (roundingMethod) {
    1541            case BigDecimal.ROUND_CEILING :
    1542                if (sign == -1) {
    1543                    unscaled = FastMath.floor(nextAfter(unscaled, Double.NEGATIVE_INFINITY));
    1544                } else {
    1545                    unscaled = FastMath.ceil(nextAfter(unscaled, Double.POSITIVE_INFINITY));
    1546                }
    1547                break;
    1548            case BigDecimal.ROUND_DOWN :
    1549                unscaled = FastMath.floor(nextAfter(unscaled, Double.NEGATIVE_INFINITY));
    1550                break;
    1551            case BigDecimal.ROUND_FLOOR :
    1552                if (sign == -1) {
    1553                    unscaled = FastMath.ceil(nextAfter(unscaled, Double.POSITIVE_INFINITY));
    1554                } else {
    1555                    unscaled = FastMath.floor(nextAfter(unscaled, Double.NEGATIVE_INFINITY));
    1556                }
    1557                break;
    1558            case BigDecimal.ROUND_HALF_DOWN : {
    1559                unscaled = nextAfter(unscaled, Double.NEGATIVE_INFINITY);
    1560                double fraction = unscaled - FastMath.floor(unscaled);
    1561                if (fraction > 0.5) {
    1562                    unscaled = FastMath.ceil(unscaled);
    1563                } else {
    1564                    unscaled = FastMath.floor(unscaled);
    1565                }
    1566                break;
    1567            }
    1568            case BigDecimal.ROUND_HALF_EVEN : {
    1569                double fraction = unscaled - FastMath.floor(unscaled);
    1570                if (fraction > 0.5) {
    1571                    unscaled = FastMath.ceil(unscaled);
    1572                } else if (fraction < 0.5) {
    1573                    unscaled = FastMath.floor(unscaled);
    1574                } else {
    1575                    // The following equality test is intentional and needed for rounding purposes
    1576                    if (FastMath.floor(unscaled) / 2.0 == FastMath.floor(Math
    1577                        .floor(unscaled) / 2.0)) { // even
    1578                        unscaled = FastMath.floor(unscaled);
    1579                    } else { // odd
    1580                        unscaled = FastMath.ceil(unscaled);
    1581                    }
    1582                }
    1583                break;
    1584            }
    1585            case BigDecimal.ROUND_HALF_UP : {
    1586                unscaled = nextAfter(unscaled, Double.POSITIVE_INFINITY);
    1587                double fraction = unscaled - FastMath.floor(unscaled);
    1588                if (fraction >= 0.5) {
    1589                    unscaled = FastMath.ceil(unscaled);
    1590                } else {
    1591                    unscaled = FastMath.floor(unscaled);
    1592                }
    1593                break;
    1594            }
    1595            case BigDecimal.ROUND_UNNECESSARY :
    1596                if (unscaled != FastMath.floor(unscaled)) {
    1597                    throw new ArithmeticException("Inexact result from rounding");
    1598                }
    1599                break;
    1600            case BigDecimal.ROUND_UP :
    1601                unscaled = FastMath.ceil(nextAfter(unscaled,  Double.POSITIVE_INFINITY));
    1602                break;
    1603            default :
    1604                throw MathRuntimeException.createIllegalArgumentException(
    1605                      LocalizedFormats.INVALID_ROUNDING_METHOD,
    1606                      roundingMethod,
    1607                      "ROUND_CEILING",     BigDecimal.ROUND_CEILING,
    1608                      "ROUND_DOWN",        BigDecimal.ROUND_DOWN,
    1609                      "ROUND_FLOOR",       BigDecimal.ROUND_FLOOR,
    1610                      "ROUND_HALF_DOWN",   BigDecimal.ROUND_HALF_DOWN,
    1611                      "ROUND_HALF_EVEN",   BigDecimal.ROUND_HALF_EVEN,
    1612                      "ROUND_HALF_UP",     BigDecimal.ROUND_HALF_UP,
    1613                      "ROUND_UNNECESSARY", BigDecimal.ROUND_UNNECESSARY,
    1614                      "ROUND_UP",          BigDecimal.ROUND_UP);
    1615            }
    1616            return unscaled;
    1617        }
    1618
    1619        /**
    1620         * Returns the <a href="http://mathworld.wolfram.com/Sign.html"> sign</a>
    1621         * for byte value <code>x</code>.
    1622         * <p>
    1623         * For a byte value x, this method returns (byte)(+1) if x > 0, (byte)(0) if
    1624         * x = 0, and (byte)(-1) if x < 0.</p>
    1625         *
    1626         * @param x the value, a byte
    1627         * @return (byte)(+1), (byte)(0), or (byte)(-1), depending on the sign of x
    1628         */
    1629        public static byte sign(final byte x) {
    1630            return (x == ZB) ? ZB : (x > ZB) ? PB : NB;
    1631        }
    1632
    1633        /**
    1634         * Returns the <a href="http://mathworld.wolfram.com/Sign.html"> sign</a>
    1635         * for double precision <code>x</code>.
    1636         * <p>
    1637         * For a double value <code>x</code>, this method returns
    1638         * <code>+1.0</code> if <code>x > 0</code>, <code>0.0</code> if
    1639         * <code>x = 0.0</code>, and <code>-1.0</code> if <code>x < 0</code>.
    1640         * Returns <code>NaN</code> if <code>x</code> is <code>NaN</code>.</p>
    1641         *
    1642         * @param x the value, a double
    1643         * @return +1.0, 0.0, or -1.0, depending on the sign of x
    1644         */
    1645        public static double sign(final double x) {
    1646            if (Double.isNaN(x)) {
    1647                return Double.NaN;
    1648            }
    1649            return (x == 0.0) ? 0.0 : (x > 0.0) ? 1.0 : -1.0;
    1650        }
    1651
    1652        /**
    1653         * Returns the <a href="http://mathworld.wolfram.com/Sign.html"> sign</a>
    1654         * for float value <code>x</code>.
    1655         * <p>
    1656         * For a float value x, this method returns +1.0F if x > 0, 0.0F if x =
    1657         * 0.0F, and -1.0F if x < 0. Returns <code>NaN</code> if <code>x</code>
    1658         * is <code>NaN</code>.</p>
    1659         *
    1660         * @param x the value, a float
    1661         * @return +1.0F, 0.0F, or -1.0F, depending on the sign of x
    1662         */
    1663        public static float sign(final float x) {
    1664            if (Float.isNaN(x)) {
    1665                return Float.NaN;
    1666            }
    1667            return (x == 0.0F) ? 0.0F : (x > 0.0F) ? 1.0F : -1.0F;
    1668        }
    1669
    1670        /**
    1671         * Returns the <a href="http://mathworld.wolfram.com/Sign.html"> sign</a>
    1672         * for int value <code>x</code>.
    1673         * <p>
    1674         * For an int value x, this method returns +1 if x > 0, 0 if x = 0, and -1
    1675         * if x < 0.</p>
    1676         *
    1677         * @param x the value, an int
    1678         * @return +1, 0, or -1, depending on the sign of x
    1679         */
    1680        public static int sign(final int x) {
    1681            return (x == 0) ? 0 : (x > 0) ? 1 : -1;
    1682        }
    1683
    1684        /**
    1685         * Returns the <a href="http://mathworld.wolfram.com/Sign.html"> sign</a>
    1686         * for long value <code>x</code>.
    1687         * <p>
    1688         * For a long value x, this method returns +1L if x > 0, 0L if x = 0, and
    1689         * -1L if x < 0.</p>
    1690         *
    1691         * @param x the value, a long
    1692         * @return +1L, 0L, or -1L, depending on the sign of x
    1693         */
    1694        public static long sign(final long x) {
    1695            return (x == 0L) ? 0L : (x > 0L) ? 1L : -1L;
    1696        }
    1697
    1698        /**
    1699         * Returns the <a href="http://mathworld.wolfram.com/Sign.html"> sign</a>
    1700         * for short value <code>x</code>.
    1701         * <p>
    1702         * For a short value x, this method returns (short)(+1) if x > 0, (short)(0)
    1703         * if x = 0, and (short)(-1) if x < 0.</p>
    1704         *
    1705         * @param x the value, a short
    1706         * @return (short)(+1), (short)(0), or (short)(-1), depending on the sign of
    1707         *         x
    1708         */
    1709        public static short sign(final short x) {
    1710            return (x == ZS) ? ZS : (x > ZS) ? PS : NS;
    1711        }
    1712
    1713        /**
    1714         * Returns the <a href="http://mathworld.wolfram.com/HyperbolicSine.html">
    1715         * hyperbolic sine</a> of x.
    1716         *
    1717         * @param x double value for which to find the hyperbolic sine
    1718         * @return hyperbolic sine of x
    1719         */
    1720        public static double sinh(double x) {
    1721            return (FastMath.exp(x) - FastMath.exp(-x)) / 2.0;
    1722        }
    1723
    1724        /**
    1725         * Subtract two integers, checking for overflow.
    1726         *
    1727         * @param x the minuend
    1728         * @param y the subtrahend
    1729         * @return the difference <code>x-y</code>
    1730         * @throws ArithmeticException if the result can not be represented as an
    1731         *         int
    1732         * @since 1.1
    1733         */
    1734        public static int subAndCheck(int x, int y) {
    1735            long s = (long)x - (long)y;
    1736            if (s < Integer.MIN_VALUE || s > Integer.MAX_VALUE) {
    1737                throw MathRuntimeException.createArithmeticException(LocalizedFormats.OVERFLOW_IN_SUBTRACTION, x, y);
    1738            }
    1739            return (int)s;
    1740        }
    1741
    1742        /**
    1743         * Subtract two long integers, checking for overflow.
    1744         *
    1745         * @param a first value
    1746         * @param b second value
    1747         * @return the difference <code>a-b</code>
    1748         * @throws ArithmeticException if the result can not be represented as an
    1749         *         long
    1750         * @since 1.2
    1751         */
    1752        public static long subAndCheck(long a, long b) {
    1753            long ret;
    1754            String msg = "overflow: subtract";
    1755            if (b == Long.MIN_VALUE) {
    1756                if (a < 0) {
    1757                    ret = a - b;
    1758                } else {
    1759                    throw new ArithmeticException(msg);
    1760                }
    1761            } else {
    1762                // use additive inverse
    1763                ret = addAndCheck(a, -b, LocalizedFormats.OVERFLOW_IN_ADDITION);
    1764            }
    1765            return ret;
    1766        }
    1767
    1768        /**
    1769         * Raise an int to an int power.
    1770         * @param k number to raise
    1771         * @param e exponent (must be positive or null)
    1772         * @return k<sup>e</sup>
    1773         * @exception IllegalArgumentException if e is negative
    1774         */
    1775        public static int pow(final int k, int e)
    1776            throws IllegalArgumentException {
    1777
    1778            if (e < 0) {
    1779                throw MathRuntimeException.createIllegalArgumentException(
    1780                    LocalizedFormats.POWER_NEGATIVE_PARAMETERS,
    1781                    k, e);
    1782            }
    1783
    1784            int result = 1;
    1785            int k2p    = k;
    1786            while (e != 0) {
    1787                if ((e & 0x1) != 0) {
    1788                    result *= k2p;
    1789                }
    1790                k2p *= k2p;
    1791                e = e >> 1;
    1792            }
    1793
    1794            return result;
    1795
    1796        }
    1797
    1798        /**
    1799         * Raise an int to a long power.
    1800         * @param k number to raise
    1801         * @param e exponent (must be positive or null)
    1802         * @return k<sup>e</sup>
    1803         * @exception IllegalArgumentException if e is negative
    1804         */
    1805        public static int pow(final int k, long e)
    1806            throws IllegalArgumentException {
    1807
    1808            if (e < 0) {
    1809                throw MathRuntimeException.createIllegalArgumentException(
    1810                    LocalizedFormats.POWER_NEGATIVE_PARAMETERS,
    1811                    k, e);
    1812            }
    1813
    1814            int result = 1;
    1815            int k2p    = k;
    1816            while (e != 0) {
    1817                if ((e & 0x1) != 0) {
    1818                    result *= k2p;
    1819                }
    1820                k2p *= k2p;
    1821                e = e >> 1;
    1822            }
    1823
    1824            return result;
    1825
    1826        }
    1827
    1828        /**
    1829         * Raise a long to an int power.
    1830         * @param k number to raise
    1831         * @param e exponent (must be positive or null)
    1832         * @return k<sup>e</sup>
    1833         * @exception IllegalArgumentException if e is negative
    1834         */
    1835        public static long pow(final long k, int e)
    1836            throws IllegalArgumentException {
    1837
    1838            if (e < 0) {
    1839                throw MathRuntimeException.createIllegalArgumentException(
    1840                    LocalizedFormats.POWER_NEGATIVE_PARAMETERS,
    1841                    k, e);
    1842            }
    1843
    1844            long result = 1l;
    1845            long k2p    = k;
    1846            while (e != 0) {
    1847                if ((e & 0x1) != 0) {
    1848                    result *= k2p;
    1849                }
    1850                k2p *= k2p;
    1851                e = e >> 1;
    1852            }
    1853
    1854            return result;
    1855
    1856        }
    1857
    1858        /**
    1859         * Raise a long to a long power.
    1860         * @param k number to raise
    1861         * @param e exponent (must be positive or null)
    1862         * @return k<sup>e</sup>
    1863         * @exception IllegalArgumentException if e is negative
    1864         */
    1865        public static long pow(final long k, long e)
    1866            throws IllegalArgumentException {
    1867
    1868            if (e < 0) {
    1869                throw MathRuntimeException.createIllegalArgumentException(
    1870                    LocalizedFormats.POWER_NEGATIVE_PARAMETERS,
    1871                    k, e);
    1872            }
    1873
    1874            long result = 1l;
    1875            long k2p    = k;
    1876            while (e != 0) {
    1877                if ((e & 0x1) != 0) {
    1878                    result *= k2p;
    1879                }
    1880                k2p *= k2p;
    1881                e = e >> 1;
    1882            }
    1883
    1884            return result;
    1885
    1886        }
    1887
    1888        /**
    1889         * Raise a BigInteger to an int power.
    1890         * @param k number to raise
    1891         * @param e exponent (must be positive or null)
    1892         * @return k<sup>e</sup>
    1893         * @exception IllegalArgumentException if e is negative
    1894         */
    1895        public static BigInteger pow(final BigInteger k, int e)
    1896            throws IllegalArgumentException {
    1897
    1898            if (e < 0) {
    1899                throw MathRuntimeException.createIllegalArgumentException(
    1900                    LocalizedFormats.POWER_NEGATIVE_PARAMETERS,
    1901                    k, e);
    1902            }
    1903
    1904            return k.pow(e);
    1905
    1906        }
    1907
    1908        /**
    1909         * Raise a BigInteger to a long power.
    1910         * @param k number to raise
    1911         * @param e exponent (must be positive or null)
    1912         * @return k<sup>e</sup>
    1913         * @exception IllegalArgumentException if e is negative
    1914         */
    1915        public static BigInteger pow(final BigInteger k, long e)
    1916            throws IllegalArgumentException {
    1917
    1918            if (e < 0) {
    1919                throw MathRuntimeException.createIllegalArgumentException(
    1920                    LocalizedFormats.POWER_NEGATIVE_PARAMETERS,
    1921                    k, e);
    1922            }
    1923
    1924            BigInteger result = BigInteger.ONE;
    1925            BigInteger k2p    = k;
    1926            while (e != 0) {
    1927                if ((e & 0x1) != 0) {
    1928                    result = result.multiply(k2p);
    1929                }
    1930                k2p = k2p.multiply(k2p);
    1931                e = e >> 1;
    1932            }
    1933
    1934            return result;
    1935
    1936        }
    1937
    1938        /**
    1939         * Raise a BigInteger to a BigInteger power.
    1940         * @param k number to raise
    1941         * @param e exponent (must be positive or null)
    1942         * @return k<sup>e</sup>
    1943         * @exception IllegalArgumentException if e is negative
    1944         */
    1945        public static BigInteger pow(final BigInteger k, BigInteger e)
    1946            throws IllegalArgumentException {
    1947
    1948            if (e.compareTo(BigInteger.ZERO) < 0) {
    1949                throw MathRuntimeException.createIllegalArgumentException(
    1950                    LocalizedFormats.POWER_NEGATIVE_PARAMETERS,
    1951                    k, e);
    1952            }
    1953
    1954            BigInteger result = BigInteger.ONE;
    1955            BigInteger k2p    = k;
    1956            while (!BigInteger.ZERO.equals(e)) {
    1957                if (e.testBit(0)) {
    1958                    result = result.multiply(k2p);
    1959                }
    1960                k2p = k2p.multiply(k2p);
    1961                e = e.shiftRight(1);
    1962            }
    1963
    1964            return result;
    1965
    1966        }
    1967
    1968        /**
    1969         * Calculates the L<sub>1</sub> (sum of abs) distance between two points.
    1970         *
    1971         * @param p1 the first point
    1972         * @param p2 the second point
    1973         * @return the L<sub>1</sub> distance between the two points
    1974         */
    1975        public static double distance1(double[] p1, double[] p2) {
    1976            double sum = 0;
    1977            for (int i = 0; i < p1.length; i++) {
    1978                sum += FastMath.abs(p1[i] - p2[i]);
    1979            }
    1980            return sum;
    1981        }
    1982
    1983        /**
    1984         * Calculates the L<sub>1</sub> (sum of abs) distance between two points.
    1985         *
    1986         * @param p1 the first point
    1987         * @param p2 the second point
    1988         * @return the L<sub>1</sub> distance between the two points
    1989         */
    1990        public static int distance1(int[] p1, int[] p2) {
    1991          int sum = 0;
    1992          for (int i = 0; i < p1.length; i++) {
    1993              sum += FastMath.abs(p1[i] - p2[i]);
    1994          }
    1995          return sum;
    1996        }
    1997
    1998        /**
    1999         * Calculates the L<sub>2</sub> (Euclidean) distance between two points.
    2000         *
    2001         * @param p1 the first point
    2002         * @param p2 the second point
    2003         * @return the L<sub>2</sub> distance between the two points
    2004         */
    2005        public static double distance(double[] p1, double[] p2) {
    2006            double sum = 0;
    2007            for (int i = 0; i < p1.length; i++) {
    2008                final double dp = p1[i] - p2[i];
    2009                sum += dp * dp;
    2010            }
    2011            return FastMath.sqrt(sum);
    2012        }
    2013
    2014        /**
    2015         * Calculates the L<sub>2</sub> (Euclidean) distance between two points.
    2016         *
    2017         * @param p1 the first point
    2018         * @param p2 the second point
    2019         * @return the L<sub>2</sub> distance between the two points
    2020         */
    2021        public static double distance(int[] p1, int[] p2) {
    2022          double sum = 0;
    2023          for (int i = 0; i < p1.length; i++) {
    2024              final double dp = p1[i] - p2[i];
    2025              sum += dp * dp;
    2026          }
    2027          return FastMath.sqrt(sum);
    2028        }
    2029
    2030        /**
    2031         * Calculates the L<sub>&infin;</sub> (max of abs) distance between two points.
    2032         *
    2033         * @param p1 the first point
    2034         * @param p2 the second point
    2035         * @return the L<sub>&infin;</sub> distance between the two points
    2036         */
    2037        public static double distanceInf(double[] p1, double[] p2) {
    2038            double max = 0;
    2039            for (int i = 0; i < p1.length; i++) {
    2040                max = FastMath.max(max, FastMath.abs(p1[i] - p2[i]));
    2041            }
    2042            return max;
    2043        }
    2044
    2045        /**
    2046         * Calculates the L<sub>&infin;</sub> (max of abs) distance between two points.
    2047         *
    2048         * @param p1 the first point
    2049         * @param p2 the second point
    2050         * @return the L<sub>&infin;</sub> distance between the two points
    2051         */
    2052        public static int distanceInf(int[] p1, int[] p2) {
    2053            int max = 0;
    2054            for (int i = 0; i < p1.length; i++) {
    2055                max = FastMath.max(max, FastMath.abs(p1[i] - p2[i]));
    2056            }
    2057            return max;
    2058        }
    2059
    2060        /**
    2061         * Specification of ordering direction.
    2062         */
    2063        public static enum OrderDirection {
    2064            /** Constant for increasing direction. */
    2065            INCREASING,
    2066            /** Constant for decreasing direction. */
    2067            DECREASING
    2068        }
    2069
    2070        /**
    2071         * Checks that the given array is sorted.
    2072         *
    2073         * @param val Values.
    2074         * @param dir Ordering direction.
    2075         * @param strict Whether the order should be strict.
    2076         * @throws NonMonotonousSequenceException if the array is not sorted.
    2077         * @since 2.2
    2078         */
    2079        public static void checkOrder(double[] val, OrderDirection dir, boolean strict) {
    2080            double previous = val[0];
    2081            boolean ok = true;
    2082
    2083            int max = val.length;
    2084            for (int i = 1; i < max; i++) {
    2085                switch (dir) {
    2086                case INCREASING:
    2087                    if (strict) {
    2088                        if (val[i] <= previous) {
    2089                            ok = false;
    2090                        }
    2091                    } else {
    2092                        if (val[i] < previous) {
    2093                            ok = false;
    2094                        }
    2095                    }
    2096                    break;
    2097                case DECREASING:
    2098                    if (strict) {
    2099                        if (val[i] >= previous) {
    2100                            ok = false;
    2101                        }
    2102                    } else {
    2103                        if (val[i] > previous) {
    2104                            ok = false;
    2105                        }
    2106                    }
    2107                    break;
    2108                default:
    2109                    // Should never happen.
    2110                    throw new IllegalArgumentException();
    2111                }
    2112
    2113                if (!ok) {
    2114                    throw new NonMonotonousSequenceException(val[i], previous, i, dir, strict);
    2115                }
    2116                previous = val[i];
    2117            }
    2118        }
    2119
    2120        /**
    2121         * Checks that the given array is sorted in strictly increasing order.
    2122         *
    2123         * @param val Values.
    2124         * @throws NonMonotonousSequenceException if the array is not sorted.
    2125         * @since 2.2
    2126         */
    2127        public static void checkOrder(double[] val) {
    2128            checkOrder(val, OrderDirection.INCREASING, true);
    2129        }
    2130
    2131        /**
    2132         * Checks that the given array is sorted.
    2133         *
    2134         * @param val Values
    2135         * @param dir Order direction (-1 for decreasing, 1 for increasing)
    2136         * @param strict Whether the order should be strict
    2137         * @throws NonMonotonousSequenceException if the array is not sorted.
    2138         * @deprecated as of 2.2 (please use the new {@link #checkOrder(double[],OrderDirection,boolean)
    2139         * checkOrder} method). To be removed in 3.0.
    2140         */
    2141        @Deprecated
    2142        public static void checkOrder(double[] val, int dir, boolean strict) {
    2143            if (dir > 0) {
    2144                checkOrder(val, OrderDirection.INCREASING, strict);
    2145            } else {
    2146                checkOrder(val, OrderDirection.DECREASING, strict);
    2147            }
    2148        }
    2149
    2150        /**
    2151         * Returns the Cartesian norm (2-norm), handling both overflow and underflow.
    2152         * Translation of the minpack enorm subroutine.
    2153         *
    2154         * The redistribution policy for MINPACK is available <a
    2155         * href="http://www.netlib.org/minpack/disclaimer">here</a>, for convenience, it
    2156         * is reproduced below.</p>
    2157         *
    2158         * <table border="0" width="80%" cellpadding="10" align="center" bgcolor="#E0E0E0">
    2159         * <tr><td>
    2160         *    Minpack Copyright Notice (1999) University of Chicago.
    2161         *    All rights reserved
    2162         * </td></tr>
    2163         * <tr><td>
    2164         * Redistribution and use in source and binary forms, with or without
    2165         * modification, are permitted provided that the following conditions
    2166         * are met:
    2167         * <ol>
    2168         *  <li>Redistributions of source code must retain the above copyright
    2169         *      notice, this list of conditions and the following disclaimer.</li>
    2170         * <li>Redistributions in binary form must reproduce the above
    2171         *     copyright notice, this list of conditions and the following
    2172         *     disclaimer in the documentation and/or other materials provided
    2173         *     with the distribution.</li>
    2174         * <li>The end-user documentation included with the redistribution, if any,
    2175         *     must include the following acknowledgment:
    2176         *     <code>This product includes software developed by the University of
    2177         *           Chicago, as Operator of Argonne National Laboratory.</code>
    2178         *     Alternately, this acknowledgment may appear in the software itself,
    2179         *     if and wherever such third-party acknowledgments normally appear.</li>
    2180         * <li><strong>WARRANTY DISCLAIMER. THE SOFTWARE IS SUPPLIED "AS IS"
    2181         *     WITHOUT WARRANTY OF ANY KIND. THE COPYRIGHT HOLDER, THE
    2182         *     UNITED STATES, THE UNITED STATES DEPARTMENT OF ENERGY, AND
    2183         *     THEIR EMPLOYEES: (1) DISCLAIM ANY WARRANTIES, EXPRESS OR
    2184         *     IMPLIED, INCLUDING BUT NOT LIMITED TO ANY IMPLIED WARRANTIES
    2185         *     OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, TITLE
    2186         *     OR NON-INFRINGEMENT, (2) DO NOT ASSUME ANY LEGAL LIABILITY
    2187         *     OR RESPONSIBILITY FOR THE ACCURACY, COMPLETENESS, OR
    2188         *     USEFULNESS OF THE SOFTWARE, (3) DO NOT REPRESENT THAT USE OF
    2189         *     THE SOFTWARE WOULD NOT INFRINGE PRIVATELY OWNED RIGHTS, (4)
    2190         *     DO NOT WARRANT THAT THE SOFTWARE WILL FUNCTION
    2191         *     UNINTERRUPTED, THAT IT IS ERROR-FREE OR THAT ANY ERRORS WILL
    2192         *     BE CORRECTED.</strong></li>
    2193         * <li><strong>LIMITATION OF LIABILITY. IN NO EVENT WILL THE COPYRIGHT
    2194         *     HOLDER, THE UNITED STATES, THE UNITED STATES DEPARTMENT OF
    2195         *     ENERGY, OR THEIR EMPLOYEES: BE LIABLE FOR ANY INDIRECT,
    2196         *     INCIDENTAL, CONSEQUENTIAL, SPECIAL OR PUNITIVE DAMAGES OF
    2197         *     ANY KIND OR NATURE, INCLUDING BUT NOT LIMITED TO LOSS OF
    2198         *     PROFITS OR LOSS OF DATA, FOR ANY REASON WHATSOEVER, WHETHER
    2199         *     SUCH LIABILITY IS ASSERTED ON THE BASIS OF CONTRACT, TORT
    2200         *     (INCLUDING NEGLIGENCE OR STRICT LIABILITY), OR OTHERWISE,
    2201         *     EVEN IF ANY OF SAID PARTIES HAS BEEN WARNED OF THE
    2202         *     POSSIBILITY OF SUCH LOSS OR DAMAGES.</strong></li>
    2203         * <ol></td></tr>
    2204         * </table>
    2205         *
    2206         * @param v vector of doubles
    2207         * @return the 2-norm of the vector
    2208         * @since 2.2
    2209         */
    2210        public static double safeNorm(double[] v) {
    2211        double rdwarf = 3.834e-20;
    2212        double rgiant = 1.304e+19;
    2213        double s1=0.0;
    2214        double s2=0.0;
    2215        double s3=0.0;
    2216        double x1max = 0.0;
    2217        double x3max = 0.0;
    2218        double floatn = (double)v.length;
    2219        double agiant = rgiant/floatn;
    2220        for (int i=0;i<v.length;i++) {
    2221            double xabs = Math.abs(v[i]);
    2222            if (xabs<rdwarf || xabs>agiant) {
    2223                if (xabs>rdwarf) {
    2224                    if (xabs>x1max) {
    2225                        double r=x1max/xabs;
    2226                        s1=1.0+s1*r*r;
    2227                        x1max=xabs;
    2228                    } else {
    2229                        double r=xabs/x1max;
    2230                        s1+=r*r;
    2231                    }
    2232                } else {
    2233                    if (xabs>x3max) {
    2234                     double r=x3max/xabs;
    2235                     s3=1.0+s3*r*r;
    2236                     x3max=xabs;
    2237                    } else {
    2238                        if (xabs!=0.0) {
    2239                            double r=xabs/x3max;
    2240                            s3+=r*r;
    2241                        }
    2242                    }
    2243                }
    2244            } else {
    2245             s2+=xabs*xabs;
    2246            }
    2247        }
    2248        double norm;
    2249        if (s1!=0.0) {
    2250            norm = x1max*Math.sqrt(s1+(s2/x1max)/x1max);
    2251        } else {
    2252            if (s2==0.0) {
    2253                norm = x3max*Math.sqrt(s3);
    2254            } else {
    2255                if (s2>=x3max) {
    2256                    norm = Math.sqrt(s2*(1.0+(x3max/s2)*(x3max*s3)));
    2257                } else {
    2258                    norm = Math.sqrt(x3max*((s2/x3max)+(x3max*s3)));
    2259                }
    2260            }
    2261        }
    2262        return norm;
    2263    }
    2264
    2265    }