diff --git a/src/main/java/org/apache/commons/math3/util/ArithmeticUtils.java b/src/main/java/org/apache/commons/math3/util/ArithmeticUtils.java index 64d35a78e..2bdb5c41b 100644 --- a/src/main/java/org/apache/commons/math3/util/ArithmeticUtils.java +++ b/src/main/java/org/apache/commons/math3/util/ArithmeticUtils.java @@ -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 STIRLING_S2 = new AtomicReference (null); + /** Private constructor. */ private ArithmeticUtils() { super(); @@ -882,6 +887,91 @@ public final class ArithmeticUtils { return result; } + /** + * Returns the + * Stirling number of the second kind, "{@code S(n,k)}", the number of + * ways of partitioning an {@code n}-element set into {@code k} non-empty + * subsets. + *

+ * The preconditions are {@code 0 <= k <= n } (otherwise + * {@code NotPositiveException} is thrown) + *

+ * @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. * diff --git a/src/site/xdoc/userguide/index.xml b/src/site/xdoc/userguide/index.xml index 5cd2771bd..46a565149 100644 --- a/src/site/xdoc/userguide/index.xml +++ b/src/site/xdoc/userguide/index.xml @@ -89,7 +89,7 @@
  • 6.2 Double array utilities
  • 6.3 int/double hash map
  • 6.4 Continued Fractions
    • -
  • 6.5 Binomial coefficients, factorials and other common math functions
    • +
  • 6.5 Binomial coefficients, factorials, Stirling numbers and other common math functions
  • 6.6 Fast mathematical functions
  • 6.7 Miscellaneous
    • diff --git a/src/site/xdoc/userguide/utilities.xml b/src/site/xdoc/userguide/utilities.xml index b1a7bc0b9..7d19bd33d 100644 --- a/src/site/xdoc/userguide/utilities.xml +++ b/src/site/xdoc/userguide/utilities.xml @@ -146,11 +146,11 @@

    - +

    A collection of reusable math functions is provided in the - MathUtils - utility class. MathUtils currently includes methods to compute the following:

      + ArithmeticUtils + utility class. ArithmeticUtils currently includes methods to compute the following:
      • Binomial coefficients -- "n choose k" available as an (exact) long value, binomialCoefficient(int, int) for small n, k; as a double, @@ -158,24 +158,13 @@ a "super-sized" version, binomialCoefficientLog(int, int) that returns the natural logarithm of the value.
      • + Stirling numbers of the second kind -- S(n,k) as an exact long value + stirlingS2(int, int) for small n, k.
        • +
      • Factorials -- like binomial coefficients, these are available as exact long values, factorial(int); doubles, factorialDouble(int); or logs, factorialLog(int).
      • - Hyperbolic sine and cosine functions -- - cosh(double), sinh(double)
        • -
      • - 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.
        • -
      • - a hash function, hash(double), returning a long-valued - hash code for a double value. -
        • -
      • - Convience methods to round floating-point number to arbitrary precision. -
        • -
      • Least common multiple and greatest common denominator functions.
      @@ -226,6 +215,7 @@
    • asinh(double)
    • acosh(double)
    • atanh(double)
      • +
    • pow(double,int)
    The following methods are found in Math/StrictMath since 1.6 only, they are provided by FastMath even in 1.5 Java virtual machines diff --git a/src/test/java/org/apache/commons/math3/util/ArithmeticUtilsTest.java b/src/test/java/org/apache/commons/math3/util/ArithmeticUtilsTest.java index e609e3fd6..969000318 100644 --- a/src/test/java/org/apache/commons/math3/util/ArithmeticUtilsTest.java +++ b/src/test/java/org/apache/commons/math3/util/ArithmeticUtilsTest.java @@ -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 */