/*
* LingPipe v. 3.9
* Copyright (C) 2003-2010 Alias-i
*
* This program is licensed under the Alias-i Royalty Free License
* Version 1 WITHOUT ANY WARRANTY, without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the Alias-i
* Royalty Free License Version 1 for more details.
*
* You should have received a copy of the Alias-i Royalty Free License
* Version 1 along with this program; if not, visit
* http://alias-i.com/lingpipe/licenses/lingpipe-license-1.txt or contact
* Alias-i, Inc. at 181 North 11th Street, Suite 401, Brooklyn, NY 11211,
* +1 (718) 290-9170.
*/
//package com.aliasi.util;
/**
* A math utility class with static methods.
*
* @author Bob Carpenter
* @version 4.0.0
* @since LingPipe1.0
*/
public class Math {
// forbid instances
private Math() {
/* no instances */
}
/**
* The value of the golden ratio. The golden ratio is defined to
* be the value φ such that:
*
*
* φ = (φ + 1) / φ
*
*
* Note that this is a quadratic equation (multiply both sides by
* φ) with the solution roughly 1.61803399.
*
* See the following for a fascinating tour of the properties
* of the golden ratio:
*
*
*
*
*/
public static final double GOLDEN_RATIO = (1.0 + java.lang.Math.sqrt(5))/2.0;
/**
* The natural logarithm of 2.
*/
public static final double LN_2 = java.lang.Math.log(2.0);
static final double INV_LN_2 = 1.0/LN_2;
/**
* An array of the Fibonacci sequence up the maximum value
* representable as a long integer. The array is defined as
* follows:
*
*
* FIBONACCI_SEQUENCE[0] = 1
* FIBONACCI_SEQUENCE[1] = 2
* FIBONACCI_SEQUENCE[n+2] = FIBONACCI_SEQUENCE[n+1] + FIBONACCI_SEQUENCE[n]
*
*
* So FIBONACCI_SEQUENCE[0] represents the second
* Fibonacci number in the traditional numbering. The inital entries
* are:
*
*
* 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597,
* 2584, ...
*
*
* The length of the array is 91, and the largest value is:
*
*
* FIBONACCI_SEQUENCE[90] = 7540113804746346429
*
*
*
* See the following references for more information on
* the fascinating properties of Fibonacci numbers:
*
*
* -
Wikipedia: Fibonacci Number
* -
Mathworld: Fibonacci Number
*/
public static final long[] FIBONACCI_SEQUENCE = new long[] {
1l,
2l,
3l,
5l,
8l,
13l,
21l,
34l,
55l,
89l,
144l,
233l,
377l,
610l,
987l,
1597l,
2584l,
4181l,
6765l,
10946l,
17711l,
28657l,
46368l,
75025l,
121393l,
196418l,
317811l,
514229l,
832040l,
1346269l,
2178309l,
3524578l,
5702887l,
9227465l,
14930352l,
24157817l,
39088169l,
63245986l,
102334155l,
165580141l,
267914296l,
433494437l,
701408733l,
1134903170l,
1836311903l,
2971215073l,
4807526976l,
7778742049l,
12586269025l,
20365011074l,
32951280099l,
53316291173l,
86267571272l,
139583862445l,
225851433717l,
365435296162l,
591286729879l,
956722026041l,
1548008755920l,
2504730781961l,
4052739537881l,
6557470319842l,
10610209857723l,
17167680177565l,
27777890035288l,
44945570212853l,
72723460248141l,
117669030460994l,
190392490709135l,
308061521170129l,
498454011879264l,
806515533049393l,
1304969544928657l,
2111485077978050l,
3416454622906707l,
5527939700884757l,
8944394323791464l,
14472334024676221l,
23416728348467685l,
37889062373143906l,
61305790721611591l,
99194853094755497l,
160500643816367088l,
259695496911122585l,
420196140727489673l,
679891637638612258l,
1100087778366101931l,
1779979416004714189l,
2880067194370816120l,
4660046610375530309l,
7540113804746346429l
};
/**
* Returns true if the specified number is prime. A
* prime is a positive number greater than 1 with no
* divisors other than 1 and itself, thus
* {2,3,5,7,11,13,...}.
*
* @param num Number to test for primality.
* @return true if the specified number is prime.
*/
public static boolean isPrime(int num) {
if (num < 2) return false;
for (int i = 2; i <= num/2; ++i)
if (num % i == 0) return false;
return true;
}
/**
* Returns the smallest prime number that is strictly larger than
* the specified integer. See {@link #isPrime(int)} for the
* definition of primality.
*
* @param num Base from which to look for the next prime.
* @return Smallest prime number strictly larget than specified
* number.
*/
public static int nextPrime(int num) {
if (num < 2) return 2;
for (int i = num + 1; ; ++i)
if (isPrime(i)) return i;
}
/**
* Converts a natural logarithm to a base 2 logarithm.
* This inverts the operation of {@link #logBase2ToNaturalLog(double)}.
* If the input is x = ln z, then
* the return value is log2 z.
* Recall that log2 z = ln z / ln 2.
*
* @param x Natural log of value.
* @return Log base 2 of value.
*/
public static double naturalLogToBase2Log(double x) {
return x * INV_LN_2;
}
/**
* Returns the log base 2 of the specivied value.
*
* @param x Value whose log is taken.
* @return Log of specified value.
*/
public static double log2(double x) {
return naturalLogToBase2Log(java.lang.Math.log(x));
}
/**
* Returns the integer value of reading the specified byte as an
* unsigned value. The computation is carried out by subtracting
* the minimum value, as defined by the constant {@link
* Byte#MIN_VALUE}.
*
* @param b Byte to convert.
* @return Unsigned value of specified byte.
*/
public static int byteAsUnsigned(byte b) {
return (b >= 0) ? (int)b : (256+(int)b);
}
/**
* Returns the log (base 2) of the factorial of the specified long
* integer. The factorial of n is defined for
* n > 0 by:
*
*
* n!
* = Πi < 0 <= n i
*
*
* Taking logs of both sides gives:
*
*
* log2 n!
* = Σi < 0 <= n
* log2 i
*
*
* By convention, 0! is taken to be 1, and hence ln 0! = 0.
*
* @param n Specified long integer.
* @return Log of factorial of specified integer.
* @throws IllegalArgumentException If the argument is negative.
*/
public static double log2Factorial(long n) {
if (n < 0) {
String msg = "Factorials only defined for non-negative arguments."
+ " Found argument=" + n;
throw new IllegalArgumentException(msg);
}
double sum = 0.0;
for (long i = 1; i <= n; ++i)
sum += log2(i);
return sum;
}
/**
* Returns the sum of the specified array of double values.
*
* @param xs Variable length list of values, or an array of values.
* @return The sum of the values.
*/
public static double sum(double... xs) {
double sum = 0.0;
for (int i = 0; i < xs.length; ++i)
sum += xs[i];
return sum;
}
/**
* Returns the minimum of the specified array of double values.
* If the length of the array is zero, the result is {@link
* Double#NaN}.
*
* @param xs Variable length list of values, or an array.
* @return Minimum value in array.
*/
public static double minimum(double... xs) {
if (xs.length == 0) return Double.NaN;
double min = xs[0];
for (int i = 1; i < xs.length; ++i)
if (xs[i] < min) min = xs[i];
return min;
}
/**
* Returns the maximum of the specified array of double values.
* If the length of the array is zero, the result is {@link
* Double#NaN}.
*
* @param xs Variable length list of values, or an array.
* @return Maximum value in array.
*/
public static double maximum(double... xs) {
if (xs.length == 0) return Double.NaN;
double max = xs[0];
for (int i = 1; i < xs.length; ++i)
if (xs[i] > max) max = xs[i];
return max;
}
/**
* Returns the log (base 2) of the binomial coefficient of the
* specified arguments. The binomial coefficient is equal to the
* number of ways to choose a subset of size m from a
* set of n objects, which is pronounced "n choose
* m", and is given by:
*
*
* choose(n,m) = n! / ( m! * (n-m)!)
*
* log2 choose(n,m)
* = log2 n - log2 m
* - log2 (n-m)
*
*
* @return The log (base 2) of the binomial coefficient of the
* specified arguments.
*/
public static double log2BinomialCoefficient(long n, long m) {
return log2(n) - log2(m) - log2(n-m);
}
static double[] LANCZOS_COEFFS = new double[] {
0.99999999999980993,
676.5203681218851,
-1259.1392167224028,
771.32342877765313,
-176.61502916214059,
12.507343278686905,
-0.13857109526572012,
9.9843695780195716e-6,
1.5056327351493116e-7
};
static double SQRT_2_PI = java.lang.Math.sqrt(2.0 * java.lang.Math.PI);
// assumes input in [0.5,1.5] inclusive
static double lanczosGamma(double z) {
double zMinus1 = z - 1;
double x = LANCZOS_COEFFS[0];
for (int i = 1; i < LANCZOS_COEFFS.length - 2; ++i)
x += LANCZOS_COEFFS[i] / (zMinus1 + i);
double t = zMinus1 + (LANCZOS_COEFFS.length - 2) + 0.5;
return SQRT_2_PI
* java.lang.Math.pow(t, zMinus1 + 0.5)
* java.lang.Math.exp(-t) * x;
}
/**
* Returns the value of the digamma function for the specified
* value. The returned values are accurate to at least 13
* decimal places.
*
* The digamma function is the derivative of the log of the
* gamma function; see the method documentation for {@link
* #log2Gamma(double)} for more information on the gamma function
* itself.
*
*
* Ψ(z)
* = d log Γ(z) / dz
* = Γ'(z) / Γ(z)
*
*
* The numerical approximation is derived from:
*
*
* - Richard J. Mathar. 2005.
* Chebyshev Series Expansion of Inverse Polynomials.
*
* - Richard J. Mathar. 2005.
* digamma.c.
* (C Program implementing algorithm.)
*
*
*
* Implementation Note: The recursive calls in the C
* implementation have been transformed into loops and
* accumulators, and the recursion for values greater than three
* replaced with a simpler reduction. The number of loops
* required before the fixed length expansion is approximately
* integer value of the absolute value of the input. Each loop
* requires a floating point division, two additions and a local
* variable assignment. The fixed portion of the algorithm is
* roughly 30 steps requiring four multiplications, three
* additions, one static final array lookup, and four assignments per
* loop iteration.
*
* @param x Value at which to evaluate the digamma function.
* @return The value of the digamma function at the specified
* value.
*/
public static double digamma(double x)
{
if (x <= 0.0 && (x == (double)((long) x)))
return Double.NaN;
double accum = 0.0;
if (x < 0.0) {
accum += java.lang.Math.PI
/ java.lang.Math.tan(java.lang.Math.PI * (1.0 - x));
x = 1.0 - x;
}
if (x < 1.0 ) {
while (x < 1.0)
accum -= 1.0 / x++;
}
if (x == 1.0)
return accum - NEGATIVE_DIGAMMA_1;
if (x == 2.0)
return accum + 1.0 - NEGATIVE_DIGAMMA_1;
if (x == 3.0)
return accum + 1.5 - NEGATIVE_DIGAMMA_1;
// simpler recursion than Mahar to reduce recursion
if (x > 3.0) {
while (x > 3.0)
accum += 1.0 / --x;
return accum + digamma(x);
}
x -= 2.0;
double tNMinus1 = 1.0;
double tN = x;
double digamma = DIGAMMA_COEFFS[0] + DIGAMMA_COEFFS[1] * tN;
for (int n = 2; n < DIGAMMA_COEFFS.length; n++) {
double tN1 = 2.0 * x * tN - tNMinus1;
digamma += DIGAMMA_COEFFS[n] * tN1;
tNMinus1 = tN;
tN = tN1;
}
return accum + digamma;
}
/**
* The γ constant for computing the digamma function.
*
* The value is defined as the negative of the digamma funtion
* evaluated at 1:
*
*
* γ = - Ψ(1)
*
*/
static double NEGATIVE_DIGAMMA_1 = 0.5772156649015328606065120900824024;
private static final double DIGAMMA_COEFFS[]
= {
.30459198558715155634315638246624251,
.72037977439182833573548891941219706,
-.12454959243861367729528855995001087,
.27769457331927827002810119567456810e-1,
-.67762371439822456447373550186163070e-2,
.17238755142247705209823876688592170e-2,
-.44817699064252933515310345718960928e-3,
.11793660000155572716272710617753373e-3,
-.31253894280980134452125172274246963e-4,
.83173997012173283398932708991137488e-5,
-.22191427643780045431149221890172210e-5,
.59302266729329346291029599913617915e-6,
-.15863051191470655433559920279603632e-6,
.42459203983193603241777510648681429e-7,
-.11369129616951114238848106591780146e-7,
.304502217295931698401459168423403510e-8,
-.81568455080753152802915013641723686e-9,
.21852324749975455125936715817306383e-9,
-.58546491441689515680751900276454407e-10,
.15686348450871204869813586459513648e-10,
-.42029496273143231373796179302482033e-11,
.11261435719264907097227520956710754e-11,
-.30174353636860279765375177200637590e-12,
.80850955256389526647406571868193768e-13,
-.21663779809421233144009565199997351e-13,
.58047634271339391495076374966835526e-14,
-.15553767189204733561108869588173845e-14,
.41676108598040807753707828039353330e-15,
-.11167065064221317094734023242188463e-15 };
/**
* Returns the relative absolute difference between the specified
* values, defined to be:
*
*
* relAbsDiff(x,y) = abs(x-y) / (abs(x) + abs(y))
*
* @param x First value.
* @param y Second value.
* @return The absolute relative difference between the values.
*/
public static double relativeAbsoluteDifference(double x, double y) {
return (Double.isInfinite(x) || Double.isInfinite(y))
? Double.POSITIVE_INFINITY
: (java.lang.Math.abs(x - y)
/ (java.lang.Math.abs(x) + java.lang.Math.abs(y)));
}
/**
* This method returns the log of the sum of the natural
* exponentiated values in the specified array. Mathematically,
* the result is
*
*
* logSumOfExponentials(xs) = log ( Σi exp(xs[i]) )
*
* But the result is not calculated directly. Instead, the
* calculation performed is:
*
*
* logSumOfExponentials(xs) = max(xs) + log ( Σi exp(xs[i] - max(xs)) )
*
* which produces the same result, but is much more arithmetically
* stable, because the largest value for which exp()
* is calculated is 0.0.
*
* Values of {@code Double.NEGATIVE_INFINITY} are treated as
* having exponentials of 0 and logs of negative infinity.
* That is, they are ignored for the purposes of this computation.
*
* @param xs Array of values.
* @return The log of the sum of the exponentiated values in the
* array.
*/
public static double logSumOfExponentials(double[] xs) {
if (xs.length == 1) return xs[0];
double max = maximum(xs);
double sum = 0.0;
for (int i = 0; i < xs.length; ++i)
if (xs[i] != Double.NEGATIVE_INFINITY)
sum += java.lang.Math.exp(xs[i] - max);
return max + java.lang.Math.log(sum);
}
/**
* Returns the maximum value of an element in xs. If any of the
* values are {@code Double.NaN}, or if the input array is empty,
* the result is {@code Double.NaN}.
*
* @param xs Array in which to find maximum.
* @return Maximum value in array.
*/
public static double max(double... xs) {
if (xs.length == 0)
return Double.NaN;
double max = xs[0];
for (int i = 1; i < xs.length; ++i)
max = java.lang.Math.max(max,xs[i]);
return max;
}
/**
* Returns the maximum value of an element in the specified array.
*
* @param xs Array in which to find maximum.
* @return Maximum value in the array.
* @throws ArrayIndexOutOfBoundsException If the specified array does
* not contai at least one element.
*/
public static int max(int... xs) {
int max = xs[0];
for (int i = 1; i < xs.length; ++i)
if (xs[i] > max)
max = xs[i];
return max;
}
/**
* Returns the sum of the specified integer array. Note that
* there is no check for overflow. If the array is of length 0,
* the sum is defined to be 0.
*
* @param xs Array of integers to sum.
* @return Sum of the array.
*/
public static int sum(int... xs) {
int sum = 0;
for (int i = 0; i < xs.length; ++i)
sum += xs[i];
return sum;
}
}