1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.apache.commons.numbers.combinatorics;
18
19 import java.math.BigInteger;
20 import org.apache.commons.numbers.core.Precision;
21 import org.junit.jupiter.api.Assertions;
22 import org.junit.jupiter.api.Test;
23 import org.junit.jupiter.params.ParameterizedTest;
24 import org.junit.jupiter.params.provider.CsvSource;
25
26
27
28
29 class BinomialCoefficientDoubleTest {
30 @ParameterizedTest
31 @CsvSource({
32 "4, 5",
33 "-1, 1",
34 "10, -1",
35 "-1, -1",
36 "-1, -2",
37 })
38 void testBinomialCoefficientIllegalArguments(int n, int k) {
39 Assertions.assertThrows(CombinatoricsException.class, () -> BinomialCoefficientDouble.value(n, k),
40 () -> n + " choose " + k);
41 }
42
43 @ParameterizedTest
44 @CsvSource({
45
46
47 "0, 0, 1, 0",
48 "5, 0, 1, 0",
49 "5, 1, 5, 0",
50 "5, 2, 10, 0",
51 "6, 0, 1, 0",
52 "6, 1, 6, 0",
53 "6, 2, 15, 0",
54 "6, 3, 20, 0",
55 "34, 17, 2333606220, 0",
56 "66, 33, 7219428434016265740, 0",
57 "100, 10, 17310309456440, 0",
58 "1500, 4, 210094780875, 0",
59 "300, 3, 4455100, 0",
60 "700, 697, 56921900, 0",
61 "10000, 3, 166616670000, 0",
62 "412, 9, 863177604710689620, 0",
63 "678, 7, 12667255449994080, 0",
64 "66, 33, 7219428434016265740, 0",
65
66 "67, 30, 9989690752182277136, 1",
67 "67, 33, 14226520737620288370, 0",
68 "68, 34, 28453041475240576740, 0",
69
70
71 "1040, 450, 2.3101613255412135615e307, 11",
72 "1029, 514, 1.4298206864989040819e308, 5",
73 "1786388282, 38, 7.187239013254065384599502085053593e306, 0",
74 "1914878305, 38, 100.6570419073661447979173868523364e306, 1",
75 "1179067476, 39, 30.22890249420109200962786203300876e306, 2",
76 "2147483647, 37, 1.388890512412231479281222156415993e302, 4",
77 "20000, 116, 1.75293130532995289393810309132e308, 8",
78 "20000, 117, 2.97908427992998148231326853571e310, 0",
79 "1028, 514, 7.156051054877897008430135897e307, 8",
80 "1030, 496, 1.41941785031194251722295917039e308, 0",
81 "1030, 497, 1.52508879691464246317315935007e308, 32",
82 "1030, 498, 1.63227375252109323869737737668e308, 0",
83 "1030, 499, 1.74021971210665651901203359598e308, 12",
84 "1030, 500, 1.84811333425726922319077967894e308, 0",
85 "1020, 510, 2.80626776829962271039414307883e305, 8",
86 "1022, 511, 1.12140876377061244121816833013e306, 14",
87 "1024, 512, 4.48125455209897081002416485048e306, 3",
88 })
89 void testBinomialCoefficient(int n, int k, double nCk, int ulp) {
90 assertBinomial(n, k, nCk, ulp);
91 }
92
93 @ParameterizedTest
94 @CsvSource({
95 "1030, 515",
96 "10000000, 10000",
97 })
98 void testBinomialCoefficientOverflow(int n, int k) {
99 Assertions.assertEquals(Double.POSITIVE_INFINITY, BinomialCoefficientDouble.value(n, k),
100 () -> n + " choose " + k);
101 }
102
103
104
105
106
107 @Test
108 void testBinomialCoefficientLarge() {
109
110 final int size = 200;
111 final BigInteger[] factorials = new BigInteger[size + 1];
112 factorials[0] = BigInteger.ONE;
113 for (int n = 1; n <= size; n++) {
114 factorials[n] = factorials[n - 1].multiply(BigInteger.valueOf(n));
115 }
116
117 for (int n = 0; n <= size; n++) {
118 int ulp;
119 if (n <= 66) {
120 ulp = 0;
121 } else if (n <= 100) {
122 ulp = 5;
123 } else if (n <= 150) {
124 ulp = 10;
125 } else {
126 ulp = 15;
127 }
128 for (int k = 0; k <= n / 2; k++) {
129 final BigInteger nCk = factorials[n].divide(factorials[n - k]).divide(factorials[k]);
130 final double expected = nCk.doubleValue();
131 assertBinomial(n, k, expected, ulp);
132 }
133 }
134 }
135
136 private static void assertBinomial(int n, int k, double expected, int ulp) {
137 final double actual = BinomialCoefficientDouble.value(n, k);
138 if (expected == Double.POSITIVE_INFINITY) {
139 Assertions.assertEquals(expected, actual, () -> n + " choose " + k);
140 } else {
141 Assertions.assertTrue(Precision.equals(expected, actual, ulp),
142 () -> String.format("C(%d, %d) = %s : actual %s : ULP error = %d", n, k,
143 expected, actual,
144 Double.doubleToRawLongBits(actual) - Double.doubleToRawLongBits(expected)));
145 }
146
147 Assertions.assertEquals(actual, BinomialCoefficientDouble.value(n, n - k), () -> n + " choose " + k);
148 }
149 }