Added Stirling numbers of the second kind in ArithmeticUtils.
git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1371680 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe
parent
f2fe724505
commit
add45005fc
|
|
@ -17,6 +17,8 @@
|
|||
package org.apache.commons.math3.util;
|
||||
|
||||
import java.math.BigInteger;
|
||||
import java.util.concurrent.atomic.AtomicReference;
|
||||
|
||||
import org.apache.commons.math3.exception.MathArithmeticException;
|
||||
import org.apache.commons.math3.exception.NotPositiveException;
|
||||
import org.apache.commons.math3.exception.NumberIsTooLargeException;
|
||||
|
|
@ -41,6 +43,9 @@ public final class ArithmeticUtils {
|
|||
1307674368000l, 20922789888000l, 355687428096000l,
|
||||
6402373705728000l, 121645100408832000l, 2432902008176640000l };
|
||||
|
||||
/** Stirling numbers of the second kind. */
|
||||
static final AtomicReference<long[][]> STIRLING_S2 = new AtomicReference<long[][]> (null);
|
||||
|
||||
/** Private constructor. */
|
||||
private ArithmeticUtils() {
|
||||
super();
|
||||
|
|
@ -882,6 +887,91 @@ public final class ArithmeticUtils {
|
|||
return result;
|
||||
}
|
||||
|
||||
/**
|
||||
* Returns the <a
|
||||
* href="http://mathworld.wolfram.com/StirlingNumberoftheSecondKind.html">
|
||||
* Stirling number of the second kind</a>, "{@code S(n,k)}", the number of
|
||||
* ways of partitioning an {@code n}-element set into {@code k} non-empty
|
||||
* subsets.
|
||||
* <p>
|
||||
* The preconditions are {@code 0 <= k <= n } (otherwise
|
||||
* {@code NotPositiveException} is thrown)
|
||||
* </p>
|
||||
* @param n the size of the set
|
||||
* @param k the number of non-empty subsets
|
||||
* @return {@code S(n,k)}
|
||||
* @throws NotPositiveException if {@code k < 0}.
|
||||
* @throws NumberIsTooLargeException if {@code k > n}.
|
||||
* @throws MathArithmeticException if some overflow happens, typically for n exceeding 25 and
|
||||
* k between 20 and n-2 (S(n,n-1) is handled specifically and does not overflow)
|
||||
*/
|
||||
public static long stirlingS2(final int n, final int k)
|
||||
throws NotPositiveException, NumberIsTooLargeException, MathArithmeticException {
|
||||
if (k < 0) {
|
||||
throw new NotPositiveException(k);
|
||||
}
|
||||
if (k > n) {
|
||||
throw new NumberIsTooLargeException(k, n, true);
|
||||
}
|
||||
|
||||
long[][] stirlingS2 = STIRLING_S2.get();
|
||||
|
||||
if (stirlingS2 == null) {
|
||||
// the cache has never been initialized, compute the first numbers
|
||||
// by direct recurrence relation
|
||||
|
||||
// as S(26,9) = 11201516780955125625 is larger than Long.MAX_VALUE
|
||||
// we must stop computation at row 26
|
||||
final int maxIndex = 26;
|
||||
stirlingS2 = new long[maxIndex][];
|
||||
stirlingS2[0] = new long[] { 1l };
|
||||
for (int i = 1; i < stirlingS2.length; ++i) {
|
||||
stirlingS2[i] = new long[i + 1];
|
||||
stirlingS2[i][0] = 0;
|
||||
stirlingS2[i][1] = 1;
|
||||
stirlingS2[i][i] = 1;
|
||||
for (int j = 2; j < i; ++j) {
|
||||
stirlingS2[i][j] = j * stirlingS2[i - 1][j] + stirlingS2[i - 1][j - 1];
|
||||
}
|
||||
}
|
||||
|
||||
// atomically save the cache
|
||||
STIRLING_S2.compareAndSet(null, stirlingS2);
|
||||
|
||||
}
|
||||
|
||||
if (n < stirlingS2.length) {
|
||||
// the number is in the small cache
|
||||
return stirlingS2[n][k];
|
||||
} else {
|
||||
// use explicit formula to compute the number without caching it
|
||||
if (k == 0) {
|
||||
return 0;
|
||||
} else if (k == 1 || k == n) {
|
||||
return 1;
|
||||
} else if (k == 2) {
|
||||
return (1l << (n - 1)) - 1l;
|
||||
} else if (k == n - 1) {
|
||||
return binomialCoefficient(n, 2);
|
||||
} else {
|
||||
// definition formula: note that this may trigger some overflow
|
||||
long sum = 0;
|
||||
long sign = ((k & 0x1) == 0) ? 1 : -1;
|
||||
for (int j = 1; j <= k; ++j) {
|
||||
sign = -sign;
|
||||
sum += sign * binomialCoefficient(k, j) * pow(j, n);
|
||||
if (sum < 0) {
|
||||
// there was an overflow somewhere
|
||||
throw new MathArithmeticException(LocalizedFormats.ARGUMENT_OUTSIDE_DOMAIN,
|
||||
n, 0, stirlingS2.length - 1);
|
||||
}
|
||||
}
|
||||
return sum / factorial(k);
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
/**
|
||||
* Add two long integers, checking for overflow.
|
||||
*
|
||||
|
|
|
|||
|
|
@ -89,7 +89,7 @@
|
|||
<li><a href="utilities.html#a6.2_Double_array_utilities">6.2 Double array utilities</a></li>
|
||||
<li><a href="utilities.html#a6.3_intdouble_hash_map">6.3 int/double hash map</a></li>
|
||||
<li><a href="utilities.html#a6.4_Continued_Fractions">6.4 Continued Fractions</a></li>
|
||||
<li><a href="utilities.html#a6.5_binomial_coefficients_factorials_and_other_common_math_functions">6.5 Binomial coefficients, factorials and other common math functions</a></li>
|
||||
<li><a href="utilities.html#a6.5_binomial_coefficients_factorials_Stirling_numbers_and_other_common_math_functions">6.5 Binomial coefficients, factorials, Stirling numbers and other common math functions</a></li>
|
||||
<li><a href="utilities.html#a6.6_fast_math">6.6 Fast mathematical functions</a></li>
|
||||
<li><a href="utilities.html#a6.7_miscellaneous">6.7 Miscellaneous</a></li>
|
||||
</ul></li>
|
||||
|
|
|
|||
|
|
@ -146,11 +146,11 @@
|
|||
</p>
|
||||
</subsection>
|
||||
|
||||
<subsection name="6.5 Binomial coefficients, factorials and other common math functions" href="binomial_coefficients_factorials_and_other_common_math_functions">
|
||||
<subsection name="6.5 Binomial coefficients, factorials, Stirling numbers and other common math functions" href="binomial_coefficients_factorials_and_other_common_math_functions">
|
||||
<p>
|
||||
A collection of reusable math functions is provided in the
|
||||
<a href="../apidocs/org/apache/commons/math3/util/MathUtils.html">MathUtils</a>
|
||||
utility class. MathUtils currently includes methods to compute the following: <ul>
|
||||
<a href="../apidocs/org/apache/commons/math3/util/ArithmeticUtils.html">ArithmeticUtils</a>
|
||||
utility class. ArithmeticUtils currently includes methods to compute the following: <ul>
|
||||
<li>
|
||||
Binomial coefficients -- "n choose k" available as an (exact) long value,
|
||||
<code>binomialCoefficient(int, int)</code> for small n, k; as a double,
|
||||
|
|
@ -158,24 +158,13 @@
|
|||
a "super-sized" version, <code>binomialCoefficientLog(int, int)</code>
|
||||
that returns the natural logarithm of the value.</li>
|
||||
<li>
|
||||
Stirling numbers of the second kind -- S(n,k) as an exact long value
|
||||
<code>stirlingS2(int, int)</code> for small n, k.</li>
|
||||
<li>
|
||||
Factorials -- like binomial coefficients, these are available as exact long
|
||||
values, <code>factorial(int)</code>; doubles,
|
||||
<code>factorialDouble(int)</code>; or logs, <code>factorialLog(int)</code>. </li>
|
||||
<li>
|
||||
Hyperbolic sine and cosine functions --
|
||||
<code>cosh(double), sinh(double)</code></li>
|
||||
<li>
|
||||
sign (+1 if argument > 0, 0 if x = 0, and -1 if x < 0) and
|
||||
indicator (+1.0 if argument >= 0 and -1.0 if argument < 0) functions
|
||||
for variables of all primitive numeric types.</li>
|
||||
<li>
|
||||
a hash function, <code>hash(double),</code> returning a long-valued
|
||||
hash code for a double value.
|
||||
</li>
|
||||
<li>
|
||||
Convience methods to round floating-point number to arbitrary precision.
|
||||
</li>
|
||||
<li>
|
||||
Least common multiple and greatest common denominator functions.
|
||||
</li>
|
||||
</ul>
|
||||
|
|
@ -226,6 +215,7 @@
|
|||
<li>asinh(double)</li>
|
||||
<li>acosh(double)</li>
|
||||
<li>atanh(double)</li>
|
||||
<li>pow(double,int)</li>
|
||||
</ul>
|
||||
The following methods are found in Math/StrictMath since 1.6 only, they are
|
||||
provided by FastMath even in 1.5 Java virtual machines
|
||||
|
|
|
|||
|
|
@ -25,6 +25,8 @@ import java.math.BigInteger;
|
|||
|
||||
import org.apache.commons.math3.exception.MathArithmeticException;
|
||||
import org.apache.commons.math3.exception.MathIllegalArgumentException;
|
||||
import org.apache.commons.math3.exception.NotPositiveException;
|
||||
import org.apache.commons.math3.exception.NumberIsTooLargeException;
|
||||
import org.apache.commons.math3.random.RandomDataImpl;
|
||||
import org.junit.Assert;
|
||||
import org.junit.Test;
|
||||
|
|
@ -693,6 +695,63 @@ public class ArithmeticUtilsTest {
|
|||
}
|
||||
}
|
||||
|
||||
@Test
|
||||
public void testStirlingS2() {
|
||||
|
||||
Assert.assertEquals(1, ArithmeticUtils.stirlingS2(0, 0));
|
||||
|
||||
for (int n = 1; n < 30; ++n) {
|
||||
Assert.assertEquals(0, ArithmeticUtils.stirlingS2(n, 0));
|
||||
Assert.assertEquals(1, ArithmeticUtils.stirlingS2(n, 1));
|
||||
if (n > 2) {
|
||||
Assert.assertEquals((1l << (n - 1)) - 1l, ArithmeticUtils.stirlingS2(n, 2));
|
||||
Assert.assertEquals(ArithmeticUtils.binomialCoefficient(n, 2),
|
||||
ArithmeticUtils.stirlingS2(n, n - 1));
|
||||
}
|
||||
Assert.assertEquals(1, ArithmeticUtils.stirlingS2(n, n));
|
||||
}
|
||||
Assert.assertEquals(536870911l, ArithmeticUtils.stirlingS2(30, 2));
|
||||
Assert.assertEquals(576460752303423487l, ArithmeticUtils.stirlingS2(60, 2));
|
||||
|
||||
Assert.assertEquals( 25, ArithmeticUtils.stirlingS2( 5, 3));
|
||||
Assert.assertEquals( 90, ArithmeticUtils.stirlingS2( 6, 3));
|
||||
Assert.assertEquals( 65, ArithmeticUtils.stirlingS2( 6, 4));
|
||||
Assert.assertEquals( 301, ArithmeticUtils.stirlingS2( 7, 3));
|
||||
Assert.assertEquals( 350, ArithmeticUtils.stirlingS2( 7, 4));
|
||||
Assert.assertEquals( 140, ArithmeticUtils.stirlingS2( 7, 5));
|
||||
Assert.assertEquals( 966, ArithmeticUtils.stirlingS2( 8, 3));
|
||||
Assert.assertEquals( 1701, ArithmeticUtils.stirlingS2( 8, 4));
|
||||
Assert.assertEquals( 1050, ArithmeticUtils.stirlingS2( 8, 5));
|
||||
Assert.assertEquals( 266, ArithmeticUtils.stirlingS2( 8, 6));
|
||||
Assert.assertEquals( 3025, ArithmeticUtils.stirlingS2( 9, 3));
|
||||
Assert.assertEquals( 7770, ArithmeticUtils.stirlingS2( 9, 4));
|
||||
Assert.assertEquals( 6951, ArithmeticUtils.stirlingS2( 9, 5));
|
||||
Assert.assertEquals( 2646, ArithmeticUtils.stirlingS2( 9, 6));
|
||||
Assert.assertEquals( 462, ArithmeticUtils.stirlingS2( 9, 7));
|
||||
Assert.assertEquals( 9330, ArithmeticUtils.stirlingS2(10, 3));
|
||||
Assert.assertEquals(34105, ArithmeticUtils.stirlingS2(10, 4));
|
||||
Assert.assertEquals(42525, ArithmeticUtils.stirlingS2(10, 5));
|
||||
Assert.assertEquals(22827, ArithmeticUtils.stirlingS2(10, 6));
|
||||
Assert.assertEquals( 5880, ArithmeticUtils.stirlingS2(10, 7));
|
||||
Assert.assertEquals( 750, ArithmeticUtils.stirlingS2(10, 8));
|
||||
|
||||
}
|
||||
|
||||
@Test(expected=NotPositiveException.class)
|
||||
public void testStirlingS2NegativeN() {
|
||||
ArithmeticUtils.stirlingS2(3, -1);
|
||||
}
|
||||
|
||||
@Test(expected=NumberIsTooLargeException.class)
|
||||
public void testStirlingS2LargeK() {
|
||||
ArithmeticUtils.stirlingS2(3, 4);
|
||||
}
|
||||
|
||||
@Test(expected=MathArithmeticException.class)
|
||||
public void testStirlingS2Overflow() {
|
||||
ArithmeticUtils.stirlingS2(26, 9);
|
||||
}
|
||||
|
||||
/**
|
||||
* Exact (caching) recursive implementation to test against
|
||||
*/
|
||||
|
|
|
|||
Loading…
Reference in New Issue