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