special_functions/detail/lgamma_small.hpp

0001 //  (C) Copyright John Maddock 2006.
0002 //  Use, modification and distribution are subject to the
0003 //  Boost Software License, Version 1.0. (See accompanying file
0004 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0005
0006 #ifndef BOOST_MATH_SPECIAL_FUNCTIONS_DETAIL_LGAMMA_SMALL
0007 #define BOOST_MATH_SPECIAL_FUNCTIONS_DETAIL_LGAMMA_SMALL
0008
0009 #ifdef _MSC_VER
0010 #pragma once
0011 #endif
0012
0013 #include <boost/math/tools/big_constant.hpp>
0014
0015 #if defined(__GNUC__) && defined(BOOST_MATH_USE_FLOAT128)
0016 //
0017 // This is the only way we can avoid
0018 // warning: non-standard suffix on floating constant [-Wpedantic]
0019 // when building with -Wall -pedantic.  Neither __extension__
0020 // nor #pragma diagnostic ignored work :(
0021 //
0022 #pragma GCC system_header
0023 #endif
0024
0025 namespace boost{ namespace math{ namespace detail{
0026
0027 //
0028 // These need forward declaring to keep GCC happy:
0029 //
0030 template <class T, class Policy, class Lanczos>
0031 T gamma_imp(T z, const Policy& pol, const Lanczos& l);
0032 template <class T, class Policy>
0033 T gamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos& l);
0034
0035 //
0036 // lgamma for small arguments:
0037 //
0038 template <class T, class Policy, class Lanczos>
0039 T lgamma_small_imp(T z, T zm1, T zm2, const std::integral_constant<int, 64>&, const Policy& /* l */, const Lanczos&)
0040 {
0041    // This version uses rational approximations for small
0042    // values of z accurate enough for 64-bit mantissas
0043    // (80-bit long doubles), works well for 53-bit doubles as well.
0044    // Lanczos is only used to select the Lanczos function.
0045
0046    BOOST_MATH_STD_USING  // for ADL of std names
0047    T result = 0;
0048    if(z < tools::epsilon<T>())
0049    {
0050       result = -log(z);
0051    }
0052    else if((zm1 == 0) || (zm2 == 0))
0053    {
0054       // nothing to do, result is zero....
0055    }
0056    else if(z > 2)
0057    {
0058       //
0059       // Begin by performing argument reduction until
0060       // z is in [2,3):
0061       //
0062       if(z >= 3)
0063       {
0064          do
0065          {
0066             z -= 1;
0067             zm2 -= 1;
0068             result += log(z);
0069          }while(z >= 3);
0070          // Update zm2, we need it below:
0071          zm2 = z - 2;
0072       }
0073
0074       //
0075       // Use the following form:
0076       //
0077       // lgamma(z) = (z-2)(z+1)(Y + R(z-2))
0078       //
0079       // where R(z-2) is a rational approximation optimised for
0080       // low absolute error - as long as it's absolute error
0081       // is small compared to the constant Y - then any rounding
0082       // error in it's computation will get wiped out.
0083       //
0084       // R(z-2) has the following properties:
0085       //
0086       // At double: Max error found:                    4.231e-18
0087       // At long double: Max error found:               1.987e-21
0088       // Maximum Deviation Found (approximation error): 5.900e-24
0089       //
0090       static const T P[] = {
0091          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -0.180355685678449379109e-1)),
0092          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.25126649619989678683e-1)),
0093          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.494103151567532234274e-1)),
0094          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.172491608709613993966e-1)),
0095          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -0.259453563205438108893e-3)),
0096          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -0.541009869215204396339e-3)),
0097          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -0.324588649825948492091e-4))
0098       };
0099       static const T Q[] = {
0100          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.1e1)),
0101          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.196202987197795200688e1)),
0102          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.148019669424231326694e1)),
0103          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.541391432071720958364e0)),
0104          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.988504251128010129477e-1)),
0105          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.82130967464889339326e-2)),
0106          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.224936291922115757597e-3)),
0107          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -0.223352763208617092964e-6))
0108       };
0109
0110       static const float Y = 0.158963680267333984375e0f;
0111
0112       T r = zm2 * (z + 1);
0113       T R = tools::evaluate_polynomial(P, zm2);
0114       R /= tools::evaluate_polynomial(Q, zm2);
0115
0116       result +=  r * Y + r * R;
0117    }
0118    else
0119    {
0120       //
0121       // If z is less than 1 use recurrence to shift to
0122       // z in the interval [1,2]:
0123       //
0124       if(z < 1)
0125       {
0126          result += -log(z);
0127          zm2 = zm1;
0128          zm1 = z;
0129          z += 1;
0130       }
0131       //
0132       // Two approximations, on for z in [1,1.5] and
0133       // one for z in [1.5,2]:
0134       //
0135       if(z <= T(1.5))
0136       {
0137          //
0138          // Use the following form:
0139          //
0140          // lgamma(z) = (z-1)(z-2)(Y + R(z-1))
0141          //
0142          // where R(z-1) is a rational approximation optimised for
0143          // low absolute error - as long as it's absolute error
0144          // is small compared to the constant Y - then any rounding
0145          // error in it's computation will get wiped out.
0146          //
0147          // R(z-1) has the following properties:
0148          //
0149          // At double precision: Max error found:                1.230011e-17
0150          // At 80-bit long double precision:   Max error found:  5.631355e-21
0151          // Maximum Deviation Found:                             3.139e-021
0152          // Expected Error Term:                                 3.139e-021
0153
0154          //
0155          static const float Y = 0.52815341949462890625f;
0156
0157          static const T P[] = {
0158             static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.490622454069039543534e-1)),
0159             static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -0.969117530159521214579e-1)),
0160             static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -0.414983358359495381969e0)),
0161             static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -0.406567124211938417342e0)),
0162             static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -0.158413586390692192217e0)),
0163             static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -0.240149820648571559892e-1)),
0164             static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -0.100346687696279557415e-2))
0165          };
0166          static const T Q[] = {
0167             static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.1e1)),
0168             static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.302349829846463038743e1)),
0169             static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.348739585360723852576e1)),
0170             static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.191415588274426679201e1)),
0171             static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.507137738614363510846e0)),
0172             static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.577039722690451849648e-1)),
0173             static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.195768102601107189171e-2))
0174          };
0175
0176          T r = tools::evaluate_polynomial(P, zm1) / tools::evaluate_polynomial(Q, zm1);
0177          T prefix = zm1 * zm2;
0178
0179          result += prefix * Y + prefix * r;
0180       }
0181       else
0182       {
0183          //
0184          // Use the following form:
0185          //
0186          // lgamma(z) = (2-z)(1-z)(Y + R(2-z))
0187          //
0188          // where R(2-z) is a rational approximation optimised for
0189          // low absolute error - as long as it's absolute error
0190          // is small compared to the constant Y - then any rounding
0191          // error in it's computation will get wiped out.
0192          //
0193          // R(2-z) has the following properties:
0194          //
0195          // At double precision, max error found:              1.797565e-17
0196          // At 80-bit long double precision, max error found:  9.306419e-21
0197          // Maximum Deviation Found:                           2.151e-021
0198          // Expected Error Term:                               2.150e-021
0199          //
0200          static const float Y = 0.452017307281494140625f;
0201
0202          static const T P[] = {
0203             static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -0.292329721830270012337e-1)),
0204             static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.144216267757192309184e0)),
0205             static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -0.142440390738631274135e0)),
0206             static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.542809694055053558157e-1)),
0207             static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -0.850535976868336437746e-2)),
0208             static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.431171342679297331241e-3))
0209          };
0210          static const T Q[] = {
0211             static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.1e1)),
0212             static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -0.150169356054485044494e1)),
0213             static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.846973248876495016101e0)),
0214             static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -0.220095151814995745555e0)),
0215             static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.25582797155975869989e-1)),
0216             static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -0.100666795539143372762e-2)),
0217             static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -0.827193521891290553639e-6))
0218          };
0219          T r = zm2 * zm1;
0220          T R = tools::evaluate_polynomial(P, T(-zm2)) / tools::evaluate_polynomial(Q, T(-zm2));
0221
0222          result += r * Y + r * R;
0223       }
0224    }
0225    return result;
0226 }
0227 template <class T, class Policy, class Lanczos>
0228 T lgamma_small_imp(T z, T zm1, T zm2, const std::integral_constant<int, 113>&, const Policy& /* l */, const Lanczos&)
0229 {
0230    //
0231    // This version uses rational approximations for small
0232    // values of z accurate enough for 113-bit mantissas
0233    // (128-bit long doubles).
0234    //
0235    BOOST_MATH_STD_USING  // for ADL of std names
0236    T result = 0;
0237    if(z < tools::epsilon<T>())
0238    {
0239       result = -log(z);
0240       BOOST_MATH_INSTRUMENT_CODE(result);
0241    }
0242    else if((zm1 == 0) || (zm2 == 0))
0243    {
0244       // nothing to do, result is zero....
0245    }
0246    else if(z > 2)
0247    {
0248       //
0249       // Begin by performing argument reduction until
0250       // z is in [2,3):
0251       //
0252       if(z >= 3)
0253       {
0254          do
0255          {
0256             z -= 1;
0257             result += log(z);
0258          }while(z >= 3);
0259          zm2 = z - 2;
0260       }
0261       BOOST_MATH_INSTRUMENT_CODE(zm2);
0262       BOOST_MATH_INSTRUMENT_CODE(z);
0263       BOOST_MATH_INSTRUMENT_CODE(result);
0264
0265       //
0266       // Use the following form:
0267       //
0268       // lgamma(z) = (z-2)(z+1)(Y + R(z-2))
0269       //
0270       // where R(z-2) is a rational approximation optimised for
0271       // low absolute error - as long as it's absolute error
0272       // is small compared to the constant Y - then any rounding
0273       // error in it's computation will get wiped out.
0274       //
0275       // Maximum Deviation Found (approximation error)      3.73e-37
0276
0277       static const T P[] = {
0278          BOOST_MATH_BIG_CONSTANT(T, 113, -0.018035568567844937910504030027467476655),
0279          BOOST_MATH_BIG_CONSTANT(T, 113, 0.013841458273109517271750705401202404195),
0280          BOOST_MATH_BIG_CONSTANT(T, 113, 0.062031842739486600078866923383017722399),
0281          BOOST_MATH_BIG_CONSTANT(T, 113, 0.052518418329052161202007865149435256093),
0282          BOOST_MATH_BIG_CONSTANT(T, 113, 0.01881718142472784129191838493267755758),
0283          BOOST_MATH_BIG_CONSTANT(T, 113, 0.0025104830367021839316463675028524702846),
0284          BOOST_MATH_BIG_CONSTANT(T, 113, -0.00021043176101831873281848891452678568311),
0285          BOOST_MATH_BIG_CONSTANT(T, 113, -0.00010249622350908722793327719494037981166),
0286          BOOST_MATH_BIG_CONSTANT(T, 113, -0.11381479670982006841716879074288176994e-4),
0287          BOOST_MATH_BIG_CONSTANT(T, 113, -0.49999811718089980992888533630523892389e-6),
0288          BOOST_MATH_BIG_CONSTANT(T, 113, -0.70529798686542184668416911331718963364e-8)
0289       };
0290       static const T Q[] = {
0291          BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),
0292          BOOST_MATH_BIG_CONSTANT(T, 113, 2.5877485070422317542808137697939233685),
0293          BOOST_MATH_BIG_CONSTANT(T, 113, 2.8797959228352591788629602533153837126),
0294          BOOST_MATH_BIG_CONSTANT(T, 113, 1.8030885955284082026405495275461180977),
0295          BOOST_MATH_BIG_CONSTANT(T, 113, 0.69774331297747390169238306148355428436),
0296          BOOST_MATH_BIG_CONSTANT(T, 113, 0.17261566063277623942044077039756583802),
0297          BOOST_MATH_BIG_CONSTANT(T, 113, 0.02729301254544230229429621192443000121),
0298          BOOST_MATH_BIG_CONSTANT(T, 113, 0.0026776425891195270663133581960016620433),
0299          BOOST_MATH_BIG_CONSTANT(T, 113, 0.00015244249160486584591370355730402168106),
0300          BOOST_MATH_BIG_CONSTANT(T, 113, 0.43997034032479866020546814475414346627e-5),
0301          BOOST_MATH_BIG_CONSTANT(T, 113, 0.46295080708455613044541885534408170934e-7),
0302          BOOST_MATH_BIG_CONSTANT(T, 113, -0.93326638207459533682980757982834180952e-11),
0303          BOOST_MATH_BIG_CONSTANT(T, 113, 0.42316456553164995177177407325292867513e-13)
0304       };
0305
0306       T R = tools::evaluate_polynomial(P, zm2);
0307       R /= tools::evaluate_polynomial(Q, zm2);
0308
0309       static const float Y = 0.158963680267333984375F;
0310
0311       T r = zm2 * (z + 1);
0312
0313       result +=  r * Y + r * R;
0314       BOOST_MATH_INSTRUMENT_CODE(result);
0315    }
0316    else
0317    {
0318       //
0319       // If z is less than 1 use recurrence to shift to
0320       // z in the interval [1,2]:
0321       //
0322       if(z < 1)
0323       {
0324          result += -log(z);
0325          zm2 = zm1;
0326          zm1 = z;
0327          z += 1;
0328       }
0329       BOOST_MATH_INSTRUMENT_CODE(result);
0330       BOOST_MATH_INSTRUMENT_CODE(z);
0331       BOOST_MATH_INSTRUMENT_CODE(zm2);
0332       //
0333       // Three approximations, on for z in [1,1.35], [1.35,1.625] and [1.625,1]
0334       //
0335       if(z <= 1.35)
0336       {
0337          //
0338          // Use the following form:
0339          //
0340          // lgamma(z) = (z-1)(z-2)(Y + R(z-1))
0341          //
0342          // where R(z-1) is a rational approximation optimised for
0343          // low absolute error - as long as it's absolute error
0344          // is small compared to the constant Y - then any rounding
0345          // error in it's computation will get wiped out.
0346          //
0347          // R(z-1) has the following properties:
0348          //
0349          // Maximum Deviation Found (approximation error)            1.659e-36
0350          // Expected Error Term (theoretical error)                  1.343e-36
0351          // Max error found at 128-bit long double precision         1.007e-35
0352          //
0353          static const float Y = 0.54076099395751953125f;
0354
0355          static const T P[] = {
0356             BOOST_MATH_BIG_CONSTANT(T, 113, 0.036454670944013329356512090082402429697),
0357             BOOST_MATH_BIG_CONSTANT(T, 113, -0.066235835556476033710068679907798799959),
0358             BOOST_MATH_BIG_CONSTANT(T, 113, -0.67492399795577182387312206593595565371),
0359             BOOST_MATH_BIG_CONSTANT(T, 113, -1.4345555263962411429855341651960000166),
0360             BOOST_MATH_BIG_CONSTANT(T, 113, -1.4894319559821365820516771951249649563),
0361             BOOST_MATH_BIG_CONSTANT(T, 113, -0.87210277668067964629483299712322411566),
0362             BOOST_MATH_BIG_CONSTANT(T, 113, -0.29602090537771744401524080430529369136),
0363             BOOST_MATH_BIG_CONSTANT(T, 113, -0.0561832587517836908929331992218879676),
0364             BOOST_MATH_BIG_CONSTANT(T, 113, -0.0053236785487328044334381502530383140443),
0365             BOOST_MATH_BIG_CONSTANT(T, 113, -0.00018629360291358130461736386077971890789),
0366             BOOST_MATH_BIG_CONSTANT(T, 113, -0.10164985672213178500790406939467614498e-6),
0367             BOOST_MATH_BIG_CONSTANT(T, 113, 0.13680157145361387405588201461036338274e-8)
0368          };
0369          static const T Q[] = {
0370             BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),
0371             BOOST_MATH_BIG_CONSTANT(T, 113, 4.9106336261005990534095838574132225599),
0372             BOOST_MATH_BIG_CONSTANT(T, 113, 10.258804800866438510889341082793078432),
0373             BOOST_MATH_BIG_CONSTANT(T, 113, 11.88588976846826108836629960537466889),
0374             BOOST_MATH_BIG_CONSTANT(T, 113, 8.3455000546999704314454891036700998428),
0375             BOOST_MATH_BIG_CONSTANT(T, 113, 3.6428823682421746343233362007194282703),
0376             BOOST_MATH_BIG_CONSTANT(T, 113, 0.97465989807254572142266753052776132252),
0377             BOOST_MATH_BIG_CONSTANT(T, 113, 0.15121052897097822172763084966793352524),
0378             BOOST_MATH_BIG_CONSTANT(T, 113, 0.012017363555383555123769849654484594893),
0379             BOOST_MATH_BIG_CONSTANT(T, 113, 0.0003583032812720649835431669893011257277)
0380          };
0381
0382          T r = tools::evaluate_polynomial(P, zm1) / tools::evaluate_polynomial(Q, zm1);
0383          T prefix = zm1 * zm2;
0384
0385          result += prefix * Y + prefix * r;
0386          BOOST_MATH_INSTRUMENT_CODE(result);
0387       }
0388       else if(z <= 1.625)
0389       {
0390          //
0391          // Use the following form:
0392          //
0393          // lgamma(z) = (2-z)(1-z)(Y + R(2-z))
0394          //
0395          // where R(2-z) is a rational approximation optimised for
0396          // low absolute error - as long as it's absolute error
0397          // is small compared to the constant Y - then any rounding
0398          // error in it's computation will get wiped out.
0399          //
0400          // R(2-z) has the following properties:
0401          //
0402          // Max error found at 128-bit long double precision  9.634e-36
0403          // Maximum Deviation Found (approximation error)     1.538e-37
0404          // Expected Error Term (theoretical error)           2.350e-38
0405          //
0406          static const float Y = 0.483787059783935546875f;
0407
0408          static const T P[] = {
0409             BOOST_MATH_BIG_CONSTANT(T, 113, -0.017977422421608624353488126610933005432),
0410             BOOST_MATH_BIG_CONSTANT(T, 113, 0.18484528905298309555089509029244135703),
0411             BOOST_MATH_BIG_CONSTANT(T, 113, -0.40401251514859546989565001431430884082),
0412             BOOST_MATH_BIG_CONSTANT(T, 113, 0.40277179799147356461954182877921388182),
0413             BOOST_MATH_BIG_CONSTANT(T, 113, -0.21993421441282936476709677700477598816),
0414             BOOST_MATH_BIG_CONSTANT(T, 113, 0.069595742223850248095697771331107571011),
0415             BOOST_MATH_BIG_CONSTANT(T, 113, -0.012681481427699686635516772923547347328),
0416             BOOST_MATH_BIG_CONSTANT(T, 113, 0.0012489322866834830413292771335113136034),
0417             BOOST_MATH_BIG_CONSTANT(T, 113, -0.57058739515423112045108068834668269608e-4),
0418             BOOST_MATH_BIG_CONSTANT(T, 113, 0.8207548771933585614380644961342925976e-6)
0419          };
0420          static const T Q[] = {
0421             BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),
0422             BOOST_MATH_BIG_CONSTANT(T, 113, -2.9629552288944259229543137757200262073),
0423             BOOST_MATH_BIG_CONSTANT(T, 113, 3.7118380799042118987185957298964772755),
0424             BOOST_MATH_BIG_CONSTANT(T, 113, -2.5569815272165399297600586376727357187),
0425             BOOST_MATH_BIG_CONSTANT(T, 113, 1.0546764918220835097855665680632153367),
0426             BOOST_MATH_BIG_CONSTANT(T, 113, -0.26574021300894401276478730940980810831),
0427             BOOST_MATH_BIG_CONSTANT(T, 113, 0.03996289731752081380552901986471233462),
0428             BOOST_MATH_BIG_CONSTANT(T, 113, -0.0033398680924544836817826046380586480873),
0429             BOOST_MATH_BIG_CONSTANT(T, 113, 0.00013288854760548251757651556792598235735),
0430             BOOST_MATH_BIG_CONSTANT(T, 113, -0.17194794958274081373243161848194745111e-5)
0431          };
0432          T r = zm2 * zm1;
0433          T R = tools::evaluate_polynomial(P, T(0.625 - zm1)) / tools::evaluate_polynomial(Q, T(0.625 - zm1));
0434
0435          result += r * Y + r * R;
0436          BOOST_MATH_INSTRUMENT_CODE(result);
0437       }
0438       else
0439       {
0440          //
0441          // Same form as above.
0442          //
0443          // Max error found (at 128-bit long double precision) 1.831e-35
0444          // Maximum Deviation Found (approximation error)      8.588e-36
0445          // Expected Error Term (theoretical error)            1.458e-36
0446          //
0447          static const float Y = 0.443811893463134765625f;
0448
0449          static const T P[] = {
0450             BOOST_MATH_BIG_CONSTANT(T, 113, -0.021027558364667626231512090082402429494),
0451             BOOST_MATH_BIG_CONSTANT(T, 113, 0.15128811104498736604523586803722368377),
0452             BOOST_MATH_BIG_CONSTANT(T, 113, -0.26249631480066246699388544451126410278),
0453             BOOST_MATH_BIG_CONSTANT(T, 113, 0.21148748610533489823742352180628489742),
0454             BOOST_MATH_BIG_CONSTANT(T, 113, -0.093964130697489071999873506148104370633),
0455             BOOST_MATH_BIG_CONSTANT(T, 113, 0.024292059227009051652542804957550866827),
0456             BOOST_MATH_BIG_CONSTANT(T, 113, -0.0036284453226534839926304745756906117066),
0457             BOOST_MATH_BIG_CONSTANT(T, 113, 0.0002939230129315195346843036254392485984),
0458             BOOST_MATH_BIG_CONSTANT(T, 113, -0.11088589183158123733132268042570710338e-4),
0459             BOOST_MATH_BIG_CONSTANT(T, 113, 0.13240510580220763969511741896361984162e-6)
0460          };
0461          static const T Q[] = {
0462             BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),
0463             BOOST_MATH_BIG_CONSTANT(T, 113, -2.4240003754444040525462170802796471996),
0464             BOOST_MATH_BIG_CONSTANT(T, 113, 2.4868383476933178722203278602342786002),
0465             BOOST_MATH_BIG_CONSTANT(T, 113, -1.4047068395206343375520721509193698547),
0466             BOOST_MATH_BIG_CONSTANT(T, 113, 0.47583809087867443858344765659065773369),
0467             BOOST_MATH_BIG_CONSTANT(T, 113, -0.09865724264554556400463655444270700132),
0468             BOOST_MATH_BIG_CONSTANT(T, 113, 0.012238223514176587501074150988445109735),
0469             BOOST_MATH_BIG_CONSTANT(T, 113, -0.00084625068418239194670614419707491797097),
0470             BOOST_MATH_BIG_CONSTANT(T, 113, 0.2796574430456237061420839429225710602e-4),
0471             BOOST_MATH_BIG_CONSTANT(T, 113, -0.30202973883316730694433702165188835331e-6)
0472          };
0473          // (2 - x) * (1 - x) * (c + R(2 - x))
0474          T r = zm2 * zm1;
0475          T R = tools::evaluate_polynomial(P, T(-zm2)) / tools::evaluate_polynomial(Q, T(-zm2));
0476
0477          result += r * Y + r * R;
0478          BOOST_MATH_INSTRUMENT_CODE(result);
0479       }
0480    }
0481    BOOST_MATH_INSTRUMENT_CODE(result);
0482    return result;
0483 }
0484 template <class T, class Policy, class Lanczos>
0485 T lgamma_small_imp(T z, T zm1, T zm2, const std::integral_constant<int, 0>&, const Policy& pol, const Lanczos& l)
0486 {
0487    //
0488    // No rational approximations are available because either
0489    // T has no numeric_limits support (so we can't tell how
0490    // many digits it has), or T has more digits than we know
0491    // what to do with.... we do have a Lanczos approximation
0492    // though, and that can be used to keep errors under control.
0493    //
0494    BOOST_MATH_STD_USING  // for ADL of std names
0495    T result = 0;
0496    if(z < tools::epsilon<T>())
0497    {
0498       result = -log(z);
0499    }
0500    else if(z < 0.5)
0501    {
0502       // taking the log of tgamma reduces the error, no danger of overflow here:
0503       result = log(gamma_imp(z, pol, Lanczos()));
0504    }
0505    else if(z >= 3)
0506    {
0507       // taking the log of tgamma reduces the error, no danger of overflow here:
0508       result = log(gamma_imp(z, pol, Lanczos()));
0509    }
0510    else if(z >= 1.5)
0511    {
0512       // special case near 2:
0513       T dz = zm2;
0514       result = dz * log((z + lanczos_g_near_1_and_2(l) - T(0.5)) / boost::math::constants::e<T>());
0515       result += boost::math::log1p(dz / (lanczos_g_near_1_and_2(l) + T(1.5)), pol) * T(1.5);
0516       result += boost::math::log1p(Lanczos::lanczos_sum_near_2(dz), pol);
0517    }
0518    else
0519    {
0520       // special case near 1:
0521       T dz = zm1;
0522       result = dz * log((z + lanczos_g_near_1_and_2(l) - T(0.5)) / boost::math::constants::e<T>());
0523       result += boost::math::log1p(dz / (lanczos_g_near_1_and_2(l) + T(0.5)), pol) / 2;
0524       result += boost::math::log1p(Lanczos::lanczos_sum_near_1(dz), pol);
0525    }
0526    return result;
0527 }
0528
0529 }}} // namespaces
0530
0531 #endif // BOOST_MATH_SPECIAL_FUNCTIONS_DETAIL_LGAMMA_SMALL
0532