You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by ah...@apache.org on 2022/11/02 16:52:19 UTC
[commons-numbers] 01/03: NUMBERS-29: Add computation of Stirling number of the second kind
This is an automated email from the ASF dual-hosted git repository.
aherbert pushed a commit to branch master
in repository https://gitbox.apache.org/repos/asf/commons-numbers.git
commit 916d8ad09f1c9da2cbe00798281c0d179f33f8e9
Author: aherbert <ah...@apache.org>
AuthorDate: Tue Nov 1 16:08:11 2022 +0000
NUMBERS-29: Add computation of Stirling number of the second kind
---
.../commons/numbers/combinatorics/Stirling.java | 157 +++++++++++++++++++++
.../numbers/combinatorics/StirlingTest.java | 149 +++++++++++++++++++
2 files changed, 306 insertions(+)
diff --git a/commons-numbers-combinatorics/src/main/java/org/apache/commons/numbers/combinatorics/Stirling.java b/commons-numbers-combinatorics/src/main/java/org/apache/commons/numbers/combinatorics/Stirling.java
new file mode 100644
index 00000000..d5250fc4
--- /dev/null
+++ b/commons-numbers-combinatorics/src/main/java/org/apache/commons/numbers/combinatorics/Stirling.java
@@ -0,0 +1,157 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+package org.apache.commons.numbers.combinatorics;
+
+/**
+ * Computation of <a href="https://en.wikipedia.org/wiki/Stirling_number">Stirling numbers</a>.
+ *
+ * @since 1.2
+ */
+public final class Stirling {
+ /** Stirling S2 error message. */
+ private static final String S2_ERROR_FORMAT = "S(n=%d, k=%d)";
+ /** Overflow threshold for n when computing S(n, n-2). */
+ private static final int S2_OVERFLOW_K_EQUALS_NM2 = 92683;
+
+ /**
+ * Precomputed Stirling numbers of the second kind.
+ * Provides a thread-safe lazy initialization of the cache.
+ */
+ private static class StirlingS2Cache {
+ /** Maximum n to compute (exclusive).
+ * As S(26,9) = 11201516780955125625 is larger than Long.MAX_VALUE
+ * we must stop computation at row 26. */
+ static final int MAX_N = 26;
+ /** Stirling numbers of the second kind. */
+ static final long[][] STIRLING_S2;
+
+ static {
+ STIRLING_S2 = new long[MAX_N][];
+ STIRLING_S2[0] = new long[] {1};
+ for (int n = 1; n < STIRLING_S2.length; n++) {
+ STIRLING_S2[n] = new long[n + 1];
+ STIRLING_S2[n][0] = 0;
+ STIRLING_S2[n][1] = 1;
+ STIRLING_S2[n][n] = 1;
+ for (int k = 2; k < n; k++) {
+ STIRLING_S2[n][k] = k * STIRLING_S2[n - 1][k] + STIRLING_S2[n - 1][k - 1];
+ }
+ }
+ }
+ }
+
+ /** Private constructor. */
+ private Stirling() {
+ // intentionally empty.
+ }
+
+ /**
+ * Returns the <a
+ * href="https://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.
+ *
+ * @param n Size of the set
+ * @param k Number of non-empty subsets ({@code 0 <= k <= n})
+ * @return {@code S(n,k)}
+ * @throws IllegalArgumentException if {@code k < 0} or {@code k > n}.
+ * @throws ArithmeticException 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(int n, int k) {
+ if (k < 0) {
+ throw new CombinatoricsException(CombinatoricsException.NEGATIVE, k);
+ }
+ if (k > n) {
+ throw new CombinatoricsException(CombinatoricsException.OUT_OF_RANGE, k, 0, n);
+ }
+
+ if (n < StirlingS2Cache.MAX_N) {
+ // The number is in the small cache
+ return StirlingS2Cache.STIRLING_S2[n][k];
+ }
+
+ // Simple cases
+ if (k == 0) {
+ return 0;
+ } else if (k == 1 || k == n) {
+ return 1;
+ } else if (k == 2) {
+ checkN(n, k, 64);
+ return (1L << (n - 1)) - 1L;
+ } else if (k == n - 1) {
+ return BinomialCoefficient.value(n, 2);
+ }
+
+ // Compute using: S(n, k) = k * S(n - 1, k) + S(n - 1, k - 1)
+
+ if (k == n - 2) {
+ // Given:
+ // k * S(n - 1, k) == (n-2) * S(n-1, n-2)) == (n-2) * binom(n-1, 2))
+ // the recursion reduces to a sum of binomial coefficients:
+ // for i in [1, k]:
+ // sum (i * binom(i+1, 2))
+ // Avoid overflow checks using the known limit for n when k=n-2
+ checkN(n, k, S2_OVERFLOW_K_EQUALS_NM2);
+ long binom = BinomialCoefficient.value(k + 1, 2);
+ long sum = 0;
+ for (int i = k; i > 0; i--) {
+ sum += i * binom;
+ // update binomial coefficient:
+ // binom(i, 2) = binom(i+1, 2) - i
+ binom -= i;
+ }
+ return sum;
+ }
+
+ // n >= 26 (MAX_N)
+ // 3 <= k <= n-3
+
+ // Start at the largest easily computed value: n < MAX_N or k < 3
+ final int reduction = Math.min(n - StirlingS2Cache.MAX_N, k - 3) + 1;
+ int n0 = n - reduction;
+ int k0 = k - reduction;
+
+ long sum = stirlingS2(n0, k0);
+ while (n0 < n) {
+ n0++;
+ k0++;
+ sum = Math.addExact(
+ Math.multiplyExact(k0, stirlingS2(n0 - 1, k0)),
+ sum
+ );
+ }
+
+ return sum;
+ }
+
+ /**
+ * Check {@code n <= threshold}, or else throw an {@link ArithmeticException}.
+ *
+ * @param n N
+ * @param k K
+ * @param threshold Threshold for {@code n}
+ * @throws ArithmeticException if overflow is expected to happen
+ */
+ private static void checkN(int n, int k, int threshold) {
+ if (n > threshold) {
+ throw new ArithmeticException(String.format(S2_ERROR_FORMAT, n, k));
+ }
+ }
+}
diff --git a/commons-numbers-combinatorics/src/test/java/org/apache/commons/numbers/combinatorics/StirlingTest.java b/commons-numbers-combinatorics/src/test/java/org/apache/commons/numbers/combinatorics/StirlingTest.java
new file mode 100644
index 00000000..0264d200
--- /dev/null
+++ b/commons-numbers-combinatorics/src/test/java/org/apache/commons/numbers/combinatorics/StirlingTest.java
@@ -0,0 +1,149 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.commons.numbers.combinatorics;
+
+import org.junit.jupiter.api.Assertions;
+import org.junit.jupiter.api.Test;
+import org.junit.jupiter.params.ParameterizedTest;
+import org.junit.jupiter.params.provider.CsvSource;
+
+/**
+ * Test cases for the {@link Stirling} class.
+ */
+class StirlingTest {
+
+ @ParameterizedTest
+ @CsvSource({
+ "1, -1",
+ "-1, -1",
+ "-1, 1",
+ "10, 15",
+ })
+ void testStirlingS2IllegalArgument(int n, int k) {
+ Assertions.assertThrows(IllegalArgumentException.class, () -> Stirling.stirlingS2(n, k));
+ }
+
+ @Test
+ void testStirlingS2StandardCases() {
+ Assertions.assertEquals(1, Stirling.stirlingS2(0, 0));
+
+ for (int n = 1; n < 64; ++n) {
+ Assertions.assertEquals(0, Stirling.stirlingS2(n, 0));
+ Assertions.assertEquals(1, Stirling.stirlingS2(n, 1));
+ if (n > 2) {
+ Assertions.assertEquals((1L << (n - 1)) - 1L, Stirling.stirlingS2(n, 2));
+ Assertions.assertEquals(BinomialCoefficient.value(n, 2),
+ Stirling.stirlingS2(n, n - 1));
+ }
+ Assertions.assertEquals(1, Stirling.stirlingS2(n, n));
+ }
+ }
+
+ @ParameterizedTest
+ @CsvSource({
+ // Data verified using Mathematica StirlingS2[n, k]
+ "5, 3, 25",
+ "6, 3, 90",
+ "6, 4, 65",
+ "7, 3, 301",
+ "7, 4, 350",
+ "7, 5, 140",
+ "8, 3, 966",
+ "8, 4, 1701",
+ "8, 5, 1050",
+ "8, 6, 266",
+ "9, 3, 3025",
+ "9, 4, 7770",
+ "9, 5, 6951",
+ "9, 6, 2646",
+ "9, 7, 462",
+ "10, 3, 9330",
+ "10, 4, 34105",
+ "10, 5, 42525",
+ "10, 6, 22827",
+ "10, 7, 5880",
+ "10, 8, 750",
+ "30, 2, 536870911",
+ // Limits for k in [n-1, n]
+ "2147483647, 2147483646, 2305843005992468481",
+ "2147483647, 2147483647, 1",
+ })
+ void testStirlingS2(int n, int k, long expected) {
+ Assertions.assertEquals(expected, Stirling.stirlingS2(n, k));
+ }
+
+ @ParameterizedTest
+ @CsvSource({
+ // Upper limits for n with k in [2, 22]
+ "64, 2, 9223372036854775807",
+ "41, 3, 6078831630016836625",
+ "33, 4, 3073530837671316330",
+ "30, 5, 7713000216608565075",
+ "28, 6, 8220146115188676396",
+ "26, 7, 1631853797991016600",
+ "26, 8, 5749622251945664950",
+ "25, 9, 1167921451092973005",
+ "25, 10, 1203163392175387500",
+ "25, 11, 802355904438462660",
+ "26, 12, 5149507353856958820",
+ "26, 13, 1850568574253550060",
+ "27, 14, 8541149231801585700",
+ "27, 15, 1834634071262848260",
+ "28, 16, 6539643128396047620",
+ "28, 17, 898741468057510350",
+ "29, 18, 2598531274376323650",
+ "30, 19, 7145845579888333500",
+ "30, 20, 581535955088511150",
+ "31, 21, 1359760239259935240",
+ "32, 22, 3069483578649883980",
+ // Upper limits for n with with k in [n-10, n-2]
+ "33, 23, 6708404338089491700",
+ "38, 29, 6766081393022256030",
+ "47, 39, 8248929419122431611",
+ "63, 56, 8426132708490143472",
+ "96, 90, 8130394568857873780",
+ "183, 178, 9208213764702344301",
+ "496, 492, 9161200151995742994",
+ "2762, 2759, 9212349555946145400",
+ "92683, 92681, 9223345488487980291",
+ })
+ void testStirlingS2LimitsN(int n, int k, long expected) {
+ Assertions.assertEquals(expected, Stirling.stirlingS2(n, k));
+ Assertions.assertThrows(ArithmeticException.class, () -> Stirling.stirlingS2(n + 1, k));
+ Assertions.assertThrows(ArithmeticException.class, () -> Stirling.stirlingS2(n + 100, k));
+ Assertions.assertThrows(ArithmeticException.class, () -> Stirling.stirlingS2(n + 10000, k));
+ }
+
+ @ParameterizedTest
+ @CsvSource({
+ // Large numbers that should easily overflow. Verifies the exception is correct
+ // (e.g. no StackOverflowError occurs due to recursion)
+ "123, 32",
+ "612534, 56123",
+ "261388631, 213",
+ "678688997, 213879",
+ "1000000002, 1000000000",
+ "1000000003, 1000000000",
+ "1000000004, 1000000000",
+ "1000000005, 1000000000",
+ "1000000010, 1000000000",
+ "1000000100, 1000000000",
+ })
+ void testStirlingS2Overflow(int n, int k) {
+ Assertions.assertThrows(ArithmeticException.class, () -> Stirling.stirlingS2(n, k));
+ }
+}