special_functions/detail/lgamma_small.hpp

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