ceph/src/boost/libs/math/include/boost/math/special_functions/trigamma.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_SF_TRIGAMMA_HPP
   7 #define BOOST_MATH_SF_TRIGAMMA_HPP
   8
   9 #ifdef _MSC_VER
  10 #pragma once
  11 #endif
  12
  13 #include <boost/math/special_functions/math_fwd.hpp>
  14 #include <boost/math/tools/rational.hpp>
  15 #include <boost/math/tools/series.hpp>
  16 #include <boost/math/tools/promotion.hpp>
  17 #include <boost/math/policies/error_handling.hpp>
  18 #include <boost/math/constants/constants.hpp>
  19 #include <boost/mpl/comparison.hpp>
  20 #include <boost/math/tools/big_constant.hpp>
  21 #include <boost/math/special_functions/polygamma.hpp>
  22
  23 namespace boost{
  24 namespace math{
  25 namespace detail{
  26
  27 template<class T, class Policy>
  28 T polygamma_imp(const int n, T x, const Policy &pol);
  29
  30 template <class T, class Policy>
  31 T trigamma_prec(T x, const mpl::int_<53>*, const Policy&)
  32 {
  33    // Max error in interpolated form: 3.736e-017
  34    static const T offset = BOOST_MATH_BIG_CONSTANT(T, 53, 2.1093254089355469);
  35    static const T P_1_2[] = {
  36       BOOST_MATH_BIG_CONSTANT(T, 53, -1.1093280605946045),
  37       BOOST_MATH_BIG_CONSTANT(T, 53, -3.8310674472619321),
  38       BOOST_MATH_BIG_CONSTANT(T, 53, -3.3703848401898283),
  39       BOOST_MATH_BIG_CONSTANT(T, 53, 0.28080574467981213),
  40       BOOST_MATH_BIG_CONSTANT(T, 53, 1.6638069578676164),
  41       BOOST_MATH_BIG_CONSTANT(T, 53, 0.64468386819102836),
  42    };
  43    static const T Q_1_2[] = {
  44       BOOST_MATH_BIG_CONSTANT(T, 53, 1.0),
  45       BOOST_MATH_BIG_CONSTANT(T, 53, 3.4535389668541151),
  46       BOOST_MATH_BIG_CONSTANT(T, 53, 4.5208926987851437),
  47       BOOST_MATH_BIG_CONSTANT(T, 53, 2.7012734178351534),
  48       BOOST_MATH_BIG_CONSTANT(T, 53, 0.64468798399785611),
  49       BOOST_MATH_BIG_CONSTANT(T, 53, -0.20314516859987728e-6),
  50    };
  51    // Max error in interpolated form: 1.159e-017
  52    static const T P_2_4[] = {
  53       BOOST_MATH_BIG_CONSTANT(T, 53, -0.13803835004508849e-7),
  54       BOOST_MATH_BIG_CONSTANT(T, 53, 0.50000049158540261),
  55       BOOST_MATH_BIG_CONSTANT(T, 53, 1.6077979838469348),
  56       BOOST_MATH_BIG_CONSTANT(T, 53, 2.5645435828098254),
  57       BOOST_MATH_BIG_CONSTANT(T, 53, 2.0534873203680393),
  58       BOOST_MATH_BIG_CONSTANT(T, 53, 0.74566981111565923),
  59    };
  60    static const T Q_2_4[] = {
  61       BOOST_MATH_BIG_CONSTANT(T, 53, 1.0),
  62       BOOST_MATH_BIG_CONSTANT(T, 53, 2.8822787662376169),
  63       BOOST_MATH_BIG_CONSTANT(T, 53, 4.1681660554090917),
  64       BOOST_MATH_BIG_CONSTANT(T, 53, 2.7853527819234466),
  65       BOOST_MATH_BIG_CONSTANT(T, 53, 0.74967671848044792),
  66       BOOST_MATH_BIG_CONSTANT(T, 53, -0.00057069112416246805),
  67    };
  68    // Maximum Deviation Found:                     6.896e-018
  69    // Expected Error Term :                       -6.895e-018
  70    // Maximum Relative Change in Control Points :  8.497e-004
  71    static const T P_4_inf[] = {
  72       static_cast<T>(0.68947581948701249e-17L),
  73       static_cast<T>(0.49999999999998975L),
  74       static_cast<T>(1.0177274392923795L),
  75       static_cast<T>(2.498208511343429L),
  76       static_cast<T>(2.1921221359427595L),
  77       static_cast<T>(1.5897035272532764L),
  78       static_cast<T>(0.40154388356961734L),
  79    };
  80    static const T Q_4_inf[] = {
  81       static_cast<T>(1.0L),
  82       static_cast<T>(1.7021215452463932L),
  83       static_cast<T>(4.4290431747556469L),
  84       static_cast<T>(2.9745631894384922L),
  85       static_cast<T>(2.3013614809773616L),
  86       static_cast<T>(0.28360399799075752L),
  87       static_cast<T>(0.022892987908906897L),
  88    };
  89
  90    if(x <= 2)
  91    {
  92       return (offset + boost::math::tools::evaluate_polynomial(P_1_2, x) / tools::evaluate_polynomial(Q_1_2, x)) / (x * x);
  93    }
  94    else if(x <= 4)
  95    {
  96       T y = 1 / x;
  97       return (1 + tools::evaluate_polynomial(P_2_4, y) / tools::evaluate_polynomial(Q_2_4, y)) / x;
  98    }
  99    T y = 1 / x;
 100    return (1 + tools::evaluate_polynomial(P_4_inf, y) / tools::evaluate_polynomial(Q_4_inf, y)) / x;
 101 }
 102
 103 template <class T, class Policy>
 104 T trigamma_prec(T x, const mpl::int_<64>*, const Policy&)
 105 {
 106    // Max error in interpolated form: 1.178e-020
 107    static const T offset_1_2 = BOOST_MATH_BIG_CONSTANT(T, 64, 2.109325408935546875);
 108    static const T P_1_2[] = {
 109       BOOST_MATH_BIG_CONSTANT(T, 64, -1.10932535608960258341),
 110       BOOST_MATH_BIG_CONSTANT(T, 64, -4.18793841543017129052),
 111       BOOST_MATH_BIG_CONSTANT(T, 64, -4.63865531898487734531),
 112       BOOST_MATH_BIG_CONSTANT(T, 64, -0.919832884430500908047),
 113       BOOST_MATH_BIG_CONSTANT(T, 64, 1.68074038333180423012),
 114       BOOST_MATH_BIG_CONSTANT(T, 64, 1.21172611429185622377),
 115       BOOST_MATH_BIG_CONSTANT(T, 64, 0.259635673503366427284),
 116    };
 117    static const T Q_1_2[] = {
 118       BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
 119       BOOST_MATH_BIG_CONSTANT(T, 64, 3.77521119359546982995),
 120       BOOST_MATH_BIG_CONSTANT(T, 64, 5.664338024578956321),
 121       BOOST_MATH_BIG_CONSTANT(T, 64, 4.25995134879278028361),
 122       BOOST_MATH_BIG_CONSTANT(T, 64, 1.62956638448940402182),
 123       BOOST_MATH_BIG_CONSTANT(T, 64, 0.259635512844691089868),
 124       BOOST_MATH_BIG_CONSTANT(T, 64, 0.629642219810618032207e-8),
 125    };
 126    // Max error in interpolated form: 3.912e-020
 127    static const T P_2_8[] = {
 128       BOOST_MATH_BIG_CONSTANT(T, 64, -0.387540035162952880976e-11),
 129       BOOST_MATH_BIG_CONSTANT(T, 64, 0.500000000276430504),
 130       BOOST_MATH_BIG_CONSTANT(T, 64, 3.21926880986360957306),
 131       BOOST_MATH_BIG_CONSTANT(T, 64, 10.2550347708483445775),
 132       BOOST_MATH_BIG_CONSTANT(T, 64, 18.9002075150709144043),
 133       BOOST_MATH_BIG_CONSTANT(T, 64, 21.0357215832399705625),
 134       BOOST_MATH_BIG_CONSTANT(T, 64, 13.4346512182925923978),
 135       BOOST_MATH_BIG_CONSTANT(T, 64, 3.98656291026448279118),
 136    };
 137    static const T Q_2_8[] = {
 138       BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
 139       BOOST_MATH_BIG_CONSTANT(T, 64, 6.10520430478613667724),
 140       BOOST_MATH_BIG_CONSTANT(T, 64, 18.475001060603645512),
 141       BOOST_MATH_BIG_CONSTANT(T, 64, 31.7087534567758405638),
 142       BOOST_MATH_BIG_CONSTANT(T, 64, 31.908814523890465398),
 143       BOOST_MATH_BIG_CONSTANT(T, 64, 17.4175479039227084798),
 144       BOOST_MATH_BIG_CONSTANT(T, 64, 3.98749106958394941276),
 145       BOOST_MATH_BIG_CONSTANT(T, 64, -0.000115917322224411128566),
 146    };
 147    // Maximum Deviation Found:                     2.635e-020
 148    // Expected Error Term :                        2.635e-020
 149    // Maximum Relative Change in Control Points :  1.791e-003
 150    static const T P_8_inf[] = {
 151       BOOST_MATH_BIG_CONSTANT(T, 64, -0.263527875092466899848e-19),
 152       BOOST_MATH_BIG_CONSTANT(T, 64, 0.500000000000000058145),
 153       BOOST_MATH_BIG_CONSTANT(T, 64, 0.0730121433777364138677),
 154       BOOST_MATH_BIG_CONSTANT(T, 64, 1.94505878379957149534),
 155       BOOST_MATH_BIG_CONSTANT(T, 64, 0.0517092358874932620529),
 156       BOOST_MATH_BIG_CONSTANT(T, 64, 1.07995383547483921121),
 157    };
 158    static const T Q_8_inf[] = {
 159       BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
 160       BOOST_MATH_BIG_CONSTANT(T, 64, -0.187309046577818095504),
 161       BOOST_MATH_BIG_CONSTANT(T, 64, 3.95255391645238842975),
 162       BOOST_MATH_BIG_CONSTANT(T, 64, -1.14743283327078949087),
 163       BOOST_MATH_BIG_CONSTANT(T, 64, 2.52989799376344914499),
 164       BOOST_MATH_BIG_CONSTANT(T, 64, -0.627414303172402506396),
 165       BOOST_MATH_BIG_CONSTANT(T, 64, 0.141554248216425512536),
 166    };
 167
 168    if(x <= 2)
 169    {
 170       return (offset_1_2 + boost::math::tools::evaluate_polynomial(P_1_2, x) / tools::evaluate_polynomial(Q_1_2, x)) / (x * x);
 171    }
 172    else if(x <= 8)
 173    {
 174       T y = 1 / x;
 175       return (1 + tools::evaluate_polynomial(P_2_8, y) / tools::evaluate_polynomial(Q_2_8, y)) / x;
 176    }
 177    T y = 1 / x;
 178    return (1 + tools::evaluate_polynomial(P_8_inf, y) / tools::evaluate_polynomial(Q_8_inf, y)) / x;
 179 }
 180
 181 template <class T, class Policy>
 182 T trigamma_prec(T x, const mpl::int_<113>*, const Policy&)
 183 {
 184    // Max error in interpolated form: 1.916e-035
 185
 186    static const T P_1_2[] = {
 187       BOOST_MATH_BIG_CONSTANT(T, 113, -0.999999999999999082554457936871832533),
 188       BOOST_MATH_BIG_CONSTANT(T, 113, -4.71237311120865266379041700054847734),
 189       BOOST_MATH_BIG_CONSTANT(T, 113, -7.94125711970499027763789342500817316),
 190       BOOST_MATH_BIG_CONSTANT(T, 113, -5.74657746697664735258222071695644535),
 191       BOOST_MATH_BIG_CONSTANT(T, 113, -0.404213349456398905981223965160595687),
 192       BOOST_MATH_BIG_CONSTANT(T, 113, 2.47877781178642876561595890095758896),
 193       BOOST_MATH_BIG_CONSTANT(T, 113, 2.07714151702455125992166949812126433),
 194       BOOST_MATH_BIG_CONSTANT(T, 113, 0.858877899162360138844032265418028567),
 195       BOOST_MATH_BIG_CONSTANT(T, 113, 0.20499222604410032375789018837922397),
 196       BOOST_MATH_BIG_CONSTANT(T, 113, 0.0272103140348194747360175268778415049),
 197       BOOST_MATH_BIG_CONSTANT(T, 113, 0.0015764849020876949848954081173520686),
 198    };
 199    static const T Q_1_2[] = {
 200       BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),
 201       BOOST_MATH_BIG_CONSTANT(T, 113, 4.71237311120863419878375031457715223),
 202       BOOST_MATH_BIG_CONSTANT(T, 113, 9.58619118655339853449127952145877467),
 203       BOOST_MATH_BIG_CONSTANT(T, 113, 11.0940067269829372437561421279054968),
 204       BOOST_MATH_BIG_CONSTANT(T, 113, 8.09075424749327792073276309969037885),
 205       BOOST_MATH_BIG_CONSTANT(T, 113, 3.87705890159891405185343806884451286),
 206       BOOST_MATH_BIG_CONSTANT(T, 113, 1.22758678701914477836330837816976782),
 207       BOOST_MATH_BIG_CONSTANT(T, 113, 0.249092040606385004109672077814668716),
 208       BOOST_MATH_BIG_CONSTANT(T, 113, 0.0295750413900655597027079600025569048),
 209       BOOST_MATH_BIG_CONSTANT(T, 113, 0.00157648490200498142247694709728858139),
 210       BOOST_MATH_BIG_CONSTANT(T, 113, 0.161264050344059471721062360645432809e-14),
 211    };
 212
 213    // Max error in interpolated form: 8.958e-035
 214    static const T P_2_4[] = {
 215       BOOST_MATH_BIG_CONSTANT(T, 113, -2.55843734739907925764326773972215085),
 216       BOOST_MATH_BIG_CONSTANT(T, 113, -12.2830208240542011967952466273455887),
 217       BOOST_MATH_BIG_CONSTANT(T, 113, -23.9195022162767993526575786066414403),
 218       BOOST_MATH_BIG_CONSTANT(T, 113, -24.9256431504823483094158828285470862),
 219       BOOST_MATH_BIG_CONSTANT(T, 113, -14.7979122765478779075108064826412285),
 220       BOOST_MATH_BIG_CONSTANT(T, 113, -4.46654453928610666393276765059122272),
 221       BOOST_MATH_BIG_CONSTANT(T, 113, -0.0191439033405649675717082465687845002),
 222       BOOST_MATH_BIG_CONSTANT(T, 113, 0.515412052554351265708917209749037352),
 223       BOOST_MATH_BIG_CONSTANT(T, 113, 0.195378348786064304378247325360320038),
 224       BOOST_MATH_BIG_CONSTANT(T, 113, 0.0334761282624174313035014426794245393),
 225       BOOST_MATH_BIG_CONSTANT(T, 113, 0.002373665205942206348500250056602687),
 226    };
 227    static const T Q_2_4[] = {
 228       BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),
 229       BOOST_MATH_BIG_CONSTANT(T, 113, 4.80098558454419907830670928248659245),
 230       BOOST_MATH_BIG_CONSTANT(T, 113, 9.99220727843170133895059300223445265),
 231       BOOST_MATH_BIG_CONSTANT(T, 113, 11.8896146167631330735386697123464976),
 232       BOOST_MATH_BIG_CONSTANT(T, 113, 8.96613256683809091593793565879092581),
 233       BOOST_MATH_BIG_CONSTANT(T, 113, 4.47254136149624110878909334574485751),
 234       BOOST_MATH_BIG_CONSTANT(T, 113, 1.48600982028196527372434773913633152),
 235       BOOST_MATH_BIG_CONSTANT(T, 113, 0.319570735766764237068541501137990078),
 236       BOOST_MATH_BIG_CONSTANT(T, 113, 0.0407358345787680953107374215319322066),
 237       BOOST_MATH_BIG_CONSTANT(T, 113, 0.00237366520593271641375755486420859837),
 238       BOOST_MATH_BIG_CONSTANT(T, 113, 0.239554887903526152679337256236302116e-15),
 239       BOOST_MATH_BIG_CONSTANT(T, 113, -0.294749244740618656265237072002026314e-17),
 240    };
 241
 242    static const T y_offset_2_4 = BOOST_MATH_BIG_CONSTANT(T, 113, 3.558437347412109375);
 243
 244    // Max error in interpolated form: 4.319e-035
 245    static const T P_4_8[] = {
 246       BOOST_MATH_BIG_CONSTANT(T, 113, 0.166626112697021464248967707021688845e-16),
 247       BOOST_MATH_BIG_CONSTANT(T, 113, 0.499999999999997739552090249208808197),
 248       BOOST_MATH_BIG_CONSTANT(T, 113, 6.40270945019053817915772473771553187),
 249       BOOST_MATH_BIG_CONSTANT(T, 113, 41.3833374155000608013677627389343329),
 250       BOOST_MATH_BIG_CONSTANT(T, 113, 166.803341854562809335667241074035245),
 251       BOOST_MATH_BIG_CONSTANT(T, 113, 453.39964786925369319960722793414521),
 252       BOOST_MATH_BIG_CONSTANT(T, 113, 851.153712317697055375935433362983944),
 253       BOOST_MATH_BIG_CONSTANT(T, 113, 1097.70657567285059133109286478004458),
 254       BOOST_MATH_BIG_CONSTANT(T, 113, 938.431232478455316020076349367632922),
 255       BOOST_MATH_BIG_CONSTANT(T, 113, 487.268001604651932322080970189930074),
 256       BOOST_MATH_BIG_CONSTANT(T, 113, 119.953445242335730062471193124820659),
 257    };
 258    static const T Q_4_8[] = {
 259       BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),
 260       BOOST_MATH_BIG_CONSTANT(T, 113, 12.4720855670474488978638945855932398),
 261       BOOST_MATH_BIG_CONSTANT(T, 113, 78.6093129753298570701376952709727391),
 262       BOOST_MATH_BIG_CONSTANT(T, 113, 307.470246050318322489781182863190127),
 263       BOOST_MATH_BIG_CONSTANT(T, 113, 805.140686101151538537565264188630079),
 264       BOOST_MATH_BIG_CONSTANT(T, 113, 1439.12019760292146454787601409644413),
 265       BOOST_MATH_BIG_CONSTANT(T, 113, 1735.6105285756048831268586001383127),
 266       BOOST_MATH_BIG_CONSTANT(T, 113, 1348.32500712856328019355198611280536),
 267       BOOST_MATH_BIG_CONSTANT(T, 113, 607.225985860570846699704222144650563),
 268       BOOST_MATH_BIG_CONSTANT(T, 113, 119.952317857277045332558673164517227),
 269       BOOST_MATH_BIG_CONSTANT(T, 113, 0.000140165918355036060868680809129436084),
 270    };
 271
 272    // Maximum Deviation Found:                     2.867e-035
 273    // Expected Error Term :                        2.866e-035
 274    // Maximum Relative Change in Control Points :  2.662e-004
 275    static const T P_8_16[] = {
 276       BOOST_MATH_BIG_CONSTANT(T, 113, -0.184828315274146610610872315609837439e-19),
 277       BOOST_MATH_BIG_CONSTANT(T, 113, 0.500000000000000004122475157735807738),
 278       BOOST_MATH_BIG_CONSTANT(T, 113, 3.02533865247313349284875558880415875),
 279       BOOST_MATH_BIG_CONSTANT(T, 113, 13.5995927517457371243039532492642734),
 280       BOOST_MATH_BIG_CONSTANT(T, 113, 35.3132224283087906757037999452941588),
 281       BOOST_MATH_BIG_CONSTANT(T, 113, 67.1639424550714159157603179911505619),
 282       BOOST_MATH_BIG_CONSTANT(T, 113, 83.5767733658513967581959839367419891),
 283       BOOST_MATH_BIG_CONSTANT(T, 113, 71.073491212235705900866411319363501),
 284       BOOST_MATH_BIG_CONSTANT(T, 113, 35.8621515614725564575893663483998663),
 285       BOOST_MATH_BIG_CONSTANT(T, 113, 8.72152231639983491987779743154333318),
 286    };
 287    static const T Q_8_16[] = {
 288       BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),
 289       BOOST_MATH_BIG_CONSTANT(T, 113, 5.71734397161293452310624822415866372),
 290       BOOST_MATH_BIG_CONSTANT(T, 113, 25.293404179620438179337103263274815),
 291       BOOST_MATH_BIG_CONSTANT(T, 113, 62.2619767967468199111077640625328469),
 292       BOOST_MATH_BIG_CONSTANT(T, 113, 113.955048909238993473389714972250235),
 293       BOOST_MATH_BIG_CONSTANT(T, 113, 130.807138328938966981862203944329408),
 294       BOOST_MATH_BIG_CONSTANT(T, 113, 102.423146902337654110717764213057753),
 295       BOOST_MATH_BIG_CONSTANT(T, 113, 44.0424772805245202514468199602123565),
 296       BOOST_MATH_BIG_CONSTANT(T, 113, 8.89898032477904072082994913461386099),
 297       BOOST_MATH_BIG_CONSTANT(T, 113, -0.0296627336872039988632793863671456398),
 298    };
 299    // Maximum Deviation Found:                     1.079e-035
 300    // Expected Error Term :                       -1.079e-035
 301    // Maximum Relative Change in Control Points :  7.884e-003
 302    static const T P_16_inf[] = {
 303       BOOST_MATH_BIG_CONSTANT(T, 113, 0.0),
 304       BOOST_MATH_BIG_CONSTANT(T, 113, 0.500000000000000000000000000000087317),
 305       BOOST_MATH_BIG_CONSTANT(T, 113, 0.345625669885456215194494735902663968),
 306       BOOST_MATH_BIG_CONSTANT(T, 113, 9.62895499360842232127552650044647769),
 307       BOOST_MATH_BIG_CONSTANT(T, 113, 3.5936085382439026269301003761320812),
 308       BOOST_MATH_BIG_CONSTANT(T, 113, 49.459599118438883265036646019410669),
 309       BOOST_MATH_BIG_CONSTANT(T, 113, 7.77519237321893917784735690560496607),
 310       BOOST_MATH_BIG_CONSTANT(T, 113, 74.4536074488178075948642351179304121),
 311       BOOST_MATH_BIG_CONSTANT(T, 113, 2.75209340397069050436806159297952699),
 312       BOOST_MATH_BIG_CONSTANT(T, 113, 23.9292359711471667884504840186561598),
 313    };
 314    static const T Q_16_inf[] = {
 315       BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),
 316       BOOST_MATH_BIG_CONSTANT(T, 113, 0.357918006437579097055656138920742037),
 317       BOOST_MATH_BIG_CONSTANT(T, 113, 19.1386039850709849435325005484512944),
 318       BOOST_MATH_BIG_CONSTANT(T, 113, 0.874349081464143606016221431763364517),
 319       BOOST_MATH_BIG_CONSTANT(T, 113, 98.6516097434855572678195488061432509),
 320       BOOST_MATH_BIG_CONSTANT(T, 113, -16.1051972833382893468655223662534306),
 321       BOOST_MATH_BIG_CONSTANT(T, 113, 154.316860216253720989145047141653727),
 322       BOOST_MATH_BIG_CONSTANT(T, 113, -40.2026880424378986053105969312264534),
 323       BOOST_MATH_BIG_CONSTANT(T, 113, 60.1679136674264778074736441126810223),
 324       BOOST_MATH_BIG_CONSTANT(T, 113, -13.3414844622256422644504472438320114),
 325       BOOST_MATH_BIG_CONSTANT(T, 113, 2.53795636200649908779512969030363442),
 326    };
 327
 328    if(x <= 2)
 329    {
 330       return (2 + boost::math::tools::evaluate_polynomial(P_1_2, x) / tools::evaluate_polynomial(Q_1_2, x)) / (x * x);
 331    }
 332    else if(x <= 4)
 333    {
 334       return (y_offset_2_4 + boost::math::tools::evaluate_polynomial(P_2_4, x) / tools::evaluate_polynomial(Q_2_4, x)) / (x * x);
 335    }
 336    else if(x <= 8)
 337    {
 338       T y = 1 / x;
 339       return (1 + tools::evaluate_polynomial(P_4_8, y) / tools::evaluate_polynomial(Q_4_8, y)) / x;
 340    }
 341    else if(x <= 16)
 342    {
 343       T y = 1 / x;
 344       return (1 + tools::evaluate_polynomial(P_8_16, y) / tools::evaluate_polynomial(Q_8_16, y)) / x;
 345    }
 346    T y = 1 / x;
 347    return (1 + tools::evaluate_polynomial(P_16_inf, y) / tools::evaluate_polynomial(Q_16_inf, y)) / x;
 348 }
 349
 350 template <class T, class Tag, class Policy>
 351 T trigamma_imp(T x, const Tag* t, const Policy& pol)
 352 {
 353    //
 354    // This handles reflection of negative arguments, and all our
 355    // error handling, then forwards to the T-specific approximation.
 356    //
 357    BOOST_MATH_STD_USING // ADL of std functions.
 358
 359    T result = 0;
 360    //
 361    // Check for negative arguments and use reflection:
 362    //
 363    if(x <= 0)
 364    {
 365       // Reflect:
 366       T z = 1 - x;
 367       // Argument reduction for tan:
 368       if(floor(x) == x)
 369       {
 370          return policies::raise_pole_error<T>("boost::math::trigamma<%1%>(%1%)", 0, (1-x), pol);
 371       }
 372       T s = fabs(x) < fabs(z) ? boost::math::sin_pi(x, pol) : boost::math::sin_pi(z, pol);
 373       return -trigamma_imp(z, t, pol) + boost::math::pow<2>(constants::pi<T>()) / (s * s);
 374    }
 375    if(x < 1)
 376    {
 377       result = 1 / (x * x);
 378       x += 1;
 379    }
 380    return result + trigamma_prec(x, t, pol);
 381 }
 382
 383 template <class T, class Policy>
 384 T trigamma_imp(T x, const mpl::int_<0>*, const Policy& pol)
 385 {
 386    return polygamma_imp(1, x, pol);
 387 }
 388 //
 389 // Initializer: ensure all our constants are initialized prior to the first call of main:
 390 //
 391 template <class T, class Policy>
 392 struct trigamma_initializer
 393 {
 394    struct init
 395    {
 396       init()
 397       {
 398          typedef typename policies::precision<T, Policy>::type precision_type;
 399          do_init(mpl::bool_<precision_type::value && (precision_type::value <= 113)>());
 400       }
 401       void do_init(const mpl::true_&)
 402       {
 403          boost::math::trigamma(T(2.5), Policy());
 404       }
 405       void do_init(const mpl::false_&){}
 406       void force_instantiate()const{}
 407    };
 408    static const init initializer;
 409    static void force_instantiate()
 410    {
 411       initializer.force_instantiate();
 412    }
 413 };
 414
 415 template <class T, class Policy>
 416 const typename trigamma_initializer<T, Policy>::init trigamma_initializer<T, Policy>::initializer;
 417
 418 } // namespace detail
 419
 420 template <class T, class Policy>
 421 inline typename tools::promote_args<T>::type
 422    trigamma(T x, const Policy&)
 423 {
 424    typedef typename tools::promote_args<T>::type result_type;
 425    typedef typename policies::evaluation<result_type, Policy>::type value_type;
 426    typedef typename policies::precision<T, Policy>::type precision_type;
 427    typedef typename mpl::if_<
 428       mpl::or_<
 429          mpl::less_equal<precision_type, mpl::int_<0> >,
 430          mpl::greater<precision_type, mpl::int_<114> >
 431       >,
 432       mpl::int_<0>,
 433       typename mpl::if_<
 434          mpl::less<precision_type, mpl::int_<54> >,
 435          mpl::int_<53>,
 436          typename mpl::if_<
 437             mpl::less<precision_type, mpl::int_<65> >,
 438             mpl::int_<64>,
 439             mpl::int_<113>
 440          >::type
 441       >::type
 442    >::type tag_type;
 443
 444    typedef typename policies::normalise<
 445       Policy,
 446       policies::promote_float<false>,
 447       policies::promote_double<false>,
 448       policies::discrete_quantile<>,
 449       policies::assert_undefined<> >::type forwarding_policy;
 450
 451    // Force initialization of constants:
 452    detail::trigamma_initializer<value_type, forwarding_policy>::force_instantiate();
 453
 454    return policies::checked_narrowing_cast<result_type, Policy>(detail::trigamma_imp(
 455       static_cast<value_type>(x),
 456       static_cast<const tag_type*>(0), forwarding_policy()), "boost::math::trigamma<%1%>(%1%)");
 457 }
 458
 459 template <class T>
 460 inline typename tools::promote_args<T>::type
 461    trigamma(T x)
 462 {
 463    return trigamma(x, policies::policy<>());
 464 }
 465
 466 } // namespace math
 467 } // namespace boost
 468 #endif
 469