boost/math/special_functions/detail/lgamma_small.hpp

   1 //  (C) 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 #ifndef BOOST_MATH_SPECIAL_FUNCTIONS_DETAIL_LGAMMA_SMALL
   7 #define BOOST_MATH_SPECIAL_FUNCTIONS_DETAIL_LGAMMA_SMALL
   8
   9 #ifdef _MSC_VER
  10 #pragma once
  11 #endif
  12
  13 namespace boost{ namespace math{ namespace detail{
  14
  15 //
  16 // lgamma for small arguments:
  17 //
  18 template <class T, class Policy, class L>
  19 T lgamma_small_imp(T z, T zm1, T zm2, const mpl::int_<64>&, const Policy& /* l */, const L&)
  20 {
  21    // This version uses rational approximations for small
  22    // values of z accurate enough for 64-bit mantissas
  23    // (80-bit long doubles), works well for 53-bit doubles as well.
  24    // L is only used to select the Lanczos function.
  25
  26    BOOST_MATH_STD_USING  // for ADL of std names
  27    T result = 0;
  28    if(z < tools::epsilon<T>())
  29    {
  30       result = -log(z);
  31    }
  32    else if((zm1 == 0) || (zm2 == 0))
  33    {
  34       // nothing to do, result is zero....
  35    }
  36    else if(z > 2)
  37    {
  38       //
  39       // Begin by performing argument reduction until
  40       // z is in [2,3):
  41       //
  42       if(z >= 3)
  43       {
  44          do
  45          {
  46             z -= 1;
  47             zm2 -= 1;
  48             result += log(z);
  49          }while(z >= 3);
  50          // Update zm2, we need it below:
  51          zm2 = z - 2;
  52       }
  53
  54       //
  55       // Use the following form:
  56       //
  57       // lgamma(z) = (z-2)(z+1)(Y + R(z-2))
  58       //
  59       // where R(z-2) is a rational approximation optimised for
  60       // low absolute error - as long as it's absolute error
  61       // is small compared to the constant Y - then any rounding
  62       // error in it's computation will get wiped out.
  63       //
  64       // R(z-2) has the following properties:
  65       //
  66       // At double: Max error found:                    4.231e-18
  67       // At long double: Max error found:               1.987e-21
  68       // Maximum Deviation Found (approximation error): 5.900e-24
  69       //
  70       static const T P[] = {
  71          static_cast<T>(-0.180355685678449379109e-1L),
  72          static_cast<T>(0.25126649619989678683e-1L),
  73          static_cast<T>(0.494103151567532234274e-1L),
  74          static_cast<T>(0.172491608709613993966e-1L),
  75          static_cast<T>(-0.259453563205438108893e-3L),
  76          static_cast<T>(-0.541009869215204396339e-3L),
  77          static_cast<T>(-0.324588649825948492091e-4L)
  78       };
  79       static const T Q[] = {
  80          static_cast<T>(0.1e1),
  81          static_cast<T>(0.196202987197795200688e1L),
  82          static_cast<T>(0.148019669424231326694e1L),
  83          static_cast<T>(0.541391432071720958364e0L),
  84          static_cast<T>(0.988504251128010129477e-1L),
  85          static_cast<T>(0.82130967464889339326e-2L),
  86          static_cast<T>(0.224936291922115757597e-3L),
  87          static_cast<T>(-0.223352763208617092964e-6L)
  88       };
  89
  90       static const float Y = 0.158963680267333984375e0f;
  91
  92       T r = zm2 * (z + 1);
  93       T R = tools::evaluate_polynomial(P, zm2);
  94       R /= tools::evaluate_polynomial(Q, zm2);
  95
  96       result +=  r * Y + r * R;
  97    }
  98    else
  99    {
 100       //
 101       // If z is less than 1 use recurrance to shift to
 102       // z in the interval [1,2]:
 103       //
 104       if(z < 1)
 105       {
 106          result += -log(z);
 107          zm2 = zm1;
 108          zm1 = z;
 109          z += 1;
 110       }
 111       //
 112       // Two approximations, on for z in [1,1.5] and
 113       // one for z in [1.5,2]:
 114       //
 115       if(z <= 1.5)
 116       {
 117          //
 118          // Use the following form:
 119          //
 120          // lgamma(z) = (z-1)(z-2)(Y + R(z-1))
 121          //
 122          // where R(z-1) is a rational approximation optimised for
 123          // low absolute error - as long as it's absolute error
 124          // is small compared to the constant Y - then any rounding
 125          // error in it's computation will get wiped out.
 126          //
 127          // R(z-1) has the following properties:
 128          //
 129          // At double precision: Max error found:                1.230011e-17
 130          // At 80-bit long double precision:   Max error found:  5.631355e-21
 131          // Maximum Deviation Found:                             3.139e-021
 132          // Expected Error Term:                                 3.139e-021
 133
 134          //
 135          static const float Y = 0.52815341949462890625f;
 136
 137          static const T P[] = {
 138             static_cast<T>(0.490622454069039543534e-1L),
 139             static_cast<T>(-0.969117530159521214579e-1L),
 140             static_cast<T>(-0.414983358359495381969e0L),
 141             static_cast<T>(-0.406567124211938417342e0L),
 142             static_cast<T>(-0.158413586390692192217e0L),
 143             static_cast<T>(-0.240149820648571559892e-1L),
 144             static_cast<T>(-0.100346687696279557415e-2L)
 145          };
 146          static const T Q[] = {
 147             static_cast<T>(0.1e1L),
 148             static_cast<T>(0.302349829846463038743e1L),
 149             static_cast<T>(0.348739585360723852576e1L),
 150             static_cast<T>(0.191415588274426679201e1L),
 151             static_cast<T>(0.507137738614363510846e0L),
 152             static_cast<T>(0.577039722690451849648e-1L),
 153             static_cast<T>(0.195768102601107189171e-2L)
 154          };
 155
 156          T r = tools::evaluate_polynomial(P, zm1) / tools::evaluate_polynomial(Q, zm1);
 157          T prefix = zm1 * zm2;
 158
 159          result += prefix * Y + prefix * r;
 160       }
 161       else
 162       {
 163          //
 164          // Use the following form:
 165          //
 166          // lgamma(z) = (2-z)(1-z)(Y + R(2-z))
 167          //
 168          // where R(2-z) is a rational approximation optimised for
 169          // low absolute error - as long as it's absolute error
 170          // is small compared to the constant Y - then any rounding
 171          // error in it's computation will get wiped out.
 172          //
 173          // R(2-z) has the following properties:
 174          //
 175          // At double precision, max error found:              1.797565e-17
 176          // At 80-bit long double precision, max error found:  9.306419e-21
 177          // Maximum Deviation Found:                           2.151e-021
 178          // Expected Error Term:                               2.150e-021
 179          //
 180          static const float Y = 0.452017307281494140625f;
 181
 182          static const T P[] = {
 183             static_cast<T>(-0.292329721830270012337e-1L),
 184             static_cast<T>(0.144216267757192309184e0L),
 185             static_cast<T>(-0.142440390738631274135e0L),
 186             static_cast<T>(0.542809694055053558157e-1L),
 187             static_cast<T>(-0.850535976868336437746e-2L),
 188             static_cast<T>(0.431171342679297331241e-3L)
 189          };
 190          static const T Q[] = {
 191             static_cast<T>(0.1e1),
 192             static_cast<T>(-0.150169356054485044494e1L),
 193             static_cast<T>(0.846973248876495016101e0L),
 194             static_cast<T>(-0.220095151814995745555e0L),
 195             static_cast<T>(0.25582797155975869989e-1L),
 196             static_cast<T>(-0.100666795539143372762e-2L),
 197             static_cast<T>(-0.827193521891290553639e-6L)
 198          };
 199          T r = zm2 * zm1;
 200          T R = tools::evaluate_polynomial(P, -zm2) / tools::evaluate_polynomial(Q, -zm2);
 201
 202          result += r * Y + r * R;
 203       }
 204    }
 205    return result;
 206 }
 207 template <class T, class Policy, class L>
 208 T lgamma_small_imp(T z, T zm1, T zm2, const mpl::int_<113>&, const Policy& /* l */, const L&)
 209 {
 210    //
 211    // This version uses rational approximations for small
 212    // values of z accurate enough for 113-bit mantissas
 213    // (128-bit long doubles).
 214    //
 215    BOOST_MATH_STD_USING  // for ADL of std names
 216    T result = 0;
 217    if(z < tools::epsilon<T>())
 218    {
 219       result = -log(z);
 220       BOOST_MATH_INSTRUMENT_CODE(result);
 221    }
 222    else if((zm1 == 0) || (zm2 == 0))
 223    {
 224       // nothing to do, result is zero....
 225    }
 226    else if(z > 2)
 227    {
 228       //
 229       // Begin by performing argument reduction until
 230       // z is in [2,3):
 231       //
 232       if(z >= 3)
 233       {
 234          do
 235          {
 236             z -= 1;
 237             result += log(z);
 238          }while(z >= 3);
 239          zm2 = z - 2;
 240       }
 241       BOOST_MATH_INSTRUMENT_CODE(zm2);
 242       BOOST_MATH_INSTRUMENT_CODE(z);
 243       BOOST_MATH_INSTRUMENT_CODE(result);
 244
 245       //
 246       // Use the following form:
 247       //
 248       // lgamma(z) = (z-2)(z+1)(Y + R(z-2))
 249       //
 250       // where R(z-2) is a rational approximation optimised for
 251       // low absolute error - as long as it's absolute error
 252       // is small compared to the constant Y - then any rounding
 253       // error in it's computation will get wiped out.
 254       //
 255       // Maximum Deviation Found (approximation error)      3.73e-37
 256
 257       static const T P[] = {
 258          -0.018035568567844937910504030027467476655L,
 259          0.013841458273109517271750705401202404195L,
 260          0.062031842739486600078866923383017722399L,
 261          0.052518418329052161202007865149435256093L,
 262          0.01881718142472784129191838493267755758L,
 263          0.0025104830367021839316463675028524702846L,
 264          -0.00021043176101831873281848891452678568311L,
 265          -0.00010249622350908722793327719494037981166L,
 266          -0.11381479670982006841716879074288176994e-4L,
 267          -0.49999811718089980992888533630523892389e-6L,
 268          -0.70529798686542184668416911331718963364e-8L
 269       };
 270       static const T Q[] = {
 271          1L,
 272          2.5877485070422317542808137697939233685L,
 273          2.8797959228352591788629602533153837126L,
 274          1.8030885955284082026405495275461180977L,
 275          0.69774331297747390169238306148355428436L,
 276          0.17261566063277623942044077039756583802L,
 277          0.02729301254544230229429621192443000121L,
 278          0.0026776425891195270663133581960016620433L,
 279          0.00015244249160486584591370355730402168106L,
 280          0.43997034032479866020546814475414346627e-5L,
 281          0.46295080708455613044541885534408170934e-7L,
 282          -0.93326638207459533682980757982834180952e-11L,
 283          0.42316456553164995177177407325292867513e-13L
 284       };
 285
 286       T R = tools::evaluate_polynomial(P, zm2);
 287       R /= tools::evaluate_polynomial(Q, zm2);
 288
 289       static const float Y = 0.158963680267333984375F;
 290
 291       T r = zm2 * (z + 1);
 292
 293       result +=  r * Y + r * R;
 294       BOOST_MATH_INSTRUMENT_CODE(result);
 295    }
 296    else
 297    {
 298       //
 299       // If z is less than 1 use recurrance to shift to
 300       // z in the interval [1,2]:
 301       //
 302       if(z < 1)
 303       {
 304          result += -log(z);
 305          zm2 = zm1;
 306          zm1 = z;
 307          z += 1;
 308       }
 309       BOOST_MATH_INSTRUMENT_CODE(result);
 310       BOOST_MATH_INSTRUMENT_CODE(z);
 311       BOOST_MATH_INSTRUMENT_CODE(zm2);
 312       //
 313       // Three approximations, on for z in [1,1.35], [1.35,1.625] and [1.625,1]
 314       //
 315       if(z <= 1.35)
 316       {
 317          //
 318          // Use the following form:
 319          //
 320          // lgamma(z) = (z-1)(z-2)(Y + R(z-1))
 321          //
 322          // where R(z-1) is a rational approximation optimised for
 323          // low absolute error - as long as it's absolute error
 324          // is small compared to the constant Y - then any rounding
 325          // error in it's computation will get wiped out.
 326          //
 327          // R(z-1) has the following properties:
 328          //
 329          // Maximum Deviation Found (approximation error)            1.659e-36
 330          // Expected Error Term (theoretical error)                  1.343e-36
 331          // Max error found at 128-bit long double precision         1.007e-35
 332          //
 333          static const float Y = 0.54076099395751953125f;
 334
 335          static const T P[] = {
 336             0.036454670944013329356512090082402429697L,
 337             -0.066235835556476033710068679907798799959L,
 338             -0.67492399795577182387312206593595565371L,
 339             -1.4345555263962411429855341651960000166L,
 340             -1.4894319559821365820516771951249649563L,
 341             -0.87210277668067964629483299712322411566L,
 342             -0.29602090537771744401524080430529369136L,
 343             -0.0561832587517836908929331992218879676L,
 344             -0.0053236785487328044334381502530383140443L,
 345             -0.00018629360291358130461736386077971890789L,
 346             -0.10164985672213178500790406939467614498e-6L,
 347             0.13680157145361387405588201461036338274e-8L
 348          };
 349          static const T Q[] = {
 350             1,
 351             4.9106336261005990534095838574132225599L,
 352             10.258804800866438510889341082793078432L,
 353             11.88588976846826108836629960537466889L,
 354             8.3455000546999704314454891036700998428L,
 355             3.6428823682421746343233362007194282703L,
 356             0.97465989807254572142266753052776132252L,
 357             0.15121052897097822172763084966793352524L,
 358             0.012017363555383555123769849654484594893L,
 359             0.0003583032812720649835431669893011257277L
 360          };
 361
 362          T r = tools::evaluate_polynomial(P, zm1) / tools::evaluate_polynomial(Q, zm1);
 363          T prefix = zm1 * zm2;
 364
 365          result += prefix * Y + prefix * r;
 366          BOOST_MATH_INSTRUMENT_CODE(result);
 367       }
 368       else if(z <= 1.625)
 369       {
 370          //
 371          // Use the following form:
 372          //
 373          // lgamma(z) = (2-z)(1-z)(Y + R(2-z))
 374          //
 375          // where R(2-z) is a rational approximation optimised for
 376          // low absolute error - as long as it's absolute error
 377          // is small compared to the constant Y - then any rounding
 378          // error in it's computation will get wiped out.
 379          //
 380          // R(2-z) has the following properties:
 381          //
 382          // Max error found at 128-bit long double precision  9.634e-36
 383          // Maximum Deviation Found (approximation error)     1.538e-37
 384          // Expected Error Term (theoretical error)           2.350e-38
 385          //
 386          static const float Y = 0.483787059783935546875f;
 387
 388          static const T P[] = {
 389             -0.017977422421608624353488126610933005432L,
 390             0.18484528905298309555089509029244135703L,
 391             -0.40401251514859546989565001431430884082L,
 392             0.40277179799147356461954182877921388182L,
 393             -0.21993421441282936476709677700477598816L,
 394             0.069595742223850248095697771331107571011L,
 395             -0.012681481427699686635516772923547347328L,
 396             0.0012489322866834830413292771335113136034L,
 397             -0.57058739515423112045108068834668269608e-4L,
 398             0.8207548771933585614380644961342925976e-6L
 399          };
 400          static const T Q[] = {
 401             1,
 402             -2.9629552288944259229543137757200262073L,
 403             3.7118380799042118987185957298964772755L,
 404             -2.5569815272165399297600586376727357187L,
 405             1.0546764918220835097855665680632153367L,
 406             -0.26574021300894401276478730940980810831L,
 407             0.03996289731752081380552901986471233462L,
 408             -0.0033398680924544836817826046380586480873L,
 409             0.00013288854760548251757651556792598235735L,
 410             -0.17194794958274081373243161848194745111e-5L
 411          };
 412          T r = zm2 * zm1;
 413          T R = tools::evaluate_polynomial(P, 0.625 - zm1) / tools::evaluate_polynomial(Q, 0.625 - zm1);
 414
 415          result += r * Y + r * R;
 416          BOOST_MATH_INSTRUMENT_CODE(result);
 417       }
 418       else
 419       {
 420          //
 421          // Same form as above.
 422          //
 423          // Max error found (at 128-bit long double precision) 1.831e-35
 424          // Maximum Deviation Found (approximation error)      8.588e-36
 425          // Expected Error Term (theoretical error)            1.458e-36
 426          //
 427          static const float Y = 0.443811893463134765625f;
 428
 429          static const T P[] = {
 430             -0.021027558364667626231512090082402429494L,
 431             0.15128811104498736604523586803722368377L,
 432             -0.26249631480066246699388544451126410278L,
 433             0.21148748610533489823742352180628489742L,
 434             -0.093964130697489071999873506148104370633L,
 435             0.024292059227009051652542804957550866827L,
 436             -0.0036284453226534839926304745756906117066L,
 437             0.0002939230129315195346843036254392485984L,
 438             -0.11088589183158123733132268042570710338e-4L,
 439             0.13240510580220763969511741896361984162e-6L
 440          };
 441          static const T Q[] = {
 442             1,
 443             -2.4240003754444040525462170802796471996L,
 444             2.4868383476933178722203278602342786002L,
 445             -1.4047068395206343375520721509193698547L,
 446             0.47583809087867443858344765659065773369L,
 447             -0.09865724264554556400463655444270700132L,
 448             0.012238223514176587501074150988445109735L,
 449             -0.00084625068418239194670614419707491797097L,
 450             0.2796574430456237061420839429225710602e-4L,
 451             -0.30202973883316730694433702165188835331e-6L
 452          };
 453          // (2 - x) * (1 - x) * (c + R(2 - x))
 454          T r = zm2 * zm1;
 455          T R = tools::evaluate_polynomial(P, -zm2) / tools::evaluate_polynomial(Q, -zm2);
 456
 457          result += r * Y + r * R;
 458          BOOST_MATH_INSTRUMENT_CODE(result);
 459       }
 460    }
 461    BOOST_MATH_INSTRUMENT_CODE(result);
 462    return result;
 463 }
 464 template <class T, class Policy, class L>
 465 T lgamma_small_imp(T z, T zm1, T zm2, const mpl::int_<0>&, const Policy& pol, const L&)
 466 {
 467    //
 468    // No rational approximations are available because either
 469    // T has no numeric_limits support (so we can't tell how
 470    // many digits it has), or T has more digits than we know
 471    // what to do with.... we do have a Lanczos approximation
 472    // though, and that can be used to keep errors under control.
 473    //
 474    BOOST_MATH_STD_USING  // for ADL of std names
 475    T result = 0;
 476    if(z < tools::epsilon<T>())
 477    {
 478       result = -log(z);
 479    }
 480    else if(z < 0.5)
 481    {
 482       // taking the log of tgamma reduces the error, no danger of overflow here:
 483       result = log(gamma_imp(z, pol, L()));
 484    }
 485    else if(z >= 3)
 486    {
 487       // taking the log of tgamma reduces the error, no danger of overflow here:
 488       result = log(gamma_imp(z, pol, L()));
 489    }
 490    else if(z >= 1.5)
 491    {
 492       // special case near 2:
 493       T dz = zm2;
 494       result = dz * log((z + L::g() - T(0.5)) / boost::math::constants::e<T>());
 495       result += boost::math::log1p(dz / (L::g() + T(1.5)), pol) * T(1.5);
 496       result += boost::math::log1p(L::lanczos_sum_near_2(dz), pol);
 497    }
 498    else
 499    {
 500       // special case near 1:
 501       T dz = zm1;
 502       result = dz * log((z + L::g() - T(0.5)) / boost::math::constants::e<T>());
 503       result += boost::math::log1p(dz / (L::g() + T(0.5)), pol) / 2;
 504       result += boost::math::log1p(L::lanczos_sum_near_1(dz), pol);
 505    }
 506    return result;
 507 }
 508
 509 }}} // namespaces
 510
 511 #endif // BOOST_MATH_SPECIAL_FUNCTIONS_DETAIL_LGAMMA_SMALL
 512