libs/math/test/float128/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 #ifdef _MSC_VER
   7 #  pragma warning(disable: 4127) // conditional expression is constant.
   8 #  pragma warning(disable: 4245) // int/unsigned int conversion
   9 #endif
  10
  11 // Return infinities not exceptions:
  12 #define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error
  13
  14 #include <boost/cstdfloat.hpp>
  15 #define BOOST_TEST_MAIN
  16 #include <boost/test/unit_test.hpp>
  17 #include <boost/test/tools/floating_point_comparison.hpp>
  18 #include <boost/math/special_functions/factorials.hpp>
  19 #include <boost/math/special_functions/gamma.hpp>
  20 #include <boost/math/tools/stats.hpp>
  21 #include <boost/math/tools/test.hpp>
  22
  23 #include <iostream>
  24   using std::cout;
  25   using std::endl;
  26
  27 template <class T>
  28 T naive_falling_factorial(T x, unsigned n)
  29 {
  30    if(n == 0)
  31       return 1;
  32    T result = x;
  33    while(--n)
  34    {
  35       x -= 1;
  36       result *= x;
  37    }
  38    return result;
  39 }
  40
  41 template <class T>
  42 void test_spots(T)
  43 {
  44    //
  45    // Basic sanity checks.
  46    //
  47    T tolerance = boost::math::tools::epsilon<T>() * 100 * 2;  // 2 eps as a percent.
  48    BOOST_CHECK_CLOSE(
  49       ::boost::math::factorial<T>(0),
  50       static_cast<T>(1), tolerance);
  51    BOOST_CHECK_CLOSE(
  52       ::boost::math::factorial<T>(1),
  53       static_cast<T>(1), tolerance);
  54    BOOST_CHECK_CLOSE(
  55       ::boost::math::factorial<T>(10),
  56       static_cast<T>(3628800L), tolerance);
  57    BOOST_CHECK_CLOSE(
  58       ::boost::math::unchecked_factorial<T>(0),
  59       static_cast<T>(1), tolerance);
  60    BOOST_CHECK_CLOSE(
  61       ::boost::math::unchecked_factorial<T>(1),
  62       static_cast<T>(1), tolerance);
  63    BOOST_CHECK_CLOSE(
  64       ::boost::math::unchecked_factorial<T>(10),
  65       static_cast<T>(3628800L), tolerance);
  66
  67    //
  68    // Try some double factorials:
  69    //
  70    BOOST_CHECK_CLOSE(
  71       ::boost::math::double_factorial<T>(0),
  72       static_cast<T>(1), tolerance);
  73    BOOST_CHECK_CLOSE(
  74       ::boost::math::double_factorial<T>(1),
  75       static_cast<T>(1), tolerance);
  76    BOOST_CHECK_CLOSE(
  77       ::boost::math::double_factorial<T>(2),
  78       static_cast<T>(2), tolerance);
  79    BOOST_CHECK_CLOSE(
  80       ::boost::math::double_factorial<T>(5),
  81       static_cast<T>(15), tolerance);
  82    BOOST_CHECK_CLOSE(
  83       ::boost::math::double_factorial<T>(10),
  84       static_cast<T>(3840), tolerance);
  85    BOOST_CHECK_CLOSE(
  86       ::boost::math::double_factorial<T>(19),
  87       static_cast<T>(6.547290750e8Q), tolerance);
  88    BOOST_CHECK_CLOSE(
  89       ::boost::math::double_factorial<T>(24),
  90       static_cast<T>(1.961990553600000e12Q), tolerance);
  91    BOOST_CHECK_CLOSE(
  92       ::boost::math::double_factorial<T>(33),
  93       static_cast<T>(6.33265987076285062500000e18Q), tolerance);
  94    BOOST_CHECK_CLOSE(
  95       ::boost::math::double_factorial<T>(42),
  96       static_cast<T>(1.0714547155728479551488000000e26Q), tolerance);
  97    BOOST_CHECK_CLOSE(
  98       ::boost::math::double_factorial<T>(47),
  99       static_cast<T>(1.19256819277443412353990764062500000e30Q), tolerance);
 100
 101    if((std::numeric_limits<T>::has_infinity) && (std::numeric_limits<T>::max_exponent <= 1024))
 102    {
 103       BOOST_CHECK_EQUAL(
 104          ::boost::math::double_factorial<T>(320),
 105          std::numeric_limits<T>::infinity());
 106       BOOST_CHECK_EQUAL(
 107          ::boost::math::double_factorial<T>(301),
 108          std::numeric_limits<T>::infinity());
 109    }
 110    //
 111    // Rising factorials:
 112    //
 113    tolerance = boost::math::tools::epsilon<T>() * 100 * 20;  // 20 eps as a percent.
 114    if(std::numeric_limits<T>::is_specialized == 0)
 115       tolerance *= 5;  // higher error rates without Lanczos support
 116    BOOST_CHECK_CLOSE(
 117       ::boost::math::rising_factorial(static_cast<T>(3), 4),
 118       static_cast<T>(360), tolerance);
 119    BOOST_CHECK_CLOSE(
 120       ::boost::math::rising_factorial(static_cast<T>(7), -4),
 121       static_cast<T>(0.00277777777777777777777777777777777777777777777777777777777778Q), tolerance);
 122    BOOST_CHECK_CLOSE(
 123       ::boost::math::rising_factorial(static_cast<T>(120.5f), 8),
 124       static_cast<T>(5.58187566784927180664062500e16Q), tolerance);
 125    BOOST_CHECK_CLOSE(
 126       ::boost::math::rising_factorial(static_cast<T>(120.5f), -4),
 127       static_cast<T>(5.15881498170104646868208445266116850161120996179812063177241e-9Q), tolerance);
 128    BOOST_CHECK_CLOSE(
 129       ::boost::math::rising_factorial(static_cast<T>(5000.25f), 8),
 130       static_cast<T>(3.92974581976666067544013393509103775024414062500000e29Q), tolerance);
 131    BOOST_CHECK_CLOSE(
 132       ::boost::math::rising_factorial(static_cast<T>(5000.25f), -7),
 133       static_cast<T>(1.28674092710208810281923019294164707555099052561945725535047e-26Q), tolerance);
 134    BOOST_CHECK_CLOSE(
 135       ::boost::math::rising_factorial(static_cast<T>(30.25), 21),
 136       static_cast<T>(3.93286957998925490693364184100209193343633629069699964020401e33Q), tolerance * 2);
 137    BOOST_CHECK_CLOSE(
 138       ::boost::math::rising_factorial(static_cast<T>(30.25), -21),
 139       static_cast<T>(3.35010902064291983728782493133164809108646650368560147505884e-27Q), tolerance);
 140    BOOST_CHECK_CLOSE(
 141       ::boost::math::rising_factorial(static_cast<T>(-30.25), 21),
 142       static_cast<T>(-9.76168312768123676601980433377916854311706629232503473758698e26Q), tolerance * 2);
 143    BOOST_CHECK_CLOSE(
 144       ::boost::math::rising_factorial(static_cast<T>(-30.25), -21),
 145       static_cast<T>(-1.50079704000923674318934280259377728203516775215430875839823e-34Q), 2 * tolerance);
 146    BOOST_CHECK_CLOSE(
 147       ::boost::math::rising_factorial(static_cast<T>(-30.25), 5),
 148       static_cast<T>(-1.78799177197265625000000e7Q), tolerance);
 149    BOOST_CHECK_CLOSE(
 150       ::boost::math::rising_factorial(static_cast<T>(-30.25), -5),
 151       static_cast<T>(-2.47177487004482195012362027432181137141899692171397467859150e-8Q), tolerance);
 152    BOOST_CHECK_CLOSE(
 153       ::boost::math::rising_factorial(static_cast<T>(-30.25), 6),
 154       static_cast<T>(4.5146792242309570312500000e8Q), tolerance);
 155    BOOST_CHECK_CLOSE(
 156       ::boost::math::rising_factorial(static_cast<T>(-30.25), -6),
 157       static_cast<T>(6.81868929667537089689274558433603136943171564610751635473516e-10Q), tolerance);
 158    BOOST_CHECK_CLOSE(
 159       ::boost::math::rising_factorial(static_cast<T>(-3), 6),
 160       static_cast<T>(0), tolerance);
 161    BOOST_CHECK_CLOSE(
 162       ::boost::math::rising_factorial(static_cast<T>(-3.25), 6),
 163       static_cast<T>(2.99926757812500Q), tolerance);
 164    BOOST_CHECK_CLOSE(
 165       ::boost::math::rising_factorial(static_cast<T>(-5.25), 6),
 166       static_cast<T>(50.987548828125000000000000Q), tolerance);
 167    BOOST_CHECK_CLOSE(
 168       ::boost::math::rising_factorial(static_cast<T>(-5.25), 13),
 169       static_cast<T>(127230.91046623885631561279296875000Q), tolerance);
 170    BOOST_CHECK_CLOSE(
 171       ::boost::math::rising_factorial(static_cast<T>(-3.25), -6),
 172       static_cast<T>(0.0000129609865918182348202632178291407500332449622510474437452125Q), tolerance);
 173    BOOST_CHECK_CLOSE(
 174       ::boost::math::rising_factorial(static_cast<T>(-5.25), -6),
 175       static_cast<T>(2.50789821857946332294524052303699065683926911849535903362649e-6Q), tolerance);
 176    BOOST_CHECK_CLOSE(
 177       ::boost::math::rising_factorial(static_cast<T>(-5.25), -13),
 178       static_cast<T>(-1.38984989447269128946284683518361786049649013886981662962096e-14Q), tolerance);
 179
 180    //
 181    // Falling factorials:
 182    //
 183    BOOST_CHECK_CLOSE(
 184       ::boost::math::falling_factorial(static_cast<T>(30.25), 0),
 185       static_cast<T>(naive_falling_factorial(30.25Q, 0)),
 186       tolerance);
 187    BOOST_CHECK_CLOSE(
 188       ::boost::math::falling_factorial(static_cast<T>(30.25), 1),
 189       static_cast<T>(naive_falling_factorial(30.25Q, 1)),
 190       tolerance);
 191    BOOST_CHECK_CLOSE(
 192       ::boost::math::falling_factorial(static_cast<T>(30.25), 2),
 193       static_cast<T>(naive_falling_factorial(30.25Q, 2)),
 194       tolerance);
 195    BOOST_CHECK_CLOSE(
 196       ::boost::math::falling_factorial(static_cast<T>(30.25), 5),
 197       static_cast<T>(naive_falling_factorial(30.25Q, 5)),
 198       tolerance);
 199    BOOST_CHECK_CLOSE(
 200       ::boost::math::falling_factorial(static_cast<T>(30.25), 22),
 201       static_cast<T>(naive_falling_factorial(30.25Q, 22)),
 202       tolerance);
 203    BOOST_CHECK_CLOSE(
 204       ::boost::math::falling_factorial(static_cast<T>(100.5), 6),
 205       static_cast<T>(naive_falling_factorial(100.5Q, 6)),
 206       tolerance);
 207    BOOST_CHECK_CLOSE(
 208       ::boost::math::falling_factorial(static_cast<T>(30.75), 30),
 209       static_cast<T>(naive_falling_factorial(30.75Q, 30)),
 210       tolerance * 3);
 211    if(boost::math::policies::digits<T, boost::math::policies::policy<> >() > 50)
 212    {
 213       BOOST_CHECK_CLOSE(
 214          ::boost::math::falling_factorial(static_cast<T>(-30.75Q), 30),
 215          static_cast<T>(naive_falling_factorial(-30.75Q, 30)),
 216          tolerance * 3);
 217       BOOST_CHECK_CLOSE(
 218          ::boost::math::falling_factorial(static_cast<T>(-30.75Q), 27),
 219          static_cast<T>(naive_falling_factorial(-30.75Q, 27)),
 220          tolerance * 3);
 221    }
 222    BOOST_CHECK_CLOSE(
 223       ::boost::math::falling_factorial(static_cast<T>(-12.0), 6),
 224       static_cast<T>(naive_falling_factorial(-12.0Q, 6)),
 225       tolerance);
 226    BOOST_CHECK_CLOSE(
 227       ::boost::math::falling_factorial(static_cast<T>(-12), 5),
 228       static_cast<T>(naive_falling_factorial(-12.0Q, 5)),
 229       tolerance);
 230    BOOST_CHECK_CLOSE(
 231       ::boost::math::falling_factorial(static_cast<T>(-3.0), 6),
 232       static_cast<T>(naive_falling_factorial(-3.0Q, 6)),
 233       tolerance);
 234    BOOST_CHECK_CLOSE(
 235       ::boost::math::falling_factorial(static_cast<T>(-3), 5),
 236       static_cast<T>(naive_falling_factorial(-3.0Q, 5)),
 237       tolerance);
 238    BOOST_CHECK_CLOSE(
 239       ::boost::math::falling_factorial(static_cast<T>(3.0), 6),
 240       static_cast<T>(naive_falling_factorial(3.0Q, 6)),
 241       tolerance);
 242    BOOST_CHECK_CLOSE(
 243       ::boost::math::falling_factorial(static_cast<T>(3), 5),
 244       static_cast<T>(naive_falling_factorial(3.0Q, 5)),
 245       tolerance);
 246    BOOST_CHECK_CLOSE(
 247       ::boost::math::falling_factorial(static_cast<T>(3.25), 4),
 248       static_cast<T>(naive_falling_factorial(3.25Q, 4)),
 249       tolerance);
 250    BOOST_CHECK_CLOSE(
 251       ::boost::math::falling_factorial(static_cast<T>(3.25), 5),
 252       static_cast<T>(naive_falling_factorial(3.25Q, 5)),
 253       tolerance);
 254    BOOST_CHECK_CLOSE(
 255       ::boost::math::falling_factorial(static_cast<T>(3.25), 6),
 256       static_cast<T>(naive_falling_factorial(3.25Q, 6)),
 257       tolerance);
 258    BOOST_CHECK_CLOSE(
 259       ::boost::math::falling_factorial(static_cast<T>(3.25), 7),
 260       static_cast<T>(naive_falling_factorial(3.25Q, 7)),
 261       tolerance);
 262    BOOST_CHECK_CLOSE(
 263       ::boost::math::falling_factorial(static_cast<T>(8.25), 12),
 264       static_cast<T>(naive_falling_factorial(8.25Q, 12)),
 265       tolerance);
 266
 267
 268    tolerance = boost::math::tools::epsilon<T>() * 100 * 20;  // 20 eps as a percent.
 269    unsigned i = boost::math::max_factorial<T>::value;
 270    if((boost::is_floating_point<T>::value) && (sizeof(T) <= sizeof(double)))
 271    {
 272       // Without Lanczos support, tgamma isn't accurate enough for this test:
 273       BOOST_CHECK_CLOSE(
 274          ::boost::math::unchecked_factorial<T>(i),
 275          boost::math::tgamma(static_cast<T>(i+1)), tolerance);
 276    }
 277
 278    i += 10;
 279    while(boost::math::lgamma(static_cast<T>(i+1)) < boost::math::tools::log_max_value<T>())
 280    {
 281       BOOST_CHECK_CLOSE(
 282          ::boost::math::factorial<T>(i),
 283          boost::math::tgamma(static_cast<T>(i+1)), tolerance);
 284       i += 10;
 285    }
 286 } // template <class T> void test_spots(T)
 287
 288 BOOST_AUTO_TEST_CASE( test_main )
 289 {
 290    BOOST_MATH_CONTROL_FP;
 291    test_spots(0.0Q);
 292    cout << "max factorial for __float128"  << boost::math::max_factorial<boost::floatmax_t>::value  << endl;
 293 }
 294
 295
 296