math/special_functions/trigamma.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_SF_TRIGAMMA_HPP
0008 #define BOOST_MATH_SF_TRIGAMMA_HPP
0009
0010 #ifdef _MSC_VER
0011 #pragma once
0012 #endif
0013
0014 #include <boost/math/tools/config.hpp>
0015 #include <boost/math/tools/rational.hpp>
0016 #include <boost/math/tools/promotion.hpp>
0017 #include <boost/math/tools/big_constant.hpp>
0018 #include <boost/math/tools/type_traits.hpp>
0019 #include <boost/math/policies/policy.hpp>
0020 #include <boost/math/policies/error_handling.hpp>
0021 #include <boost/math/constants/constants.hpp>
0022 #include <boost/math/special_functions/sin_pi.hpp>
0023 #include <boost/math/special_functions/pow.hpp>
0024
0025 #ifndef BOOST_MATH_HAS_NVRTC
0026 #include <boost/math/special_functions/math_fwd.hpp>
0027 #include <boost/math/special_functions/polygamma.hpp>
0028 #include <boost/math/tools/series.hpp>
0029 #endif
0030
0031 #if defined(__GNUC__) && defined(BOOST_MATH_USE_FLOAT128)
0032 //
0033 // This is the only way we can avoid
0034 // warning: non-standard suffix on floating constant [-Wpedantic]
0035 // when building with -Wall -pedantic.  Neither __extension__
0036 // nor #pragma diagnostic ignored work :(
0037 //
0038 #pragma GCC system_header
0039 #endif
0040
0041 namespace boost{
0042 namespace math{
0043 namespace detail{
0044
0045 // TODO(mborland): Temporary for NVRTC
0046 #ifndef BOOST_MATH_HAS_NVRTC
0047 template<class T, class Policy>
0048 T polygamma_imp(const int n, T x, const Policy &pol);
0049
0050 template <class T, class Policy>
0051 T trigamma_prec(T x, const Policy& pol, const boost::math::integral_constant<int, 0>&)
0052 {
0053    return polygamma_imp(1, x, pol);
0054 }
0055 #endif
0056
0057 template <class T, class Policy>
0058 BOOST_MATH_GPU_ENABLED T trigamma_prec(T x, const Policy&, const boost::math::integral_constant<int, 53>&)
0059 {
0060    // Max error in interpolated form: 3.736e-017
0061    BOOST_MATH_STATIC const T offset = BOOST_MATH_BIG_CONSTANT(T, 53, 2.1093254089355469);
0062    BOOST_MATH_STATIC const T P_1_2[] = {
0063       BOOST_MATH_BIG_CONSTANT(T, 53, -1.1093280605946045),
0064       BOOST_MATH_BIG_CONSTANT(T, 53, -3.8310674472619321),
0065       BOOST_MATH_BIG_CONSTANT(T, 53, -3.3703848401898283),
0066       BOOST_MATH_BIG_CONSTANT(T, 53, 0.28080574467981213),
0067       BOOST_MATH_BIG_CONSTANT(T, 53, 1.6638069578676164),
0068       BOOST_MATH_BIG_CONSTANT(T, 53, 0.64468386819102836),
0069    };
0070    BOOST_MATH_STATIC const T Q_1_2[] = {
0071       BOOST_MATH_BIG_CONSTANT(T, 53, 1.0),
0072       BOOST_MATH_BIG_CONSTANT(T, 53, 3.4535389668541151),
0073       BOOST_MATH_BIG_CONSTANT(T, 53, 4.5208926987851437),
0074       BOOST_MATH_BIG_CONSTANT(T, 53, 2.7012734178351534),
0075       BOOST_MATH_BIG_CONSTANT(T, 53, 0.64468798399785611),
0076       BOOST_MATH_BIG_CONSTANT(T, 53, -0.20314516859987728e-6),
0077    };
0078    // Max error in interpolated form: 1.159e-017
0079    BOOST_MATH_STATIC const T P_2_4[] = {
0080       BOOST_MATH_BIG_CONSTANT(T, 53, -0.13803835004508849e-7),
0081       BOOST_MATH_BIG_CONSTANT(T, 53, 0.50000049158540261),
0082       BOOST_MATH_BIG_CONSTANT(T, 53, 1.6077979838469348),
0083       BOOST_MATH_BIG_CONSTANT(T, 53, 2.5645435828098254),
0084       BOOST_MATH_BIG_CONSTANT(T, 53, 2.0534873203680393),
0085       BOOST_MATH_BIG_CONSTANT(T, 53, 0.74566981111565923),
0086    };
0087    BOOST_MATH_STATIC const T Q_2_4[] = {
0088       BOOST_MATH_BIG_CONSTANT(T, 53, 1.0),
0089       BOOST_MATH_BIG_CONSTANT(T, 53, 2.8822787662376169),
0090       BOOST_MATH_BIG_CONSTANT(T, 53, 4.1681660554090917),
0091       BOOST_MATH_BIG_CONSTANT(T, 53, 2.7853527819234466),
0092       BOOST_MATH_BIG_CONSTANT(T, 53, 0.74967671848044792),
0093       BOOST_MATH_BIG_CONSTANT(T, 53, -0.00057069112416246805),
0094    };
0095    // Maximum Deviation Found:                     6.896e-018
0096    // Expected Error Term :                       -6.895e-018
0097    // Maximum Relative Change in Control Points :  8.497e-004
0098    BOOST_MATH_STATIC const T P_4_inf[] = {
0099       static_cast<T>(0.68947581948701249e-17L),
0100       static_cast<T>(0.49999999999998975L),
0101       static_cast<T>(1.0177274392923795L),
0102       static_cast<T>(2.498208511343429L),
0103       static_cast<T>(2.1921221359427595L),
0104       static_cast<T>(1.5897035272532764L),
0105       static_cast<T>(0.40154388356961734L),
0106    };
0107    BOOST_MATH_STATIC const T Q_4_inf[] = {
0108       static_cast<T>(1.0L),
0109       static_cast<T>(1.7021215452463932L),
0110       static_cast<T>(4.4290431747556469L),
0111       static_cast<T>(2.9745631894384922L),
0112       static_cast<T>(2.3013614809773616L),
0113       static_cast<T>(0.28360399799075752L),
0114       static_cast<T>(0.022892987908906897L),
0115    };
0116
0117    if(x <= 2)
0118    {
0119       return (offset + boost::math::tools::evaluate_polynomial(P_1_2, x) / tools::evaluate_polynomial(Q_1_2, x)) / (x * x);
0120    }
0121    else if(x <= 4)
0122    {
0123       T y = 1 / x;
0124       return (1 + tools::evaluate_polynomial(P_2_4, y) / tools::evaluate_polynomial(Q_2_4, y)) / x;
0125    }
0126    T y = 1 / x;
0127    return (1 + tools::evaluate_polynomial(P_4_inf, y) / tools::evaluate_polynomial(Q_4_inf, y)) / x;
0128 }
0129
0130 template <class T, class Policy>
0131 BOOST_MATH_GPU_ENABLED T trigamma_prec(T x, const Policy&, const boost::math::integral_constant<int, 64>&)
0132 {
0133    // Max error in interpolated form: 1.178e-020
0134    BOOST_MATH_STATIC const T offset_1_2 = BOOST_MATH_BIG_CONSTANT(T, 64, 2.109325408935546875);
0135    BOOST_MATH_STATIC const T P_1_2[] = {
0136       BOOST_MATH_BIG_CONSTANT(T, 64, -1.10932535608960258341),
0137       BOOST_MATH_BIG_CONSTANT(T, 64, -4.18793841543017129052),
0138       BOOST_MATH_BIG_CONSTANT(T, 64, -4.63865531898487734531),
0139       BOOST_MATH_BIG_CONSTANT(T, 64, -0.919832884430500908047),
0140       BOOST_MATH_BIG_CONSTANT(T, 64, 1.68074038333180423012),
0141       BOOST_MATH_BIG_CONSTANT(T, 64, 1.21172611429185622377),
0142       BOOST_MATH_BIG_CONSTANT(T, 64, 0.259635673503366427284),
0143    };
0144    BOOST_MATH_STATIC const T Q_1_2[] = {
0145       BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
0146       BOOST_MATH_BIG_CONSTANT(T, 64, 3.77521119359546982995),
0147       BOOST_MATH_BIG_CONSTANT(T, 64, 5.664338024578956321),
0148       BOOST_MATH_BIG_CONSTANT(T, 64, 4.25995134879278028361),
0149       BOOST_MATH_BIG_CONSTANT(T, 64, 1.62956638448940402182),
0150       BOOST_MATH_BIG_CONSTANT(T, 64, 0.259635512844691089868),
0151       BOOST_MATH_BIG_CONSTANT(T, 64, 0.629642219810618032207e-8),
0152    };
0153    // Max error in interpolated form: 3.912e-020
0154    BOOST_MATH_STATIC const T P_2_8[] = {
0155       BOOST_MATH_BIG_CONSTANT(T, 64, -0.387540035162952880976e-11),
0156       BOOST_MATH_BIG_CONSTANT(T, 64, 0.500000000276430504),
0157       BOOST_MATH_BIG_CONSTANT(T, 64, 3.21926880986360957306),
0158       BOOST_MATH_BIG_CONSTANT(T, 64, 10.2550347708483445775),
0159       BOOST_MATH_BIG_CONSTANT(T, 64, 18.9002075150709144043),
0160       BOOST_MATH_BIG_CONSTANT(T, 64, 21.0357215832399705625),
0161       BOOST_MATH_BIG_CONSTANT(T, 64, 13.4346512182925923978),
0162       BOOST_MATH_BIG_CONSTANT(T, 64, 3.98656291026448279118),
0163    };
0164    BOOST_MATH_STATIC const T Q_2_8[] = {
0165       BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
0166       BOOST_MATH_BIG_CONSTANT(T, 64, 6.10520430478613667724),
0167       BOOST_MATH_BIG_CONSTANT(T, 64, 18.475001060603645512),
0168       BOOST_MATH_BIG_CONSTANT(T, 64, 31.7087534567758405638),
0169       BOOST_MATH_BIG_CONSTANT(T, 64, 31.908814523890465398),
0170       BOOST_MATH_BIG_CONSTANT(T, 64, 17.4175479039227084798),
0171       BOOST_MATH_BIG_CONSTANT(T, 64, 3.98749106958394941276),
0172       BOOST_MATH_BIG_CONSTANT(T, 64, -0.000115917322224411128566),
0173    };
0174    // Maximum Deviation Found:                     2.635e-020
0175    // Expected Error Term :                        2.635e-020
0176    // Maximum Relative Change in Control Points :  1.791e-003
0177    BOOST_MATH_STATIC const T P_8_inf[] = {
0178       BOOST_MATH_BIG_CONSTANT(T, 64, -0.263527875092466899848e-19),
0179       BOOST_MATH_BIG_CONSTANT(T, 64, 0.500000000000000058145),
0180       BOOST_MATH_BIG_CONSTANT(T, 64, 0.0730121433777364138677),
0181       BOOST_MATH_BIG_CONSTANT(T, 64, 1.94505878379957149534),
0182       BOOST_MATH_BIG_CONSTANT(T, 64, 0.0517092358874932620529),
0183       BOOST_MATH_BIG_CONSTANT(T, 64, 1.07995383547483921121),
0184    };
0185    BOOST_MATH_STATIC const T Q_8_inf[] = {
0186       BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
0187       BOOST_MATH_BIG_CONSTANT(T, 64, -0.187309046577818095504),
0188       BOOST_MATH_BIG_CONSTANT(T, 64, 3.95255391645238842975),
0189       BOOST_MATH_BIG_CONSTANT(T, 64, -1.14743283327078949087),
0190       BOOST_MATH_BIG_CONSTANT(T, 64, 2.52989799376344914499),
0191       BOOST_MATH_BIG_CONSTANT(T, 64, -0.627414303172402506396),
0192       BOOST_MATH_BIG_CONSTANT(T, 64, 0.141554248216425512536),
0193    };
0194
0195    if(x <= 2)
0196    {
0197       return (offset_1_2 + boost::math::tools::evaluate_polynomial(P_1_2, x) / tools::evaluate_polynomial(Q_1_2, x)) / (x * x);
0198    }
0199    else if(x <= 8)
0200    {
0201       T y = 1 / x;
0202       return (1 + tools::evaluate_polynomial(P_2_8, y) / tools::evaluate_polynomial(Q_2_8, y)) / x;
0203    }
0204    T y = 1 / x;
0205    return (1 + tools::evaluate_polynomial(P_8_inf, y) / tools::evaluate_polynomial(Q_8_inf, y)) / x;
0206 }
0207
0208 template <class T, class Policy>
0209 BOOST_MATH_GPU_ENABLED T trigamma_prec(T x, const Policy&, const boost::math::integral_constant<int, 113>&)
0210 {
0211    // Max error in interpolated form: 1.916e-035
0212
0213    static const T P_1_2[] = {
0214       BOOST_MATH_BIG_CONSTANT(T, 113, -0.999999999999999082554457936871832533),
0215       BOOST_MATH_BIG_CONSTANT(T, 113, -4.71237311120865266379041700054847734),
0216       BOOST_MATH_BIG_CONSTANT(T, 113, -7.94125711970499027763789342500817316),
0217       BOOST_MATH_BIG_CONSTANT(T, 113, -5.74657746697664735258222071695644535),
0218       BOOST_MATH_BIG_CONSTANT(T, 113, -0.404213349456398905981223965160595687),
0219       BOOST_MATH_BIG_CONSTANT(T, 113, 2.47877781178642876561595890095758896),
0220       BOOST_MATH_BIG_CONSTANT(T, 113, 2.07714151702455125992166949812126433),
0221       BOOST_MATH_BIG_CONSTANT(T, 113, 0.858877899162360138844032265418028567),
0222       BOOST_MATH_BIG_CONSTANT(T, 113, 0.20499222604410032375789018837922397),
0223       BOOST_MATH_BIG_CONSTANT(T, 113, 0.0272103140348194747360175268778415049),
0224       BOOST_MATH_BIG_CONSTANT(T, 113, 0.0015764849020876949848954081173520686),
0225    };
0226    static const T Q_1_2[] = {
0227       BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),
0228       BOOST_MATH_BIG_CONSTANT(T, 113, 4.71237311120863419878375031457715223),
0229       BOOST_MATH_BIG_CONSTANT(T, 113, 9.58619118655339853449127952145877467),
0230       BOOST_MATH_BIG_CONSTANT(T, 113, 11.0940067269829372437561421279054968),
0231       BOOST_MATH_BIG_CONSTANT(T, 113, 8.09075424749327792073276309969037885),
0232       BOOST_MATH_BIG_CONSTANT(T, 113, 3.87705890159891405185343806884451286),
0233       BOOST_MATH_BIG_CONSTANT(T, 113, 1.22758678701914477836330837816976782),
0234       BOOST_MATH_BIG_CONSTANT(T, 113, 0.249092040606385004109672077814668716),
0235       BOOST_MATH_BIG_CONSTANT(T, 113, 0.0295750413900655597027079600025569048),
0236       BOOST_MATH_BIG_CONSTANT(T, 113, 0.00157648490200498142247694709728858139),
0237       BOOST_MATH_BIG_CONSTANT(T, 113, 0.161264050344059471721062360645432809e-14),
0238    };
0239
0240    // Max error in interpolated form: 8.958e-035
0241    static const T P_2_4[] = {
0242       BOOST_MATH_BIG_CONSTANT(T, 113, -2.55843734739907925764326773972215085),
0243       BOOST_MATH_BIG_CONSTANT(T, 113, -12.2830208240542011967952466273455887),
0244       BOOST_MATH_BIG_CONSTANT(T, 113, -23.9195022162767993526575786066414403),
0245       BOOST_MATH_BIG_CONSTANT(T, 113, -24.9256431504823483094158828285470862),
0246       BOOST_MATH_BIG_CONSTANT(T, 113, -14.7979122765478779075108064826412285),
0247       BOOST_MATH_BIG_CONSTANT(T, 113, -4.46654453928610666393276765059122272),
0248       BOOST_MATH_BIG_CONSTANT(T, 113, -0.0191439033405649675717082465687845002),
0249       BOOST_MATH_BIG_CONSTANT(T, 113, 0.515412052554351265708917209749037352),
0250       BOOST_MATH_BIG_CONSTANT(T, 113, 0.195378348786064304378247325360320038),
0251       BOOST_MATH_BIG_CONSTANT(T, 113, 0.0334761282624174313035014426794245393),
0252       BOOST_MATH_BIG_CONSTANT(T, 113, 0.002373665205942206348500250056602687),
0253    };
0254    static const T Q_2_4[] = {
0255       BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),
0256       BOOST_MATH_BIG_CONSTANT(T, 113, 4.80098558454419907830670928248659245),
0257       BOOST_MATH_BIG_CONSTANT(T, 113, 9.99220727843170133895059300223445265),
0258       BOOST_MATH_BIG_CONSTANT(T, 113, 11.8896146167631330735386697123464976),
0259       BOOST_MATH_BIG_CONSTANT(T, 113, 8.96613256683809091593793565879092581),
0260       BOOST_MATH_BIG_CONSTANT(T, 113, 4.47254136149624110878909334574485751),
0261       BOOST_MATH_BIG_CONSTANT(T, 113, 1.48600982028196527372434773913633152),
0262       BOOST_MATH_BIG_CONSTANT(T, 113, 0.319570735766764237068541501137990078),
0263       BOOST_MATH_BIG_CONSTANT(T, 113, 0.0407358345787680953107374215319322066),
0264       BOOST_MATH_BIG_CONSTANT(T, 113, 0.00237366520593271641375755486420859837),
0265       BOOST_MATH_BIG_CONSTANT(T, 113, 0.239554887903526152679337256236302116e-15),
0266       BOOST_MATH_BIG_CONSTANT(T, 113, -0.294749244740618656265237072002026314e-17),
0267    };
0268
0269    static const T y_offset_2_4 = BOOST_MATH_BIG_CONSTANT(T, 113, 3.558437347412109375);
0270
0271    // Max error in interpolated form: 4.319e-035
0272    static const T P_4_8[] = {
0273       BOOST_MATH_BIG_CONSTANT(T, 113, 0.166626112697021464248967707021688845e-16),
0274       BOOST_MATH_BIG_CONSTANT(T, 113, 0.499999999999997739552090249208808197),
0275       BOOST_MATH_BIG_CONSTANT(T, 113, 6.40270945019053817915772473771553187),
0276       BOOST_MATH_BIG_CONSTANT(T, 113, 41.3833374155000608013677627389343329),
0277       BOOST_MATH_BIG_CONSTANT(T, 113, 166.803341854562809335667241074035245),
0278       BOOST_MATH_BIG_CONSTANT(T, 113, 453.39964786925369319960722793414521),
0279       BOOST_MATH_BIG_CONSTANT(T, 113, 851.153712317697055375935433362983944),
0280       BOOST_MATH_BIG_CONSTANT(T, 113, 1097.70657567285059133109286478004458),
0281       BOOST_MATH_BIG_CONSTANT(T, 113, 938.431232478455316020076349367632922),
0282       BOOST_MATH_BIG_CONSTANT(T, 113, 487.268001604651932322080970189930074),
0283       BOOST_MATH_BIG_CONSTANT(T, 113, 119.953445242335730062471193124820659),
0284    };
0285    static const T Q_4_8[] = {
0286       BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),
0287       BOOST_MATH_BIG_CONSTANT(T, 113, 12.4720855670474488978638945855932398),
0288       BOOST_MATH_BIG_CONSTANT(T, 113, 78.6093129753298570701376952709727391),
0289       BOOST_MATH_BIG_CONSTANT(T, 113, 307.470246050318322489781182863190127),
0290       BOOST_MATH_BIG_CONSTANT(T, 113, 805.140686101151538537565264188630079),
0291       BOOST_MATH_BIG_CONSTANT(T, 113, 1439.12019760292146454787601409644413),
0292       BOOST_MATH_BIG_CONSTANT(T, 113, 1735.6105285756048831268586001383127),
0293       BOOST_MATH_BIG_CONSTANT(T, 113, 1348.32500712856328019355198611280536),
0294       BOOST_MATH_BIG_CONSTANT(T, 113, 607.225985860570846699704222144650563),
0295       BOOST_MATH_BIG_CONSTANT(T, 113, 119.952317857277045332558673164517227),
0296       BOOST_MATH_BIG_CONSTANT(T, 113, 0.000140165918355036060868680809129436084),
0297    };
0298
0299    // Maximum Deviation Found:                     2.867e-035
0300    // Expected Error Term :                        2.866e-035
0301    // Maximum Relative Change in Control Points :  2.662e-004
0302    static const T P_8_16[] = {
0303       BOOST_MATH_BIG_CONSTANT(T, 113, -0.184828315274146610610872315609837439e-19),
0304       BOOST_MATH_BIG_CONSTANT(T, 113, 0.500000000000000004122475157735807738),
0305       BOOST_MATH_BIG_CONSTANT(T, 113, 3.02533865247313349284875558880415875),
0306       BOOST_MATH_BIG_CONSTANT(T, 113, 13.5995927517457371243039532492642734),
0307       BOOST_MATH_BIG_CONSTANT(T, 113, 35.3132224283087906757037999452941588),
0308       BOOST_MATH_BIG_CONSTANT(T, 113, 67.1639424550714159157603179911505619),
0309       BOOST_MATH_BIG_CONSTANT(T, 113, 83.5767733658513967581959839367419891),
0310       BOOST_MATH_BIG_CONSTANT(T, 113, 71.073491212235705900866411319363501),
0311       BOOST_MATH_BIG_CONSTANT(T, 113, 35.8621515614725564575893663483998663),
0312       BOOST_MATH_BIG_CONSTANT(T, 113, 8.72152231639983491987779743154333318),
0313    };
0314    static const T Q_8_16[] = {
0315       BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),
0316       BOOST_MATH_BIG_CONSTANT(T, 113, 5.71734397161293452310624822415866372),
0317       BOOST_MATH_BIG_CONSTANT(T, 113, 25.293404179620438179337103263274815),
0318       BOOST_MATH_BIG_CONSTANT(T, 113, 62.2619767967468199111077640625328469),
0319       BOOST_MATH_BIG_CONSTANT(T, 113, 113.955048909238993473389714972250235),
0320       BOOST_MATH_BIG_CONSTANT(T, 113, 130.807138328938966981862203944329408),
0321       BOOST_MATH_BIG_CONSTANT(T, 113, 102.423146902337654110717764213057753),
0322       BOOST_MATH_BIG_CONSTANT(T, 113, 44.0424772805245202514468199602123565),
0323       BOOST_MATH_BIG_CONSTANT(T, 113, 8.89898032477904072082994913461386099),
0324       BOOST_MATH_BIG_CONSTANT(T, 113, -0.0296627336872039988632793863671456398),
0325    };
0326    // Maximum Deviation Found:                     1.079e-035
0327    // Expected Error Term :                       -1.079e-035
0328    // Maximum Relative Change in Control Points :  7.884e-003
0329    static const T P_16_inf[] = {
0330       BOOST_MATH_BIG_CONSTANT(T, 113, 0.0),
0331       BOOST_MATH_BIG_CONSTANT(T, 113, 0.500000000000000000000000000000087317),
0332       BOOST_MATH_BIG_CONSTANT(T, 113, 0.345625669885456215194494735902663968),
0333       BOOST_MATH_BIG_CONSTANT(T, 113, 9.62895499360842232127552650044647769),
0334       BOOST_MATH_BIG_CONSTANT(T, 113, 3.5936085382439026269301003761320812),
0335       BOOST_MATH_BIG_CONSTANT(T, 113, 49.459599118438883265036646019410669),
0336       BOOST_MATH_BIG_CONSTANT(T, 113, 7.77519237321893917784735690560496607),
0337       BOOST_MATH_BIG_CONSTANT(T, 113, 74.4536074488178075948642351179304121),
0338       BOOST_MATH_BIG_CONSTANT(T, 113, 2.75209340397069050436806159297952699),
0339       BOOST_MATH_BIG_CONSTANT(T, 113, 23.9292359711471667884504840186561598),
0340    };
0341    static const T Q_16_inf[] = {
0342       BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),
0343       BOOST_MATH_BIG_CONSTANT(T, 113, 0.357918006437579097055656138920742037),
0344       BOOST_MATH_BIG_CONSTANT(T, 113, 19.1386039850709849435325005484512944),
0345       BOOST_MATH_BIG_CONSTANT(T, 113, 0.874349081464143606016221431763364517),
0346       BOOST_MATH_BIG_CONSTANT(T, 113, 98.6516097434855572678195488061432509),
0347       BOOST_MATH_BIG_CONSTANT(T, 113, -16.1051972833382893468655223662534306),
0348       BOOST_MATH_BIG_CONSTANT(T, 113, 154.316860216253720989145047141653727),
0349       BOOST_MATH_BIG_CONSTANT(T, 113, -40.2026880424378986053105969312264534),
0350       BOOST_MATH_BIG_CONSTANT(T, 113, 60.1679136674264778074736441126810223),
0351       BOOST_MATH_BIG_CONSTANT(T, 113, -13.3414844622256422644504472438320114),
0352       BOOST_MATH_BIG_CONSTANT(T, 113, 2.53795636200649908779512969030363442),
0353    };
0354
0355    if(x <= 2)
0356    {
0357       return (2 + boost::math::tools::evaluate_polynomial(P_1_2, x) / tools::evaluate_polynomial(Q_1_2, x)) / (x * x);
0358    }
0359    else if(x <= 4)
0360    {
0361       return (y_offset_2_4 + boost::math::tools::evaluate_polynomial(P_2_4, x) / tools::evaluate_polynomial(Q_2_4, x)) / (x * x);
0362    }
0363    else if(x <= 8)
0364    {
0365       T y = 1 / x;
0366       return (1 + tools::evaluate_polynomial(P_4_8, y) / tools::evaluate_polynomial(Q_4_8, y)) / x;
0367    }
0368    else if(x <= 16)
0369    {
0370       T y = 1 / x;
0371       return (1 + tools::evaluate_polynomial(P_8_16, y) / tools::evaluate_polynomial(Q_8_16, y)) / x;
0372    }
0373    T y = 1 / x;
0374    return (1 + tools::evaluate_polynomial(P_16_inf, y) / tools::evaluate_polynomial(Q_16_inf, y)) / x;
0375 }
0376
0377 template <class T, class Policy, class Tag>
0378 BOOST_MATH_GPU_ENABLED T trigamma_dispatch(T x, const Policy& pol, const Tag& tag)
0379 {
0380    //
0381    // This handles reflection of negative arguments, and all our
0382    // error handling, then forwards to the T-specific approximation.
0383    //
0384    BOOST_MATH_STD_USING // ADL of std functions.
0385
0386    T result = 0;
0387    //
0388    // Check for negative arguments and use reflection:
0389    //
0390    if(x <= 0)
0391    {
0392       // Reflect:
0393       T z = 1 - x;
0394
0395       if(z < 1)
0396       {
0397          result = 1 / (z * z);
0398          z += 1;
0399       }
0400
0401       // Argument reduction for tan:
0402       if(floor(x) == x)
0403       {
0404          return policies::raise_pole_error<T>("boost::math::trigamma<%1%>(%1%)", nullptr, (1-x), pol);
0405       }
0406       T s = fabs(x) < fabs(z) ? boost::math::sin_pi(x, pol) : boost::math::sin_pi(z, pol);
0407       return result - trigamma_prec(T(z), pol, tag) + boost::math::pow<2>(constants::pi<T>()) / (s * s);
0408    }
0409    if(x < 1)
0410    {
0411       result = 1 / (x * x);
0412       x += 1;
0413    }
0414    return result + trigamma_prec(x, pol, tag);
0415 }
0416
0417 //
0418 // Initializer: ensure all our constants are initialized prior to the first call of main:
0419 //
0420 template <class T, class Policy>
0421 struct trigamma_initializer
0422 {
0423    struct init
0424    {
0425       BOOST_MATH_GPU_ENABLED init()
0426       {
0427          typedef typename policies::precision<T, Policy>::type precision_type;
0428          do_init(boost::math::integral_constant<bool, precision_type::value && (precision_type::value <= 113)>());
0429       }
0430       BOOST_MATH_GPU_ENABLED void do_init(const boost::math::true_type&)
0431       {
0432          boost::math::trigamma(T(2.5), Policy());
0433       }
0434       BOOST_MATH_GPU_ENABLED void do_init(const boost::math::false_type&){}
0435       BOOST_MATH_GPU_ENABLED void force_instantiate()const{}
0436    };
0437    static const init initializer;
0438    BOOST_MATH_GPU_ENABLED static void force_instantiate()
0439    {
0440       #ifndef BOOST_MATH_HAS_GPU_SUPPORT
0441       initializer.force_instantiate();
0442       #endif
0443    }
0444 };
0445
0446 template <class T, class Policy>
0447 const typename trigamma_initializer<T, Policy>::init trigamma_initializer<T, Policy>::initializer;
0448
0449 } // namespace detail
0450
0451 template <class T, class Policy>
0452 BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T>::type
0453    trigamma(T x, const Policy&)
0454 {
0455    typedef typename tools::promote_args<T>::type result_type;
0456    typedef typename policies::evaluation<result_type, Policy>::type value_type;
0457    typedef typename policies::precision<T, Policy>::type precision_type;
0458    typedef boost::math::integral_constant<int,
0459       precision_type::value <= 0 ? 0 :
0460       precision_type::value <= 53 ? 53 :
0461       precision_type::value <= 64 ? 64 :
0462       precision_type::value <= 113 ? 113 : 0
0463    > tag_type;
0464    typedef typename policies::normalise<
0465       Policy,
0466       policies::promote_float<false>,
0467       policies::promote_double<false>,
0468       policies::discrete_quantile<>,
0469       policies::assert_undefined<> >::type forwarding_policy;
0470
0471    // Force initialization of constants:
0472    detail::trigamma_initializer<value_type, forwarding_policy>::force_instantiate();
0473
0474    return policies::checked_narrowing_cast<result_type, Policy>(detail::trigamma_dispatch(
0475       static_cast<value_type>(x),
0476       forwarding_policy(),
0477       tag_type()), "boost::math::trigamma<%1%>(%1%)");
0478 }
0479
0480 template <class T>
0481 BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T>::type
0482    trigamma(T x)
0483 {
0484    return trigamma(x, policies::policy<>());
0485 }
0486
0487 } // namespace math
0488 } // namespace boost
0489 #endif
0490