commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From l..@apache.org
Subject svn commit: r1371680 - in /commons/proper/math/trunk/src: main/java/org/apache/commons/math3/util/ArithmeticUtils.java site/xdoc/userguide/index.xml site/xdoc/userguide/utilities.xml test/java/org/apache/commons/math3/util/ArithmeticUtilsTest.java
Date Fri, 10 Aug 2012 12:15:24 GMT
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 &gt; 0, 0 if x = 0, and -1 if x &lt; 0) and 
-    indicator (+1.0 if argument  &gt;= 0 and -1.0 if argument &lt; 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
      */



Mime
View raw message