igamma_large.hpp 40 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824
  1. // Copyright John Maddock 2006.
  2. // Copyright Matt Borland 2024.
  3. // Use, modification and distribution are subject to the
  4. // Boost Software License, Version 1.0. (See accompanying file
  5. // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
  6. //
  7. // This file implements the asymptotic expansions of the incomplete
  8. // gamma functions P(a, x) and Q(a, x), used when a is large and
  9. // x ~ a.
  10. //
  11. // The primary reference is:
  12. //
  13. // "The Asymptotic Expansion of the Incomplete Gamma Functions"
  14. // N. M. Temme.
  15. // Siam J. Math Anal. Vol 10 No 4, July 1979, p757.
  16. //
  17. // A different way of evaluating these expansions,
  18. // plus a lot of very useful background information is in:
  19. //
  20. // "A Set of Algorithms For the Incomplete Gamma Functions."
  21. // N. M. Temme.
  22. // Probability in the Engineering and Informational Sciences,
  23. // 8, 1994, 291.
  24. //
  25. // An alternative implementation is in:
  26. //
  27. // "Computation of the Incomplete Gamma Function Ratios and their Inverse."
  28. // A. R. Didonato and A. H. Morris.
  29. // ACM TOMS, Vol 12, No 4, Dec 1986, p377.
  30. //
  31. // There are various versions of the same code below, each accurate
  32. // to a different precision. To understand the code, refer to Didonato
  33. // and Morris, from Eq 17 and 18 onwards.
  34. //
  35. // The coefficients used here are not taken from Didonato and Morris:
  36. // the domain over which these expansions are used is slightly different
  37. // to theirs, and their constants are not quite accurate enough for
  38. // 128-bit long double's. Instead the coefficients were calculated
  39. // using the methods described by Temme p762 from Eq 3.8 onwards.
  40. // The values obtained agree with those obtained by Didonato and Morris
  41. // (at least to the first 30 digits that they provide).
  42. // At double precision the degrees of polynomial required for full
  43. // machine precision are close to those recommended to Didonato and Morris,
  44. // but of course many more terms are needed for larger types.
  45. //
  46. #ifndef BOOST_MATH_DETAIL_IGAMMA_LARGE
  47. #define BOOST_MATH_DETAIL_IGAMMA_LARGE
  48. #ifdef _MSC_VER
  49. #pragma once
  50. #endif
  51. #if defined(__GNUC__) && defined(BOOST_MATH_USE_FLOAT128)
  52. //
  53. // This is the only way we can avoid
  54. // warning: non-standard suffix on floating constant [-Wpedantic]
  55. // when building with -Wall -pedantic. Neither __extension__
  56. // nor #pragma diagnostic ignored work :(
  57. //
  58. #pragma GCC system_header
  59. #endif
  60. #include <boost/math/tools/config.hpp>
  61. #include <boost/math/tools/type_traits.hpp>
  62. namespace boost{ namespace math{ namespace detail{
  63. // This version will never be called (at runtime), it's a stub used
  64. // when T is unsuitable to be passed to these routines:
  65. //
  66. template <class T, class Policy>
  67. BOOST_MATH_GPU_ENABLED inline T igamma_temme_large(T, T, const Policy& /* pol */, const boost::math::integral_constant<int, 0>&)
  68. {
  69. // stub function, should never actually be called
  70. BOOST_MATH_ASSERT(0);
  71. return 0;
  72. }
  73. //
  74. // This version is accurate for up to 64-bit mantissa's,
  75. // (80-bit long double, or 10^-20).
  76. //
  77. #ifndef BOOST_MATH_HAS_GPU_SUPPORT
  78. template <class T, class Policy>
  79. BOOST_MATH_GPU_ENABLED T igamma_temme_large(T a, T x, const Policy& pol, const boost::math::integral_constant<int, 64>&)
  80. {
  81. BOOST_MATH_STD_USING // ADL of std functions
  82. T sigma = (x - a) / a;
  83. T phi = -boost::math::log1pmx(sigma, pol);
  84. T y = a * phi;
  85. T z = sqrt(2 * phi);
  86. if(x < a)
  87. z = -z;
  88. T workspace[13];
  89. // LCOV_EXCL_START
  90. BOOST_MATH_STATIC const T C0[] = {
  91. BOOST_MATH_BIG_CONSTANT(T, 64, -0.333333333333333333333),
  92. BOOST_MATH_BIG_CONSTANT(T, 64, 0.0833333333333333333333),
  93. BOOST_MATH_BIG_CONSTANT(T, 64, -0.0148148148148148148148),
  94. BOOST_MATH_BIG_CONSTANT(T, 64, 0.00115740740740740740741),
  95. BOOST_MATH_BIG_CONSTANT(T, 64, 0.000352733686067019400353),
  96. BOOST_MATH_BIG_CONSTANT(T, 64, -0.0001787551440329218107),
  97. BOOST_MATH_BIG_CONSTANT(T, 64, 0.39192631785224377817e-4),
  98. BOOST_MATH_BIG_CONSTANT(T, 64, -0.218544851067999216147e-5),
  99. BOOST_MATH_BIG_CONSTANT(T, 64, -0.18540622107151599607e-5),
  100. BOOST_MATH_BIG_CONSTANT(T, 64, 0.829671134095308600502e-6),
  101. BOOST_MATH_BIG_CONSTANT(T, 64, -0.176659527368260793044e-6),
  102. BOOST_MATH_BIG_CONSTANT(T, 64, 0.670785354340149858037e-8),
  103. BOOST_MATH_BIG_CONSTANT(T, 64, 0.102618097842403080426e-7),
  104. BOOST_MATH_BIG_CONSTANT(T, 64, -0.438203601845335318655e-8),
  105. BOOST_MATH_BIG_CONSTANT(T, 64, 0.914769958223679023418e-9),
  106. BOOST_MATH_BIG_CONSTANT(T, 64, -0.255141939949462497669e-10),
  107. BOOST_MATH_BIG_CONSTANT(T, 64, -0.583077213255042506746e-10),
  108. BOOST_MATH_BIG_CONSTANT(T, 64, 0.243619480206674162437e-10),
  109. BOOST_MATH_BIG_CONSTANT(T, 64, -0.502766928011417558909e-11),
  110. };
  111. workspace[0] = tools::evaluate_polynomial(C0, z);
  112. BOOST_MATH_STATIC const T C1[] = {
  113. BOOST_MATH_BIG_CONSTANT(T, 64, -0.00185185185185185185185),
  114. BOOST_MATH_BIG_CONSTANT(T, 64, -0.00347222222222222222222),
  115. BOOST_MATH_BIG_CONSTANT(T, 64, 0.00264550264550264550265),
  116. BOOST_MATH_BIG_CONSTANT(T, 64, -0.000990226337448559670782),
  117. BOOST_MATH_BIG_CONSTANT(T, 64, 0.000205761316872427983539),
  118. BOOST_MATH_BIG_CONSTANT(T, 64, -0.40187757201646090535e-6),
  119. BOOST_MATH_BIG_CONSTANT(T, 64, -0.18098550334489977837e-4),
  120. BOOST_MATH_BIG_CONSTANT(T, 64, 0.764916091608111008464e-5),
  121. BOOST_MATH_BIG_CONSTANT(T, 64, -0.161209008945634460038e-5),
  122. BOOST_MATH_BIG_CONSTANT(T, 64, 0.464712780280743434226e-8),
  123. BOOST_MATH_BIG_CONSTANT(T, 64, 0.137863344691572095931e-6),
  124. BOOST_MATH_BIG_CONSTANT(T, 64, -0.575254560351770496402e-7),
  125. BOOST_MATH_BIG_CONSTANT(T, 64, 0.119516285997781473243e-7),
  126. BOOST_MATH_BIG_CONSTANT(T, 64, -0.175432417197476476238e-10),
  127. BOOST_MATH_BIG_CONSTANT(T, 64, -0.100915437106004126275e-8),
  128. BOOST_MATH_BIG_CONSTANT(T, 64, 0.416279299184258263623e-9),
  129. BOOST_MATH_BIG_CONSTANT(T, 64, -0.856390702649298063807e-10),
  130. };
  131. workspace[1] = tools::evaluate_polynomial(C1, z);
  132. BOOST_MATH_STATIC const T C2[] = {
  133. BOOST_MATH_BIG_CONSTANT(T, 64, 0.00413359788359788359788),
  134. BOOST_MATH_BIG_CONSTANT(T, 64, -0.00268132716049382716049),
  135. BOOST_MATH_BIG_CONSTANT(T, 64, 0.000771604938271604938272),
  136. BOOST_MATH_BIG_CONSTANT(T, 64, 0.200938786008230452675e-5),
  137. BOOST_MATH_BIG_CONSTANT(T, 64, -0.000107366532263651605215),
  138. BOOST_MATH_BIG_CONSTANT(T, 64, 0.529234488291201254164e-4),
  139. BOOST_MATH_BIG_CONSTANT(T, 64, -0.127606351886187277134e-4),
  140. BOOST_MATH_BIG_CONSTANT(T, 64, 0.342357873409613807419e-7),
  141. BOOST_MATH_BIG_CONSTANT(T, 64, 0.137219573090629332056e-5),
  142. BOOST_MATH_BIG_CONSTANT(T, 64, -0.629899213838005502291e-6),
  143. BOOST_MATH_BIG_CONSTANT(T, 64, 0.142806142060642417916e-6),
  144. BOOST_MATH_BIG_CONSTANT(T, 64, -0.204770984219908660149e-9),
  145. BOOST_MATH_BIG_CONSTANT(T, 64, -0.140925299108675210533e-7),
  146. BOOST_MATH_BIG_CONSTANT(T, 64, 0.622897408492202203356e-8),
  147. BOOST_MATH_BIG_CONSTANT(T, 64, -0.136704883966171134993e-8),
  148. };
  149. workspace[2] = tools::evaluate_polynomial(C2, z);
  150. BOOST_MATH_STATIC const T C3[] = {
  151. BOOST_MATH_BIG_CONSTANT(T, 64, 0.000649434156378600823045),
  152. BOOST_MATH_BIG_CONSTANT(T, 64, 0.000229472093621399176955),
  153. BOOST_MATH_BIG_CONSTANT(T, 64, -0.000469189494395255712128),
  154. BOOST_MATH_BIG_CONSTANT(T, 64, 0.000267720632062838852962),
  155. BOOST_MATH_BIG_CONSTANT(T, 64, -0.756180167188397641073e-4),
  156. BOOST_MATH_BIG_CONSTANT(T, 64, -0.239650511386729665193e-6),
  157. BOOST_MATH_BIG_CONSTANT(T, 64, 0.110826541153473023615e-4),
  158. BOOST_MATH_BIG_CONSTANT(T, 64, -0.56749528269915965675e-5),
  159. BOOST_MATH_BIG_CONSTANT(T, 64, 0.142309007324358839146e-5),
  160. BOOST_MATH_BIG_CONSTANT(T, 64, -0.278610802915281422406e-10),
  161. BOOST_MATH_BIG_CONSTANT(T, 64, -0.169584040919302772899e-6),
  162. BOOST_MATH_BIG_CONSTANT(T, 64, 0.809946490538808236335e-7),
  163. BOOST_MATH_BIG_CONSTANT(T, 64, -0.191111684859736540607e-7),
  164. };
  165. workspace[3] = tools::evaluate_polynomial(C3, z);
  166. BOOST_MATH_STATIC const T C4[] = {
  167. BOOST_MATH_BIG_CONSTANT(T, 64, -0.000861888290916711698605),
  168. BOOST_MATH_BIG_CONSTANT(T, 64, 0.000784039221720066627474),
  169. BOOST_MATH_BIG_CONSTANT(T, 64, -0.000299072480303190179733),
  170. BOOST_MATH_BIG_CONSTANT(T, 64, -0.146384525788434181781e-5),
  171. BOOST_MATH_BIG_CONSTANT(T, 64, 0.664149821546512218666e-4),
  172. BOOST_MATH_BIG_CONSTANT(T, 64, -0.396836504717943466443e-4),
  173. BOOST_MATH_BIG_CONSTANT(T, 64, 0.113757269706784190981e-4),
  174. BOOST_MATH_BIG_CONSTANT(T, 64, 0.250749722623753280165e-9),
  175. BOOST_MATH_BIG_CONSTANT(T, 64, -0.169541495365583060147e-5),
  176. BOOST_MATH_BIG_CONSTANT(T, 64, 0.890750753220530968883e-6),
  177. BOOST_MATH_BIG_CONSTANT(T, 64, -0.229293483400080487057e-6),
  178. };
  179. workspace[4] = tools::evaluate_polynomial(C4, z);
  180. BOOST_MATH_STATIC const T C5[] = {
  181. BOOST_MATH_BIG_CONSTANT(T, 64, -0.000336798553366358150309),
  182. BOOST_MATH_BIG_CONSTANT(T, 64, -0.697281375836585777429e-4),
  183. BOOST_MATH_BIG_CONSTANT(T, 64, 0.000277275324495939207873),
  184. BOOST_MATH_BIG_CONSTANT(T, 64, -0.000199325705161888477003),
  185. BOOST_MATH_BIG_CONSTANT(T, 64, 0.679778047793720783882e-4),
  186. BOOST_MATH_BIG_CONSTANT(T, 64, 0.141906292064396701483e-6),
  187. BOOST_MATH_BIG_CONSTANT(T, 64, -0.135940481897686932785e-4),
  188. BOOST_MATH_BIG_CONSTANT(T, 64, 0.801847025633420153972e-5),
  189. BOOST_MATH_BIG_CONSTANT(T, 64, -0.229148117650809517038e-5),
  190. };
  191. workspace[5] = tools::evaluate_polynomial(C5, z);
  192. BOOST_MATH_STATIC const T C6[] = {
  193. BOOST_MATH_BIG_CONSTANT(T, 64, 0.000531307936463992223166),
  194. BOOST_MATH_BIG_CONSTANT(T, 64, -0.000592166437353693882865),
  195. BOOST_MATH_BIG_CONSTANT(T, 64, 0.000270878209671804482771),
  196. BOOST_MATH_BIG_CONSTANT(T, 64, 0.790235323266032787212e-6),
  197. BOOST_MATH_BIG_CONSTANT(T, 64, -0.815396936756196875093e-4),
  198. BOOST_MATH_BIG_CONSTANT(T, 64, 0.561168275310624965004e-4),
  199. BOOST_MATH_BIG_CONSTANT(T, 64, -0.183291165828433755673e-4),
  200. BOOST_MATH_BIG_CONSTANT(T, 64, -0.307961345060330478256e-8),
  201. BOOST_MATH_BIG_CONSTANT(T, 64, 0.346515536880360908674e-5),
  202. BOOST_MATH_BIG_CONSTANT(T, 64, -0.20291327396058603727e-5),
  203. BOOST_MATH_BIG_CONSTANT(T, 64, 0.57887928631490037089e-6),
  204. };
  205. workspace[6] = tools::evaluate_polynomial(C6, z);
  206. BOOST_MATH_STATIC const T C7[] = {
  207. BOOST_MATH_BIG_CONSTANT(T, 64, 0.000344367606892377671254),
  208. BOOST_MATH_BIG_CONSTANT(T, 64, 0.517179090826059219337e-4),
  209. BOOST_MATH_BIG_CONSTANT(T, 64, -0.000334931610811422363117),
  210. BOOST_MATH_BIG_CONSTANT(T, 64, 0.000281269515476323702274),
  211. BOOST_MATH_BIG_CONSTANT(T, 64, -0.000109765822446847310235),
  212. BOOST_MATH_BIG_CONSTANT(T, 64, -0.127410090954844853795e-6),
  213. BOOST_MATH_BIG_CONSTANT(T, 64, 0.277444515115636441571e-4),
  214. BOOST_MATH_BIG_CONSTANT(T, 64, -0.182634888057113326614e-4),
  215. BOOST_MATH_BIG_CONSTANT(T, 64, 0.578769494973505239894e-5),
  216. };
  217. workspace[7] = tools::evaluate_polynomial(C7, z);
  218. BOOST_MATH_STATIC const T C8[] = {
  219. BOOST_MATH_BIG_CONSTANT(T, 64, -0.000652623918595309418922),
  220. BOOST_MATH_BIG_CONSTANT(T, 64, 0.000839498720672087279993),
  221. BOOST_MATH_BIG_CONSTANT(T, 64, -0.000438297098541721005061),
  222. BOOST_MATH_BIG_CONSTANT(T, 64, -0.696909145842055197137e-6),
  223. BOOST_MATH_BIG_CONSTANT(T, 64, 0.000166448466420675478374),
  224. BOOST_MATH_BIG_CONSTANT(T, 64, -0.000127835176797692185853),
  225. BOOST_MATH_BIG_CONSTANT(T, 64, 0.462995326369130429061e-4),
  226. };
  227. workspace[8] = tools::evaluate_polynomial(C8, z);
  228. BOOST_MATH_STATIC const T C9[] = {
  229. BOOST_MATH_BIG_CONSTANT(T, 64, -0.000596761290192746250124),
  230. BOOST_MATH_BIG_CONSTANT(T, 64, -0.720489541602001055909e-4),
  231. BOOST_MATH_BIG_CONSTANT(T, 64, 0.000678230883766732836162),
  232. BOOST_MATH_BIG_CONSTANT(T, 64, -0.0006401475260262758451),
  233. BOOST_MATH_BIG_CONSTANT(T, 64, 0.000277501076343287044992),
  234. };
  235. workspace[9] = tools::evaluate_polynomial(C9, z);
  236. BOOST_MATH_STATIC const T C10[] = {
  237. BOOST_MATH_BIG_CONSTANT(T, 64, 0.00133244544948006563713),
  238. BOOST_MATH_BIG_CONSTANT(T, 64, -0.0019144384985654775265),
  239. BOOST_MATH_BIG_CONSTANT(T, 64, 0.00110893691345966373396),
  240. };
  241. workspace[10] = tools::evaluate_polynomial(C10, z);
  242. BOOST_MATH_STATIC const T C11[] = {
  243. BOOST_MATH_BIG_CONSTANT(T, 64, 0.00157972766073083495909),
  244. BOOST_MATH_BIG_CONSTANT(T, 64, 0.000162516262783915816899),
  245. BOOST_MATH_BIG_CONSTANT(T, 64, -0.00206334210355432762645),
  246. BOOST_MATH_BIG_CONSTANT(T, 64, 0.00213896861856890981541),
  247. BOOST_MATH_BIG_CONSTANT(T, 64, -0.00101085593912630031708),
  248. };
  249. workspace[11] = tools::evaluate_polynomial(C11, z);
  250. BOOST_MATH_STATIC const T C12[] = {
  251. BOOST_MATH_BIG_CONSTANT(T, 64, -0.00407251211951401664727),
  252. BOOST_MATH_BIG_CONSTANT(T, 64, 0.00640336283380806979482),
  253. BOOST_MATH_BIG_CONSTANT(T, 64, -0.00404101610816766177474),
  254. };
  255. // LCOV_EXCL_STOP
  256. workspace[12] = tools::evaluate_polynomial(C12, z);
  257. T result = tools::evaluate_polynomial<13, T, T>(workspace, 1/a);
  258. result *= exp(-y) / sqrt(2 * constants::pi<T>() * a);
  259. if(x < a)
  260. result = -result;
  261. result += boost::math::erfc(sqrt(y), pol) / 2;
  262. return result;
  263. }
  264. #endif
  265. //
  266. // This one is accurate for 53-bit mantissa's
  267. // (IEEE double precision or 10^-17).
  268. //
  269. template <class T, class Policy>
  270. BOOST_MATH_GPU_ENABLED T igamma_temme_large(T a, T x, const Policy& pol, const boost::math::integral_constant<int, 53>&)
  271. {
  272. BOOST_MATH_STD_USING // ADL of std functions
  273. T sigma = (x - a) / a;
  274. T phi = -boost::math::log1pmx(sigma, pol);
  275. T y = a * phi;
  276. T z = sqrt(2 * phi);
  277. if(x < a)
  278. z = -z;
  279. T workspace[10];
  280. // LCOV_EXCL_START
  281. BOOST_MATH_STATIC const T C0[] = {
  282. static_cast<T>(-0.33333333333333333L),
  283. static_cast<T>(0.083333333333333333L),
  284. static_cast<T>(-0.014814814814814815L),
  285. static_cast<T>(0.0011574074074074074L),
  286. static_cast<T>(0.0003527336860670194L),
  287. static_cast<T>(-0.00017875514403292181L),
  288. static_cast<T>(0.39192631785224378e-4L),
  289. static_cast<T>(-0.21854485106799922e-5L),
  290. static_cast<T>(-0.185406221071516e-5L),
  291. static_cast<T>(0.8296711340953086e-6L),
  292. static_cast<T>(-0.17665952736826079e-6L),
  293. static_cast<T>(0.67078535434014986e-8L),
  294. static_cast<T>(0.10261809784240308e-7L),
  295. static_cast<T>(-0.43820360184533532e-8L),
  296. static_cast<T>(0.91476995822367902e-9L),
  297. };
  298. workspace[0] = tools::evaluate_polynomial(C0, z);
  299. BOOST_MATH_STATIC const T C1[] = {
  300. static_cast<T>(-0.0018518518518518519L),
  301. static_cast<T>(-0.0034722222222222222L),
  302. static_cast<T>(0.0026455026455026455L),
  303. static_cast<T>(-0.00099022633744855967L),
  304. static_cast<T>(0.00020576131687242798L),
  305. static_cast<T>(-0.40187757201646091e-6L),
  306. static_cast<T>(-0.18098550334489978e-4L),
  307. static_cast<T>(0.76491609160811101e-5L),
  308. static_cast<T>(-0.16120900894563446e-5L),
  309. static_cast<T>(0.46471278028074343e-8L),
  310. static_cast<T>(0.1378633446915721e-6L),
  311. static_cast<T>(-0.5752545603517705e-7L),
  312. static_cast<T>(0.11951628599778147e-7L),
  313. };
  314. workspace[1] = tools::evaluate_polynomial(C1, z);
  315. BOOST_MATH_STATIC const T C2[] = {
  316. static_cast<T>(0.0041335978835978836L),
  317. static_cast<T>(-0.0026813271604938272L),
  318. static_cast<T>(0.00077160493827160494L),
  319. static_cast<T>(0.20093878600823045e-5L),
  320. static_cast<T>(-0.00010736653226365161L),
  321. static_cast<T>(0.52923448829120125e-4L),
  322. static_cast<T>(-0.12760635188618728e-4L),
  323. static_cast<T>(0.34235787340961381e-7L),
  324. static_cast<T>(0.13721957309062933e-5L),
  325. static_cast<T>(-0.6298992138380055e-6L),
  326. static_cast<T>(0.14280614206064242e-6L),
  327. };
  328. workspace[2] = tools::evaluate_polynomial(C2, z);
  329. BOOST_MATH_STATIC const T C3[] = {
  330. static_cast<T>(0.00064943415637860082L),
  331. static_cast<T>(0.00022947209362139918L),
  332. static_cast<T>(-0.00046918949439525571L),
  333. static_cast<T>(0.00026772063206283885L),
  334. static_cast<T>(-0.75618016718839764e-4L),
  335. static_cast<T>(-0.23965051138672967e-6L),
  336. static_cast<T>(0.11082654115347302e-4L),
  337. static_cast<T>(-0.56749528269915966e-5L),
  338. static_cast<T>(0.14230900732435884e-5L),
  339. };
  340. workspace[3] = tools::evaluate_polynomial(C3, z);
  341. BOOST_MATH_STATIC const T C4[] = {
  342. static_cast<T>(-0.0008618882909167117L),
  343. static_cast<T>(0.00078403922172006663L),
  344. static_cast<T>(-0.00029907248030319018L),
  345. static_cast<T>(-0.14638452578843418e-5L),
  346. static_cast<T>(0.66414982154651222e-4L),
  347. static_cast<T>(-0.39683650471794347e-4L),
  348. static_cast<T>(0.11375726970678419e-4L),
  349. };
  350. workspace[4] = tools::evaluate_polynomial(C4, z);
  351. BOOST_MATH_STATIC const T C5[] = {
  352. static_cast<T>(-0.00033679855336635815L),
  353. static_cast<T>(-0.69728137583658578e-4L),
  354. static_cast<T>(0.00027727532449593921L),
  355. static_cast<T>(-0.00019932570516188848L),
  356. static_cast<T>(0.67977804779372078e-4L),
  357. static_cast<T>(0.1419062920643967e-6L),
  358. static_cast<T>(-0.13594048189768693e-4L),
  359. static_cast<T>(0.80184702563342015e-5L),
  360. static_cast<T>(-0.22914811765080952e-5L),
  361. };
  362. workspace[5] = tools::evaluate_polynomial(C5, z);
  363. BOOST_MATH_STATIC const T C6[] = {
  364. static_cast<T>(0.00053130793646399222L),
  365. static_cast<T>(-0.00059216643735369388L),
  366. static_cast<T>(0.00027087820967180448L),
  367. static_cast<T>(0.79023532326603279e-6L),
  368. static_cast<T>(-0.81539693675619688e-4L),
  369. static_cast<T>(0.56116827531062497e-4L),
  370. static_cast<T>(-0.18329116582843376e-4L),
  371. };
  372. workspace[6] = tools::evaluate_polynomial(C6, z);
  373. BOOST_MATH_STATIC const T C7[] = {
  374. static_cast<T>(0.00034436760689237767L),
  375. static_cast<T>(0.51717909082605922e-4L),
  376. static_cast<T>(-0.00033493161081142236L),
  377. static_cast<T>(0.0002812695154763237L),
  378. static_cast<T>(-0.00010976582244684731L),
  379. };
  380. workspace[7] = tools::evaluate_polynomial(C7, z);
  381. BOOST_MATH_STATIC const T C8[] = {
  382. static_cast<T>(-0.00065262391859530942L),
  383. static_cast<T>(0.00083949872067208728L),
  384. static_cast<T>(-0.00043829709854172101L),
  385. };
  386. // LCOV_EXCL_STOP
  387. workspace[8] = tools::evaluate_polynomial(C8, z);
  388. workspace[9] = static_cast<T>(-0.00059676129019274625L);
  389. T result = tools::evaluate_polynomial<10, T, T>(workspace, 1/a);
  390. result *= exp(-y) / sqrt(2 * constants::pi<T>() * a);
  391. if(x < a)
  392. result = -result;
  393. #ifdef BOOST_MATH_HAS_NVRTC
  394. if (boost::math::is_same_v<T, float>)
  395. {
  396. result += ::erfcf(::sqrtf(y)) / 2;
  397. }
  398. else
  399. {
  400. result += ::erfc(::sqrt(y)) / 2;
  401. }
  402. #else
  403. result += boost::math::erfc(sqrt(y), pol) / 2;
  404. #endif
  405. return result;
  406. }
  407. //
  408. // This one is accurate for 24-bit mantissa's
  409. // (IEEE float precision, or 10^-8)
  410. //
  411. template <class T, class Policy>
  412. BOOST_MATH_GPU_ENABLED T igamma_temme_large(T a, T x, const Policy& pol, const boost::math::integral_constant<int, 24>&)
  413. {
  414. BOOST_MATH_STD_USING // ADL of std functions
  415. T sigma = (x - a) / a;
  416. T phi = -boost::math::log1pmx(sigma, pol);
  417. T y = a * phi;
  418. T z = sqrt(2 * phi);
  419. if(x < a)
  420. z = -z;
  421. T workspace[3];
  422. // LCOV_EXCL_START
  423. BOOST_MATH_STATIC const T C0[] = {
  424. static_cast<T>(-0.333333333L),
  425. static_cast<T>(0.0833333333L),
  426. static_cast<T>(-0.0148148148L),
  427. static_cast<T>(0.00115740741L),
  428. static_cast<T>(0.000352733686L),
  429. static_cast<T>(-0.000178755144L),
  430. static_cast<T>(0.391926318e-4L),
  431. };
  432. workspace[0] = tools::evaluate_polynomial(C0, z);
  433. BOOST_MATH_STATIC const T C1[] = {
  434. static_cast<T>(-0.00185185185L),
  435. static_cast<T>(-0.00347222222L),
  436. static_cast<T>(0.00264550265L),
  437. static_cast<T>(-0.000990226337L),
  438. static_cast<T>(0.000205761317L),
  439. };
  440. workspace[1] = tools::evaluate_polynomial(C1, z);
  441. BOOST_MATH_STATIC const T C2[] = {
  442. static_cast<T>(0.00413359788L),
  443. static_cast<T>(-0.00268132716L),
  444. static_cast<T>(0.000771604938L),
  445. };
  446. workspace[2] = tools::evaluate_polynomial(C2, z);
  447. // LCOV_EXCL_STOP
  448. T result = tools::evaluate_polynomial(workspace, 1/a);
  449. result *= exp(-y) / sqrt(2 * constants::pi<T>() * a);
  450. if(x < a)
  451. result = -result;
  452. #ifdef BOOST_MATH_HAS_NVRTC
  453. if (boost::math::is_same_v<T, float>)
  454. {
  455. result += ::erfcf(::sqrtf(y)) / 2;
  456. }
  457. else
  458. {
  459. result += ::erfc(::sqrt(y)) / 2;
  460. }
  461. #else
  462. result += boost::math::erfc(sqrt(y), pol) / 2;
  463. #endif
  464. return result;
  465. }
  466. //
  467. // And finally, a version for 113-bit mantissa's
  468. // (128-bit long doubles, or 10^-34).
  469. // Note this one has been optimised for a > 200
  470. // It's use for a < 200 is not recommended, that would
  471. // require many more terms in the polynomials.
  472. //
  473. // LCOV_EXCL_START: 128-bit floats not deliberately tested in our coverage tests (takes too long)
  474. #ifndef BOOST_MATH_HAS_GPU_SUPPORT
  475. template <class T, class Policy>
  476. BOOST_MATH_GPU_ENABLED T igamma_temme_large(T a, T x, const Policy& pol, const boost::math::integral_constant<int, 113>&)
  477. {
  478. BOOST_MATH_STD_USING // ADL of std functions
  479. T sigma = (x - a) / a;
  480. T phi = -boost::math::log1pmx(sigma, pol);
  481. T y = a * phi;
  482. T z = sqrt(2 * phi);
  483. if(x < a)
  484. z = -z;
  485. T workspace[14];
  486. BOOST_MATH_STATIC const T C0[] = {
  487. BOOST_MATH_BIG_CONSTANT(T, 113, -0.333333333333333333333333333333333333),
  488. BOOST_MATH_BIG_CONSTANT(T, 113, 0.0833333333333333333333333333333333333),
  489. BOOST_MATH_BIG_CONSTANT(T, 113, -0.0148148148148148148148148148148148148),
  490. BOOST_MATH_BIG_CONSTANT(T, 113, 0.00115740740740740740740740740740740741),
  491. BOOST_MATH_BIG_CONSTANT(T, 113, 0.0003527336860670194003527336860670194),
  492. BOOST_MATH_BIG_CONSTANT(T, 113, -0.000178755144032921810699588477366255144),
  493. BOOST_MATH_BIG_CONSTANT(T, 113, 0.391926317852243778169704095630021556e-4),
  494. BOOST_MATH_BIG_CONSTANT(T, 113, -0.218544851067999216147364295512443661e-5),
  495. BOOST_MATH_BIG_CONSTANT(T, 113, -0.185406221071515996070179883622956325e-5),
  496. BOOST_MATH_BIG_CONSTANT(T, 113, 0.829671134095308600501624213166443227e-6),
  497. BOOST_MATH_BIG_CONSTANT(T, 113, -0.17665952736826079304360054245742403e-6),
  498. BOOST_MATH_BIG_CONSTANT(T, 113, 0.670785354340149858036939710029613572e-8),
  499. BOOST_MATH_BIG_CONSTANT(T, 113, 0.102618097842403080425739573227252951e-7),
  500. BOOST_MATH_BIG_CONSTANT(T, 113, -0.438203601845335318655297462244719123e-8),
  501. BOOST_MATH_BIG_CONSTANT(T, 113, 0.914769958223679023418248817633113681e-9),
  502. BOOST_MATH_BIG_CONSTANT(T, 113, -0.255141939949462497668779537993887013e-10),
  503. BOOST_MATH_BIG_CONSTANT(T, 113, -0.583077213255042506746408945040035798e-10),
  504. BOOST_MATH_BIG_CONSTANT(T, 113, 0.243619480206674162436940696707789943e-10),
  505. BOOST_MATH_BIG_CONSTANT(T, 113, -0.502766928011417558909054985925744366e-11),
  506. BOOST_MATH_BIG_CONSTANT(T, 113, 0.110043920319561347708374174497293411e-12),
  507. BOOST_MATH_BIG_CONSTANT(T, 113, 0.337176326240098537882769884169200185e-12),
  508. BOOST_MATH_BIG_CONSTANT(T, 113, -0.13923887224181620659193661848957998e-12),
  509. BOOST_MATH_BIG_CONSTANT(T, 113, 0.285348938070474432039669099052828299e-13),
  510. BOOST_MATH_BIG_CONSTANT(T, 113, -0.513911183424257261899064580300494205e-15),
  511. BOOST_MATH_BIG_CONSTANT(T, 113, -0.197522882943494428353962401580710912e-14),
  512. BOOST_MATH_BIG_CONSTANT(T, 113, 0.809952115670456133407115668702575255e-15),
  513. BOOST_MATH_BIG_CONSTANT(T, 113, -0.165225312163981618191514820265351162e-15),
  514. BOOST_MATH_BIG_CONSTANT(T, 113, 0.253054300974788842327061090060267385e-17),
  515. BOOST_MATH_BIG_CONSTANT(T, 113, 0.116869397385595765888230876507793475e-16),
  516. BOOST_MATH_BIG_CONSTANT(T, 113, -0.477003704982048475822167804084816597e-17),
  517. BOOST_MATH_BIG_CONSTANT(T, 113, 0.969912605905623712420709685898585354e-18),
  518. };
  519. workspace[0] = tools::evaluate_polynomial(C0, z);
  520. BOOST_MATH_STATIC const T C1[] = {
  521. BOOST_MATH_BIG_CONSTANT(T, 113, -0.00185185185185185185185185185185185185),
  522. BOOST_MATH_BIG_CONSTANT(T, 113, -0.00347222222222222222222222222222222222),
  523. BOOST_MATH_BIG_CONSTANT(T, 113, 0.0026455026455026455026455026455026455),
  524. BOOST_MATH_BIG_CONSTANT(T, 113, -0.000990226337448559670781893004115226337),
  525. BOOST_MATH_BIG_CONSTANT(T, 113, 0.000205761316872427983539094650205761317),
  526. BOOST_MATH_BIG_CONSTANT(T, 113, -0.401877572016460905349794238683127572e-6),
  527. BOOST_MATH_BIG_CONSTANT(T, 113, -0.180985503344899778370285914867533523e-4),
  528. BOOST_MATH_BIG_CONSTANT(T, 113, 0.76491609160811100846374214980916921e-5),
  529. BOOST_MATH_BIG_CONSTANT(T, 113, -0.16120900894563446003775221882217767e-5),
  530. BOOST_MATH_BIG_CONSTANT(T, 113, 0.464712780280743434226135033938722401e-8),
  531. BOOST_MATH_BIG_CONSTANT(T, 113, 0.137863344691572095931187533077488877e-6),
  532. BOOST_MATH_BIG_CONSTANT(T, 113, -0.575254560351770496402194531835048307e-7),
  533. BOOST_MATH_BIG_CONSTANT(T, 113, 0.119516285997781473243076536699698169e-7),
  534. BOOST_MATH_BIG_CONSTANT(T, 113, -0.175432417197476476237547551202312502e-10),
  535. BOOST_MATH_BIG_CONSTANT(T, 113, -0.100915437106004126274577504686681675e-8),
  536. BOOST_MATH_BIG_CONSTANT(T, 113, 0.416279299184258263623372347219858628e-9),
  537. BOOST_MATH_BIG_CONSTANT(T, 113, -0.856390702649298063807431562579670208e-10),
  538. BOOST_MATH_BIG_CONSTANT(T, 113, 0.606721510160475861512701762169919581e-13),
  539. BOOST_MATH_BIG_CONSTANT(T, 113, 0.716249896481148539007961017165545733e-11),
  540. BOOST_MATH_BIG_CONSTANT(T, 113, -0.293318664377143711740636683615595403e-11),
  541. BOOST_MATH_BIG_CONSTANT(T, 113, 0.599669636568368872330374527568788909e-12),
  542. BOOST_MATH_BIG_CONSTANT(T, 113, -0.216717865273233141017100472779701734e-15),
  543. BOOST_MATH_BIG_CONSTANT(T, 113, -0.497833997236926164052815522048108548e-13),
  544. BOOST_MATH_BIG_CONSTANT(T, 113, 0.202916288237134247736694804325894226e-13),
  545. BOOST_MATH_BIG_CONSTANT(T, 113, -0.413125571381061004935108332558187111e-14),
  546. BOOST_MATH_BIG_CONSTANT(T, 113, 0.828651623988309644380188591057589316e-18),
  547. BOOST_MATH_BIG_CONSTANT(T, 113, 0.341003088693333279336339355910600992e-15),
  548. BOOST_MATH_BIG_CONSTANT(T, 113, -0.138541953028939715357034547426313703e-15),
  549. BOOST_MATH_BIG_CONSTANT(T, 113, 0.281234665322887466568860332727259483e-16),
  550. };
  551. workspace[1] = tools::evaluate_polynomial(C1, z);
  552. BOOST_MATH_STATIC const T C2[] = {
  553. BOOST_MATH_BIG_CONSTANT(T, 113, 0.0041335978835978835978835978835978836),
  554. BOOST_MATH_BIG_CONSTANT(T, 113, -0.00268132716049382716049382716049382716),
  555. BOOST_MATH_BIG_CONSTANT(T, 113, 0.000771604938271604938271604938271604938),
  556. BOOST_MATH_BIG_CONSTANT(T, 113, 0.200938786008230452674897119341563786e-5),
  557. BOOST_MATH_BIG_CONSTANT(T, 113, -0.000107366532263651605215391223621676297),
  558. BOOST_MATH_BIG_CONSTANT(T, 113, 0.529234488291201254164217127180090143e-4),
  559. BOOST_MATH_BIG_CONSTANT(T, 113, -0.127606351886187277133779191392360117e-4),
  560. BOOST_MATH_BIG_CONSTANT(T, 113, 0.34235787340961380741902003904747389e-7),
  561. BOOST_MATH_BIG_CONSTANT(T, 113, 0.137219573090629332055943852926020279e-5),
  562. BOOST_MATH_BIG_CONSTANT(T, 113, -0.629899213838005502290672234278391876e-6),
  563. BOOST_MATH_BIG_CONSTANT(T, 113, 0.142806142060642417915846008822771748e-6),
  564. BOOST_MATH_BIG_CONSTANT(T, 113, -0.204770984219908660149195854409200226e-9),
  565. BOOST_MATH_BIG_CONSTANT(T, 113, -0.140925299108675210532930244154315272e-7),
  566. BOOST_MATH_BIG_CONSTANT(T, 113, 0.622897408492202203356394293530327112e-8),
  567. BOOST_MATH_BIG_CONSTANT(T, 113, -0.136704883966171134992724380284402402e-8),
  568. BOOST_MATH_BIG_CONSTANT(T, 113, 0.942835615901467819547711211663208075e-12),
  569. BOOST_MATH_BIG_CONSTANT(T, 113, 0.128722524000893180595479368872770442e-9),
  570. BOOST_MATH_BIG_CONSTANT(T, 113, -0.556459561343633211465414765894951439e-10),
  571. BOOST_MATH_BIG_CONSTANT(T, 113, 0.119759355463669810035898150310311343e-10),
  572. BOOST_MATH_BIG_CONSTANT(T, 113, -0.416897822518386350403836626692480096e-14),
  573. BOOST_MATH_BIG_CONSTANT(T, 113, -0.109406404278845944099299008640802908e-11),
  574. BOOST_MATH_BIG_CONSTANT(T, 113, 0.4662239946390135746326204922464679e-12),
  575. BOOST_MATH_BIG_CONSTANT(T, 113, -0.990510576390690597844122258212382301e-13),
  576. BOOST_MATH_BIG_CONSTANT(T, 113, 0.189318767683735145056885183170630169e-16),
  577. BOOST_MATH_BIG_CONSTANT(T, 113, 0.885922187259112726176031067028740667e-14),
  578. BOOST_MATH_BIG_CONSTANT(T, 113, -0.373782039804640545306560251777191937e-14),
  579. BOOST_MATH_BIG_CONSTANT(T, 113, 0.786883363903515525774088394065960751e-15),
  580. };
  581. workspace[2] = tools::evaluate_polynomial(C2, z);
  582. BOOST_MATH_STATIC const T C3[] = {
  583. BOOST_MATH_BIG_CONSTANT(T, 113, 0.000649434156378600823045267489711934156),
  584. BOOST_MATH_BIG_CONSTANT(T, 113, 0.000229472093621399176954732510288065844),
  585. BOOST_MATH_BIG_CONSTANT(T, 113, -0.000469189494395255712128140111679206329),
  586. BOOST_MATH_BIG_CONSTANT(T, 113, 0.000267720632062838852962309752433209223),
  587. BOOST_MATH_BIG_CONSTANT(T, 113, -0.756180167188397641072538191879755666e-4),
  588. BOOST_MATH_BIG_CONSTANT(T, 113, -0.239650511386729665193314027333231723e-6),
  589. BOOST_MATH_BIG_CONSTANT(T, 113, 0.110826541153473023614770299726861227e-4),
  590. BOOST_MATH_BIG_CONSTANT(T, 113, -0.567495282699159656749963105701560205e-5),
  591. BOOST_MATH_BIG_CONSTANT(T, 113, 0.14230900732435883914551894470580433e-5),
  592. BOOST_MATH_BIG_CONSTANT(T, 113, -0.278610802915281422405802158211174452e-10),
  593. BOOST_MATH_BIG_CONSTANT(T, 113, -0.16958404091930277289864168795820267e-6),
  594. BOOST_MATH_BIG_CONSTANT(T, 113, 0.809946490538808236335278504852724081e-7),
  595. BOOST_MATH_BIG_CONSTANT(T, 113, -0.191111684859736540606728140872727635e-7),
  596. BOOST_MATH_BIG_CONSTANT(T, 113, 0.239286204398081179686413514022282056e-11),
  597. BOOST_MATH_BIG_CONSTANT(T, 113, 0.206201318154887984369925818486654549e-8),
  598. BOOST_MATH_BIG_CONSTANT(T, 113, -0.946049666185513217375417988510192814e-9),
  599. BOOST_MATH_BIG_CONSTANT(T, 113, 0.215410497757749078380130268468744512e-9),
  600. BOOST_MATH_BIG_CONSTANT(T, 113, -0.138882333681390304603424682490735291e-13),
  601. BOOST_MATH_BIG_CONSTANT(T, 113, -0.218947616819639394064123400466489455e-10),
  602. BOOST_MATH_BIG_CONSTANT(T, 113, 0.979099895117168512568262802255883368e-11),
  603. BOOST_MATH_BIG_CONSTANT(T, 113, -0.217821918801809621153859472011393244e-11),
  604. BOOST_MATH_BIG_CONSTANT(T, 113, 0.62088195734079014258166361684972205e-16),
  605. BOOST_MATH_BIG_CONSTANT(T, 113, 0.212697836327973697696702537114614471e-12),
  606. BOOST_MATH_BIG_CONSTANT(T, 113, -0.934468879151743333127396765626749473e-13),
  607. BOOST_MATH_BIG_CONSTANT(T, 113, 0.204536712267828493249215913063207436e-13),
  608. };
  609. workspace[3] = tools::evaluate_polynomial(C3, z);
  610. BOOST_MATH_STATIC const T C4[] = {
  611. BOOST_MATH_BIG_CONSTANT(T, 113, -0.000861888290916711698604702719929057378),
  612. BOOST_MATH_BIG_CONSTANT(T, 113, 0.00078403922172006662747403488144228885),
  613. BOOST_MATH_BIG_CONSTANT(T, 113, -0.000299072480303190179733389609932819809),
  614. BOOST_MATH_BIG_CONSTANT(T, 113, -0.146384525788434181781232535690697556e-5),
  615. BOOST_MATH_BIG_CONSTANT(T, 113, 0.664149821546512218665853782451862013e-4),
  616. BOOST_MATH_BIG_CONSTANT(T, 113, -0.396836504717943466443123507595386882e-4),
  617. BOOST_MATH_BIG_CONSTANT(T, 113, 0.113757269706784190980552042885831759e-4),
  618. BOOST_MATH_BIG_CONSTANT(T, 113, 0.250749722623753280165221942390057007e-9),
  619. BOOST_MATH_BIG_CONSTANT(T, 113, -0.169541495365583060147164356781525752e-5),
  620. BOOST_MATH_BIG_CONSTANT(T, 113, 0.890750753220530968882898422505515924e-6),
  621. BOOST_MATH_BIG_CONSTANT(T, 113, -0.229293483400080487057216364891158518e-6),
  622. BOOST_MATH_BIG_CONSTANT(T, 113, 0.295679413754404904696572852500004588e-10),
  623. BOOST_MATH_BIG_CONSTANT(T, 113, 0.288658297427087836297341274604184504e-7),
  624. BOOST_MATH_BIG_CONSTANT(T, 113, -0.141897394378032193894774303903982717e-7),
  625. BOOST_MATH_BIG_CONSTANT(T, 113, 0.344635804994648970659527720474194356e-8),
  626. BOOST_MATH_BIG_CONSTANT(T, 113, -0.230245171745280671320192735850147087e-12),
  627. BOOST_MATH_BIG_CONSTANT(T, 113, -0.394092330280464052750697640085291799e-9),
  628. BOOST_MATH_BIG_CONSTANT(T, 113, 0.186023389685045019134258533045185639e-9),
  629. BOOST_MATH_BIG_CONSTANT(T, 113, -0.435632300505661804380678327446262424e-10),
  630. BOOST_MATH_BIG_CONSTANT(T, 113, 0.127860010162962312660550463349930726e-14),
  631. BOOST_MATH_BIG_CONSTANT(T, 113, 0.467927502665791946200382739991760062e-11),
  632. BOOST_MATH_BIG_CONSTANT(T, 113, -0.214924647061348285410535341910721086e-11),
  633. BOOST_MATH_BIG_CONSTANT(T, 113, 0.490881561480965216323649688463984082e-12),
  634. };
  635. workspace[4] = tools::evaluate_polynomial(C4, z);
  636. BOOST_MATH_STATIC const T C5[] = {
  637. BOOST_MATH_BIG_CONSTANT(T, 113, -0.000336798553366358150308767592718210002),
  638. BOOST_MATH_BIG_CONSTANT(T, 113, -0.697281375836585777429398828575783308e-4),
  639. BOOST_MATH_BIG_CONSTANT(T, 113, 0.00027727532449593920787336425196507501),
  640. BOOST_MATH_BIG_CONSTANT(T, 113, -0.000199325705161888477003360405280844238),
  641. BOOST_MATH_BIG_CONSTANT(T, 113, 0.679778047793720783881640176604435742e-4),
  642. BOOST_MATH_BIG_CONSTANT(T, 113, 0.141906292064396701483392727105575757e-6),
  643. BOOST_MATH_BIG_CONSTANT(T, 113, -0.135940481897686932784583938837504469e-4),
  644. BOOST_MATH_BIG_CONSTANT(T, 113, 0.80184702563342015397192571980419684e-5),
  645. BOOST_MATH_BIG_CONSTANT(T, 113, -0.229148117650809517038048790128781806e-5),
  646. BOOST_MATH_BIG_CONSTANT(T, 113, -0.325247355129845395166230137750005047e-9),
  647. BOOST_MATH_BIG_CONSTANT(T, 113, 0.346528464910852649559195496827579815e-6),
  648. BOOST_MATH_BIG_CONSTANT(T, 113, -0.184471871911713432765322367374920978e-6),
  649. BOOST_MATH_BIG_CONSTANT(T, 113, 0.482409670378941807563762631738989002e-7),
  650. BOOST_MATH_BIG_CONSTANT(T, 113, -0.179894667217435153025754291716644314e-13),
  651. BOOST_MATH_BIG_CONSTANT(T, 113, -0.630619450001352343517516981425944698e-8),
  652. BOOST_MATH_BIG_CONSTANT(T, 113, 0.316241762877456793773762181540969623e-8),
  653. BOOST_MATH_BIG_CONSTANT(T, 113, -0.784092425369742929000839303523267545e-9),
  654. };
  655. workspace[5] = tools::evaluate_polynomial(C5, z);
  656. BOOST_MATH_STATIC const T C6[] = {
  657. BOOST_MATH_BIG_CONSTANT(T, 113, 0.00053130793646399222316574854297762391),
  658. BOOST_MATH_BIG_CONSTANT(T, 113, -0.000592166437353693882864836225604401187),
  659. BOOST_MATH_BIG_CONSTANT(T, 113, 0.000270878209671804482771279183488328692),
  660. BOOST_MATH_BIG_CONSTANT(T, 113, 0.790235323266032787212032944390816666e-6),
  661. BOOST_MATH_BIG_CONSTANT(T, 113, -0.815396936756196875092890088464682624e-4),
  662. BOOST_MATH_BIG_CONSTANT(T, 113, 0.561168275310624965003775619041471695e-4),
  663. BOOST_MATH_BIG_CONSTANT(T, 113, -0.183291165828433755673259749374098313e-4),
  664. BOOST_MATH_BIG_CONSTANT(T, 113, -0.307961345060330478256414192546677006e-8),
  665. BOOST_MATH_BIG_CONSTANT(T, 113, 0.346515536880360908673728529745376913e-5),
  666. BOOST_MATH_BIG_CONSTANT(T, 113, -0.202913273960586037269527254582695285e-5),
  667. BOOST_MATH_BIG_CONSTANT(T, 113, 0.578879286314900370889997586203187687e-6),
  668. BOOST_MATH_BIG_CONSTANT(T, 113, 0.233863067382665698933480579231637609e-12),
  669. BOOST_MATH_BIG_CONSTANT(T, 113, -0.88286007463304835250508524317926246e-7),
  670. BOOST_MATH_BIG_CONSTANT(T, 113, 0.474359588804081278032150770595852426e-7),
  671. BOOST_MATH_BIG_CONSTANT(T, 113, -0.125454150207103824457130611214783073e-7),
  672. };
  673. workspace[6] = tools::evaluate_polynomial(C6, z);
  674. BOOST_MATH_STATIC const T C7[] = {
  675. BOOST_MATH_BIG_CONSTANT(T, 113, 0.000344367606892377671254279625108523655),
  676. BOOST_MATH_BIG_CONSTANT(T, 113, 0.517179090826059219337057843002058823e-4),
  677. BOOST_MATH_BIG_CONSTANT(T, 113, -0.000334931610811422363116635090580012327),
  678. BOOST_MATH_BIG_CONSTANT(T, 113, 0.000281269515476323702273722110707777978),
  679. BOOST_MATH_BIG_CONSTANT(T, 113, -0.000109765822446847310235396824500789005),
  680. BOOST_MATH_BIG_CONSTANT(T, 113, -0.127410090954844853794579954588107623e-6),
  681. BOOST_MATH_BIG_CONSTANT(T, 113, 0.277444515115636441570715073933712622e-4),
  682. BOOST_MATH_BIG_CONSTANT(T, 113, -0.182634888057113326614324442681892723e-4),
  683. BOOST_MATH_BIG_CONSTANT(T, 113, 0.578769494973505239894178121070843383e-5),
  684. BOOST_MATH_BIG_CONSTANT(T, 113, 0.493875893393627039981813418398565502e-9),
  685. BOOST_MATH_BIG_CONSTANT(T, 113, -0.105953670140260427338098566209633945e-5),
  686. BOOST_MATH_BIG_CONSTANT(T, 113, 0.616671437611040747858836254004890765e-6),
  687. BOOST_MATH_BIG_CONSTANT(T, 113, -0.175629733590604619378669693914265388e-6),
  688. };
  689. workspace[7] = tools::evaluate_polynomial(C7, z);
  690. BOOST_MATH_STATIC const T C8[] = {
  691. BOOST_MATH_BIG_CONSTANT(T, 113, -0.000652623918595309418922034919726622692),
  692. BOOST_MATH_BIG_CONSTANT(T, 113, 0.000839498720672087279993357516764983445),
  693. BOOST_MATH_BIG_CONSTANT(T, 113, -0.000438297098541721005061087953050560377),
  694. BOOST_MATH_BIG_CONSTANT(T, 113, -0.696909145842055197136911097362072702e-6),
  695. BOOST_MATH_BIG_CONSTANT(T, 113, 0.00016644846642067547837384572662326101),
  696. BOOST_MATH_BIG_CONSTANT(T, 113, -0.000127835176797692185853344001461664247),
  697. BOOST_MATH_BIG_CONSTANT(T, 113, 0.462995326369130429061361032704489636e-4),
  698. BOOST_MATH_BIG_CONSTANT(T, 113, 0.455790986792270771162749294232219616e-8),
  699. BOOST_MATH_BIG_CONSTANT(T, 113, -0.105952711258051954718238500312872328e-4),
  700. BOOST_MATH_BIG_CONSTANT(T, 113, 0.678334290486516662273073740749269432e-5),
  701. BOOST_MATH_BIG_CONSTANT(T, 113, -0.210754766662588042469972680229376445e-5),
  702. };
  703. workspace[8] = tools::evaluate_polynomial(C8, z);
  704. BOOST_MATH_STATIC const T C9[] = {
  705. BOOST_MATH_BIG_CONSTANT(T, 113, -0.000596761290192746250124390067179459605),
  706. BOOST_MATH_BIG_CONSTANT(T, 113, -0.720489541602001055908571930225015052e-4),
  707. BOOST_MATH_BIG_CONSTANT(T, 113, 0.000678230883766732836161951166000673426),
  708. BOOST_MATH_BIG_CONSTANT(T, 113, -0.000640147526026275845100045652582354779),
  709. BOOST_MATH_BIG_CONSTANT(T, 113, 0.000277501076343287044992374518205845463),
  710. BOOST_MATH_BIG_CONSTANT(T, 113, 0.181970083804651510461686554030325202e-6),
  711. BOOST_MATH_BIG_CONSTANT(T, 113, -0.847950711706850318239732559632810086e-4),
  712. BOOST_MATH_BIG_CONSTANT(T, 113, 0.610519208250153101764709122740859458e-4),
  713. BOOST_MATH_BIG_CONSTANT(T, 113, -0.210739201834048624082975255893773306e-4),
  714. };
  715. workspace[9] = tools::evaluate_polynomial(C9, z);
  716. BOOST_MATH_STATIC const T C10[] = {
  717. BOOST_MATH_BIG_CONSTANT(T, 113, 0.00133244544948006563712694993432717968),
  718. BOOST_MATH_BIG_CONSTANT(T, 113, -0.00191443849856547752650089885832852254),
  719. BOOST_MATH_BIG_CONSTANT(T, 113, 0.0011089369134596637339607446329267522),
  720. BOOST_MATH_BIG_CONSTANT(T, 113, 0.993240412264229896742295262075817566e-6),
  721. BOOST_MATH_BIG_CONSTANT(T, 113, -0.000508745012930931989848393025305956774),
  722. BOOST_MATH_BIG_CONSTANT(T, 113, 0.00042735056665392884328432271160040444),
  723. BOOST_MATH_BIG_CONSTANT(T, 113, -0.000168588537679107988033552814662382059),
  724. };
  725. workspace[10] = tools::evaluate_polynomial(C10, z);
  726. BOOST_MATH_STATIC const T C11[] = {
  727. BOOST_MATH_BIG_CONSTANT(T, 113, 0.00157972766073083495908785631307733022),
  728. BOOST_MATH_BIG_CONSTANT(T, 113, 0.000162516262783915816898635123980270998),
  729. BOOST_MATH_BIG_CONSTANT(T, 113, -0.00206334210355432762645284467690276817),
  730. BOOST_MATH_BIG_CONSTANT(T, 113, 0.00213896861856890981541061922797693947),
  731. BOOST_MATH_BIG_CONSTANT(T, 113, -0.00101085593912630031708085801712479376),
  732. };
  733. workspace[11] = tools::evaluate_polynomial(C11, z);
  734. BOOST_MATH_STATIC const T C12[] = {
  735. BOOST_MATH_BIG_CONSTANT(T, 113, -0.00407251211951401664727281097914544601),
  736. BOOST_MATH_BIG_CONSTANT(T, 113, 0.00640336283380806979482363809026579583),
  737. BOOST_MATH_BIG_CONSTANT(T, 113, -0.00404101610816766177473974858518094879),
  738. };
  739. workspace[12] = tools::evaluate_polynomial(C12, z);
  740. workspace[13] = -0.0059475779383993002845382844736066323L;
  741. T result = tools::evaluate_polynomial(workspace, T(1/a));
  742. result *= exp(-y) / sqrt(2 * constants::pi<T>() * a);
  743. if(x < a)
  744. result = -result;
  745. result += boost::math::erfc(sqrt(y), pol) / 2;
  746. return result;
  747. }
  748. // LCOV_EXCL_STOP
  749. #endif
  750. } // namespace detail
  751. } // namespace math
  752. } // namespace math
  753. #endif // BOOST_MATH_DETAIL_IGAMMA_LARGE