449 lines · plain
1// (C) Copyright John Maddock 2006.2// (C) Copyright Matt Borland 2024.3// Use, modification and distribution are subject to the4// Boost Software License, Version 1.0. (See accompanying file5// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)6 7#ifndef BOOST_MATH_SF_TRIGAMMA_HPP8#define BOOST_MATH_SF_TRIGAMMA_HPP9 10#ifdef _MSC_VER11#pragma once12#endif13 14#include <boost/math/tools/config.hpp>15#include <boost/math/tools/rational.hpp>16#include <boost/math/tools/promotion.hpp>17#include <boost/math/tools/big_constant.hpp>18#include <boost/math/tools/type_traits.hpp>19#include <boost/math/policies/policy.hpp>20#include <boost/math/policies/error_handling.hpp>21#include <boost/math/constants/constants.hpp>22#include <boost/math/special_functions/sin_pi.hpp>23#include <boost/math/special_functions/pow.hpp>24 25#ifndef BOOST_MATH_HAS_NVRTC26#include <boost/math/special_functions/math_fwd.hpp>27#include <boost/math/special_functions/polygamma.hpp>28#include <boost/math/tools/series.hpp>29#endif30 31#if defined(__GNUC__) && defined(BOOST_MATH_USE_FLOAT128)32//33// This is the only way we can avoid34// warning: non-standard suffix on floating constant [-Wpedantic]35// when building with -Wall -pedantic. Neither __extension__36// nor #pragma diagnostic ignored work :(37//38#pragma GCC system_header39#endif40 41namespace boost{42namespace math{43namespace detail{44 45// TODO(mborland): Temporary for NVRTC46#ifndef BOOST_MATH_HAS_NVRTC47template<class T, class Policy>48T polygamma_imp(const int n, T x, const Policy &pol);49 50template <class T, class Policy>51T trigamma_prec(T x, const Policy& pol, const boost::math::integral_constant<int, 0>&)52{53 return polygamma_imp(1, x, pol);54}55#endif56 57template <class T, class Policy>58BOOST_MATH_GPU_ENABLED T trigamma_prec(T x, const Policy&, const boost::math::integral_constant<int, 53>&)59{60 // Max error in interpolated form: 3.736e-01761 BOOST_MATH_STATIC const T offset = BOOST_MATH_BIG_CONSTANT(T, 53, 2.1093254089355469);62 BOOST_MATH_STATIC const T P_1_2[] = {63 BOOST_MATH_BIG_CONSTANT(T, 53, -1.1093280605946045),64 BOOST_MATH_BIG_CONSTANT(T, 53, -3.8310674472619321),65 BOOST_MATH_BIG_CONSTANT(T, 53, -3.3703848401898283),66 BOOST_MATH_BIG_CONSTANT(T, 53, 0.28080574467981213),67 BOOST_MATH_BIG_CONSTANT(T, 53, 1.6638069578676164),68 BOOST_MATH_BIG_CONSTANT(T, 53, 0.64468386819102836),69 };70 BOOST_MATH_STATIC const T Q_1_2[] = {71 BOOST_MATH_BIG_CONSTANT(T, 53, 1.0),72 BOOST_MATH_BIG_CONSTANT(T, 53, 3.4535389668541151),73 BOOST_MATH_BIG_CONSTANT(T, 53, 4.5208926987851437),74 BOOST_MATH_BIG_CONSTANT(T, 53, 2.7012734178351534),75 BOOST_MATH_BIG_CONSTANT(T, 53, 0.64468798399785611),76 BOOST_MATH_BIG_CONSTANT(T, 53, -0.20314516859987728e-6),77 };78 // Max error in interpolated form: 1.159e-01779 BOOST_MATH_STATIC const T P_2_4[] = {80 BOOST_MATH_BIG_CONSTANT(T, 53, -0.13803835004508849e-7),81 BOOST_MATH_BIG_CONSTANT(T, 53, 0.50000049158540261),82 BOOST_MATH_BIG_CONSTANT(T, 53, 1.6077979838469348),83 BOOST_MATH_BIG_CONSTANT(T, 53, 2.5645435828098254),84 BOOST_MATH_BIG_CONSTANT(T, 53, 2.0534873203680393),85 BOOST_MATH_BIG_CONSTANT(T, 53, 0.74566981111565923),86 };87 BOOST_MATH_STATIC const T Q_2_4[] = {88 BOOST_MATH_BIG_CONSTANT(T, 53, 1.0),89 BOOST_MATH_BIG_CONSTANT(T, 53, 2.8822787662376169),90 BOOST_MATH_BIG_CONSTANT(T, 53, 4.1681660554090917),91 BOOST_MATH_BIG_CONSTANT(T, 53, 2.7853527819234466),92 BOOST_MATH_BIG_CONSTANT(T, 53, 0.74967671848044792),93 BOOST_MATH_BIG_CONSTANT(T, 53, -0.00057069112416246805),94 };95 // Maximum Deviation Found: 6.896e-01896 // Expected Error Term : -6.895e-01897 // Maximum Relative Change in Control Points : 8.497e-00498 BOOST_MATH_STATIC const T P_4_inf[] = {99 static_cast<T>(0.68947581948701249e-17L),100 static_cast<T>(0.49999999999998975L),101 static_cast<T>(1.0177274392923795L),102 static_cast<T>(2.498208511343429L),103 static_cast<T>(2.1921221359427595L),104 static_cast<T>(1.5897035272532764L),105 static_cast<T>(0.40154388356961734L),106 };107 BOOST_MATH_STATIC const T Q_4_inf[] = {108 static_cast<T>(1.0L),109 static_cast<T>(1.7021215452463932L),110 static_cast<T>(4.4290431747556469L),111 static_cast<T>(2.9745631894384922L),112 static_cast<T>(2.3013614809773616L),113 static_cast<T>(0.28360399799075752L),114 static_cast<T>(0.022892987908906897L),115 };116 117 if(x <= 2)118 {119 return (offset + boost::math::tools::evaluate_polynomial(P_1_2, x) / tools::evaluate_polynomial(Q_1_2, x)) / (x * x);120 }121 else if(x <= 4)122 {123 T y = 1 / x;124 return (1 + tools::evaluate_polynomial(P_2_4, y) / tools::evaluate_polynomial(Q_2_4, y)) / x;125 }126 T y = 1 / x;127 return (1 + tools::evaluate_polynomial(P_4_inf, y) / tools::evaluate_polynomial(Q_4_inf, y)) / x;128}129 130template <class T, class Policy>131BOOST_MATH_GPU_ENABLED T trigamma_prec(T x, const Policy&, const boost::math::integral_constant<int, 64>&)132{133 // Max error in interpolated form: 1.178e-020134 BOOST_MATH_STATIC const T offset_1_2 = BOOST_MATH_BIG_CONSTANT(T, 64, 2.109325408935546875);135 BOOST_MATH_STATIC const T P_1_2[] = {136 BOOST_MATH_BIG_CONSTANT(T, 64, -1.10932535608960258341),137 BOOST_MATH_BIG_CONSTANT(T, 64, -4.18793841543017129052),138 BOOST_MATH_BIG_CONSTANT(T, 64, -4.63865531898487734531),139 BOOST_MATH_BIG_CONSTANT(T, 64, -0.919832884430500908047),140 BOOST_MATH_BIG_CONSTANT(T, 64, 1.68074038333180423012),141 BOOST_MATH_BIG_CONSTANT(T, 64, 1.21172611429185622377),142 BOOST_MATH_BIG_CONSTANT(T, 64, 0.259635673503366427284),143 };144 BOOST_MATH_STATIC const T Q_1_2[] = {145 BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),146 BOOST_MATH_BIG_CONSTANT(T, 64, 3.77521119359546982995),147 BOOST_MATH_BIG_CONSTANT(T, 64, 5.664338024578956321),148 BOOST_MATH_BIG_CONSTANT(T, 64, 4.25995134879278028361),149 BOOST_MATH_BIG_CONSTANT(T, 64, 1.62956638448940402182),150 BOOST_MATH_BIG_CONSTANT(T, 64, 0.259635512844691089868),151 BOOST_MATH_BIG_CONSTANT(T, 64, 0.629642219810618032207e-8),152 };153 // Max error in interpolated form: 3.912e-020154 BOOST_MATH_STATIC const T P_2_8[] = {155 BOOST_MATH_BIG_CONSTANT(T, 64, -0.387540035162952880976e-11),156 BOOST_MATH_BIG_CONSTANT(T, 64, 0.500000000276430504),157 BOOST_MATH_BIG_CONSTANT(T, 64, 3.21926880986360957306),158 BOOST_MATH_BIG_CONSTANT(T, 64, 10.2550347708483445775),159 BOOST_MATH_BIG_CONSTANT(T, 64, 18.9002075150709144043),160 BOOST_MATH_BIG_CONSTANT(T, 64, 21.0357215832399705625),161 BOOST_MATH_BIG_CONSTANT(T, 64, 13.4346512182925923978),162 BOOST_MATH_BIG_CONSTANT(T, 64, 3.98656291026448279118),163 };164 BOOST_MATH_STATIC const T Q_2_8[] = {165 BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),166 BOOST_MATH_BIG_CONSTANT(T, 64, 6.10520430478613667724),167 BOOST_MATH_BIG_CONSTANT(T, 64, 18.475001060603645512),168 BOOST_MATH_BIG_CONSTANT(T, 64, 31.7087534567758405638),169 BOOST_MATH_BIG_CONSTANT(T, 64, 31.908814523890465398),170 BOOST_MATH_BIG_CONSTANT(T, 64, 17.4175479039227084798),171 BOOST_MATH_BIG_CONSTANT(T, 64, 3.98749106958394941276),172 BOOST_MATH_BIG_CONSTANT(T, 64, -0.000115917322224411128566),173 };174 // Maximum Deviation Found: 2.635e-020175 // Expected Error Term : 2.635e-020176 // Maximum Relative Change in Control Points : 1.791e-003177 BOOST_MATH_STATIC const T P_8_inf[] = {178 BOOST_MATH_BIG_CONSTANT(T, 64, -0.263527875092466899848e-19),179 BOOST_MATH_BIG_CONSTANT(T, 64, 0.500000000000000058145),180 BOOST_MATH_BIG_CONSTANT(T, 64, 0.0730121433777364138677),181 BOOST_MATH_BIG_CONSTANT(T, 64, 1.94505878379957149534),182 BOOST_MATH_BIG_CONSTANT(T, 64, 0.0517092358874932620529),183 BOOST_MATH_BIG_CONSTANT(T, 64, 1.07995383547483921121),184 };185 BOOST_MATH_STATIC const T Q_8_inf[] = {186 BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),187 BOOST_MATH_BIG_CONSTANT(T, 64, -0.187309046577818095504),188 BOOST_MATH_BIG_CONSTANT(T, 64, 3.95255391645238842975),189 BOOST_MATH_BIG_CONSTANT(T, 64, -1.14743283327078949087),190 BOOST_MATH_BIG_CONSTANT(T, 64, 2.52989799376344914499),191 BOOST_MATH_BIG_CONSTANT(T, 64, -0.627414303172402506396),192 BOOST_MATH_BIG_CONSTANT(T, 64, 0.141554248216425512536),193 };194 195 if(x <= 2)196 {197 return (offset_1_2 + boost::math::tools::evaluate_polynomial(P_1_2, x) / tools::evaluate_polynomial(Q_1_2, x)) / (x * x);198 }199 else if(x <= 8)200 {201 T y = 1 / x;202 return (1 + tools::evaluate_polynomial(P_2_8, y) / tools::evaluate_polynomial(Q_2_8, y)) / x;203 }204 T y = 1 / x;205 return (1 + tools::evaluate_polynomial(P_8_inf, y) / tools::evaluate_polynomial(Q_8_inf, y)) / x;206}207 208template <class T, class Policy>209BOOST_MATH_GPU_ENABLED T trigamma_prec(T x, const Policy&, const boost::math::integral_constant<int, 113>&)210{211 // Max error in interpolated form: 1.916e-035212 213 static const T P_1_2[] = {214 BOOST_MATH_BIG_CONSTANT(T, 113, -0.999999999999999082554457936871832533),215 BOOST_MATH_BIG_CONSTANT(T, 113, -4.71237311120865266379041700054847734),216 BOOST_MATH_BIG_CONSTANT(T, 113, -7.94125711970499027763789342500817316),217 BOOST_MATH_BIG_CONSTANT(T, 113, -5.74657746697664735258222071695644535),218 BOOST_MATH_BIG_CONSTANT(T, 113, -0.404213349456398905981223965160595687),219 BOOST_MATH_BIG_CONSTANT(T, 113, 2.47877781178642876561595890095758896),220 BOOST_MATH_BIG_CONSTANT(T, 113, 2.07714151702455125992166949812126433),221 BOOST_MATH_BIG_CONSTANT(T, 113, 0.858877899162360138844032265418028567),222 BOOST_MATH_BIG_CONSTANT(T, 113, 0.20499222604410032375789018837922397),223 BOOST_MATH_BIG_CONSTANT(T, 113, 0.0272103140348194747360175268778415049),224 BOOST_MATH_BIG_CONSTANT(T, 113, 0.0015764849020876949848954081173520686),225 };226 static const T Q_1_2[] = {227 BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),228 BOOST_MATH_BIG_CONSTANT(T, 113, 4.71237311120863419878375031457715223),229 BOOST_MATH_BIG_CONSTANT(T, 113, 9.58619118655339853449127952145877467),230 BOOST_MATH_BIG_CONSTANT(T, 113, 11.0940067269829372437561421279054968),231 BOOST_MATH_BIG_CONSTANT(T, 113, 8.09075424749327792073276309969037885),232 BOOST_MATH_BIG_CONSTANT(T, 113, 3.87705890159891405185343806884451286),233 BOOST_MATH_BIG_CONSTANT(T, 113, 1.22758678701914477836330837816976782),234 BOOST_MATH_BIG_CONSTANT(T, 113, 0.249092040606385004109672077814668716),235 BOOST_MATH_BIG_CONSTANT(T, 113, 0.0295750413900655597027079600025569048),236 BOOST_MATH_BIG_CONSTANT(T, 113, 0.00157648490200498142247694709728858139),237 BOOST_MATH_BIG_CONSTANT(T, 113, 0.161264050344059471721062360645432809e-14),238 };239 240 // Max error in interpolated form: 8.958e-035241 static const T P_2_4[] = {242 BOOST_MATH_BIG_CONSTANT(T, 113, -2.55843734739907925764326773972215085),243 BOOST_MATH_BIG_CONSTANT(T, 113, -12.2830208240542011967952466273455887),244 BOOST_MATH_BIG_CONSTANT(T, 113, -23.9195022162767993526575786066414403),245 BOOST_MATH_BIG_CONSTANT(T, 113, -24.9256431504823483094158828285470862),246 BOOST_MATH_BIG_CONSTANT(T, 113, -14.7979122765478779075108064826412285),247 BOOST_MATH_BIG_CONSTANT(T, 113, -4.46654453928610666393276765059122272),248 BOOST_MATH_BIG_CONSTANT(T, 113, -0.0191439033405649675717082465687845002),249 BOOST_MATH_BIG_CONSTANT(T, 113, 0.515412052554351265708917209749037352),250 BOOST_MATH_BIG_CONSTANT(T, 113, 0.195378348786064304378247325360320038),251 BOOST_MATH_BIG_CONSTANT(T, 113, 0.0334761282624174313035014426794245393),252 BOOST_MATH_BIG_CONSTANT(T, 113, 0.002373665205942206348500250056602687),253 };254 static const T Q_2_4[] = {255 BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),256 BOOST_MATH_BIG_CONSTANT(T, 113, 4.80098558454419907830670928248659245),257 BOOST_MATH_BIG_CONSTANT(T, 113, 9.99220727843170133895059300223445265),258 BOOST_MATH_BIG_CONSTANT(T, 113, 11.8896146167631330735386697123464976),259 BOOST_MATH_BIG_CONSTANT(T, 113, 8.96613256683809091593793565879092581),260 BOOST_MATH_BIG_CONSTANT(T, 113, 4.47254136149624110878909334574485751),261 BOOST_MATH_BIG_CONSTANT(T, 113, 1.48600982028196527372434773913633152),262 BOOST_MATH_BIG_CONSTANT(T, 113, 0.319570735766764237068541501137990078),263 BOOST_MATH_BIG_CONSTANT(T, 113, 0.0407358345787680953107374215319322066),264 BOOST_MATH_BIG_CONSTANT(T, 113, 0.00237366520593271641375755486420859837),265 BOOST_MATH_BIG_CONSTANT(T, 113, 0.239554887903526152679337256236302116e-15),266 BOOST_MATH_BIG_CONSTANT(T, 113, -0.294749244740618656265237072002026314e-17),267 };268 269 static const T y_offset_2_4 = BOOST_MATH_BIG_CONSTANT(T, 113, 3.558437347412109375);270 271 // Max error in interpolated form: 4.319e-035272 static const T P_4_8[] = {273 BOOST_MATH_BIG_CONSTANT(T, 113, 0.166626112697021464248967707021688845e-16),274 BOOST_MATH_BIG_CONSTANT(T, 113, 0.499999999999997739552090249208808197),275 BOOST_MATH_BIG_CONSTANT(T, 113, 6.40270945019053817915772473771553187),276 BOOST_MATH_BIG_CONSTANT(T, 113, 41.3833374155000608013677627389343329),277 BOOST_MATH_BIG_CONSTANT(T, 113, 166.803341854562809335667241074035245),278 BOOST_MATH_BIG_CONSTANT(T, 113, 453.39964786925369319960722793414521),279 BOOST_MATH_BIG_CONSTANT(T, 113, 851.153712317697055375935433362983944),280 BOOST_MATH_BIG_CONSTANT(T, 113, 1097.70657567285059133109286478004458),281 BOOST_MATH_BIG_CONSTANT(T, 113, 938.431232478455316020076349367632922),282 BOOST_MATH_BIG_CONSTANT(T, 113, 487.268001604651932322080970189930074),283 BOOST_MATH_BIG_CONSTANT(T, 113, 119.953445242335730062471193124820659),284 };285 static const T Q_4_8[] = {286 BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),287 BOOST_MATH_BIG_CONSTANT(T, 113, 12.4720855670474488978638945855932398),288 BOOST_MATH_BIG_CONSTANT(T, 113, 78.6093129753298570701376952709727391),289 BOOST_MATH_BIG_CONSTANT(T, 113, 307.470246050318322489781182863190127),290 BOOST_MATH_BIG_CONSTANT(T, 113, 805.140686101151538537565264188630079),291 BOOST_MATH_BIG_CONSTANT(T, 113, 1439.12019760292146454787601409644413),292 BOOST_MATH_BIG_CONSTANT(T, 113, 1735.6105285756048831268586001383127),293 BOOST_MATH_BIG_CONSTANT(T, 113, 1348.32500712856328019355198611280536),294 BOOST_MATH_BIG_CONSTANT(T, 113, 607.225985860570846699704222144650563),295 BOOST_MATH_BIG_CONSTANT(T, 113, 119.952317857277045332558673164517227),296 BOOST_MATH_BIG_CONSTANT(T, 113, 0.000140165918355036060868680809129436084),297 };298 299 // Maximum Deviation Found: 2.867e-035300 // Expected Error Term : 2.866e-035301 // Maximum Relative Change in Control Points : 2.662e-004302 static const T P_8_16[] = {303 BOOST_MATH_BIG_CONSTANT(T, 113, -0.184828315274146610610872315609837439e-19),304 BOOST_MATH_BIG_CONSTANT(T, 113, 0.500000000000000004122475157735807738),305 BOOST_MATH_BIG_CONSTANT(T, 113, 3.02533865247313349284875558880415875),306 BOOST_MATH_BIG_CONSTANT(T, 113, 13.5995927517457371243039532492642734),307 BOOST_MATH_BIG_CONSTANT(T, 113, 35.3132224283087906757037999452941588),308 BOOST_MATH_BIG_CONSTANT(T, 113, 67.1639424550714159157603179911505619),309 BOOST_MATH_BIG_CONSTANT(T, 113, 83.5767733658513967581959839367419891),310 BOOST_MATH_BIG_CONSTANT(T, 113, 71.073491212235705900866411319363501),311 BOOST_MATH_BIG_CONSTANT(T, 113, 35.8621515614725564575893663483998663),312 BOOST_MATH_BIG_CONSTANT(T, 113, 8.72152231639983491987779743154333318),313 };314 static const T Q_8_16[] = {315 BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),316 BOOST_MATH_BIG_CONSTANT(T, 113, 5.71734397161293452310624822415866372),317 BOOST_MATH_BIG_CONSTANT(T, 113, 25.293404179620438179337103263274815),318 BOOST_MATH_BIG_CONSTANT(T, 113, 62.2619767967468199111077640625328469),319 BOOST_MATH_BIG_CONSTANT(T, 113, 113.955048909238993473389714972250235),320 BOOST_MATH_BIG_CONSTANT(T, 113, 130.807138328938966981862203944329408),321 BOOST_MATH_BIG_CONSTANT(T, 113, 102.423146902337654110717764213057753),322 BOOST_MATH_BIG_CONSTANT(T, 113, 44.0424772805245202514468199602123565),323 BOOST_MATH_BIG_CONSTANT(T, 113, 8.89898032477904072082994913461386099),324 BOOST_MATH_BIG_CONSTANT(T, 113, -0.0296627336872039988632793863671456398),325 };326 // Maximum Deviation Found: 1.079e-035327 // Expected Error Term : -1.079e-035328 // Maximum Relative Change in Control Points : 7.884e-003329 static const T P_16_inf[] = {330 BOOST_MATH_BIG_CONSTANT(T, 113, 0.0),331 BOOST_MATH_BIG_CONSTANT(T, 113, 0.500000000000000000000000000000087317),332 BOOST_MATH_BIG_CONSTANT(T, 113, 0.345625669885456215194494735902663968),333 BOOST_MATH_BIG_CONSTANT(T, 113, 9.62895499360842232127552650044647769),334 BOOST_MATH_BIG_CONSTANT(T, 113, 3.5936085382439026269301003761320812),335 BOOST_MATH_BIG_CONSTANT(T, 113, 49.459599118438883265036646019410669),336 BOOST_MATH_BIG_CONSTANT(T, 113, 7.77519237321893917784735690560496607),337 BOOST_MATH_BIG_CONSTANT(T, 113, 74.4536074488178075948642351179304121),338 BOOST_MATH_BIG_CONSTANT(T, 113, 2.75209340397069050436806159297952699),339 BOOST_MATH_BIG_CONSTANT(T, 113, 23.9292359711471667884504840186561598),340 };341 static const T Q_16_inf[] = {342 BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),343 BOOST_MATH_BIG_CONSTANT(T, 113, 0.357918006437579097055656138920742037),344 BOOST_MATH_BIG_CONSTANT(T, 113, 19.1386039850709849435325005484512944),345 BOOST_MATH_BIG_CONSTANT(T, 113, 0.874349081464143606016221431763364517),346 BOOST_MATH_BIG_CONSTANT(T, 113, 98.6516097434855572678195488061432509),347 BOOST_MATH_BIG_CONSTANT(T, 113, -16.1051972833382893468655223662534306),348 BOOST_MATH_BIG_CONSTANT(T, 113, 154.316860216253720989145047141653727),349 BOOST_MATH_BIG_CONSTANT(T, 113, -40.2026880424378986053105969312264534),350 BOOST_MATH_BIG_CONSTANT(T, 113, 60.1679136674264778074736441126810223),351 BOOST_MATH_BIG_CONSTANT(T, 113, -13.3414844622256422644504472438320114),352 BOOST_MATH_BIG_CONSTANT(T, 113, 2.53795636200649908779512969030363442),353 };354 355 if(x <= 2)356 {357 return (2 + boost::math::tools::evaluate_polynomial(P_1_2, x) / tools::evaluate_polynomial(Q_1_2, x)) / (x * x);358 }359 else if(x <= 4)360 {361 return (y_offset_2_4 + boost::math::tools::evaluate_polynomial(P_2_4, x) / tools::evaluate_polynomial(Q_2_4, x)) / (x * x);362 }363 else if(x <= 8)364 {365 T y = 1 / x;366 return (1 + tools::evaluate_polynomial(P_4_8, y) / tools::evaluate_polynomial(Q_4_8, y)) / x;367 }368 else if(x <= 16)369 {370 T y = 1 / x;371 return (1 + tools::evaluate_polynomial(P_8_16, y) / tools::evaluate_polynomial(Q_8_16, y)) / x;372 }373 T y = 1 / x;374 return (1 + tools::evaluate_polynomial(P_16_inf, y) / tools::evaluate_polynomial(Q_16_inf, y)) / x;375}376 377template <class T, class Policy, class Tag>378BOOST_MATH_GPU_ENABLED T trigamma_dispatch(T x, const Policy& pol, const Tag& tag)379{380 //381 // This handles reflection of negative arguments, and all our382 // error handling, then forwards to the T-specific approximation.383 //384 BOOST_MATH_STD_USING // ADL of std functions.385 386 T result = 0;387 //388 // Check for negative arguments and use reflection:389 //390 if(x <= 0)391 {392 // Reflect:393 T z = 1 - x;394 395 BOOST_MATH_ASSERT(z >= 1);396 397 // Argument reduction for tan:398 if(floor(x) == x)399 {400 return policies::raise_pole_error<T>("boost::math::trigamma<%1%>(%1%)", nullptr, (1-x), pol);401 }402 T s = fabs(x) < fabs(z) ? boost::math::sin_pi(x, pol) : boost::math::sin_pi(z, pol);403 return result - trigamma_prec(T(z), pol, tag) + boost::math::pow<2>(constants::pi<T>()) / (s * s);404 }405 if(x < 1)406 {407 result = 1 / (x * x);408 x += 1;409 }410 return result + trigamma_prec(x, pol, tag);411}412 413} // namespace detail414 415template <class T, class Policy>416BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T>::type417 trigamma(T x, const Policy&)418{419 typedef typename tools::promote_args<T>::type result_type;420 typedef typename policies::evaluation<result_type, Policy>::type value_type;421 typedef typename policies::precision<T, Policy>::type precision_type;422 typedef boost::math::integral_constant<int,423 precision_type::value <= 0 ? 0 :424 precision_type::value <= 53 ? 53 :425 precision_type::value <= 64 ? 64 :426 precision_type::value <= 113 ? 113 : 0427 > tag_type;428 typedef typename policies::normalise<429 Policy,430 policies::promote_float<false>,431 policies::promote_double<false>,432 policies::discrete_quantile<>,433 policies::assert_undefined<> >::type forwarding_policy;434 435 return policies::checked_narrowing_cast<result_type, Policy>(detail::trigamma_dispatch(static_cast<value_type>(x), forwarding_policy(), tag_type()), "boost::math::trigamma<%1%>(%1%)");436}437 438template <class T>439BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T>::type440 trigamma(T x)441{442 return trigamma(x, policies::policy<>());443}444 445} // namespace math446} // namespace boost447#endif448 449