]> git.proxmox.com Git - ceph.git/blame - ceph/src/boost/libs/math/test/test_factorials.cpp
projects / ceph.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
update ceph source to reef 18.1.2
[ceph.git] / ceph / src / boost / libs / math / test / test_factorials.cpp
CommitLineData
7c673cae
FG		 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)
5
6#include <pch.hpp>
7
8#ifdef _MSC_VER
9# pragma warning(disable: 4127) // conditional expression is constant.
10# pragma warning(disable: 4245) // int/unsigned int conversion
11#endif
12
13// Return infinities not exceptions:
14#define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error
15
16#include <boost/math/concepts/real_concept.hpp>
17#define BOOST_TEST_MAIN
18#include <boost/test/unit_test.hpp>
92f5a8d4 19#include <boost/test/tools/floating_point_comparison.hpp>
7c673cae
FG		 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>
24
25#include <iostream>
26#include <iomanip>
27using std::cout;
28using std::endl;
29
30template <class T>
31T naive_falling_factorial(T x, unsigned n)
32{
33 if(n == 0)
34 return 1;
35 T result = x;
36 while(--n)
37 {
38 x -= 1;
39 result *= x;
40 }
41 return result;
42}
43
44template <class T>
45void test_spots(T)
46{
1e59de90 47 using std::ldexp;
7c673cae
FG		 48 //
49 // Basic sanity checks.
50 //
51 T tolerance = boost::math::tools::epsilon<T>() * 100 * 2; // 2 eps as a percent.
52 BOOST_CHECK_CLOSE(
53 ::boost::math::factorial<T>(0),
54 static_cast<T>(1), tolerance);
55 BOOST_CHECK_CLOSE(
56 ::boost::math::factorial<T>(1),
57 static_cast<T>(1), tolerance);
58 BOOST_CHECK_CLOSE(
59 ::boost::math::factorial<T>(10),
60 static_cast<T>(3628800L), tolerance);
61 BOOST_CHECK_CLOSE(
62 ::boost::math::unchecked_factorial<T>(0),
63 static_cast<T>(1), tolerance);
64 BOOST_CHECK_CLOSE(
65 ::boost::math::unchecked_factorial<T>(1),
66 static_cast<T>(1), tolerance);
67 BOOST_CHECK_CLOSE(
68 ::boost::math::unchecked_factorial<T>(10),
69 static_cast<T>(3628800L), tolerance);
70
71 //
72 // Try some double factorials:
73 //
74 BOOST_CHECK_CLOSE(
75 ::boost::math::double_factorial<T>(0),
76 static_cast<T>(1), tolerance);
77 BOOST_CHECK_CLOSE(
78 ::boost::math::double_factorial<T>(1),
79 static_cast<T>(1), tolerance);
80 BOOST_CHECK_CLOSE(
81 ::boost::math::double_factorial<T>(2),
82 static_cast<T>(2), tolerance);
83 BOOST_CHECK_CLOSE(
84 ::boost::math::double_factorial<T>(5),
85 static_cast<T>(15), tolerance);
86 BOOST_CHECK_CLOSE(
87 ::boost::math::double_factorial<T>(10),
88 static_cast<T>(3840), tolerance);
89 BOOST_CHECK_CLOSE(
90 ::boost::math::double_factorial<T>(19),
91 static_cast<T>(6.547290750e8L), tolerance);
92 BOOST_CHECK_CLOSE(
93 ::boost::math::double_factorial<T>(24),
94 static_cast<T>(1.961990553600000e12L), tolerance);
95 BOOST_CHECK_CLOSE(
96 ::boost::math::double_factorial<T>(33),
97 static_cast<T>(6.33265987076285062500000e18L), tolerance);
98 BOOST_CHECK_CLOSE(
99 ::boost::math::double_factorial<T>(42),
100 static_cast<T>(1.0714547155728479551488000000e26L), tolerance);
101 BOOST_CHECK_CLOSE(
102 ::boost::math::double_factorial<T>(47),
103 static_cast<T>(1.19256819277443412353990764062500000e30L), tolerance);
104
105 if((std::numeric_limits<T>::has_infinity) && (std::numeric_limits<T>::max_exponent <= 1024))
106 {
107 BOOST_CHECK_EQUAL(
108 ::boost::math::double_factorial<T>(320),
109 std::numeric_limits<T>::infinity());
110 BOOST_CHECK_EQUAL(
111 ::boost::math::double_factorial<T>(301),
112 std::numeric_limits<T>::infinity());
113 }
114 //
115 // Rising factorials:
116 //
117 tolerance = boost::math::tools::epsilon<T>() * 100 * 20; // 20 eps as a percent.
118 if(std::numeric_limits<T>::is_specialized == 0)
119 tolerance *= 5; // higher error rates without Lanczos support
120 BOOST_CHECK_CLOSE(
121 ::boost::math::rising_factorial(static_cast<T>(3), 4),
122 static_cast<T>(360), tolerance);
123 BOOST_CHECK_CLOSE(
124 ::boost::math::rising_factorial(static_cast<T>(7), -4),
125 static_cast<T>(0.00277777777777777777777777777777777777777777777777777777777778L), tolerance);
126 BOOST_CHECK_CLOSE(
127 ::boost::math::rising_factorial(static_cast<T>(120.5f), 8),
128 static_cast<T>(5.58187566784927180664062500e16L), tolerance);
129 BOOST_CHECK_CLOSE(
130 ::boost::math::rising_factorial(static_cast<T>(120.5f), -4),
131 static_cast<T>(5.15881498170104646868208445266116850161120996179812063177241e-9L), tolerance);
132 BOOST_CHECK_CLOSE(
133 ::boost::math::rising_factorial(static_cast<T>(5000.25f), 8),
134 static_cast<T>(3.92974581976666067544013393509103775024414062500000e29L), tolerance);
135 BOOST_CHECK_CLOSE(
136 ::boost::math::rising_factorial(static_cast<T>(5000.25f), -7),
137 static_cast<T>(1.28674092710208810281923019294164707555099052561945725535047e-26L), tolerance);
138 BOOST_CHECK_CLOSE(
139 ::boost::math::rising_factorial(static_cast<T>(30.25), 21),
140 static_cast<T>(3.93286957998925490693364184100209193343633629069699964020401e33L), tolerance * 2);
141 BOOST_CHECK_CLOSE(
142 ::boost::math::rising_factorial(static_cast<T>(30.25), -21),
143 static_cast<T>(3.35010902064291983728782493133164809108646650368560147505884e-27L), tolerance);
144 BOOST_CHECK_CLOSE(
145 ::boost::math::rising_factorial(static_cast<T>(-30.25), 21),
146 static_cast<T>(-9.76168312768123676601980433377916854311706629232503473758698e26L), tolerance * 2);
147 BOOST_CHECK_CLOSE(
148 ::boost::math::rising_factorial(static_cast<T>(-30.25), -21),
149 static_cast<T>(-1.50079704000923674318934280259377728203516775215430875839823e-34L), 2 * tolerance);
150 BOOST_CHECK_CLOSE(
151 ::boost::math::rising_factorial(static_cast<T>(-30.25), 5),
152 static_cast<T>(-1.78799177197265625000000e7L), tolerance);
153 BOOST_CHECK_CLOSE(
154 ::boost::math::rising_factorial(static_cast<T>(-30.25), -5),
155 static_cast<T>(-2.47177487004482195012362027432181137141899692171397467859150e-8L), tolerance);
156 BOOST_CHECK_CLOSE(
157 ::boost::math::rising_factorial(static_cast<T>(-30.25), 6),
158 static_cast<T>(4.5146792242309570312500000e8L), tolerance);
159 BOOST_CHECK_CLOSE(
160 ::boost::math::rising_factorial(static_cast<T>(-30.25), -6),
161 static_cast<T>(6.81868929667537089689274558433603136943171564610751635473516e-10L), tolerance);
162 BOOST_CHECK_CLOSE(
163 ::boost::math::rising_factorial(static_cast<T>(-3), 6),
164 static_cast<T>(0), tolerance);
165 BOOST_CHECK_CLOSE(
166 ::boost::math::rising_factorial(static_cast<T>(-3.25), 6),
167 static_cast<T>(2.99926757812500L), tolerance);
168 BOOST_CHECK_CLOSE(
169 ::boost::math::rising_factorial(static_cast<T>(-5.25), 6),
170 static_cast<T>(50.987548828125000000000000L), tolerance);
171 BOOST_CHECK_CLOSE(
172 ::boost::math::rising_factorial(static_cast<T>(-5.25), 13),
173 static_cast<T>(127230.91046623885631561279296875000L), tolerance);
174 BOOST_CHECK_CLOSE(
175 ::boost::math::rising_factorial(static_cast<T>(-3.25), -6),
176 static_cast<T>(0.0000129609865918182348202632178291407500332449622510474437452125L), tolerance);
177 BOOST_CHECK_CLOSE(
178 ::boost::math::rising_factorial(static_cast<T>(-5.25), -6),
179 static_cast<T>(2.50789821857946332294524052303699065683926911849535903362649e-6L), tolerance);
180 BOOST_CHECK_CLOSE(
181 ::boost::math::rising_factorial(static_cast<T>(-5.25), -13),
182 static_cast<T>(-1.38984989447269128946284683518361786049649013886981662962096e-14L), tolerance);
183 //
184 // More cases reported as bugs by Rocco Romeo:
185 //
186 BOOST_CHECK_EQUAL(::boost::math::rising_factorial(static_cast<T>(0), 1), static_cast<T>(0));
187 BOOST_CHECK_EQUAL(::boost::math::rising_factorial(static_cast<T>(0), -1), static_cast<T>(-1));
188 BOOST_CHECK_CLOSE(::boost::math::rising_factorial(static_cast<T>(0.5f), -1), static_cast<T>(-2), tolerance);
189 BOOST_CHECK_CLOSE(::boost::math::rising_factorial(static_cast<T>(40.5), -41), static_cast<T>(-2.75643016796662963097096639854040835565778207128865739e-47L), tolerance);
190 BOOST_CHECK_EQUAL(::boost::math::rising_factorial(static_cast<T>(-2), 3), static_cast<T>(0));
191 BOOST_CHECK_EQUAL(::boost::math::rising_factorial(static_cast<T>(-2), 2), static_cast<T>(2));
192 BOOST_CHECK_EQUAL(::boost::math::rising_factorial(static_cast<T>(-4), 3), static_cast<T>(-24));
193 BOOST_CHECK_CLOSE(::boost::math::rising_factorial(static_cast<T>(-4), -3), static_cast<T>(-0.00476190476190476190476190476190476190476190476190476190476L), tolerance);
194 if(ldexp(T(1), -150) != 0)
195 {
196 BOOST_CHECK_CLOSE(::boost::math::rising_factorial(ldexp(T(1), -150), 0), static_cast<T>(1), tolerance);
197 BOOST_CHECK_CLOSE(::boost::math::rising_factorial(ldexp(T(1), -150), -1), static_cast<T>(-1.00000000000000000000000000000000000000000000070064923216241L), tolerance);
198 BOOST_CHECK_CLOSE(::boost::math::rising_factorial(ldexp(T(1), -150), -2), static_cast<T>(0.500000000000000000000000000000000000000000000525486924121806L), tolerance);
199 BOOST_CHECK_CLOSE(::boost::math::rising_factorial(ldexp(T(1), -150), -25), static_cast<T>(-6.44695028438447339619485321920468720529890442766578870603544e-26L), 15 * tolerance);
200 if(std::numeric_limits<T>::min_exponent10 < -50)
201 {
202 BOOST_CHECK_CLOSE(::boost::math::rising_factorial(ldexp(T(1), -150), -40), static_cast<T>(1.22561743912838584942353998493975692372556196815242899727910e-48L), tolerance);
203 }
204 }
205
206 //
207 // Falling factorials:
208 //
209 BOOST_CHECK_CLOSE(
210 ::boost::math::falling_factorial(static_cast<T>(30.25), 0),
211 static_cast<T>(naive_falling_factorial(30.25L, 0)),
212 tolerance);
213 BOOST_CHECK_CLOSE(
214 ::boost::math::falling_factorial(static_cast<T>(30.25), 1),
215 static_cast<T>(naive_falling_factorial(30.25L, 1)),
216 tolerance);
217 BOOST_CHECK_CLOSE(
218 ::boost::math::falling_factorial(static_cast<T>(30.25), 2),
219 static_cast<T>(naive_falling_factorial(30.25L, 2)),
220 tolerance);
221 BOOST_CHECK_CLOSE(
222 ::boost::math::falling_factorial(static_cast<T>(30.25), 5),
223 static_cast<T>(naive_falling_factorial(30.25L, 5)),
224 tolerance);
225 BOOST_CHECK_CLOSE(
226 ::boost::math::falling_factorial(static_cast<T>(30.25), 22),
227 static_cast<T>(naive_falling_factorial(30.25L, 22)),
228 tolerance);
229 BOOST_CHECK_CLOSE(
230 ::boost::math::falling_factorial(static_cast<T>(100.5), 6),
231 static_cast<T>(naive_falling_factorial(100.5L, 6)),
232 tolerance);
233 BOOST_CHECK_CLOSE(
234 ::boost::math::falling_factorial(static_cast<T>(30.75), 30),
235 static_cast<T>(naive_falling_factorial(30.75L, 30)),
236 tolerance * 3);
237 if(boost::math::policies::digits<T, boost::math::policies::policy<> >() > 50)
238 {
239 BOOST_CHECK_CLOSE(
240 ::boost::math::falling_factorial(static_cast<T>(-30.75L), 30),
241 static_cast<T>(naive_falling_factorial(-30.75L, 30)),
242 tolerance * 3);
243 BOOST_CHECK_CLOSE(
244 ::boost::math::falling_factorial(static_cast<T>(-30.75L), 27),
245 static_cast<T>(naive_falling_factorial(-30.75L, 27)),
246 tolerance * 3);
247 }
248 BOOST_CHECK_CLOSE(
249 ::boost::math::falling_factorial(static_cast<T>(-12.0), 6),
250 static_cast<T>(naive_falling_factorial(-12.0L, 6)),
251 tolerance);
252 BOOST_CHECK_CLOSE(
253 ::boost::math::falling_factorial(static_cast<T>(-12), 5),
254 static_cast<T>(naive_falling_factorial(-12.0L, 5)),
255 tolerance);
256 BOOST_CHECK_CLOSE(
257 ::boost::math::falling_factorial(static_cast<T>(-3.0), 6),
258 static_cast<T>(naive_falling_factorial(-3.0L, 6)),
259 tolerance);
260 BOOST_CHECK_CLOSE(
261 ::boost::math::falling_factorial(static_cast<T>(-3), 5),
262 static_cast<T>(naive_falling_factorial(-3.0L, 5)),
263 tolerance);
264 BOOST_CHECK_CLOSE(
265 ::boost::math::falling_factorial(static_cast<T>(3.0), 6),
266 static_cast<T>(naive_falling_factorial(3.0L, 6)),
267 tolerance);
268 BOOST_CHECK_CLOSE(
269 ::boost::math::falling_factorial(static_cast<T>(3), 5),
270 static_cast<T>(naive_falling_factorial(3.0L, 5)),
271 tolerance);
272 BOOST_CHECK_CLOSE(
273 ::boost::math::falling_factorial(static_cast<T>(3.25), 4),
274 static_cast<T>(naive_falling_factorial(3.25L, 4)),
275 tolerance);
276 BOOST_CHECK_CLOSE(
277 ::boost::math::falling_factorial(static_cast<T>(3.25), 5),
278 static_cast<T>(naive_falling_factorial(3.25L, 5)),
279 tolerance);
280 BOOST_CHECK_CLOSE(
281 ::boost::math::falling_factorial(static_cast<T>(3.25), 6),
282 static_cast<T>(naive_falling_factorial(3.25L, 6)),
283 tolerance);
284 BOOST_CHECK_CLOSE(
285 ::boost::math::falling_factorial(static_cast<T>(3.25), 7),
286 static_cast<T>(naive_falling_factorial(3.25L, 7)),
287 tolerance);
288 BOOST_CHECK_CLOSE(
289 ::boost::math::falling_factorial(static_cast<T>(8.25), 12),
290 static_cast<T>(naive_falling_factorial(8.25L, 12)),
291 tolerance);
292 //
293 // More cases reported as bugs by Rocco Romeo:
294 //
295 BOOST_CHECK_EQUAL(::boost::math::falling_factorial(static_cast<T>(0), 1), static_cast<T>(0));
296 BOOST_CHECK_CLOSE(::boost::math::falling_factorial(static_cast<T>(-2), 3), static_cast<T>(-24), tolerance);
297 BOOST_CHECK_CLOSE(::boost::math::falling_factorial(static_cast<T>(-2), 2), static_cast<T>(6), tolerance);
298 BOOST_CHECK_CLOSE(::boost::math::falling_factorial(static_cast<T>(-4), 3), static_cast<T>(-120), tolerance);
299 if(ldexp(static_cast<T>(1), -100) != 0)
300 {
301 BOOST_CHECK_CLOSE(::boost::math::falling_factorial(ldexp(static_cast<T>(1), -100), 1), static_cast<T>(7.888609052210118054117285652827862296732064351090230047e-31L), tolerance);
302 BOOST_CHECK_CLOSE(::boost::math::falling_factorial(ldexp(static_cast<T>(1), -100), 2), static_cast<T>(-7.88860905221011805411728565282163928145420320938308598e-31L), tolerance);
303 BOOST_CHECK_CLOSE(::boost::math::falling_factorial(ldexp(static_cast<T>(1), -100), 3), static_cast<T>(1.577721810442023610823457130563705554763054527705902790e-30L), tolerance);
304 BOOST_CHECK_CLOSE(::boost::math::falling_factorial(ldexp(static_cast<T>(1), -100), 35), static_cast<T>(2.32897613101315187323300837924676081371219290005311315885e8L), 35 * tolerance);
305 }
306 if((ldexp(static_cast<T>(1), -300) != 0) && (std::numeric_limits<T>::max_exponent10 > 290))
307 {
308 BOOST_CHECK_CLOSE(::boost::math::falling_factorial(ldexp(static_cast<T>(1), -300), 20), static_cast<T>(-5.97167167502482975928590631196751639118233432208390100e-74L), tolerance);
309 BOOST_CHECK_CLOSE(::boost::math::falling_factorial(ldexp(static_cast<T>(1), -300), 200), static_cast<T>(-1.93579759151806711025267355739174942986011285920860098569075e282L), 10 * tolerance);
310 }
311
312
313 tolerance = boost::math::tools::epsilon<T>() * 100 * 20; // 20 eps as a percent.
314 unsigned i = boost::math::max_factorial<T>::value;
315 if((boost::is_floating_point<T>::value) && (sizeof(T) <= sizeof(double)))
316 {
317 // Without Lanczos support, tgamma isn't accurate enough for this test:
318 BOOST_CHECK_CLOSE(
319 ::boost::math::unchecked_factorial<T>(i),
320 boost::math::tgamma(static_cast<T>(i+1)), tolerance);
321 }
322
323 i += 10;
324 while(boost::math::lgamma(static_cast<T>(i+1)) < boost::math::tools::log_max_value<T>())
325 {
326 BOOST_CHECK_CLOSE(
327 ::boost::math::factorial<T>(i),
328 boost::math::tgamma(static_cast<T>(i+1)), tolerance);
329 i += 10;
330 }
331} // template <class T> void test_spots(T)
332
333BOOST_AUTO_TEST_CASE( test_main )
334{
335 BOOST_MATH_CONTROL_FP;
336 test_spots(0.0F);
337 test_spots(0.0);
338#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
339 test_spots(0.0L);
340#ifndef BOOST_MATH_NO_REAL_CONCEPT_TESTS
341 test_spots(boost::math::concepts::real_concept(0.));
342#endif
343#else
344 std::cout << "<note>The long double tests have been disabled on this platform "
345 "either because the long double overloads of the usual math functions are "
346 "not available at all, or because they are too inaccurate for these tests "
347 "to pass.</note>" << std::endl;
348#endif
349 if (std::numeric_limits<double>::digits == std::numeric_limits<long double>::digits)
350 {
351 cout << "Types double and long double have the same number of floating-point significand bits ("
352 << std::numeric_limits<long double>::digits << ") on this platform." << endl;
353 }
354 if (std::numeric_limits<float>::digits == std::numeric_limits<double>::digits)
355 {
356 cout << "Types float and double have the same number of floating-point significand bits ("
357 << std::numeric_limits<double>::digits << ") on this platform." << endl;
358 }
359
360 using boost::math::max_factorial;
361 cout << "max factorial for float " << max_factorial<float>::value << endl;
362 cout << "max factorial for double " << max_factorial<double>::value << endl;
363 cout << "max factorial for long double " << max_factorial<long double>::value << endl;
364 cout << "max factorial for real_concept " << max_factorial<boost::math::concepts::real_concept>::value << endl;
365
366
367
368
369}
370
371/*
372
373Output is:
374
375Running 1 test case...
376Types double and long double have the same number of floating-point significand bits (53) on this platform.
377max factorial for float 34
378max factorial for double 170
379max factorial for long double 170
380max factorial for real_concept 100
381*** No errors detected
382
383*/
384
385
386
387
Ceph