Author: luc
Date: Fri Aug 10 12:15:23 2012
New Revision: 1371680
URL: http://svn.apache.org/viewvc?rev=1371680&view=rev
Log:
Added Stirling numbers of the second kind in ArithmeticUtils.
Modified:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/util/ArithmeticUtils.java
commons/proper/math/trunk/src/site/xdoc/userguide/index.xml
commons/proper/math/trunk/src/site/xdoc/userguide/utilities.xml
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/util/ArithmeticUtilsTest.java
Modified:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/util/ArithmeticUtils.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/util/ArithmeticUtils.java?rev=1371680&r1=1371679&r2=1371680&view=diff
==============================================================================
---
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/util/ArithmeticUtils.java
(original)
+++
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/util/ArithmeticUtils.java
Fri Aug 10 12:15:23 2012
@@ -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();
@@ -883,6 +888,91 @@ public final class ArithmeticUtils {
}
/**
+ * 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.
*
* @param a Addend.
Modified: commons/proper/math/trunk/src/site/xdoc/userguide/index.xml
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/site/xdoc/userguide/index.xml?rev=1371680&r1=1371679&r2=1371680&view=diff
==============================================================================
--- commons/proper/math/trunk/src/site/xdoc/userguide/index.xml (original)
+++ commons/proper/math/trunk/src/site/xdoc/userguide/index.xml Fri Aug 10
12:15:23 2012
@@ -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>
Modified: commons/proper/math/trunk/src/site/xdoc/userguide/utilities.xml
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/site/xdoc/userguide/utilities.xml?rev=1371680&r1=1371679&r2=1371680&view=diff
==============================================================================
--- commons/proper/math/trunk/src/site/xdoc/userguide/utilities.xml (original)
+++ commons/proper/math/trunk/src/site/xdoc/userguide/utilities.xml Fri Aug 10
12:15:23 2012
@@ -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
Modified:
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/util/ArithmeticUtilsTest.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math3/util/ArithmeticUtilsTest.java?rev=1371680&r1=1371679&r2=1371680&view=diff
==============================================================================
---
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/util/ArithmeticUtilsTest.java
(original)
+++
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/util/ArithmeticUtilsTest.java
Fri Aug 10 12:15:23 2012
@@ -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
*/
