View Javadoc
1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
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   * Test cases for the {@link BinomialCoefficientDouble} class.
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          // Data verified using maxima: bfloat(binomial(n, k)) using 30 digits of precision.
46          // Note: This test will correctly assert infinite expected values.
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          // Overflow as a long
66          "67, 30, 9989690752182277136, 1",
67          "67, 33, 14226520737620288370, 0",
68          "68, 34, 28453041475240576740, 0",
69          // Overflow without special handling
70          // See NUMBERS-183
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      * Tests correctness for large n and sharpness of upper bound in API doc.
105      * JIRA: MATH-241
106      */
107     @Test
108     void testBinomialCoefficientLarge() {
109         // This tests all values for n <= 200.
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         // Test symmetry
147         Assertions.assertEquals(actual, BinomialCoefficientDouble.value(n, n - k), () -> n + " choose " + k);
148     }
149 }