diff --git a/projects/Math/79/org/apache/commons/math/util/MathUtils.java b/projects/Math/79/org/apache/commons/math/util/MathUtils.java new file mode 100644 index 0000000..46368dc --- /dev/null +++ b/projects/Math/79/org/apache/commons/math/util/MathUtils.java @@ -0,0 +1,1663 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +package org.apache.commons.math.util; + +import java.math.BigDecimal; +import java.math.BigInteger; +import java.util.Arrays; + +import org.apache.commons.math.MathRuntimeException; + +/** + * Some useful additions to the built-in functions in {@link Math}. + * @version $Revision$ $Date$ + */ +public final class MathUtils { + + /** Smallest positive number such that 1 - EPSILON is not numerically equal to 1. */ + public static final double EPSILON = 0x1.0p-53; + + /** Safe minimum, such that 1 / SAFE_MIN does not overflow. + *

In IEEE 754 arithmetic, this is also the smallest normalized + * number 2-1022.

+ */ + public static final double SAFE_MIN = 0x1.0p-1022; + + /** 2 π. */ + public static final double TWO_PI = 2 * Math.PI; + + /** -1.0 cast as a byte. */ + private static final byte NB = (byte)-1; + + /** -1.0 cast as a short. */ + private static final short NS = (short)-1; + + /** 1.0 cast as a byte. */ + private static final byte PB = (byte)1; + + /** 1.0 cast as a short. */ + private static final short PS = (short)1; + + /** 0.0 cast as a byte. */ + private static final byte ZB = (byte)0; + + /** 0.0 cast as a short. */ + private static final short ZS = (short)0; + + /** Gap between NaN and regular numbers. */ + private static final int NAN_GAP = 4 * 1024 * 1024; + + /** Offset to order signed double numbers lexicographically. */ + private static final long SGN_MASK = 0x8000000000000000L; + + /** All long-representable factorials */ + private static final long[] FACTORIALS = new long[] { + 1l, 1l, 2l, + 6l, 24l, 120l, + 720l, 5040l, 40320l, + 362880l, 3628800l, 39916800l, + 479001600l, 6227020800l, 87178291200l, + 1307674368000l, 20922789888000l, 355687428096000l, + 6402373705728000l, 121645100408832000l, 2432902008176640000l }; + + /** + * Private Constructor + */ + private MathUtils() { + super(); + } + + /** + * Add two integers, checking for overflow. + * + * @param x an addend + * @param y an addend + * @return the sum x+y + * @throws ArithmeticException if the result can not be represented as an + * int + * @since 1.1 + */ + public static int addAndCheck(int x, int y) { + long s = (long)x + (long)y; + if (s < Integer.MIN_VALUE || s > Integer.MAX_VALUE) { + throw new ArithmeticException("overflow: add"); + } + return (int)s; + } + + /** + * Add two long integers, checking for overflow. + * + * @param a an addend + * @param b an addend + * @return the sum a+b + * @throws ArithmeticException if the result can not be represented as an + * long + * @since 1.2 + */ + public static long addAndCheck(long a, long b) { + return addAndCheck(a, b, "overflow: add"); + } + + /** + * Add two long integers, checking for overflow. + * + * @param a an addend + * @param b an addend + * @param msg the message to use for any thrown exception. + * @return the sum a+b + * @throws ArithmeticException if the result can not be represented as an + * long + * @since 1.2 + */ + private static long addAndCheck(long a, long b, String msg) { + long ret; + if (a > b) { + // use symmetry to reduce boundary cases + ret = addAndCheck(b, a, msg); + } else { + // assert a <= b + + if (a < 0) { + if (b < 0) { + // check for negative overflow + if (Long.MIN_VALUE - b <= a) { + ret = a + b; + } else { + throw new ArithmeticException(msg); + } + } else { + // opposite sign addition is always safe + ret = a + b; + } + } else { + // assert a >= 0 + // assert b >= 0 + + // check for positive overflow + if (a <= Long.MAX_VALUE - b) { + ret = a + b; + } else { + throw new ArithmeticException(msg); + } + } + } + return ret; + } + + /** + * Returns an exact representation of the Binomial + * Coefficient, "n choose k", the number of + * k-element subsets that can be selected from an + * n-element set. + *

+ * Preconditions: + *

+ * + * @param n the size of the set + * @param k the size of the subsets to be counted + * @return n choose k + * @throws IllegalArgumentException if preconditions are not met. + * @throws ArithmeticException if the result is too large to be represented + * by a long integer. + */ + public static long binomialCoefficient(final int n, final int k) { + checkBinomial(n, k); + if ((n == k) || (k == 0)) { + return 1; + } + if ((k == 1) || (k == n - 1)) { + return n; + } + // Use symmetry for large k + if (k > n / 2) + return binomialCoefficient(n, n - k); + + // We use the formula + // (n choose k) = n! / (n-k)! / k! + // (n choose k) == ((n-k+1)*...*n) / (1*...*k) + // which could be written + // (n choose k) == (n-1 choose k-1) * n / k + long result = 1; + if (n <= 61) { + // For n <= 61, the naive implementation cannot overflow. + int i = n - k + 1; + for (int j = 1; j <= k; j++) { + result = result * i / j; + i++; + } + } else if (n <= 66) { + // For n > 61 but n <= 66, the result cannot overflow, + // but we must take care not to overflow intermediate values. + int i = n - k + 1; + for (int j = 1; j <= k; j++) { + // We know that (result * i) is divisible by j, + // but (result * i) may overflow, so we split j: + // Filter out the gcd, d, so j/d and i/d are integer. + // result is divisible by (j/d) because (j/d) + // is relative prime to (i/d) and is a divisor of + // result * (i/d). + final long d = gcd(i, j); + result = (result / (j / d)) * (i / d); + i++; + } + } else { + // For n > 66, a result overflow might occur, so we check + // the multiplication, taking care to not overflow + // unnecessary. + int i = n - k + 1; + for (int j = 1; j <= k; j++) { + final long d = gcd(i, j); + result = mulAndCheck(result / (j / d), i / d); + i++; + } + } + return result; + } + + /** + * Returns a double representation of the Binomial + * Coefficient, "n choose k", the number of + * k-element subsets that can be selected from an + * n-element set. + *

+ * Preconditions: + *

+ * + * @param n the size of the set + * @param k the size of the subsets to be counted + * @return n choose k + * @throws IllegalArgumentException if preconditions are not met. + */ + public static double binomialCoefficientDouble(final int n, final int k) { + checkBinomial(n, k); + if ((n == k) || (k == 0)) { + return 1d; + } + if ((k == 1) || (k == n - 1)) { + return n; + } + if (k > n/2) { + return binomialCoefficientDouble(n, n - k); + } + if (n < 67) { + return binomialCoefficient(n,k); + } + + double result = 1d; + for (int i = 1; i <= k; i++) { + result *= (double)(n - k + i) / (double)i; + } + + return Math.floor(result + 0.5); + } + + /** + * Returns the natural log of the Binomial + * Coefficient, "n choose k", the number of + * k-element subsets that can be selected from an + * n-element set. + *

+ * Preconditions: + *

+ * + * @param n the size of the set + * @param k the size of the subsets to be counted + * @return n choose k + * @throws IllegalArgumentException if preconditions are not met. + */ + public static double binomialCoefficientLog(final int n, final int k) { + checkBinomial(n, k); + if ((n == k) || (k == 0)) { + return 0; + } + if ((k == 1) || (k == n - 1)) { + return Math.log(n); + } + + /* + * For values small enough to do exact integer computation, + * return the log of the exact value + */ + if (n < 67) { + return Math.log(binomialCoefficient(n,k)); + } + + /* + * Return the log of binomialCoefficientDouble for values that will not + * overflow binomialCoefficientDouble + */ + if (n < 1030) { + return Math.log(binomialCoefficientDouble(n, k)); + } + + if (k > n / 2) { + return binomialCoefficientLog(n, n - k); + } + + /* + * Sum logs for values that could overflow + */ + double logSum = 0; + + // n!/(n-k)! + for (int i = n - k + 1; i <= n; i++) { + logSum += Math.log(i); + } + + // divide by k! + for (int i = 2; i <= k; i++) { + logSum -= Math.log(i); + } + + return logSum; + } + + /** + * Check binomial preconditions. + * @param n the size of the set + * @param k the size of the subsets to be counted + * @exception IllegalArgumentException if preconditions are not met. + */ + private static void checkBinomial(final int n, final int k) + throws IllegalArgumentException { + if (n < k) { + throw MathRuntimeException.createIllegalArgumentException( + "must have n >= k for binomial coefficient (n,k), got n = {0}, k = {1}", + n, k); + } + if (n < 0) { + throw MathRuntimeException.createIllegalArgumentException( + "must have n >= 0 for binomial coefficient (n,k), got n = {0}", + n); + } + } + + /** + * Compares two numbers given some amount of allowed error. + * + * @param x the first number + * @param y the second number + * @param eps the amount of error to allow when checking for equality + * @return + */ + public static int compareTo(double x, double y, double eps) { + if (equals(x, y, eps)) { + return 0; + } else if (x < y) { + return -1; + } + return 1; + } + + /** + * Returns the + * hyperbolic cosine of x. + * + * @param x double value for which to find the hyperbolic cosine + * @return hyperbolic cosine of x + */ + public static double cosh(double x) { + return (Math.exp(x) + Math.exp(-x)) / 2.0; + } + + /** + * Returns true iff both arguments are NaN or neither is NaN and they are + * equal + * + * @param x first value + * @param y second value + * @return true if the values are equal or both are NaN + */ + public static boolean equals(double x, double y) { + return (Double.isNaN(x) && Double.isNaN(y)) || x == y; + } + + /** + * Returns true iff both arguments are equal or within the range of allowed + * error (inclusive). + *

+ * Two NaNs are considered equals, as are two infinities with same sign. + *

+ * + * @param x first value + * @param y second value + * @param eps the amount of absolute error to allow + * @return true if the values are equal or within range of each other + */ + public static boolean equals(double x, double y, double eps) { + return equals(x, y) || (Math.abs(y - x) <= eps); + } + + /** + * Returns true iff both arguments are equal or within the range of allowed + * error (inclusive). + * Adapted from + * Bruce Dawson + * + * @param x first value + * @param y second value + * @param maxUlps {@code (maxUlps - 1)} is the number of floating point + * values between {@code x} and {@code y}. + * @return {@code true} if there are less than {@code maxUlps} floating + * point values between {@code x} and {@code y} + */ + public static boolean equals(double x, double y, int maxUlps) { + // Check that "maxUlps" is non-negative and small enough so that the + // default NAN won't compare as equal to anything. + assert maxUlps > 0 && maxUlps < NAN_GAP; + + long xInt = Double.doubleToLongBits(x); + long yInt = Double.doubleToLongBits(y); + + // Make lexicographically ordered as a two's-complement integer. + if (xInt < 0) { + xInt = SGN_MASK - xInt; + } + if (yInt < 0) { + yInt = SGN_MASK - yInt; + } + + return Math.abs(xInt - yInt) <= maxUlps; + } + + /** + * Returns true iff both arguments are null or have same dimensions + * and all their elements are {@link #equals(double,double) equals} + * + * @param x first array + * @param y second array + * @return true if the values are both null or have same dimension + * and equal elements + * @since 1.2 + */ + public static boolean equals(double[] x, double[] y) { + if ((x == null) || (y == null)) { + return !((x == null) ^ (y == null)); + } + if (x.length != y.length) { + return false; + } + for (int i = 0; i < x.length; ++i) { + if (!equals(x[i], y[i])) { + return false; + } + } + return true; + } + + /** + * Returns n!. Shorthand for n Factorial, the + * product of the numbers 1,...,n. + *

+ * Preconditions: + *

+ *

+ * + * @param n argument + * @return n! + * @throws ArithmeticException if the result is too large to be represented + * by a long integer. + * @throws IllegalArgumentException if n < 0 + */ + public static long factorial(final int n) { + if (n < 0) { + throw MathRuntimeException.createIllegalArgumentException( + "must have n >= 0 for n!, got n = {0}", + n); + } + if (n > 20) { + throw new ArithmeticException( + "factorial value is too large to fit in a long"); + } + return FACTORIALS[n]; + } + + /** + * Returns n!. Shorthand for n Factorial, the + * product of the numbers 1,...,n as a double. + *

+ * Preconditions: + *

+ *

+ * + * @param n argument + * @return n! + * @throws IllegalArgumentException if n < 0 + */ + public static double factorialDouble(final int n) { + if (n < 0) { + throw MathRuntimeException.createIllegalArgumentException( + "must have n >= 0 for n!, got n = {0}", + n); + } + if (n < 21) { + return factorial(n); + } + return Math.floor(Math.exp(factorialLog(n)) + 0.5); + } + + /** + * Returns the natural logarithm of n!. + *

+ * Preconditions: + *

+ * + * @param n argument + * @return n! + * @throws IllegalArgumentException if preconditions are not met. + */ + public static double factorialLog(final int n) { + if (n < 0) { + throw MathRuntimeException.createIllegalArgumentException( + "must have n >= 0 for n!, got n = {0}", + n); + } + if (n < 21) { + return Math.log(factorial(n)); + } + double logSum = 0; + for (int i = 2; i <= n; i++) { + logSum += Math.log(i); + } + return logSum; + } + + /** + *

+ * Gets the greatest common divisor of the absolute value of two numbers, + * using the "binary gcd" method which avoids division and modulo + * operations. See Knuth 4.5.2 algorithm B. This algorithm is due to Josef + * Stein (1961). + *

+ * Special cases: + * + * + * @param p any number + * @param q any number + * @return the greatest common divisor, never negative + * @throws ArithmeticException + * if the result cannot be represented as a nonnegative int + * value + * @since 1.1 + */ + public static int gcd(final int p, final int q) { + int u = p; + int v = q; + if ((u == 0) || (v == 0)) { + if ((u == Integer.MIN_VALUE) || (v == Integer.MIN_VALUE)) { + throw MathRuntimeException.createArithmeticException( + "overflow: gcd({0}, {1}) is 2^31", + p, q); + } + return Math.abs(u) + Math.abs(v); + } + // keep u and v negative, as negative integers range down to + // -2^31, while positive numbers can only be as large as 2^31-1 + // (i.e. we can't necessarily negate a negative number without + // overflow) + /* assert u!=0 && v!=0; */ + if (u > 0) { + u = -u; + } // make u negative + if (v > 0) { + v = -v; + } // make v negative + // B1. [Find power of 2] + int k = 0; + while ((u & 1) == 0 && (v & 1) == 0 && k < 31) { // while u and v are + // both even... + u /= 2; + v /= 2; + k++; // cast out twos. + } + if (k == 31) { + throw MathRuntimeException.createArithmeticException( + "overflow: gcd({0}, {1}) is 2^31", + p, q); + } + // B2. Initialize: u and v have been divided by 2^k and at least + // one is odd. + int t = ((u & 1) == 1) ? v : -(u / 2)/* B3 */; + // t negative: u was odd, v may be even (t replaces v) + // t positive: u was even, v is odd (t replaces u) + do { + /* assert u<0 && v<0; */ + // B4/B3: cast out twos from t. + while ((t & 1) == 0) { // while t is even.. + t /= 2; // cast out twos + } + // B5 [reset max(u,v)] + if (t > 0) { + u = -t; + } else { + v = t; + } + // B6/B3. at this point both u and v should be odd. + t = (v - u) / 2; + // |u| larger: t positive (replace u) + // |v| larger: t negative (replace v) + } while (t != 0); + return -u * (1 << k); // gcd is u*2^k + } + + /** + * Returns an integer hash code representing the given double value. + * + * @param value the value to be hashed + * @return the hash code + */ + public static int hash(double value) { + return new Double(value).hashCode(); + } + + /** + * Returns an integer hash code representing the given double array. + * + * @param value the value to be hashed (may be null) + * @return the hash code + * @since 1.2 + */ + public static int hash(double[] value) { + return Arrays.hashCode(value); + } + + /** + * For a byte value x, this method returns (byte)(+1) if x >= 0 and + * (byte)(-1) if x < 0. + * + * @param x the value, a byte + * @return (byte)(+1) or (byte)(-1), depending on the sign of x + */ + public static byte indicator(final byte x) { + return (x >= ZB) ? PB : NB; + } + + /** + * For a double precision value x, this method returns +1.0 if x >= 0 and + * -1.0 if x < 0. Returns NaN if x is + * NaN. + * + * @param x the value, a double + * @return +1.0 or -1.0, depending on the sign of x + */ + public static double indicator(final double x) { + if (Double.isNaN(x)) { + return Double.NaN; + } + return (x >= 0.0) ? 1.0 : -1.0; + } + + /** + * For a float value x, this method returns +1.0F if x >= 0 and -1.0F if x < + * 0. Returns NaN if x is NaN. + * + * @param x the value, a float + * @return +1.0F or -1.0F, depending on the sign of x + */ + public static float indicator(final float x) { + if (Float.isNaN(x)) { + return Float.NaN; + } + return (x >= 0.0F) ? 1.0F : -1.0F; + } + + /** + * For an int value x, this method returns +1 if x >= 0 and -1 if x < 0. + * + * @param x the value, an int + * @return +1 or -1, depending on the sign of x + */ + public static int indicator(final int x) { + return (x >= 0) ? 1 : -1; + } + + /** + * For a long value x, this method returns +1L if x >= 0 and -1L if x < 0. + * + * @param x the value, a long + * @return +1L or -1L, depending on the sign of x + */ + public static long indicator(final long x) { + return (x >= 0L) ? 1L : -1L; + } + + /** + * For a short value x, this method returns (short)(+1) if x >= 0 and + * (short)(-1) if x < 0. + * + * @param x the value, a short + * @return (short)(+1) or (short)(-1), depending on the sign of x + */ + public static short indicator(final short x) { + return (x >= ZS) ? PS : NS; + } + + /** + *

+ * Returns the least common multiple of the absolute value of two numbers, + * using the formula lcm(a,b) = (a / gcd(a,b)) * b. + *

+ * Special cases: + * + * + * @param a any number + * @param b any number + * @return the least common multiple, never negative + * @throws ArithmeticException + * if the result cannot be represented as a nonnegative int + * value + * @since 1.1 + */ + public static int lcm(int a, int b) { + if (a==0 || b==0){ + return 0; + } + int lcm = Math.abs(mulAndCheck(a / gcd(a, b), b)); + if (lcm == Integer.MIN_VALUE){ + throw new ArithmeticException("overflow: lcm is 2^31"); + } + return lcm; + } + + /** + *

Returns the + * logarithm + * for base b of x. + *

+ *

Returns NaN if either argument is negative. If + * base is 0 and x is positive, 0 is returned. + * If base is positive and x is 0, + * Double.NEGATIVE_INFINITY is returned. If both arguments + * are 0, the result is NaN.

+ * + * @param base the base of the logarithm, must be greater than 0 + * @param x argument, must be greater than 0 + * @return the value of the logarithm - the number y such that base^y = x. + * @since 1.2 + */ + public static double log(double base, double x) { + return Math.log(x)/Math.log(base); + } + + /** + * Multiply two integers, checking for overflow. + * + * @param x a factor + * @param y a factor + * @return the product x*y + * @throws ArithmeticException if the result can not be represented as an + * int + * @since 1.1 + */ + public static int mulAndCheck(int x, int y) { + long m = ((long)x) * ((long)y); + if (m < Integer.MIN_VALUE || m > Integer.MAX_VALUE) { + throw new ArithmeticException("overflow: mul"); + } + return (int)m; + } + + /** + * Multiply two long integers, checking for overflow. + * + * @param a first value + * @param b second value + * @return the product a * b + * @throws ArithmeticException if the result can not be represented as an + * long + * @since 1.2 + */ + public static long mulAndCheck(long a, long b) { + long ret; + String msg = "overflow: multiply"; + if (a > b) { + // use symmetry to reduce boundary cases + ret = mulAndCheck(b, a); + } else { + if (a < 0) { + if (b < 0) { + // check for positive overflow with negative a, negative b + if (a >= Long.MAX_VALUE / b) { + ret = a * b; + } else { + throw new ArithmeticException(msg); + } + } else if (b > 0) { + // check for negative overflow with negative a, positive b + if (Long.MIN_VALUE / b <= a) { + ret = a * b; + } else { + throw new ArithmeticException(msg); + + } + } else { + // assert b == 0 + ret = 0; + } + } else if (a > 0) { + // assert a > 0 + // assert b > 0 + + // check for positive overflow with positive a, positive b + if (a <= Long.MAX_VALUE / b) { + ret = a * b; + } else { + throw new ArithmeticException(msg); + } + } else { + // assert a == 0 + ret = 0; + } + } + return ret; + } + + /** + * Get the next machine representable number after a number, moving + * in the direction of another number. + *

+ * If direction is greater than or equal tod, + * the smallest machine representable number strictly greater than + * d is returned; otherwise the largest representable number + * strictly less than d is returned.

+ *

+ * If d is NaN or Infinite, it is returned unchanged.

+ * + * @param d base number + * @param direction (the only important thing is whether + * direction is greater or smaller than d) + * @return the next machine representable number in the specified direction + * @since 1.2 + */ + public static double nextAfter(double d, double direction) { + + // handling of some important special cases + if (Double.isNaN(d) || Double.isInfinite(d)) { + return d; + } else if (d == 0) { + return (direction < 0) ? -Double.MIN_VALUE : Double.MIN_VALUE; + } + // special cases MAX_VALUE to infinity and MIN_VALUE to 0 + // are handled just as normal numbers + + // split the double in raw components + long bits = Double.doubleToLongBits(d); + long sign = bits & 0x8000000000000000L; + long exponent = bits & 0x7ff0000000000000L; + long mantissa = bits & 0x000fffffffffffffL; + + if (d * (direction - d) >= 0) { + // we should increase the mantissa + if (mantissa == 0x000fffffffffffffL) { + return Double.longBitsToDouble(sign | + (exponent + 0x0010000000000000L)); + } else { + return Double.longBitsToDouble(sign | + exponent | (mantissa + 1)); + } + } else { + // we should decrease the mantissa + if (mantissa == 0L) { + return Double.longBitsToDouble(sign | + (exponent - 0x0010000000000000L) | + 0x000fffffffffffffL); + } else { + return Double.longBitsToDouble(sign | + exponent | (mantissa - 1)); + } + } + + } + + /** + * Scale a number by 2scaleFactor. + *

If d is 0 or NaN or Infinite, it is returned unchanged.

+ * + * @param d base number + * @param scaleFactor power of two by which d sould be multiplied + * @return d × 2scaleFactor + * @since 2.0 + */ + public static double scalb(final double d, final int scaleFactor) { + + // handling of some important special cases + if ((d == 0) || Double.isNaN(d) || Double.isInfinite(d)) { + return d; + } + + // split the double in raw components + final long bits = Double.doubleToLongBits(d); + final long exponent = bits & 0x7ff0000000000000L; + final long rest = bits & 0x800fffffffffffffL; + + // shift the exponent + final long newBits = rest | (exponent + (((long) scaleFactor) << 52)); + return Double.longBitsToDouble(newBits); + + } + + /** + * Normalize an angle in a 2&pi wide interval around a center value. + *

This method has three main uses:

+ * + *

Note that due to numerical accuracy and since π cannot be represented + * exactly, the result interval is closed, it cannot be half-closed + * as would be more satisfactory in a purely mathematical view.

+ * @param a angle to normalize + * @param center center of the desired 2π interval for the result + * @return a-2kπ with integer k and center-π <= a-2kπ <= center+π + * @since 1.2 + */ + public static double normalizeAngle(double a, double center) { + return a - TWO_PI * Math.floor((a + Math.PI - center) / TWO_PI); + } + + /** + *

Normalizes an array to make it sum to a specified value. + * Returns the result of the transformation 

+      *    x |-> x * normalizedSum / sum
+      *
+ * applied to each non-NaN element x of the input array, where sum is the + * sum of the non-NaN entries in the input array.

+ * + *

Throws IllegalArgumentException if normalizedSum is infinite + * or NaN and ArithmeticException if the input array contains any infinite elements + * or sums to 0

+ * + *

Ignores (i.e., copies unchanged to the output array) NaNs in the input array.

+ * + * @param values input array to be normalized + * @param normalizedSum target sum for the normalized array + * @return normalized array + * @throws ArithmeticException if the input array contains infinite elements or sums to zero + * @throws IllegalArgumentException if the target sum is infinite or NaN + */ + public static double[] normalizeArray(double[] values, double normalizedSum) + throws ArithmeticException, IllegalArgumentException { + if (Double.isInfinite(normalizedSum)) { + throw MathRuntimeException.createIllegalArgumentException( + "Cannot normalize to an infinite value"); + } + if (Double.isNaN(normalizedSum)) { + throw MathRuntimeException.createIllegalArgumentException( + "Cannot normalize to NaN"); + } + double sum = 0d; + final int len = values.length; + double[] out = new double[len]; + for (int i = 0; i < len; i++) { + if (Double.isInfinite(values[i])) { + throw MathRuntimeException.createArithmeticException( + "Array contains an infinite element, {0} at index {1}", values[i], i); + } + if (!Double.isNaN(values[i])) { + sum += values[i]; + } + } + if (sum == 0) { + throw MathRuntimeException.createArithmeticException( + "Array sums to zero"); + } + for (int i = 0; i < len; i++) { + if (Double.isNaN(values[i])) { + out[i] = Double.NaN; + } else { + out[i] = values[i] * normalizedSum / sum; + } + } + return out; + } + + /** + * Round the given value to the specified number of decimal places. The + * value is rounded using the {@link BigDecimal#ROUND_HALF_UP} method. + * + * @param x the value to round. + * @param scale the number of digits to the right of the decimal point. + * @return the rounded value. + * @since 1.1 + */ + public static double round(double x, int scale) { + return round(x, scale, BigDecimal.ROUND_HALF_UP); + } + + /** + * Round the given value to the specified number of decimal places. The + * value is rounded using the given method which is any method defined in + * {@link BigDecimal}. + * + * @param x the value to round. + * @param scale the number of digits to the right of the decimal point. + * @param roundingMethod the rounding method as defined in + * {@link BigDecimal}. + * @return the rounded value. + * @since 1.1 + */ + public static double round(double x, int scale, int roundingMethod) { + try { + return (new BigDecimal + (Double.toString(x)) + .setScale(scale, roundingMethod)) + .doubleValue(); + } catch (NumberFormatException ex) { + if (Double.isInfinite(x)) { + return x; + } else { + return Double.NaN; + } + } + } + + /** + * Round the given value to the specified number of decimal places. The + * value is rounding using the {@link BigDecimal#ROUND_HALF_UP} method. + * + * @param x the value to round. + * @param scale the number of digits to the right of the decimal point. + * @return the rounded value. + * @since 1.1 + */ + public static float round(float x, int scale) { + return round(x, scale, BigDecimal.ROUND_HALF_UP); + } + + /** + * Round the given value to the specified number of decimal places. The + * value is rounded using the given method which is any method defined in + * {@link BigDecimal}. + * + * @param x the value to round. + * @param scale the number of digits to the right of the decimal point. + * @param roundingMethod the rounding method as defined in + * {@link BigDecimal}. + * @return the rounded value. + * @since 1.1 + */ + public static float round(float x, int scale, int roundingMethod) { + float sign = indicator(x); + float factor = (float)Math.pow(10.0f, scale) * sign; + return (float)roundUnscaled(x * factor, sign, roundingMethod) / factor; + } + + /** + * Round the given non-negative, value to the "nearest" integer. Nearest is + * determined by the rounding method specified. Rounding methods are defined + * in {@link BigDecimal}. + * + * @param unscaled the value to round. + * @param sign the sign of the original, scaled value. + * @param roundingMethod the rounding method as defined in + * {@link BigDecimal}. + * @return the rounded value. + * @since 1.1 + */ + private static double roundUnscaled(double unscaled, double sign, + int roundingMethod) { + switch (roundingMethod) { + case BigDecimal.ROUND_CEILING : + if (sign == -1) { + unscaled = Math.floor(nextAfter(unscaled, Double.NEGATIVE_INFINITY)); + } else { + unscaled = Math.ceil(nextAfter(unscaled, Double.POSITIVE_INFINITY)); + } + break; + case BigDecimal.ROUND_DOWN : + unscaled = Math.floor(nextAfter(unscaled, Double.NEGATIVE_INFINITY)); + break; + case BigDecimal.ROUND_FLOOR : + if (sign == -1) { + unscaled = Math.ceil(nextAfter(unscaled, Double.POSITIVE_INFINITY)); + } else { + unscaled = Math.floor(nextAfter(unscaled, Double.NEGATIVE_INFINITY)); + } + break; + case BigDecimal.ROUND_HALF_DOWN : { + unscaled = nextAfter(unscaled, Double.NEGATIVE_INFINITY); + double fraction = unscaled - Math.floor(unscaled); + if (fraction > 0.5) { + unscaled = Math.ceil(unscaled); + } else { + unscaled = Math.floor(unscaled); + } + break; + } + case BigDecimal.ROUND_HALF_EVEN : { + double fraction = unscaled - Math.floor(unscaled); + if (fraction > 0.5) { + unscaled = Math.ceil(unscaled); + } else if (fraction < 0.5) { + unscaled = Math.floor(unscaled); + } else { + // The following equality test is intentional and needed for rounding purposes + if (Math.floor(unscaled) / 2.0 == Math.floor(Math + .floor(unscaled) / 2.0)) { // even + unscaled = Math.floor(unscaled); + } else { // odd + unscaled = Math.ceil(unscaled); + } + } + break; + } + case BigDecimal.ROUND_HALF_UP : { + unscaled = nextAfter(unscaled, Double.POSITIVE_INFINITY); + double fraction = unscaled - Math.floor(unscaled); + if (fraction >= 0.5) { + unscaled = Math.ceil(unscaled); + } else { + unscaled = Math.floor(unscaled); + } + break; + } + case BigDecimal.ROUND_UNNECESSARY : + if (unscaled != Math.floor(unscaled)) { + throw new ArithmeticException("Inexact result from rounding"); + } + break; + case BigDecimal.ROUND_UP : + unscaled = Math.ceil(nextAfter(unscaled, Double.POSITIVE_INFINITY)); + break; + default : + throw MathRuntimeException.createIllegalArgumentException( + "invalid rounding method {0}, valid methods: {1} ({2}), {3} ({4})," + + " {5} ({6}), {7} ({8}), {9} ({10}), {11} ({12}), {13} ({14}), {15} ({16})", + roundingMethod, + "ROUND_CEILING", BigDecimal.ROUND_CEILING, + "ROUND_DOWN", BigDecimal.ROUND_DOWN, + "ROUND_FLOOR", BigDecimal.ROUND_FLOOR, + "ROUND_HALF_DOWN", BigDecimal.ROUND_HALF_DOWN, + "ROUND_HALF_EVEN", BigDecimal.ROUND_HALF_EVEN, + "ROUND_HALF_UP", BigDecimal.ROUND_HALF_UP, + "ROUND_UNNECESSARY", BigDecimal.ROUND_UNNECESSARY, + "ROUND_UP", BigDecimal.ROUND_UP); + } + return unscaled; + } + + /** + * Returns the sign + * for byte value x. + *

+ * For a byte value x, this method returns (byte)(+1) if x > 0, (byte)(0) if + * x = 0, and (byte)(-1) if x < 0.

+ * + * @param x the value, a byte + * @return (byte)(+1), (byte)(0), or (byte)(-1), depending on the sign of x + */ + public static byte sign(final byte x) { + return (x == ZB) ? ZB : (x > ZB) ? PB : NB; + } + + /** + * Returns the sign + * for double precision x. + *

+ * For a double value x, this method returns + * +1.0 if x > 0, 0.0 if + * x = 0.0, and -1.0 if x < 0. + * Returns NaN if x is NaN.

+ * + * @param x the value, a double + * @return +1.0, 0.0, or -1.0, depending on the sign of x + */ + public static double sign(final double x) { + if (Double.isNaN(x)) { + return Double.NaN; + } + return (x == 0.0) ? 0.0 : (x > 0.0) ? 1.0 : -1.0; + } + + /** + * Returns the sign + * for float value x. + *

+ * For a float value x, this method returns +1.0F if x > 0, 0.0F if x = + * 0.0F, and -1.0F if x < 0. Returns NaN if x + * is NaN.

+ * + * @param x the value, a float + * @return +1.0F, 0.0F, or -1.0F, depending on the sign of x + */ + public static float sign(final float x) { + if (Float.isNaN(x)) { + return Float.NaN; + } + return (x == 0.0F) ? 0.0F : (x > 0.0F) ? 1.0F : -1.0F; + } + + /** + * Returns the sign + * for int value x. + *

+ * For an int value x, this method returns +1 if x > 0, 0 if x = 0, and -1 + * if x < 0.

+ * + * @param x the value, an int + * @return +1, 0, or -1, depending on the sign of x + */ + public static int sign(final int x) { + return (x == 0) ? 0 : (x > 0) ? 1 : -1; + } + + /** + * Returns the sign + * for long value x. + *

+ * For a long value x, this method returns +1L if x > 0, 0L if x = 0, and + * -1L if x < 0.

+ * + * @param x the value, a long + * @return +1L, 0L, or -1L, depending on the sign of x + */ + public static long sign(final long x) { + return (x == 0L) ? 0L : (x > 0L) ? 1L : -1L; + } + + /** + * Returns the sign + * for short value x. + *

+ * For a short value x, this method returns (short)(+1) if x > 0, (short)(0) + * if x = 0, and (short)(-1) if x < 0.

+ * + * @param x the value, a short + * @return (short)(+1), (short)(0), or (short)(-1), depending on the sign of + * x + */ + public static short sign(final short x) { + return (x == ZS) ? ZS : (x > ZS) ? PS : NS; + } + + /** + * Returns the + * hyperbolic sine of x. + * + * @param x double value for which to find the hyperbolic sine + * @return hyperbolic sine of x + */ + public static double sinh(double x) { + return (Math.exp(x) - Math.exp(-x)) / 2.0; + } + + /** + * Subtract two integers, checking for overflow. + * + * @param x the minuend + * @param y the subtrahend + * @return the difference x-y + * @throws ArithmeticException if the result can not be represented as an + * int + * @since 1.1 + */ + public static int subAndCheck(int x, int y) { + long s = (long)x - (long)y; + if (s < Integer.MIN_VALUE || s > Integer.MAX_VALUE) { + throw new ArithmeticException("overflow: subtract"); + } + return (int)s; + } + + /** + * Subtract two long integers, checking for overflow. + * + * @param a first value + * @param b second value + * @return the difference a-b + * @throws ArithmeticException if the result can not be represented as an + * long + * @since 1.2 + */ + public static long subAndCheck(long a, long b) { + long ret; + String msg = "overflow: subtract"; + if (b == Long.MIN_VALUE) { + if (a < 0) { + ret = a - b; + } else { + throw new ArithmeticException(msg); + } + } else { + // use additive inverse + ret = addAndCheck(a, -b, msg); + } + return ret; + } + + /** + * Raise an int to an int power. + * @param k number to raise + * @param e exponent (must be positive or null) + * @return ke + * @exception IllegalArgumentException if e is negative + */ + public static int pow(final int k, int e) + throws IllegalArgumentException { + + if (e < 0) { + throw MathRuntimeException.createIllegalArgumentException( + "cannot raise an integral value to a negative power ({0}^{1})", + k, e); + } + + int result = 1; + int k2p = k; + while (e != 0) { + if ((e & 0x1) != 0) { + result *= k2p; + } + k2p *= k2p; + e = e >> 1; + } + + return result; + + } + + /** + * Raise an int to a long power. + * @param k number to raise + * @param e exponent (must be positive or null) + * @return ke + * @exception IllegalArgumentException if e is negative + */ + public static int pow(final int k, long e) + throws IllegalArgumentException { + + if (e < 0) { + throw MathRuntimeException.createIllegalArgumentException( + "cannot raise an integral value to a negative power ({0}^{1})", + k, e); + } + + int result = 1; + int k2p = k; + while (e != 0) { + if ((e & 0x1) != 0) { + result *= k2p; + } + k2p *= k2p; + e = e >> 1; + } + + return result; + + } + + /** + * Raise a long to an int power. + * @param k number to raise + * @param e exponent (must be positive or null) + * @return ke + * @exception IllegalArgumentException if e is negative + */ + public static long pow(final long k, int e) + throws IllegalArgumentException { + + if (e < 0) { + throw MathRuntimeException.createIllegalArgumentException( + "cannot raise an integral value to a negative power ({0}^{1})", + k, e); + } + + long result = 1l; + long k2p = k; + while (e != 0) { + if ((e & 0x1) != 0) { + result *= k2p; + } + k2p *= k2p; + e = e >> 1; + } + + return result; + + } + + /** + * Raise a long to a long power. + * @param k number to raise + * @param e exponent (must be positive or null) + * @return ke + * @exception IllegalArgumentException if e is negative + */ + public static long pow(final long k, long e) + throws IllegalArgumentException { + + if (e < 0) { + throw MathRuntimeException.createIllegalArgumentException( + "cannot raise an integral value to a negative power ({0}^{1})", + k, e); + } + + long result = 1l; + long k2p = k; + while (e != 0) { + if ((e & 0x1) != 0) { + result *= k2p; + } + k2p *= k2p; + e = e >> 1; + } + + return result; + + } + + /** + * Raise a BigInteger to an int power. + * @param k number to raise + * @param e exponent (must be positive or null) + * @return ke + * @exception IllegalArgumentException if e is negative + */ + public static BigInteger pow(final BigInteger k, int e) + throws IllegalArgumentException { + + if (e < 0) { + throw MathRuntimeException.createIllegalArgumentException( + "cannot raise an integral value to a negative power ({0}^{1})", + k, e); + } + + return k.pow(e); + + } + + /** + * Raise a BigInteger to a long power. + * @param k number to raise + * @param e exponent (must be positive or null) + * @return ke + * @exception IllegalArgumentException if e is negative + */ + public static BigInteger pow(final BigInteger k, long e) + throws IllegalArgumentException { + + if (e < 0) { + throw MathRuntimeException.createIllegalArgumentException( + "cannot raise an integral value to a negative power ({0}^{1})", + k, e); + } + + BigInteger result = BigInteger.ONE; + BigInteger k2p = k; + while (e != 0) { + if ((e & 0x1) != 0) { + result = result.multiply(k2p); + } + k2p = k2p.multiply(k2p); + e = e >> 1; + } + + return result; + + } + + /** + * Raise a BigInteger to a BigInteger power. + * @param k number to raise + * @param e exponent (must be positive or null) + * @return ke + * @exception IllegalArgumentException if e is negative + */ + public static BigInteger pow(final BigInteger k, BigInteger e) + throws IllegalArgumentException { + + if (e.compareTo(BigInteger.ZERO) < 0) { + throw MathRuntimeException.createIllegalArgumentException( + "cannot raise an integral value to a negative power ({0}^{1})", + k, e); + } + + BigInteger result = BigInteger.ONE; + BigInteger k2p = k; + while (!BigInteger.ZERO.equals(e)) { + if (e.testBit(0)) { + result = result.multiply(k2p); + } + k2p = k2p.multiply(k2p); + e = e.shiftRight(1); + } + + return result; + + } + + /** + * Calculates the L1 (sum of abs) distance between two points. + * + * @param p1 the first point + * @param p2 the second point + * @return the L1 distance between the two points + */ + public static double distance1(double[] p1, double[] p2) { + double sum = 0; + for (int i = 0; i < p1.length; i++) { + sum += Math.abs(p1[i] - p2[i]); + } + return sum; + } + + /** + * Calculates the L1 (sum of abs) distance between two points. + * + * @param p1 the first point + * @param p2 the second point + * @return the L1 distance between the two points + */ + public static int distance1(int[] p1, int[] p2) { + int sum = 0; + for (int i = 0; i < p1.length; i++) { + sum += Math.abs(p1[i] - p2[i]); + } + return sum; + } + + /** + * Calculates the L2 (Euclidean) distance between two points. + * + * @param p1 the first point + * @param p2 the second point + * @return the L2 distance between the two points + */ + public static double distance(double[] p1, double[] p2) { + double sum = 0; + for (int i = 0; i < p1.length; i++) { + final double dp = p1[i] - p2[i]; + sum += dp * dp; + } + return Math.sqrt(sum); + } + + /** + * Calculates the L2 (Euclidean) distance between two points. + * + * @param p1 the first point + * @param p2 the second point + * @return the L2 distance between the two points + */ + public static double distance(int[] p1, int[] p2) { + int sum = 0; + for (int i = 0; i < p1.length; i++) { + final int dp = p1[i] - p2[i]; + sum += dp * dp; + } + return Math.sqrt(sum); + } + + /** + * Calculates the L∞ (max of abs) distance between two points. + * + * @param p1 the first point + * @param p2 the second point + * @return the L∞ distance between the two points + */ + public static double distanceInf(double[] p1, double[] p2) { + double max = 0; + for (int i = 0; i < p1.length; i++) { + max = Math.max(max, Math.abs(p1[i] - p2[i])); + } + return max; + } + + /** + * Calculates the L∞ (max of abs) distance between two points. + * + * @param p1 the first point + * @param p2 the second point + * @return the L∞ distance between the two points + */ + public static int distanceInf(int[] p1, int[] p2) { + int max = 0; + for (int i = 0; i < p1.length; i++) { + max = Math.max(max, Math.abs(p1[i] - p2[i])); + } + return max; + } + + +}