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