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/10 14:18:01 UTC
[commons-numbers] branch master updated: NUMBERS-29: Update stirlingS2 for (n, n-2) and (n, n-3) cases
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
The following commit(s) were added to refs/heads/master by this push:
new 50bf79b2 NUMBERS-29: Update stirlingS2 for (n, n-2) and (n, n-3) cases
50bf79b2 is described below
commit 50bf79b20e84fb5946ad66ec0255266448d4b235
Author: aherbert <ah...@apache.org>
AuthorDate: Thu Nov 10 13:51:04 2022 +0000
NUMBERS-29: Update stirlingS2 for (n, n-2) and (n, n-3) cases
For S(n, n-2) the summation of binomial coefficients is replaced with a
simple identity.
A product of binomial coefficients can be used for S(n, n-3).
---
.../commons/numbers/combinatorics/Stirling.java | 66 +++++++++++++---------
.../numbers/combinatorics/StirlingTest.java | 22 ++++++++
2 files changed, 60 insertions(+), 28 deletions(-)
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
index 2d301eae..eb63797b 100644
--- 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
@@ -35,6 +35,8 @@ public final class Stirling {
private static final int S1_OVERFLOW_K_EQUALS_NM3 = 2761;
/** Overflow threshold for n when computing S(n, n-2). */
private static final int S2_OVERFLOW_K_EQUALS_NM2 = 92683;
+ /** Overflow threshold for n when computing S(n, n-3). */
+ private static final int S2_OVERFLOW_K_EQUALS_NM3 = 2762;
/**
* Precomputed Stirling numbers of the first kind.
@@ -133,14 +135,7 @@ public final class Stirling {
} else if (k == n - 2) {
checkN(n, k, S1_OVERFLOW_K_EQUALS_NM2, S1_ERROR_FORMAT);
// (3n-1) * binom(n, 3) / 4
- final long a = 3L * n - 1;
- final long b = BinomialCoefficient.value(n, 3);
- // Compute (a*b/4) without intermediate overflow.
- // The product (a*b) must be an exact multiple of 4.
- // Conditional branch on b which is typically large and even (a is 50% even)
- // If b is even: ((b/2) * a) / 2
- // If b is odd then a must be even to make a*b even: ((a/2) * b) / 2
- return (b & 1) == 0 ? ((b >>> 1) * a) >>> 1 : ((a >>> 1) * b) >>> 1;
+ return productOver4(3L * n - 1, BinomialCoefficient.value(n, 3));
} else if (k == n - 3) {
checkN(n, k, S1_OVERFLOW_K_EQUALS_NM3, S1_ERROR_FORMAT);
return -BinomialCoefficient.value(n, 2) * BinomialCoefficient.value(n, 4);
@@ -203,29 +198,18 @@ public final class Stirling {
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
+ } else if (k == n - 2) {
checkN(n, k, S2_OVERFLOW_K_EQUALS_NM2, S2_ERROR_FORMAT);
- 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;
+ // (3n-5) * binom(n, 3) / 4
+ return productOver4(3L * n - 5, BinomialCoefficient.value(n, 3));
+ } else if (k == n - 3) {
+ checkN(n, k, S2_OVERFLOW_K_EQUALS_NM3, S2_ERROR_FORMAT);
+ return BinomialCoefficient.value(n - 2, 2) * BinomialCoefficient.value(n, 4);
}
+ // Compute using:
+ // S(n, k) = k * S(n - 1, k) + S(n - 1, k - 1)
+
// n >= 26 (MAX_N)
// 3 <= k <= n-3
@@ -282,4 +266,30 @@ public final class Stirling {
throw new ArithmeticException(String.format(msgFormat, n, k));
}
}
+
+ /**
+ * Return {@code a*b/4} without intermediate overflow.
+ * It is assumed that:
+ * <ul>
+ * <li>The coefficients a and b are positive
+ * <li>The product (a*b) is an exact multiple of 4
+ * <li>The result (a*b/4) is an exact integer that does not overflow a {@code long}
+ * </ul>
+ *
+ * <p>A conditional branch is performed on the odd/even property of {@code b}.
+ * The branch is predictable if {@code b} is typically the same parity.
+ *
+ * @param a Coefficient a
+ * @param b Coefficient b
+ * @return {@code a*b/4}
+ */
+ private static long productOver4(long a, long b) {
+ // Compute (a*b/4) without intermediate overflow.
+ // The product (a*b) must be an exact multiple of 4.
+ // If b is even: ((b/2) * a) / 2
+ // If b is odd then a must be even to make a*b even: ((a/2) * b) / 2
+ return (b & 1) == 0 ?
+ ((b >>> 1) * a) >>> 1 :
+ ((a >>> 1) * b) >>> 1;
+ }
}
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
index 26587619..67af0ca2 100644
--- 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
@@ -286,6 +286,28 @@ class StirlingTest {
// Limits for k in [n-1, n] use n = Integer.MAX_VALUE
"2147483647, 2147483646, 2305843005992468481",
"2147483647, 2147483647, 1",
+ // Data for S(n, n-2)
+ "26, 24, 47450",
+ "27, 25, 55575",
+ "28, 26, 64701",
+ "29, 27, 74907",
+ "30, 28, 86275",
+ "31, 29, 98890",
+ "32, 30, 112840",
+ "92680, 92678, 9222151351858080650",
+ "92681, 92679, 9222549384516960590",
+ "92682, 92680, 9222947430060167790",
+ // Data for S(n, n-3)
+ "26, 23, 4126200",
+ "27, 24, 5265000",
+ "28, 25, 6654375",
+ "29, 26, 8336601",
+ "30, 27, 10359090",
+ "31, 28, 12774790",
+ "32, 29, 15642600",
+ "2759, 2756, 9152435640507623646",
+ "2760, 2757, 9172370757033509130",
+ "2761, 2758, 9192342044684582630",
})
void testStirlingS2(int n, int k, long expected) {
Assertions.assertEquals(expected, Stirling.stirlingS2(n, k));