* φ = (φ + 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; /** * An array of the Fibonacci sequence. 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]
* 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, ...
*
* FIBONACCI_SEQUENCE[90] = 7540113804746346429
*
* See the following references for more information on * the fascinating properties of Fibonacci numbers: * *
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.
* That is, 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 / 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));
}
public static double[] abs(double[] x){
double[] result = new double[x.length];
for(int i=0; i= 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 Array of values to sum.
* @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 sum of the specified array of double values.
*
* @param xs Array of values to sum.
* @return The sum of the values.
*/
public static double sum(ArrayList xs) {
double sum = 0.0;
for (int i = 0; i < xs.size(); ++i)
sum += xs.get(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 Array of values.
* @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 Array of values.
* @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);
}
/**
* Throws an illegal argument exception if the specified double
* value is not a finite non-negative number. The label is
* included in the generated exception's error message if
* necessary.
*
* @param label Label of variable for exception message.
* @param x Double value to check.
* @throws IllegalArgumentException If the specified double value
* is infinite, is not a number, or is negative.
*/
public static void assertFiniteNonNegative(String label, double x) {
if (x < 0.0 ||
Double.isNaN(x)
|| Double.isInfinite(x)) {
String msg = label + " must be finite and non-negative."
+ " Found " + label + "=" + x;
throw new IllegalArgumentException(msg);
}
}
/**
* Returns the log (base 2) of the Γ function. The Γ function
* is defined by:
*
*
* Γ(z) = ∫0∞ tz-1 * e-t dt
*
* The Γ function is the continuous generalization of the factorial
* function, so that for real numbers z > 0:
*
*
Γ(z+1) = z * Γ(z)
*
* In particular, integers n >= 0, we have:
*
* Γ(n+1) = n!
*
* In general, Γ satisfies:
*
*
* Γ(z) = π / (sin(π * z) * Γ(1-z))
*
*
* This method uses the Lanczos approximation which is accurate
* nearly to the full power of double-precision arithmetic. The
* Lanczos approximation is used for inputs in the range
* [0.5,1.5], converting numbers less than 0.5 using
* the above formulas, and reducing arguments greater than 1.5
* using the factorial-like expansion above.
*
*
For more information on the Γ function and its computation, see:
*
*
* - Weisstein, Eric W. Gamma Function.
* From MathWorld--A Wolfram Web Resource.
*
* -
* Weisstein, Eric W. Lanczos Approximation. From MathWorld--A Wolfram Web Resource.
*
- Wikipedia. Gamma Function.
* - Wikipedia. Lanczos Approximation.
*
*
* @param z The argument to the gamma function.
* @return The value of Γ(z).
*/
public static double log2Gamma(double z) {
if (z < 0.5) {
return Mathematics.log2(java.lang.Math.PI)
- Mathematics.log2(java.lang.Math.sin(java.lang.Math.PI * z))
- log2Gamma(1.0 - z);
}
double result = 0.0;
while (z > 1.5) {
result += Mathematics.log2(z - 1);
z -= 1.0;
}
return result + Mathematics.log2(lanczosGamma(z));
}
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;
/**
* Returns the value of the digamma function for the specified
* value and number of iterations of the series approximation.
*
* The digamma function is the derivative of the log of the
* gamma function (defined in method {@link #log2Gamma(double)}:
*
*
* Ψ(z)
* = d log Γ(z) / dz
* = Γ'(z) / Γ(z)
*
*
* The numerical approximation is derived from formula 6.3.16 of:
*
*
* - Abramowitz and Stegun.
* Handbook of Mathematical Functions.
*
*
*
* The algorithm requires on the order of 50 terms for accuracy +/- 0.1,
* 200 terms for accuracy +/- 0.01, 2000 terms for accuracy +/- 0.001,
* 20,000 terms for accuracy +/- 0.0001, and 200,000 terms for accuracy
* +/0 0.00001.
*
* @param x Value at which to evaluate the digamma function.
* @param numTerms Number of terms in the series expansion to use.
* @return The value of the digamma function at the specified
* value.
*/
public static double digamma2(double x, int numTerms) {
// not defined for 0 or negative integers??
if (x <= 0.0 && (x == (double)((long) x)))
return Double.NaN;
double digamma = -NEGATIVE_DIGAMMA_1;
// reflection if less than 0
if (x < 0.0) {
digamma += java.lang.Math.PI
/ java.lang.Math.tan(java.lang.Math.PI * (1.0 - x));
x = 1.0 - x;
}
if (x < 1.0) {
digamma -= 1.0 / x;
x += 1.0;
}
if (x == 1.0)
return digamma;
if (x == 2.0)
return digamma + 1.0;
if (x == 3.0)
return digamma + 1.5;
while (x > 3.0) {
digamma += 1.0 / (x - 1.0);
x -= 1.0;
}
x -= 1.0;
for (int n = 1; n <= numTerms; ++n)
digamma += x / (n * (n + x));
return digamma;
}
/**
* The natural logarithm of 2.
*/
public static final double LN_2 = java.lang.Math.log(2.0);
/**
* The log base 2 of the constant e, which is the base of
* the natural logarithm. The constant e is determined by
* the java constant {@link java.lang.Math#E}.
*/
public static final double LOG2_E = log2(java.lang.Math.E);
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 };
}