com.opengamma.strata.math.impl.cern.Arithmetic Maven / Gradle / Ivy
Show all versions of strata-math Show documentation
/*
* Copyright (C) 2009 - present by OpenGamma Inc. and the OpenGamma group of companies
*
* Please see distribution for license.
*/
/*
* This code is copied from the original library from the `cern.jet.math` package.
* Changes:
* - package name
* - missing Javadoc param tags
* - reformat
* - remove unused methods
* - make package scoped
*/
/*
Copyright � 1999 CERN - European Organization for Nuclear Research.
Permission to use, copy, modify, distribute and sell this software and its documentation for any purpose
is hereby granted without fee, provided that the above copyright notice appear in all copies and
that both that copyright notice and this permission notice appear in supporting documentation.
CERN makes no representations about the suitability of this software for any purpose.
It is provided "as is" without expressed or implied warranty.
*/
package com.opengamma.strata.math.impl.cern;
//CSOFF: ALL
/**
* Arithmetic functions.
*/
class Arithmetic extends Constants {
// for method stirlingCorrection(...)
private static final double[] stirlingCorrection = {
0.0,
8.106146679532726e-02, 4.134069595540929e-02,
2.767792568499834e-02, 2.079067210376509e-02,
1.664469118982119e-02, 1.387612882307075e-02,
1.189670994589177e-02, 1.041126526197209e-02,
9.255462182712733e-03, 8.330563433362871e-03,
7.573675487951841e-03, 6.942840107209530e-03,
6.408994188004207e-03, 5.951370112758848e-03,
5.554733551962801e-03, 5.207655919609640e-03,
4.901395948434738e-03, 4.629153749334029e-03,
4.385560249232324e-03, 4.166319691996922e-03,
3.967954218640860e-03, 3.787618068444430e-03,
3.622960224683090e-03, 3.472021382978770e-03,
3.333155636728090e-03, 3.204970228055040e-03,
3.086278682608780e-03, 2.976063983550410e-03,
2.873449362352470e-03, 2.777674929752690e-03,
};
// for method logFactorial(...)
// log(k!) for k = 0, ..., 29
protected static final double[] logFactorials = {
0.00000000000000000, 0.00000000000000000, 0.69314718055994531,
1.79175946922805500, 3.17805383034794562, 4.78749174278204599,
6.57925121201010100, 8.52516136106541430, 10.60460290274525023,
12.80182748008146961, 15.10441257307551530, 17.50230784587388584,
19.98721449566188615, 22.55216385312342289, 25.19122118273868150,
27.89927138384089157, 30.67186010608067280, 33.50507345013688888,
36.39544520803305358, 39.33988418719949404, 42.33561646075348503,
45.38013889847690803, 48.47118135183522388, 51.60667556776437357,
54.78472939811231919, 58.00360522298051994, 61.26170176100200198,
64.55753862700633106, 67.88974313718153498, 71.25703896716800901
};
// k! for k = 0, ..., 20
protected static final long[] longFactorials = {
1L,
1L,
2L,
6L,
24L,
120L,
720L,
5040L,
40320L,
362880L,
3628800L,
39916800L,
479001600L,
6227020800L,
87178291200L,
1307674368000L,
20922789888000L,
355687428096000L,
6402373705728000L,
121645100408832000L,
2432902008176640000L
};
// k! for k = 21, ..., 170
protected static final double[] doubleFactorials = {
5.109094217170944E19,
1.1240007277776077E21,
2.585201673888498E22,
6.204484017332394E23,
1.5511210043330984E25,
4.032914611266057E26,
1.0888869450418352E28,
3.048883446117138E29,
8.841761993739701E30,
2.652528598121911E32,
8.222838654177924E33,
2.6313083693369355E35,
8.68331761881189E36,
2.952327990396041E38,
1.0333147966386144E40,
3.719933267899013E41,
1.3763753091226346E43,
5.23022617466601E44,
2.0397882081197447E46,
8.15915283247898E47,
3.34525266131638E49,
1.4050061177528801E51,
6.041526306337384E52,
2.6582715747884495E54,
1.196222208654802E56,
5.502622159812089E57,
2.5862324151116827E59,
1.2413915592536068E61,
6.082818640342679E62,
3.0414093201713376E64,
1.5511187532873816E66,
8.06581751709439E67,
4.274883284060024E69,
2.308436973392413E71,
1.2696403353658264E73,
7.109985878048632E74,
4.052691950487723E76,
2.350561331282879E78,
1.386831185456898E80,
8.32098711274139E81,
5.075802138772246E83,
3.146997326038794E85,
1.9826083154044396E87,
1.2688693218588414E89,
8.247650592082472E90,
5.443449390774432E92,
3.6471110918188705E94,
2.48003554243683E96,
1.7112245242814127E98,
1.1978571669969892E100,
8.504785885678624E101,
6.123445837688612E103,
4.470115461512686E105,
3.307885441519387E107,
2.4809140811395404E109,
1.8854947016660506E111,
1.451830920282859E113,
1.1324281178206295E115,
8.94618213078298E116,
7.15694570462638E118,
5.797126020747369E120,
4.7536433370128435E122,
3.94552396972066E124,
3.314240134565354E126,
2.8171041143805494E128,
2.4227095383672744E130,
2.107757298379527E132,
1.854826422573984E134,
1.6507955160908465E136,
1.4857159644817605E138,
1.3520015276784033E140,
1.2438414054641305E142,
1.156772507081641E144,
1.0873661566567426E146,
1.0329978488239061E148,
9.916779348709491E149,
9.619275968248216E151,
9.426890448883248E153,
9.332621544394415E155,
9.332621544394418E157,
9.42594775983836E159,
9.614466715035125E161,
9.902900716486178E163,
1.0299016745145631E166,
1.0813967582402912E168,
1.1462805637347086E170,
1.2265202031961373E172,
1.324641819451829E174,
1.4438595832024942E176,
1.5882455415227423E178,
1.7629525510902457E180,
1.974506857221075E182,
2.2311927486598138E184,
2.543559733472186E186,
2.925093693493014E188,
3.393108684451899E190,
3.96993716080872E192,
4.6845258497542896E194,
5.574585761207606E196,
6.689502913449135E198,
8.094298525273444E200,
9.875044200833601E202,
1.2146304367025332E205,
1.506141741511141E207,
1.882677176888926E209,
2.3721732428800483E211,
3.0126600184576624E213,
3.856204823625808E215,
4.974504222477287E217,
6.466855489220473E219,
8.471580690878813E221,
1.1182486511960037E224,
1.4872707060906847E226,
1.99294274616152E228,
2.690472707318049E230,
3.6590428819525483E232,
5.0128887482749884E234,
6.917786472619482E236,
9.615723196941089E238,
1.3462012475717523E241,
1.8981437590761713E243,
2.6953641378881633E245,
3.8543707171800694E247,
5.550293832739308E249,
8.047926057471989E251,
1.1749972043909107E254,
1.72724589045464E256,
2.5563239178728637E258,
3.8089226376305687E260,
5.7133839564458575E262,
8.627209774233244E264,
1.3113358856834527E267,
2.0063439050956838E269,
3.0897696138473515E271,
4.789142901463393E273,
7.471062926282892E275,
1.1729568794264134E278,
1.8532718694937346E280,
2.946702272495036E282,
4.714723635992061E284,
7.590705053947223E286,
1.2296942187394494E289,
2.0044015765453032E291,
3.287218585534299E293,
5.423910666131583E295,
9.003691705778434E297,
1.5036165148649983E300,
2.5260757449731988E302,
4.2690680090047056E304,
7.257415615308004E306
};
/**
* Makes this class non instantiable, but still let's others inherit from it.
*/
protected Arithmetic() {
}
/**
* Efficiently returns the binomial coefficient, often also referred to as "n over k" or "n choose k".
* The binomial coefficient is defined as (n * n-1 * ... * n-k+1 ) / ( 1 * 2 * ... * k ).
*
* - k<0: 0.
*
- k==0: 1.
*
- k==1: n.
*
- else: (n * n-1 * ... * n-k+1 ) / ( 1 * 2 * ... * k ).
*
* @param n n
* @param k k
* @return the binomial coefficient.
*/
public static double binomial(double n, long k) {
if (k < 0)
return 0;
if (k == 0)
return 1;
if (k == 1)
return n;
// binomial(n,k) = (n * n-1 * ... * n-k+1 ) / ( 1 * 2 * ... * k )
double a = n - k + 1;
double b = 1;
double binomial = 1;
for (long i = k; i-- > 0;) {
binomial *= (a++) / (b++);
}
return binomial;
}
/**
* Efficiently returns the binomial coefficient, often also referred to as "n over k" or "n choose k".
* The binomial coefficient is defined as
*
* - k<0: 0.
*
- k==0 || k==n: 1.
*
- k==1 || k==n-1: n.
*
- else: (n * n-1 * ... * n-k+1 ) / ( 1 * 2 * ... * k ).
*
* @param n n
* @param inK k
* @return the binomial coefficient.
*/
public static double binomial(long n, long inK) {
long k = inK;
if (k < 0)
return 0;
if (k == 0 || k == n)
return 1;
if (k == 1 || k == n - 1)
return n;
// try quick version and see whether we get numeric overflows.
// factorial(..) is O(1); requires no loop; only a table lookup.
if (n > k) {
int max = longFactorials.length + doubleFactorials.length;
if (n < max) { // if (n! < inf && k! < inf)
double n_fac = factorial((int) n);
double k_fac = factorial((int) k);
double n_minus_k_fac = factorial((int) (n - k));
double nk = n_minus_k_fac * k_fac;
if (nk != Double.POSITIVE_INFINITY) { // no numeric overflow?
// now this is completely safe and accurate
return n_fac / nk;
}
}
if (k > n / 2)
k = n - k; // quicker
}
// binomial(n,k) = (n * n-1 * ... * n-k+1 ) / ( 1 * 2 * ... * k )
long a = n - k + 1;
long b = 1;
double binomial = 1;
for (long i = k; i-- > 0;) {
binomial *= ((double) (a++)) / (b++);
}
return binomial;
}
/**
* Returns the smallest long >= value.
* Examples: 1.0 -> 1, 1.2 -> 2, 1.9 -> 2.
* This method is safer than using (long) Math.ceil(value), because of possible rounding error.
* @param value v
* @return result
*/
public static long ceil(double value) {
return Math.round(Math.ceil(value));
}
/**
* Evaluates the series of Chebyshev polynomials Ti at argument x/2.
* The series is given by
*
* N-1
* - '
* y = > coef[i] T (x/2)
* - i
* i=0
*
* Coefficients are stored in reverse order, i.e. the zero
* order term is last in the array. Note N is the number of
* coefficients, not the order.
*
* If coefficients are for the interval a to b, x must
* have been transformed to x -> 2(2x - b - a)/(b-a) before
* entering the routine. This maps x from (a, b) to (-1, 1),
* over which the Chebyshev polynomials are defined.
*
* If the coefficients are for the inverted interval, in
* which (a, b) is mapped to (1/b, 1/a), the transformation
* required is x -> 2(2ab/x - b - a)/(b-a). If b is infinity,
* this becomes x -> 4a/x - 1.
*
* SPEED:
*
* Taking advantage of the recurrence properties of the
* Chebyshev polynomials, the routine requires one more
* addition per loop than evaluating a nested polynomial of
* the same degree.
*
* @param x argument to the polynomial.
* @param coef the coefficients of the polynomial.
* @param N the number of coefficients.
* @return result
* @throws ArithmeticException error
*/
public static double chbevl(double x, double coef[], int N) throws ArithmeticException {
double b0, b1, b2;
int p = 0;
int i;
b0 = coef[p++];
b1 = 0.0;
i = N - 1;
do {
b2 = b1;
b1 = b0;
b0 = x * b1 - b2 + coef[p++];
} while (--i > 0);
return (0.5 * (b0 - b2));
}
/**
* Instantly returns the factorial k!.
* @param k must hold k >= 0.
* @return result
*/
static public double factorial(int k) {
if (k < 0)
throw new IllegalArgumentException();
int length1 = longFactorials.length;
if (k < length1)
return longFactorials[k];
int length2 = doubleFactorials.length;
if (k < length1 + length2)
return doubleFactorials[k - length1];
else
return Double.POSITIVE_INFINITY;
}
/**
* Returns the largest long <= value.
*
Examples:
* 1.0 -> 1, 1.2 -> 1, 1.9 -> 1
* 2.0 -> 2, 2.2 -> 2, 2.9 -> 2
* This method is safer than using (long) Math.floor(value), because of possible rounding error.
* @param value v
* @return result
*/
public static long floor(double value) {
return Math.round(Math.floor(value));
}
/**
* Returns logbasevalue.
* @param base base
* @param value v
* @return result
*/
public static double log(double base, double value) {
return Math.log(value) / Math.log(base);
}
/**
* Returns log10value.
* @param value v
* @return result
*/
static public double log10(double value) {
// 1.0 / Math.log(10) == 0.43429448190325176
return Math.log(value) * 0.43429448190325176;
}
/**
* Returns log2value.
* @param value v
* @return result
*/
static public double log2(double value) {
// 1.0 / Math.log(2) == 1.4426950408889634
return Math.log(value) * 1.4426950408889634;
}
/**
* Returns log(k!).
* Tries to avoid overflows.
* For k<30 simply looks up a table in O(1).
* For k>=30 uses stirlings approximation.
* @param k must hold k >= 0.
* @return result
*/
public static double logFactorial(int k) {
if (k >= 30) {
double r, rr;
final double C0 = 9.18938533204672742e-01;
final double C1 = 8.33333333333333333e-02;
final double C3 = -2.77777777777777778e-03;
final double C5 = 7.93650793650793651e-04;
final double C7 = -5.95238095238095238e-04;
r = 1.0 / (double) k;
rr = r * r;
return (k + 0.5) * Math.log(k) - k + C0 + r * (C1 + rr * (C3 + rr * (C5 + rr * C7)));
} else
return logFactorials[k];
}
/**
* Instantly returns the factorial k!.
* @param k must hold k >= 0 && k < 21.
* @return result
* @throws IllegalArgumentException error
*/
static public long longFactorial(int k) throws IllegalArgumentException {
if (k < 0)
throw new IllegalArgumentException("Negative k");
if (k < longFactorials.length)
return longFactorials[k];
throw new IllegalArgumentException("Overflow");
}
/**
* Returns the StirlingCorrection.
*
* Correction term of the Stirling approximation for log(k!)
* (series in 1/k, or table values for small k)
* with int parameter k.
*
*
* log k! = (k + 1/2)log(k + 1) - (k + 1) + (1/2)log(2Pi) +
* stirlingCorrection(k + 1)
*
* log k! = (k + 1/2)log(k) - k + (1/2)log(2Pi) +
* stirlingCorrection(k)
*
* @param k k
* @return result
*/
public static double stirlingCorrection(int k) {
final double C1 = 8.33333333333333333e-02; // +1/12
final double C3 = -2.77777777777777778e-03; // -1/360
final double C5 = 7.93650793650793651e-04; // +1/1260
final double C7 = -5.95238095238095238e-04; // -1/1680
double r, rr;
if (k > 30) {
r = 1.0 / (double) k;
rr = r * r;
return r * (C1 + rr * (C3 + rr * (C5 + rr * C7)));
} else
return stirlingCorrection[k];
}
}