libs/math/test/test_factorials.cpp

   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>
  19 #include <boost/test/tools/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>
  24
  25 #include <iostream>
  26 #include <iomanip>
  27 using std::cout;
  28 using std::endl;
  29
  30 template <class T>
  31 T 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
  44 template <class T>
  45 void test_spots(T)
  46 {
  47    //
  48    // Basic sanity checks.
  49    //
  50    T tolerance = boost::math::tools::epsilon<T>() * 100 * 2;  // 2 eps as a percent.
  51    BOOST_CHECK_CLOSE(
  52       ::boost::math::factorial<T>(0),
  53       static_cast<T>(1), tolerance);
  54    BOOST_CHECK_CLOSE(
  55       ::boost::math::factorial<T>(1),
  56       static_cast<T>(1), tolerance);
  57    BOOST_CHECK_CLOSE(
  58       ::boost::math::factorial<T>(10),
  59       static_cast<T>(3628800L), tolerance);
  60    BOOST_CHECK_CLOSE(
  61       ::boost::math::unchecked_factorial<T>(0),
  62       static_cast<T>(1), tolerance);
  63    BOOST_CHECK_CLOSE(
  64       ::boost::math::unchecked_factorial<T>(1),
  65       static_cast<T>(1), tolerance);
  66    BOOST_CHECK_CLOSE(
  67       ::boost::math::unchecked_factorial<T>(10),
  68       static_cast<T>(3628800L), tolerance);
  69
  70    //
  71    // Try some double factorials:
  72    //
  73    BOOST_CHECK_CLOSE(
  74       ::boost::math::double_factorial<T>(0),
  75       static_cast<T>(1), tolerance);
  76    BOOST_CHECK_CLOSE(
  77       ::boost::math::double_factorial<T>(1),
  78       static_cast<T>(1), tolerance);
  79    BOOST_CHECK_CLOSE(
  80       ::boost::math::double_factorial<T>(2),
  81       static_cast<T>(2), tolerance);
  82    BOOST_CHECK_CLOSE(
  83       ::boost::math::double_factorial<T>(5),
  84       static_cast<T>(15), tolerance);
  85    BOOST_CHECK_CLOSE(
  86       ::boost::math::double_factorial<T>(10),
  87       static_cast<T>(3840), tolerance);
  88    BOOST_CHECK_CLOSE(
  89       ::boost::math::double_factorial<T>(19),
  90       static_cast<T>(6.547290750e8L), tolerance);
  91    BOOST_CHECK_CLOSE(
  92       ::boost::math::double_factorial<T>(24),
  93       static_cast<T>(1.961990553600000e12L), tolerance);
  94    BOOST_CHECK_CLOSE(
  95       ::boost::math::double_factorial<T>(33),
  96       static_cast<T>(6.33265987076285062500000e18L), tolerance);
  97    BOOST_CHECK_CLOSE(
  98       ::boost::math::double_factorial<T>(42),
  99       static_cast<T>(1.0714547155728479551488000000e26L), tolerance);
 100    BOOST_CHECK_CLOSE(
 101       ::boost::math::double_factorial<T>(47),
 102       static_cast<T>(1.19256819277443412353990764062500000e30L), tolerance);
 103
 104    if((std::numeric_limits<T>::has_infinity) && (std::numeric_limits<T>::max_exponent <= 1024))
 105    {
 106       BOOST_CHECK_EQUAL(
 107          ::boost::math::double_factorial<T>(320),
 108          std::numeric_limits<T>::infinity());
 109       BOOST_CHECK_EQUAL(
 110          ::boost::math::double_factorial<T>(301),
 111          std::numeric_limits<T>::infinity());
 112    }
 113    //
 114    // Rising factorials:
 115    //
 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
 119    BOOST_CHECK_CLOSE(
 120       ::boost::math::rising_factorial(static_cast<T>(3), 4),
 121       static_cast<T>(360), tolerance);
 122    BOOST_CHECK_CLOSE(
 123       ::boost::math::rising_factorial(static_cast<T>(7), -4),
 124       static_cast<T>(0.00277777777777777777777777777777777777777777777777777777777778L), tolerance);
 125    BOOST_CHECK_CLOSE(
 126       ::boost::math::rising_factorial(static_cast<T>(120.5f), 8),
 127       static_cast<T>(5.58187566784927180664062500e16L), tolerance);
 128    BOOST_CHECK_CLOSE(
 129       ::boost::math::rising_factorial(static_cast<T>(120.5f), -4),
 130       static_cast<T>(5.15881498170104646868208445266116850161120996179812063177241e-9L), tolerance);
 131    BOOST_CHECK_CLOSE(
 132       ::boost::math::rising_factorial(static_cast<T>(5000.25f), 8),
 133       static_cast<T>(3.92974581976666067544013393509103775024414062500000e29L), tolerance);
 134    BOOST_CHECK_CLOSE(
 135       ::boost::math::rising_factorial(static_cast<T>(5000.25f), -7),
 136       static_cast<T>(1.28674092710208810281923019294164707555099052561945725535047e-26L), tolerance);
 137    BOOST_CHECK_CLOSE(
 138       ::boost::math::rising_factorial(static_cast<T>(30.25), 21),
 139       static_cast<T>(3.93286957998925490693364184100209193343633629069699964020401e33L), tolerance * 2);
 140    BOOST_CHECK_CLOSE(
 141       ::boost::math::rising_factorial(static_cast<T>(30.25), -21),
 142       static_cast<T>(3.35010902064291983728782493133164809108646650368560147505884e-27L), tolerance);
 143    BOOST_CHECK_CLOSE(
 144       ::boost::math::rising_factorial(static_cast<T>(-30.25), 21),
 145       static_cast<T>(-9.76168312768123676601980433377916854311706629232503473758698e26L), tolerance * 2);
 146    BOOST_CHECK_CLOSE(
 147       ::boost::math::rising_factorial(static_cast<T>(-30.25), -21),
 148       static_cast<T>(-1.50079704000923674318934280259377728203516775215430875839823e-34L), 2 * tolerance);
 149    BOOST_CHECK_CLOSE(
 150       ::boost::math::rising_factorial(static_cast<T>(-30.25), 5),
 151       static_cast<T>(-1.78799177197265625000000e7L), tolerance);
 152    BOOST_CHECK_CLOSE(
 153       ::boost::math::rising_factorial(static_cast<T>(-30.25), -5),
 154       static_cast<T>(-2.47177487004482195012362027432181137141899692171397467859150e-8L), tolerance);
 155    BOOST_CHECK_CLOSE(
 156       ::boost::math::rising_factorial(static_cast<T>(-30.25), 6),
 157       static_cast<T>(4.5146792242309570312500000e8L), tolerance);
 158    BOOST_CHECK_CLOSE(
 159       ::boost::math::rising_factorial(static_cast<T>(-30.25), -6),
 160       static_cast<T>(6.81868929667537089689274558433603136943171564610751635473516e-10L), tolerance);
 161    BOOST_CHECK_CLOSE(
 162       ::boost::math::rising_factorial(static_cast<T>(-3), 6),
 163       static_cast<T>(0), tolerance);
 164    BOOST_CHECK_CLOSE(
 165       ::boost::math::rising_factorial(static_cast<T>(-3.25), 6),
 166       static_cast<T>(2.99926757812500L), tolerance);
 167    BOOST_CHECK_CLOSE(
 168       ::boost::math::rising_factorial(static_cast<T>(-5.25), 6),
 169       static_cast<T>(50.987548828125000000000000L), tolerance);
 170    BOOST_CHECK_CLOSE(
 171       ::boost::math::rising_factorial(static_cast<T>(-5.25), 13),
 172       static_cast<T>(127230.91046623885631561279296875000L), tolerance);
 173    BOOST_CHECK_CLOSE(
 174       ::boost::math::rising_factorial(static_cast<T>(-3.25), -6),
 175       static_cast<T>(0.0000129609865918182348202632178291407500332449622510474437452125L), tolerance);
 176    BOOST_CHECK_CLOSE(
 177       ::boost::math::rising_factorial(static_cast<T>(-5.25), -6),
 178       static_cast<T>(2.50789821857946332294524052303699065683926911849535903362649e-6L), tolerance);
 179    BOOST_CHECK_CLOSE(
 180       ::boost::math::rising_factorial(static_cast<T>(-5.25), -13),
 181       static_cast<T>(-1.38984989447269128946284683518361786049649013886981662962096e-14L), tolerance);
 182    //
 183    // More cases reported as bugs by Rocco Romeo:
 184    //
 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)
 194    {
 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)
 200       {
 201          BOOST_CHECK_CLOSE(::boost::math::rising_factorial(ldexp(T(1), -150), -40), static_cast<T>(1.22561743912838584942353998493975692372556196815242899727910e-48L), tolerance);
 202       }
 203    }
 204
 205    //
 206    // Falling factorials:
 207    //
 208    BOOST_CHECK_CLOSE(
 209       ::boost::math::falling_factorial(static_cast<T>(30.25), 0),
 210       static_cast<T>(naive_falling_factorial(30.25L, 0)),
 211       tolerance);
 212    BOOST_CHECK_CLOSE(
 213       ::boost::math::falling_factorial(static_cast<T>(30.25), 1),
 214       static_cast<T>(naive_falling_factorial(30.25L, 1)),
 215       tolerance);
 216    BOOST_CHECK_CLOSE(
 217       ::boost::math::falling_factorial(static_cast<T>(30.25), 2),
 218       static_cast<T>(naive_falling_factorial(30.25L, 2)),
 219       tolerance);
 220    BOOST_CHECK_CLOSE(
 221       ::boost::math::falling_factorial(static_cast<T>(30.25), 5),
 222       static_cast<T>(naive_falling_factorial(30.25L, 5)),
 223       tolerance);
 224    BOOST_CHECK_CLOSE(
 225       ::boost::math::falling_factorial(static_cast<T>(30.25), 22),
 226       static_cast<T>(naive_falling_factorial(30.25L, 22)),
 227       tolerance);
 228    BOOST_CHECK_CLOSE(
 229       ::boost::math::falling_factorial(static_cast<T>(100.5), 6),
 230       static_cast<T>(naive_falling_factorial(100.5L, 6)),
 231       tolerance);
 232    BOOST_CHECK_CLOSE(
 233       ::boost::math::falling_factorial(static_cast<T>(30.75), 30),
 234       static_cast<T>(naive_falling_factorial(30.75L, 30)),
 235       tolerance * 3);
 236    if(boost::math::policies::digits<T, boost::math::policies::policy<> >() > 50)
 237    {
 238       BOOST_CHECK_CLOSE(
 239          ::boost::math::falling_factorial(static_cast<T>(-30.75L), 30),
 240          static_cast<T>(naive_falling_factorial(-30.75L, 30)),
 241          tolerance * 3);
 242       BOOST_CHECK_CLOSE(
 243          ::boost::math::falling_factorial(static_cast<T>(-30.75L), 27),
 244          static_cast<T>(naive_falling_factorial(-30.75L, 27)),
 245          tolerance * 3);
 246    }
 247    BOOST_CHECK_CLOSE(
 248       ::boost::math::falling_factorial(static_cast<T>(-12.0), 6),
 249       static_cast<T>(naive_falling_factorial(-12.0L, 6)),
 250       tolerance);
 251    BOOST_CHECK_CLOSE(
 252       ::boost::math::falling_factorial(static_cast<T>(-12), 5),
 253       static_cast<T>(naive_falling_factorial(-12.0L, 5)),
 254       tolerance);
 255    BOOST_CHECK_CLOSE(
 256       ::boost::math::falling_factorial(static_cast<T>(-3.0), 6),
 257       static_cast<T>(naive_falling_factorial(-3.0L, 6)),
 258       tolerance);
 259    BOOST_CHECK_CLOSE(
 260       ::boost::math::falling_factorial(static_cast<T>(-3), 5),
 261       static_cast<T>(naive_falling_factorial(-3.0L, 5)),
 262       tolerance);
 263    BOOST_CHECK_CLOSE(
 264       ::boost::math::falling_factorial(static_cast<T>(3.0), 6),
 265       static_cast<T>(naive_falling_factorial(3.0L, 6)),
 266       tolerance);
 267    BOOST_CHECK_CLOSE(
 268       ::boost::math::falling_factorial(static_cast<T>(3), 5),
 269       static_cast<T>(naive_falling_factorial(3.0L, 5)),
 270       tolerance);
 271    BOOST_CHECK_CLOSE(
 272       ::boost::math::falling_factorial(static_cast<T>(3.25), 4),
 273       static_cast<T>(naive_falling_factorial(3.25L, 4)),
 274       tolerance);
 275    BOOST_CHECK_CLOSE(
 276       ::boost::math::falling_factorial(static_cast<T>(3.25), 5),
 277       static_cast<T>(naive_falling_factorial(3.25L, 5)),
 278       tolerance);
 279    BOOST_CHECK_CLOSE(
 280       ::boost::math::falling_factorial(static_cast<T>(3.25), 6),
 281       static_cast<T>(naive_falling_factorial(3.25L, 6)),
 282       tolerance);
 283    BOOST_CHECK_CLOSE(
 284       ::boost::math::falling_factorial(static_cast<T>(3.25), 7),
 285       static_cast<T>(naive_falling_factorial(3.25L, 7)),
 286       tolerance);
 287    BOOST_CHECK_CLOSE(
 288       ::boost::math::falling_factorial(static_cast<T>(8.25), 12),
 289       static_cast<T>(naive_falling_factorial(8.25L, 12)),
 290       tolerance);
 291    //
 292    // More cases reported as bugs by Rocco Romeo:
 293    //
 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)
 299    {
 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);
 304    }
 305    if((ldexp(static_cast<T>(1), -300) != 0) && (std::numeric_limits<T>::max_exponent10 > 290))
 306    {
 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);
 309    }
 310
 311
 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)))
 315    {
 316       // Without Lanczos support, tgamma isn't accurate enough for this test:
 317       BOOST_CHECK_CLOSE(
 318          ::boost::math::unchecked_factorial<T>(i),
 319          boost::math::tgamma(static_cast<T>(i+1)), tolerance);
 320    }
 321
 322    i += 10;
 323    while(boost::math::lgamma(static_cast<T>(i+1)) < boost::math::tools::log_max_value<T>())
 324    {
 325       BOOST_CHECK_CLOSE(
 326          ::boost::math::factorial<T>(i),
 327          boost::math::tgamma(static_cast<T>(i+1)), tolerance);
 328       i += 10;
 329    }
 330 } // template <class T> void test_spots(T)
 331
 332 BOOST_AUTO_TEST_CASE( test_main )
 333 {
 334    BOOST_MATH_CONTROL_FP;
 335    test_spots(0.0F);
 336    test_spots(0.0);
 337 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
 338    test_spots(0.0L);
 339 #ifndef BOOST_MATH_NO_REAL_CONCEPT_TESTS
 340    test_spots(boost::math::concepts::real_concept(0.));
 341 #endif
 342 #else
 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;
 347 #endif
 348    if (std::numeric_limits<double>::digits == std::numeric_limits<long double>::digits)
 349    {
 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;
 352    }
 353    if (std::numeric_limits<float>::digits == std::numeric_limits<double>::digits)
 354    {
 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;
 357    }
 358
 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;
 364
 365
 366
 367
 368 }
 369
 370 /*
 371
 372 Output is:
 373
 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
 381
 382 */
 383
 384
 385
 386