1 // Copyright John Maddock 2006.
2 // Use, modification and distribution are subject to the
3 // Boost Software License, Version 1.0. (See accompanying file
4 // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
9 # pragma warning(disable: 4127) // conditional expression is constant.
10 # pragma warning(disable: 4245) // int/unsigned int conversion
13 // Return infinities not exceptions:
14 #define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error
16 #include <boost/math/concepts/real_concept.hpp>
17 #define BOOST_TEST_MAIN
18 #include <boost/test/unit_test.hpp>
19 #include <boost/test/floating_point_comparison.hpp>
20 #include <boost/math/special_functions/factorials.hpp>
21 #include <boost/math/special_functions/gamma.hpp>
22 #include <boost/math/tools/stats.hpp>
23 #include <boost/math/tools/test.hpp>
31 T naive_falling_factorial(T x, unsigned n)
48 // Basic sanity checks.
50 T tolerance = boost::math::tools::epsilon<T>() * 100 * 2; // 2 eps as a percent.
52 ::boost::math::factorial<T>(0),
53 static_cast<T>(1), tolerance);
55 ::boost::math::factorial<T>(1),
56 static_cast<T>(1), tolerance);
58 ::boost::math::factorial<T>(10),
59 static_cast<T>(3628800L), tolerance);
61 ::boost::math::unchecked_factorial<T>(0),
62 static_cast<T>(1), tolerance);
64 ::boost::math::unchecked_factorial<T>(1),
65 static_cast<T>(1), tolerance);
67 ::boost::math::unchecked_factorial<T>(10),
68 static_cast<T>(3628800L), tolerance);
71 // Try some double factorials:
74 ::boost::math::double_factorial<T>(0),
75 static_cast<T>(1), tolerance);
77 ::boost::math::double_factorial<T>(1),
78 static_cast<T>(1), tolerance);
80 ::boost::math::double_factorial<T>(2),
81 static_cast<T>(2), tolerance);
83 ::boost::math::double_factorial<T>(5),
84 static_cast<T>(15), tolerance);
86 ::boost::math::double_factorial<T>(10),
87 static_cast<T>(3840), tolerance);
89 ::boost::math::double_factorial<T>(19),
90 static_cast<T>(6.547290750e8L), tolerance);
92 ::boost::math::double_factorial<T>(24),
93 static_cast<T>(1.961990553600000e12L), tolerance);
95 ::boost::math::double_factorial<T>(33),
96 static_cast<T>(6.33265987076285062500000e18L), tolerance);
98 ::boost::math::double_factorial<T>(42),
99 static_cast<T>(1.0714547155728479551488000000e26L), tolerance);
101 ::boost::math::double_factorial<T>(47),
102 static_cast<T>(1.19256819277443412353990764062500000e30L), tolerance);
104 if((std::numeric_limits<T>::has_infinity) && (std::numeric_limits<T>::max_exponent <= 1024))
107 ::boost::math::double_factorial<T>(320),
108 std::numeric_limits<T>::infinity());
110 ::boost::math::double_factorial<T>(301),
111 std::numeric_limits<T>::infinity());
114 // Rising factorials:
116 tolerance = boost::math::tools::epsilon<T>() * 100 * 20; // 20 eps as a percent.
117 if(std::numeric_limits<T>::is_specialized == 0)
118 tolerance *= 5; // higher error rates without Lanczos support
120 ::boost::math::rising_factorial(static_cast<T>(3), 4),
121 static_cast<T>(360), tolerance);
123 ::boost::math::rising_factorial(static_cast<T>(7), -4),
124 static_cast<T>(0.00277777777777777777777777777777777777777777777777777777777778L), tolerance);
126 ::boost::math::rising_factorial(static_cast<T>(120.5f), 8),
127 static_cast<T>(5.58187566784927180664062500e16L), tolerance);
129 ::boost::math::rising_factorial(static_cast<T>(120.5f), -4),
130 static_cast<T>(5.15881498170104646868208445266116850161120996179812063177241e-9L), tolerance);
132 ::boost::math::rising_factorial(static_cast<T>(5000.25f), 8),
133 static_cast<T>(3.92974581976666067544013393509103775024414062500000e29L), tolerance);
135 ::boost::math::rising_factorial(static_cast<T>(5000.25f), -7),
136 static_cast<T>(1.28674092710208810281923019294164707555099052561945725535047e-26L), tolerance);
138 ::boost::math::rising_factorial(static_cast<T>(30.25), 21),
139 static_cast<T>(3.93286957998925490693364184100209193343633629069699964020401e33L), tolerance * 2);
141 ::boost::math::rising_factorial(static_cast<T>(30.25), -21),
142 static_cast<T>(3.35010902064291983728782493133164809108646650368560147505884e-27L), tolerance);
144 ::boost::math::rising_factorial(static_cast<T>(-30.25), 21),
145 static_cast<T>(-9.76168312768123676601980433377916854311706629232503473758698e26L), tolerance * 2);
147 ::boost::math::rising_factorial(static_cast<T>(-30.25), -21),
148 static_cast<T>(-1.50079704000923674318934280259377728203516775215430875839823e-34L), 2 * tolerance);
150 ::boost::math::rising_factorial(static_cast<T>(-30.25), 5),
151 static_cast<T>(-1.78799177197265625000000e7L), tolerance);
153 ::boost::math::rising_factorial(static_cast<T>(-30.25), -5),
154 static_cast<T>(-2.47177487004482195012362027432181137141899692171397467859150e-8L), tolerance);
156 ::boost::math::rising_factorial(static_cast<T>(-30.25), 6),
157 static_cast<T>(4.5146792242309570312500000e8L), tolerance);
159 ::boost::math::rising_factorial(static_cast<T>(-30.25), -6),
160 static_cast<T>(6.81868929667537089689274558433603136943171564610751635473516e-10L), tolerance);
162 ::boost::math::rising_factorial(static_cast<T>(-3), 6),
163 static_cast<T>(0), tolerance);
165 ::boost::math::rising_factorial(static_cast<T>(-3.25), 6),
166 static_cast<T>(2.99926757812500L), tolerance);
168 ::boost::math::rising_factorial(static_cast<T>(-5.25), 6),
169 static_cast<T>(50.987548828125000000000000L), tolerance);
171 ::boost::math::rising_factorial(static_cast<T>(-5.25), 13),
172 static_cast<T>(127230.91046623885631561279296875000L), tolerance);
174 ::boost::math::rising_factorial(static_cast<T>(-3.25), -6),
175 static_cast<T>(0.0000129609865918182348202632178291407500332449622510474437452125L), tolerance);
177 ::boost::math::rising_factorial(static_cast<T>(-5.25), -6),
178 static_cast<T>(2.50789821857946332294524052303699065683926911849535903362649e-6L), tolerance);
180 ::boost::math::rising_factorial(static_cast<T>(-5.25), -13),
181 static_cast<T>(-1.38984989447269128946284683518361786049649013886981662962096e-14L), tolerance);
183 // More cases reported as bugs by Rocco Romeo:
185 BOOST_CHECK_EQUAL(::boost::math::rising_factorial(static_cast<T>(0), 1), static_cast<T>(0));
186 BOOST_CHECK_EQUAL(::boost::math::rising_factorial(static_cast<T>(0), -1), static_cast<T>(-1));
187 BOOST_CHECK_CLOSE(::boost::math::rising_factorial(static_cast<T>(0.5f), -1), static_cast<T>(-2), tolerance);
188 BOOST_CHECK_CLOSE(::boost::math::rising_factorial(static_cast<T>(40.5), -41), static_cast<T>(-2.75643016796662963097096639854040835565778207128865739e-47L), tolerance);
189 BOOST_CHECK_EQUAL(::boost::math::rising_factorial(static_cast<T>(-2), 3), static_cast<T>(0));
190 BOOST_CHECK_EQUAL(::boost::math::rising_factorial(static_cast<T>(-2), 2), static_cast<T>(2));
191 BOOST_CHECK_EQUAL(::boost::math::rising_factorial(static_cast<T>(-4), 3), static_cast<T>(-24));
192 BOOST_CHECK_CLOSE(::boost::math::rising_factorial(static_cast<T>(-4), -3), static_cast<T>(-0.00476190476190476190476190476190476190476190476190476190476L), tolerance);
193 if(ldexp(T(1), -150) != 0)
195 BOOST_CHECK_CLOSE(::boost::math::rising_factorial(ldexp(T(1), -150), 0), static_cast<T>(1), tolerance);
196 BOOST_CHECK_CLOSE(::boost::math::rising_factorial(ldexp(T(1), -150), -1), static_cast<T>(-1.00000000000000000000000000000000000000000000070064923216241L), tolerance);
197 BOOST_CHECK_CLOSE(::boost::math::rising_factorial(ldexp(T(1), -150), -2), static_cast<T>(0.500000000000000000000000000000000000000000000525486924121806L), tolerance);
198 BOOST_CHECK_CLOSE(::boost::math::rising_factorial(ldexp(T(1), -150), -25), static_cast<T>(-6.44695028438447339619485321920468720529890442766578870603544e-26L), 15 * tolerance);
199 if(std::numeric_limits<T>::min_exponent10 < -50)
201 BOOST_CHECK_CLOSE(::boost::math::rising_factorial(ldexp(T(1), -150), -40), static_cast<T>(1.22561743912838584942353998493975692372556196815242899727910e-48L), tolerance);
206 // Falling factorials:
209 ::boost::math::falling_factorial(static_cast<T>(30.25), 0),
210 static_cast<T>(naive_falling_factorial(30.25L, 0)),
213 ::boost::math::falling_factorial(static_cast<T>(30.25), 1),
214 static_cast<T>(naive_falling_factorial(30.25L, 1)),
217 ::boost::math::falling_factorial(static_cast<T>(30.25), 2),
218 static_cast<T>(naive_falling_factorial(30.25L, 2)),
221 ::boost::math::falling_factorial(static_cast<T>(30.25), 5),
222 static_cast<T>(naive_falling_factorial(30.25L, 5)),
225 ::boost::math::falling_factorial(static_cast<T>(30.25), 22),
226 static_cast<T>(naive_falling_factorial(30.25L, 22)),
229 ::boost::math::falling_factorial(static_cast<T>(100.5), 6),
230 static_cast<T>(naive_falling_factorial(100.5L, 6)),
233 ::boost::math::falling_factorial(static_cast<T>(30.75), 30),
234 static_cast<T>(naive_falling_factorial(30.75L, 30)),
236 if(boost::math::policies::digits<T, boost::math::policies::policy<> >() > 50)
239 ::boost::math::falling_factorial(static_cast<T>(-30.75L), 30),
240 static_cast<T>(naive_falling_factorial(-30.75L, 30)),
243 ::boost::math::falling_factorial(static_cast<T>(-30.75L), 27),
244 static_cast<T>(naive_falling_factorial(-30.75L, 27)),
248 ::boost::math::falling_factorial(static_cast<T>(-12.0), 6),
249 static_cast<T>(naive_falling_factorial(-12.0L, 6)),
252 ::boost::math::falling_factorial(static_cast<T>(-12), 5),
253 static_cast<T>(naive_falling_factorial(-12.0L, 5)),
256 ::boost::math::falling_factorial(static_cast<T>(-3.0), 6),
257 static_cast<T>(naive_falling_factorial(-3.0L, 6)),
260 ::boost::math::falling_factorial(static_cast<T>(-3), 5),
261 static_cast<T>(naive_falling_factorial(-3.0L, 5)),
264 ::boost::math::falling_factorial(static_cast<T>(3.0), 6),
265 static_cast<T>(naive_falling_factorial(3.0L, 6)),
268 ::boost::math::falling_factorial(static_cast<T>(3), 5),
269 static_cast<T>(naive_falling_factorial(3.0L, 5)),
272 ::boost::math::falling_factorial(static_cast<T>(3.25), 4),
273 static_cast<T>(naive_falling_factorial(3.25L, 4)),
276 ::boost::math::falling_factorial(static_cast<T>(3.25), 5),
277 static_cast<T>(naive_falling_factorial(3.25L, 5)),
280 ::boost::math::falling_factorial(static_cast<T>(3.25), 6),
281 static_cast<T>(naive_falling_factorial(3.25L, 6)),
284 ::boost::math::falling_factorial(static_cast<T>(3.25), 7),
285 static_cast<T>(naive_falling_factorial(3.25L, 7)),
288 ::boost::math::falling_factorial(static_cast<T>(8.25), 12),
289 static_cast<T>(naive_falling_factorial(8.25L, 12)),
292 // More cases reported as bugs by Rocco Romeo:
294 BOOST_CHECK_EQUAL(::boost::math::falling_factorial(static_cast<T>(0), 1), static_cast<T>(0));
295 BOOST_CHECK_CLOSE(::boost::math::falling_factorial(static_cast<T>(-2), 3), static_cast<T>(-24), tolerance);
296 BOOST_CHECK_CLOSE(::boost::math::falling_factorial(static_cast<T>(-2), 2), static_cast<T>(6), tolerance);
297 BOOST_CHECK_CLOSE(::boost::math::falling_factorial(static_cast<T>(-4), 3), static_cast<T>(-120), tolerance);
298 if(ldexp(static_cast<T>(1), -100) != 0)
300 BOOST_CHECK_CLOSE(::boost::math::falling_factorial(ldexp(static_cast<T>(1), -100), 1), static_cast<T>(7.888609052210118054117285652827862296732064351090230047e-31L), tolerance);
301 BOOST_CHECK_CLOSE(::boost::math::falling_factorial(ldexp(static_cast<T>(1), -100), 2), static_cast<T>(-7.88860905221011805411728565282163928145420320938308598e-31L), tolerance);
302 BOOST_CHECK_CLOSE(::boost::math::falling_factorial(ldexp(static_cast<T>(1), -100), 3), static_cast<T>(1.577721810442023610823457130563705554763054527705902790e-30L), tolerance);
303 BOOST_CHECK_CLOSE(::boost::math::falling_factorial(ldexp(static_cast<T>(1), -100), 35), static_cast<T>(2.32897613101315187323300837924676081371219290005311315885e8L), 35 * tolerance);
305 if((ldexp(static_cast<T>(1), -300) != 0) && (std::numeric_limits<T>::max_exponent10 > 290))
307 BOOST_CHECK_CLOSE(::boost::math::falling_factorial(ldexp(static_cast<T>(1), -300), 20), static_cast<T>(-5.97167167502482975928590631196751639118233432208390100e-74L), tolerance);
308 BOOST_CHECK_CLOSE(::boost::math::falling_factorial(ldexp(static_cast<T>(1), -300), 200), static_cast<T>(-1.93579759151806711025267355739174942986011285920860098569075e282L), 10 * tolerance);
312 tolerance = boost::math::tools::epsilon<T>() * 100 * 20; // 20 eps as a percent.
313 unsigned i = boost::math::max_factorial<T>::value;
314 if((boost::is_floating_point<T>::value) && (sizeof(T) <= sizeof(double)))
316 // Without Lanczos support, tgamma isn't accurate enough for this test:
318 ::boost::math::unchecked_factorial<T>(i),
319 boost::math::tgamma(static_cast<T>(i+1)), tolerance);
323 while(boost::math::lgamma(static_cast<T>(i+1)) < boost::math::tools::log_max_value<T>())
326 ::boost::math::factorial<T>(i),
327 boost::math::tgamma(static_cast<T>(i+1)), tolerance);
330 } // template <class T> void test_spots(T)
332 BOOST_AUTO_TEST_CASE( test_main )
334 BOOST_MATH_CONTROL_FP;
337 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
339 #ifndef BOOST_MATH_NO_REAL_CONCEPT_TESTS
340 test_spots(boost::math::concepts::real_concept(0.));
343 std::cout << "<note>The long double tests have been disabled on this platform "
344 "either because the long double overloads of the usual math functions are "
345 "not available at all, or because they are too inaccurate for these tests "
346 "to pass.</note>" << std::endl;
348 if (std::numeric_limits<double>::digits == std::numeric_limits<long double>::digits)
350 cout << "Types double and long double have the same number of floating-point significand bits ("
351 << std::numeric_limits<long double>::digits << ") on this platform." << endl;
353 if (std::numeric_limits<float>::digits == std::numeric_limits<double>::digits)
355 cout << "Types float and double have the same number of floating-point significand bits ("
356 << std::numeric_limits<double>::digits << ") on this platform." << endl;
359 using boost::math::max_factorial;
360 cout << "max factorial for float " << max_factorial<float>::value << endl;
361 cout << "max factorial for double " << max_factorial<double>::value << endl;
362 cout << "max factorial for long double " << max_factorial<long double>::value << endl;
363 cout << "max factorial for real_concept " << max_factorial<boost::math::concepts::real_concept>::value << endl;
374 Running 1 test case...
375 Types double and long double have the same number of floating-point significand bits (53) on this platform.
376 max factorial for float 34
377 max factorial for double 170
378 max factorial for long double 170
379 max factorial for real_concept 100
380 *** No errors detected