trigamma.hpp 21 KB

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