big_digamma.hpp 20 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295
  1. // (C) Copyright John Maddock 2006-8.
  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_NTL_DIGAMMA
  6. #define BOOST_MATH_NTL_DIGAMMA
  7. #include <boost/math/tools/rational.hpp>
  8. #include <boost/math/policies/error_handling.hpp>
  9. #include <boost/math/constants/constants.hpp>
  10. #include <boost/lexical_cast.hpp>
  11. namespace boost{ namespace math{ namespace detail{
  12. template <class T>
  13. T big_digamma_helper(T x)
  14. {
  15. static const T P[61] = {
  16. boost::lexical_cast<T>("0.6660133691143982067148122682345055274952e81"),
  17. boost::lexical_cast<T>("0.6365271516829242456324234577164675383137e81"),
  18. boost::lexical_cast<T>("0.2991038873096202943405966144203628966976e81"),
  19. boost::lexical_cast<T>("0.9211116495503170498076013367421231351115e80"),
  20. boost::lexical_cast<T>("0.2090792764676090716286400360584443891749e80"),
  21. boost::lexical_cast<T>("0.3730037777359591428226035156377978092809e79"),
  22. boost::lexical_cast<T>("0.5446396536956682043376492370432031543834e78"),
  23. boost::lexical_cast<T>("0.6692523966335177847425047827449069256345e77"),
  24. boost::lexical_cast<T>("0.7062543624100864681625612653756619116848e76"),
  25. boost::lexical_cast<T>("0.6499914905966283735005256964443226879158e75"),
  26. boost::lexical_cast<T>("0.5280364564853225211197557708655426736091e74"),
  27. boost::lexical_cast<T>("0.3823205608981176913075543599005095206953e73"),
  28. boost::lexical_cast<T>("0.2486733714214237704739129972671154532415e72"),
  29. boost::lexical_cast<T>("0.1462562139602039577983434547171318011675e71"),
  30. boost::lexical_cast<T>("0.7821169065036815012381267259559910324285e69"),
  31. boost::lexical_cast<T>("0.3820552182348155468636157988764435365078e68"),
  32. boost::lexical_cast<T>("0.1711618296983598244658239925535632505062e67"),
  33. boost::lexical_cast<T>("0.7056661618357643731419080738521475204245e65"),
  34. boost::lexical_cast<T>("0.2685246896473614017356264531791459936036e64"),
  35. boost::lexical_cast<T>("0.9455168125599643085283071944864977592391e62"),
  36. boost::lexical_cast<T>("0.3087541626972538362237309145177486236219e61"),
  37. boost::lexical_cast<T>("0.9367928873352980208052601301625005737407e59"),
  38. boost::lexical_cast<T>("0.2645306130689794942883818547314327466007e58"),
  39. boost::lexical_cast<T>("0.6961815141171454309161007351079576190079e56"),
  40. boost::lexical_cast<T>("0.1709637824471794552313802669803885946843e55"),
  41. boost::lexical_cast<T>("0.3921553258481531526663112728778759311158e53"),
  42. boost::lexical_cast<T>("0.8409006354449988687714450897575728228696e51"),
  43. boost::lexical_cast<T>("0.1686755204461325935742097669030363344927e50"),
  44. boost::lexical_cast<T>("0.3166653542877070999007425197729038754254e48"),
  45. boost::lexical_cast<T>("0.5566029092358215049069560272835654229637e46"),
  46. boost::lexical_cast<T>("0.9161766287916328133080586672953875116242e44"),
  47. boost::lexical_cast<T>("1412317772330871298317974693514430627922000"),
  48. boost::lexical_cast<T>("20387991986727877473732570146112459874790"),
  49. boost::lexical_cast<T>("275557928713904105182512535678580359839.3"),
  50. boost::lexical_cast<T>("3485719851040516559072031256589598330.723"),
  51. boost::lexical_cast<T>("41247046743564028399938106707656877.40859"),
  52. boost::lexical_cast<T>("456274078125709314602601667471879.0147312"),
  53. boost::lexical_cast<T>("4714450683242899367025707077155.310613012"),
  54. boost::lexical_cast<T>("45453933537925041680009544258.75073849996"),
  55. boost::lexical_cast<T>("408437900487067278846361972.302331241052"),
  56. boost::lexical_cast<T>("3415719344386166273085838.705771571751035"),
  57. boost::lexical_cast<T>("26541502879185876562320.93134691487351145"),
  58. boost::lexical_cast<T>("191261415065918713661.1571433274648417668"),
  59. boost::lexical_cast<T>("1275349770108718421.645275944284937551702"),
  60. boost::lexical_cast<T>("7849171120971773.318910987434906905704272"),
  61. boost::lexical_cast<T>("44455946386549.80866460312682983576538056"),
  62. boost::lexical_cast<T>("230920362395.3198137186361608905136598046"),
  63. boost::lexical_cast<T>("1095700096.240863858624279930600654130254"),
  64. boost::lexical_cast<T>("4727085.467506050153744334085516289728134"),
  65. boost::lexical_cast<T>("18440.75118859447173303252421991479005424"),
  66. boost::lexical_cast<T>("64.62515887799460295677071749181651317052"),
  67. boost::lexical_cast<T>("0.201851568864688406206528472883512147547"),
  68. boost::lexical_cast<T>("0.0005565091674187978029138500039504078098143"),
  69. boost::lexical_cast<T>("0.1338097668312907986354698683493366559613e-5"),
  70. boost::lexical_cast<T>("0.276308225077464312820179030238305271638e-8"),
  71. boost::lexical_cast<T>("0.4801582970473168520375942100071070575043e-11"),
  72. boost::lexical_cast<T>("0.6829184144212920949740376186058541800175e-14"),
  73. boost::lexical_cast<T>("0.7634080076358511276617829524639455399182e-17"),
  74. boost::lexical_cast<T>("0.6290035083727140966418512608156646142409e-20"),
  75. boost::lexical_cast<T>("0.339652245667538733044036638506893821352e-23"),
  76. boost::lexical_cast<T>("0.9017518064256388530773585529891677854909e-27")
  77. };
  78. static const T Q[61] = {
  79. boost::lexical_cast<T>("0"),
  80. boost::lexical_cast<T>("0.1386831185456898357379390197203894063459e81"),
  81. boost::lexical_cast<T>("0.6467076379487574703291056110838151259438e81"),
  82. boost::lexical_cast<T>("0.1394967823848615838336194279565285465161e82"),
  83. boost::lexical_cast<T>("0.1872927317344192945218570366455046340458e82"),
  84. boost::lexical_cast<T>("0.1772461045338946243584650759986310355937e82"),
  85. boost::lexical_cast<T>("0.1267294892200258648315971144069595555118e82"),
  86. boost::lexical_cast<T>("0.7157764838362416821508872117623058626589e81"),
  87. boost::lexical_cast<T>("0.329447266909948668265277828268378274513e81"),
  88. boost::lexical_cast<T>("0.1264376077317689779509250183194342571207e81"),
  89. boost::lexical_cast<T>("0.4118230304191980787640446056583623228873e80"),
  90. boost::lexical_cast<T>("0.1154393529762694616405952270558316515261e80"),
  91. boost::lexical_cast<T>("0.281655612889423906125295485693696744275e79"),
  92. boost::lexical_cast<T>("0.6037483524928743102724159846414025482077e78"),
  93. boost::lexical_cast<T>("0.1145927995397835468123576831800276999614e78"),
  94. boost::lexical_cast<T>("0.1938624296151985600348534009382865995154e77"),
  95. boost::lexical_cast<T>("0.293980925856227626211879961219188406675e76"),
  96. boost::lexical_cast<T>("0.4015574518216966910319562902099567437832e75"),
  97. boost::lexical_cast<T>("0.4961475457509727343545565970423431880907e74"),
  98. boost::lexical_cast<T>("0.5565482348278933960215521991000378896338e73"),
  99. boost::lexical_cast<T>("0.5686112924615820754631098622770303094938e72"),
  100. boost::lexical_cast<T>("0.5305988545844796293285410303747469932856e71"),
  101. boost::lexical_cast<T>("0.4533363413802585060568537458067343491358e70"),
  102. boost::lexical_cast<T>("0.3553932059473516064068322757331575565718e69"),
  103. boost::lexical_cast<T>("0.2561198565218704414618802902533972354203e68"),
  104. boost::lexical_cast<T>("0.1699519313292900324098102065697454295572e67"),
  105. boost::lexical_cast<T>("0.1039830160862334505389615281373574959236e66"),
  106. boost::lexical_cast<T>("0.5873082967977428281000961954715372504986e64"),
  107. boost::lexical_cast<T>("0.3065255179030575882202133042549783442446e63"),
  108. boost::lexical_cast<T>("0.1479494813481364701208655943688307245459e62"),
  109. boost::lexical_cast<T>("0.6608150467921598615495180659808895663164e60"),
  110. boost::lexical_cast<T>("0.2732535313770902021791888953487787496976e59"),
  111. boost::lexical_cast<T>("0.1046402297662493314531194338414508049069e58"),
  112. boost::lexical_cast<T>("0.3711375077192882936085049147920021549622e56"),
  113. boost::lexical_cast<T>("0.1219154482883895482637944309702972234576e55"),
  114. boost::lexical_cast<T>("0.3708359374149458741391374452286837880162e53"),
  115. boost::lexical_cast<T>("0.1044095509971707189716913168889769471468e52"),
  116. boost::lexical_cast<T>("0.271951506225063286130946773813524945052e50"),
  117. boost::lexical_cast<T>("0.6548016291215163843464133978454065823866e48"),
  118. boost::lexical_cast<T>("0.1456062447610542135403751730809295219344e47"),
  119. boost::lexical_cast<T>("0.2986690175077969760978388356833006028929e45"),
  120. boost::lexical_cast<T>("5643149706574013350061247429006443326844000"),
  121. boost::lexical_cast<T>("98047545414467090421964387960743688053480"),
  122. boost::lexical_cast<T>("1563378767746846395507385099301468978550"),
  123. boost::lexical_cast<T>("22823360528584500077862274918382796495"),
  124. boost::lexical_cast<T>("304215527004115213046601295970388750"),
  125. boost::lexical_cast<T>("3690289075895685793844344966820325"),
  126. boost::lexical_cast<T>("40584512015702371433911456606050"),
  127. boost::lexical_cast<T>("402834190897282802772754873905"),
  128. boost::lexical_cast<T>("3589522158493606918146495750"),
  129. boost::lexical_cast<T>("28530557707503483723634725"),
  130. boost::lexical_cast<T>("200714561335055753000730"),
  131. boost::lexical_cast<T>("1237953783437761888641"),
  132. boost::lexical_cast<T>("6614698701445762950"),
  133. boost::lexical_cast<T>("30155495647727505"),
  134. boost::lexical_cast<T>("114953256021450"),
  135. boost::lexical_cast<T>("356398020013"),
  136. boost::lexical_cast<T>("863113950"),
  137. boost::lexical_cast<T>("1531345"),
  138. boost::lexical_cast<T>("1770"),
  139. boost::lexical_cast<T>("1")
  140. };
  141. static const T PD[60] = {
  142. boost::lexical_cast<T>("0.6365271516829242456324234577164675383137e81"),
  143. 2*boost::lexical_cast<T>("0.2991038873096202943405966144203628966976e81"),
  144. 3*boost::lexical_cast<T>("0.9211116495503170498076013367421231351115e80"),
  145. 4*boost::lexical_cast<T>("0.2090792764676090716286400360584443891749e80"),
  146. 5*boost::lexical_cast<T>("0.3730037777359591428226035156377978092809e79"),
  147. 6*boost::lexical_cast<T>("0.5446396536956682043376492370432031543834e78"),
  148. 7*boost::lexical_cast<T>("0.6692523966335177847425047827449069256345e77"),
  149. 8*boost::lexical_cast<T>("0.7062543624100864681625612653756619116848e76"),
  150. 9*boost::lexical_cast<T>("0.6499914905966283735005256964443226879158e75"),
  151. 10*boost::lexical_cast<T>("0.5280364564853225211197557708655426736091e74"),
  152. 11*boost::lexical_cast<T>("0.3823205608981176913075543599005095206953e73"),
  153. 12*boost::lexical_cast<T>("0.2486733714214237704739129972671154532415e72"),
  154. 13*boost::lexical_cast<T>("0.1462562139602039577983434547171318011675e71"),
  155. 14*boost::lexical_cast<T>("0.7821169065036815012381267259559910324285e69"),
  156. 15*boost::lexical_cast<T>("0.3820552182348155468636157988764435365078e68"),
  157. 16*boost::lexical_cast<T>("0.1711618296983598244658239925535632505062e67"),
  158. 17*boost::lexical_cast<T>("0.7056661618357643731419080738521475204245e65"),
  159. 18*boost::lexical_cast<T>("0.2685246896473614017356264531791459936036e64"),
  160. 19*boost::lexical_cast<T>("0.9455168125599643085283071944864977592391e62"),
  161. 20*boost::lexical_cast<T>("0.3087541626972538362237309145177486236219e61"),
  162. 21*boost::lexical_cast<T>("0.9367928873352980208052601301625005737407e59"),
  163. 22*boost::lexical_cast<T>("0.2645306130689794942883818547314327466007e58"),
  164. 23*boost::lexical_cast<T>("0.6961815141171454309161007351079576190079e56"),
  165. 24*boost::lexical_cast<T>("0.1709637824471794552313802669803885946843e55"),
  166. 25*boost::lexical_cast<T>("0.3921553258481531526663112728778759311158e53"),
  167. 26*boost::lexical_cast<T>("0.8409006354449988687714450897575728228696e51"),
  168. 27*boost::lexical_cast<T>("0.1686755204461325935742097669030363344927e50"),
  169. 28*boost::lexical_cast<T>("0.3166653542877070999007425197729038754254e48"),
  170. 29*boost::lexical_cast<T>("0.5566029092358215049069560272835654229637e46"),
  171. 30*boost::lexical_cast<T>("0.9161766287916328133080586672953875116242e44"),
  172. 31*boost::lexical_cast<T>("1412317772330871298317974693514430627922000"),
  173. 32*boost::lexical_cast<T>("20387991986727877473732570146112459874790"),
  174. 33*boost::lexical_cast<T>("275557928713904105182512535678580359839.3"),
  175. 34*boost::lexical_cast<T>("3485719851040516559072031256589598330.723"),
  176. 35*boost::lexical_cast<T>("41247046743564028399938106707656877.40859"),
  177. 36*boost::lexical_cast<T>("456274078125709314602601667471879.0147312"),
  178. 37*boost::lexical_cast<T>("4714450683242899367025707077155.310613012"),
  179. 38*boost::lexical_cast<T>("45453933537925041680009544258.75073849996"),
  180. 39*boost::lexical_cast<T>("408437900487067278846361972.302331241052"),
  181. 40*boost::lexical_cast<T>("3415719344386166273085838.705771571751035"),
  182. 41*boost::lexical_cast<T>("26541502879185876562320.93134691487351145"),
  183. 42*boost::lexical_cast<T>("191261415065918713661.1571433274648417668"),
  184. 43*boost::lexical_cast<T>("1275349770108718421.645275944284937551702"),
  185. 44*boost::lexical_cast<T>("7849171120971773.318910987434906905704272"),
  186. 45*boost::lexical_cast<T>("44455946386549.80866460312682983576538056"),
  187. 46*boost::lexical_cast<T>("230920362395.3198137186361608905136598046"),
  188. 47*boost::lexical_cast<T>("1095700096.240863858624279930600654130254"),
  189. 48*boost::lexical_cast<T>("4727085.467506050153744334085516289728134"),
  190. 49*boost::lexical_cast<T>("18440.75118859447173303252421991479005424"),
  191. 50*boost::lexical_cast<T>("64.62515887799460295677071749181651317052"),
  192. 51*boost::lexical_cast<T>("0.201851568864688406206528472883512147547"),
  193. 52*boost::lexical_cast<T>("0.0005565091674187978029138500039504078098143"),
  194. 53*boost::lexical_cast<T>("0.1338097668312907986354698683493366559613e-5"),
  195. 54*boost::lexical_cast<T>("0.276308225077464312820179030238305271638e-8"),
  196. 55*boost::lexical_cast<T>("0.4801582970473168520375942100071070575043e-11"),
  197. 56*boost::lexical_cast<T>("0.6829184144212920949740376186058541800175e-14"),
  198. 57*boost::lexical_cast<T>("0.7634080076358511276617829524639455399182e-17"),
  199. 58*boost::lexical_cast<T>("0.6290035083727140966418512608156646142409e-20"),
  200. 59*boost::lexical_cast<T>("0.339652245667538733044036638506893821352e-23"),
  201. 60*boost::lexical_cast<T>("0.9017518064256388530773585529891677854909e-27")
  202. };
  203. static const T QD[60] = {
  204. boost::lexical_cast<T>("0.1386831185456898357379390197203894063459e81"),
  205. 2*boost::lexical_cast<T>("0.6467076379487574703291056110838151259438e81"),
  206. 3*boost::lexical_cast<T>("0.1394967823848615838336194279565285465161e82"),
  207. 4*boost::lexical_cast<T>("0.1872927317344192945218570366455046340458e82"),
  208. 5*boost::lexical_cast<T>("0.1772461045338946243584650759986310355937e82"),
  209. 6*boost::lexical_cast<T>("0.1267294892200258648315971144069595555118e82"),
  210. 7*boost::lexical_cast<T>("0.7157764838362416821508872117623058626589e81"),
  211. 8*boost::lexical_cast<T>("0.329447266909948668265277828268378274513e81"),
  212. 9*boost::lexical_cast<T>("0.1264376077317689779509250183194342571207e81"),
  213. 10*boost::lexical_cast<T>("0.4118230304191980787640446056583623228873e80"),
  214. 11*boost::lexical_cast<T>("0.1154393529762694616405952270558316515261e80"),
  215. 12*boost::lexical_cast<T>("0.281655612889423906125295485693696744275e79"),
  216. 13*boost::lexical_cast<T>("0.6037483524928743102724159846414025482077e78"),
  217. 14*boost::lexical_cast<T>("0.1145927995397835468123576831800276999614e78"),
  218. 15*boost::lexical_cast<T>("0.1938624296151985600348534009382865995154e77"),
  219. 16*boost::lexical_cast<T>("0.293980925856227626211879961219188406675e76"),
  220. 17*boost::lexical_cast<T>("0.4015574518216966910319562902099567437832e75"),
  221. 18*boost::lexical_cast<T>("0.4961475457509727343545565970423431880907e74"),
  222. 19*boost::lexical_cast<T>("0.5565482348278933960215521991000378896338e73"),
  223. 20*boost::lexical_cast<T>("0.5686112924615820754631098622770303094938e72"),
  224. 21*boost::lexical_cast<T>("0.5305988545844796293285410303747469932856e71"),
  225. 22*boost::lexical_cast<T>("0.4533363413802585060568537458067343491358e70"),
  226. 23*boost::lexical_cast<T>("0.3553932059473516064068322757331575565718e69"),
  227. 24*boost::lexical_cast<T>("0.2561198565218704414618802902533972354203e68"),
  228. 25*boost::lexical_cast<T>("0.1699519313292900324098102065697454295572e67"),
  229. 26*boost::lexical_cast<T>("0.1039830160862334505389615281373574959236e66"),
  230. 27*boost::lexical_cast<T>("0.5873082967977428281000961954715372504986e64"),
  231. 28*boost::lexical_cast<T>("0.3065255179030575882202133042549783442446e63"),
  232. 29*boost::lexical_cast<T>("0.1479494813481364701208655943688307245459e62"),
  233. 30*boost::lexical_cast<T>("0.6608150467921598615495180659808895663164e60"),
  234. 31*boost::lexical_cast<T>("0.2732535313770902021791888953487787496976e59"),
  235. 32*boost::lexical_cast<T>("0.1046402297662493314531194338414508049069e58"),
  236. 33*boost::lexical_cast<T>("0.3711375077192882936085049147920021549622e56"),
  237. 34*boost::lexical_cast<T>("0.1219154482883895482637944309702972234576e55"),
  238. 35*boost::lexical_cast<T>("0.3708359374149458741391374452286837880162e53"),
  239. 36*boost::lexical_cast<T>("0.1044095509971707189716913168889769471468e52"),
  240. 37*boost::lexical_cast<T>("0.271951506225063286130946773813524945052e50"),
  241. 38*boost::lexical_cast<T>("0.6548016291215163843464133978454065823866e48"),
  242. 39*boost::lexical_cast<T>("0.1456062447610542135403751730809295219344e47"),
  243. 40*boost::lexical_cast<T>("0.2986690175077969760978388356833006028929e45"),
  244. 41*boost::lexical_cast<T>("5643149706574013350061247429006443326844000"),
  245. 42*boost::lexical_cast<T>("98047545414467090421964387960743688053480"),
  246. 43*boost::lexical_cast<T>("1563378767746846395507385099301468978550"),
  247. 44*boost::lexical_cast<T>("22823360528584500077862274918382796495"),
  248. 45*boost::lexical_cast<T>("304215527004115213046601295970388750"),
  249. 46*boost::lexical_cast<T>("3690289075895685793844344966820325"),
  250. 47*boost::lexical_cast<T>("40584512015702371433911456606050"),
  251. 48*boost::lexical_cast<T>("402834190897282802772754873905"),
  252. 49*boost::lexical_cast<T>("3589522158493606918146495750"),
  253. 50*boost::lexical_cast<T>("28530557707503483723634725"),
  254. 51*boost::lexical_cast<T>("200714561335055753000730"),
  255. 52*boost::lexical_cast<T>("1237953783437761888641"),
  256. 53*boost::lexical_cast<T>("6614698701445762950"),
  257. 54*boost::lexical_cast<T>("30155495647727505"),
  258. 55*boost::lexical_cast<T>("114953256021450"),
  259. 56*boost::lexical_cast<T>("356398020013"),
  260. 57*boost::lexical_cast<T>("863113950"),
  261. 58*boost::lexical_cast<T>("1531345"),
  262. 59*boost::lexical_cast<T>("1770"),
  263. 60*boost::lexical_cast<T>("1")
  264. };
  265. static const double g = 63.192152;
  266. T zgh = x + g - 0.5;
  267. T result = (x - 0.5) / zgh;
  268. result += log(zgh);
  269. result += tools::evaluate_polynomial(PD, x) / tools::evaluate_polynomial(P, x);
  270. result -= tools::evaluate_polynomial(QD, x) / tools::evaluate_polynomial(Q, x);
  271. result -= 1;
  272. return result;
  273. }
  274. template <class T>
  275. T big_digamma(T x)
  276. {
  277. BOOST_MATH_STD_USING
  278. if(x < 0)
  279. {
  280. return big_digamma_helper(static_cast<T>(1-x)) + constants::pi<T>() / tan(constants::pi<T>() * (1-x));
  281. }
  282. return big_digamma_helper(x);
  283. }
  284. }}}
  285. #endif // include guard