1[section:hypergeometric Hypergeometric Functions] 2 3[section:hypergeometric_1f0 Hypergeometric [sub 1]/F/[sub 0] ] 4 5 #include <boost/math/special_functions/hypergeometric_1F0.hpp> 6 7 namespace boost { namespace math { 8 9 template <class T1, class T2> 10 ``__sf_result`` hypergeometric_1F0(T1 a, T2 z); 11 12 template <class T1, class T2, class ``__Policy``> 13 ``__sf_result`` hypergeometric_1F0(T1 a, T2 z, const ``__Policy``&); 14 15 }} // namespaces 16 17[h4 Description] 18 19The function `hypergeometric_1F0` returns the result of: 20 21[equation hyper_1f0] 22 23The return type of these functions is computed using the __arg_promotion_rules 24when `T1` and `T2` are different types. 25 26[optional_policy] 27 28The functions return the result of __domain_error whenever the result is 29undefined or complex. This occurs for `z == 1` or `1 - z < 0` and `a` not an integer. 30 31[h4 Implementation] 32 33The implementation is trivial: 34 35[expression ['[sub 1]F[sub 0](a, z) = (1-z)[super -a]]] 36 37[endsect] [/section:hyper_geometric_1f0 Hypergeometric [sub 1]/F/[sub 0] ] 38 39[section:hypergeometric_0f1 Hypergeometric [sub 0]/F/[sub 1] ] 40 41 #include <boost/math/special_functions/hypergeometric_0F1.hpp> 42 43 namespace boost { namespace math { 44 45 template <class T1, class T2> 46 ``__sf_result`` hypergeometric_0F1(T1 b, T2 z); 47 48 template <class T1, class T2, class ``__Policy``> 49 ``__sf_result`` hypergeometric_0F1(T1 b, T2 z, const ``__Policy``&); 50 51 }} // namespaces 52 53[h4 Description] 54 55The function `hypergeometric_0F1` returns the result of 56 57[equation hyper_0f1] 58 59The return type of these functions is computed using the __arg_promotion_rules 60when `T1` and `T2` are different types. 61 62[optional_policy] 63 64The functions return the result of __domain_error whenever the result is 65undefined or complex. 66This occurs only when `b` is an integer <= 0. 67 68[h4 Implementation] 69 70The function is implemented via its defining series wherever convergent. 71 72For a divergent series we use the continued fraction as long as the result is not too small: 73 74[equation hyper_0f1_cf] 75 76and one of the Bessel function relations otherwise: 77 78[equation hyper_0f1_bessel_j] 79 80[equation hyper_0f1_bessel_i] 81 82[endsect] [/section:hypergeometric_0f1 Hypergeometric [sub 0]/F/[sub 1] ] 83 84[section:hypergeometric_2f0 Hypergeometric [sub 2]/F/[sub 0]] 85 86 #include <boost/math/special_functions/hypergeometric_2F0.hpp> 87 88 namespace boost { namespace math { 89 90 template <class T1, class T2, class T3> 91 ``__sf_result`` hypergeometric_2F0(T1 a1, T2 a2, T3 z); 92 93 template <class T1, class T2, class T3, class ``__Policy``> 94 ``__sf_result`` hypergeometric_2F0(T1 a1, T2 a2, T3 z, const ``__Policy``&); 95 96 }} 97 98[h4 Description] 99 100The function `hypergeometric_2F0` returns the result of 101 102[equation hyper_2f0] 103 104The return type of these functions is computed using the __arg_promotion_rules 105when `T1` and `T2` are different types. 106 107[optional_policy] 108 109The functions return the result of __domain_error whenever the result is 110undefined or complex. The valid domain for this function occurs only when one of `a1` or 111`a2` is a negative integer: ie the polynomial case. 112 113[h4 Implementation] 114 115When `a1 == a2 - 0.5` then the function is implemented in terms of the Hermite polynomial: 116 117[equation hyper_2f0_hermite] 118 119When both `a1` and `a2` are integers then the function is implemented in terms of the associated-Laguerre polynomial: 120 121[equation hyper_2f0_laguerre] 122 123If the defining series is divergent, we use the continued fraction 124 125[equation hyper_2f0_cf] 126 127Otherwise we use the defining series. 128 129[endsect] [/section:hyper_geometric_1f0 Hypergeometric [sub 2]/F/[sub 0]] 130 131[section:hypergeometric_1f1 Hypergeometric [sub 1]/F/[sub 1]] 132 133 134 #include <boost/math/special_functions/hypergeometric_1F1.hpp> 135 136 namespace boost { namespace math { 137 138 template <class T1, class T2, class T3> 139 ``__sf_result`` hypergeometric_1F1(T1 a, T2 b, T3 z); 140 141 template <class T1, class T2, class T3, class ``__Policy``> 142 ``__sf_result`` hypergeometric_1F1(T1 a, T2 b, T3 z, const ``__Policy``&); 143 144 template <class T1, class T2, class T3> 145 ``__sf_result`` hypergeometric_1F1_regularized(T1 a, T2 b, T3 z); 146 147 template <class T1, class T2, class T3, class ``__Policy``> 148 ``__sf_result`` hypergeometric_1F1_regularized(T1 a, T2 b, T3 z, const ``__Policy``&); 149 150 template <class T1, class T2, class T3> 151 ``__sf_result`` log_hypergeometric_1F1(T1 a, T2 b, T3 z); 152 153 template <class T1, class T2, class T3, class ``__Policy``> 154 ``__sf_result`` log_hypergeometric_1F1(T1 a, T2 b, T3 z, const ``__Policy``&); 155 156 template <class T1, class T2, class T3> 157 ``__sf_result`` log_hypergeometric_1F1(T1 a, T2 b, T3 z, int* sign); 158 159 template <class T1, class T2, class T3, class ``__Policy``> 160 ``__sf_result`` log_hypergeometric_1F1(T1 a, T2 b, T3 z, int* sign, const ``__Policy``&); 161 162 }} // namespaces 163 164[h4 Description] 165 166[? __build_html '''<?dbhtml-include href="'''__base_path__'''/graphs/hypergeometric_1f1/script_include.html"?>'''] 167 168The function `hypergeometric_1F1(a, b, z)` returns the non-singular solution to 169[@https://en.wikipedia.org/wiki/Confluent_hypergeometric_function Kummer's equation] 170 171[:[/\large $$z \frac{d^2 w}{d z^2} + (b-z) \frac{dw}{dz} - aw = 0 $$] 172[$../equations/hypergeometric_1f1_2.svg]] 173 174which for |/z/| < 1 has the hypergeometric series expansion 175 176[:[$../equations/hypergeometric_1f1_1.svg]] 177 178where (/a/)[sub /n/] denotes rising factorial. 179This function has the same definition as Mathematica's `Hypergeometric1F1[a, b, z]` and Maple's `KummerM(a, b, z)`. 180 181The "regularized" versions of the function return: 182 183[:[/ \Large $$ \textbf{M}(a, b; z) = \frac{{_1F_1}(a, b; z)}{\Gamma(b)} = \sum_{n=0}^{\infty} \frac{(a)_n z^n}{\Gamma(b+n) n!} $$] 184[$../equations/hypergeometric_1f1_17.svg]] 185 186The "log" versions of the function return: 187 188[:[/ \Large $$ \ln(|_1F_1(a, b, z)|) $$ ] 189[$../equations/hypergeometric_1f1_18.svg]] 190 191When the optional `sign` argument is supplied it is set on exit to either +1 or -1 depending on the sign of [sub 1]/F/[sub 1](/a/, /b/, /z/). 192 193Both the regularized and the logarithmic versions of these functions return results without the spurious 194under/overflow that plague naive implementations. 195 196[h4 Known Issues] 197 198This function is still very much the subject of active research, 199and a full range of methods capable of calculating the function 200over its entire domain are not yet available. 201We have worked extremely hard to ensure that problem domains result in an exception being thrown 202(an __evaluation_error) rather than return a garbage result. 203Domains that are still unsolved include: 204 205[table 206[[domain][comment][example]] 207[[ [/\large $$_1F_1(-a, -b; -z)$$] [$../equations/hypergeometric_1f1_13.svg]][ /a, b/ and /z/ all large.][[sub 1]F[sub 1](-814723.75, -13586.87890625, -15.87335205078125)]] 208[[ [/\large $$_1F_1(-a, -b; z)$$] [$../equations/hypergeometric_1f1_14.svg]][ /a < b/, /b > z/, this is really the same domain as above.][ ]] 209[[ [/\large $$_1F_1(a, -b; z)$$] [$../equations/hypergeometric_1f1_15.svg]][ There is a gap in between methods where no reliable implementation is possible, the issue becomes much worse for /a/, /b/ and /z/ all large.][[sub 1]F[sub 1](9057.91796875, -1252.51318359375, 15.87335205078125)]] 210[[ [/\large $$_1F_1(a_{tiny}, b; -z)$$] [$../equations/hypergeometric_1f1_16.svg] ][This domain is mostly handled by A&S 13.3.6 (see implementation notes below), but there 211 are some values where either that series is non-convergent (most particularly for /b/ < 0) 212 or where the series becomes divergent after a few terms limiting the precision that can be achieved.][[sub 1]F[sub 1](-5.9981750131794866e-15, 0.499999999999994, -240.42092034220695)]] 213] 214 215The following graph illustrates the problem area for /b < 0/ and /a,z > 0/: 216 217[? __build_html '''<?dbhtml-include href="'''__base_path__'''/graphs/hypergeometric_1f1/negative_b_incalculable.html"?>'''] 218[?! __build_html [$../graphs/hypergeometric_1f1/negative_b_incalculable.png]] 219 220[h4 Testing] 221 222There are 3 main groups of tests: 223 224* Spot tests for special inputs with known values. 225* Sanity checks which use integrals of hypergeometric functions of known value. These are particularly useful 226since for fixed ['a] and ['b] they evaluate ['[sub 1]F[sub 1](a,b,z)] for all /z/. 227* Testing against values precomputed using arbitrary precision arithmetic and the /pFq/ series. 228 229We also have a [@../../tools/hypergeometric_1F1_error_plot.cpp small program] for testing random values over a user-specified domain of /a/, /b/ and /z/, this program 230is also used for the error rate plots below and has been extremely useful in fine-tuning the implementation. 231 232It should be noted however, that there are some domains for large negative /a/ and /b/ that have poor test coverage as we were 233simply unable to compute test values in reasonable time: it is not uncommon for the /pFq/ series to cancel many hundreds of digits 234and sometimes into the thousands of digits. 235 236[h4 Errors] 237 238We divide the domain of [sub 1]/F/[sub 1] up into several domains to explain the error rates. 239 240In the following graphs we ran 100000 random test cases over each domain, note that the scatter plots show only the very worst errors 241as otherwise the graphs are both incomprehensible and virtually unplottable (as in sudden browser death). 242 243For 0 < a,b,z the error rates are trivially small unless we are forced to take steps to avoid very-slow convergence and use a method other than direct evaluation of the series 244for performance reasons. Even so the errors rates are broadly acceptable, while the scatter graph shows where the worst errors are located: 245 246[$../graphs/hypergeometric_1f1/positive_abz_bins.svg] 247[? __build_html '''<?dbhtml-include href="'''__base_path__'''/graphs/hypergeometric_1f1/positive_abz.html"?>'''] 248[?! __build_html [$../graphs/hypergeometric_1f1/positive_abz.png]] 249 250For -1000 < a < 0 and 0 < b,z < 1000 the maximum error rates are higher, but most are still low, and the worst errors are fairly evenly distributed: 251 252[$../graphs/hypergeometric_1f1/negative_a_bins.svg] 253[? __build_html '''<?dbhtml-include href="'''__base_path__'''/graphs/hypergeometric_1f1/negative_a.html"?>'''] 254[?! __build_html [$../graphs/hypergeometric_1f1/negative_a.png]] 255 256For -1000 < /b/ < 0 and /a/,/z/ > 0 the errors are mostly low at double precision except near the "unimplementable zone". 257Note however, that the some of the methods used here fail to converge to much better than double precision. 258 259[$../graphs/hypergeometric_1f1/negative_b_bins.svg] 260[? __build_html '''<?dbhtml-include href="'''__base_path__'''/graphs/hypergeometric_1f1/negative_b.html"?>'''] 261[?! __build_html [$../graphs/hypergeometric_1f1/negative_b.png]] 262 263For both /a/ and /b/ less than zero, the errors the worst errors are clustered in a "difficult zone" with /b < a/ and /z/ [asymp] /a/: 264 265[$../graphs/hypergeometric_1f1/negative_ab_bins.svg] 266[? __build_html '''<?dbhtml-include href="'''__base_path__'''/graphs/hypergeometric_1f1/negative_ab.html"?>'''] 267[?! __build_html [$../graphs/hypergeometric_1f1/negative_ab.png]] 268 269 270[h4 Implementation] 271 272The implementation for this function is remarkably complex; 273readers will have to refer to the code for the transitions between methods, as we can only give an overview here. 274 275In almost all cases where /z < 0/ we use [@https://en.wikipedia.org/wiki/Confluent_hypergeometric_function Kummer's relation] 276to make /z/ positive: 277 278[:[/\large $$_1F_1(a, b; -z) = e^{-z} {_1F_1}(a, b; z)$$] 279[$../equations/hypergeometric_1f1_12.svg]] 280 281The main series representation 282 283[:[$../equations/hypergeometric_1f1_1.svg]] 284 285is used only when 286 287* we are sure that it is convergent and not subject to excessive cancellation, and 288* there is no other better method available. 289 290The implementation of this series is complicated by the fact that for /b/ < 0 then what appears to be a fully converged series can in fact suddenly jump back 291to life again as /b/ crosses the origin. 292 293A&S 13.3.6 gives 294 295[:[/\large $$_1F_1(a, b; z) = e^{ \frac{z}{2} } \Gamma(b-a- \frac{1}{2} ) ( \frac{z}{4})^{a-b+ \frac{1}{2}} \sum_{n=0}^{\infty} \frac{(2b-2a-1)_n(b-2a)_n(b-a-\frac{1}{2}+n)}{n!(b)_n} (-1)^n I_{b-a-\frac{1}{2}+n}(\frac{z}{2})$$] 296[$../equations/hypergeometric_1f1_3.svg]] 297 298which is particularly useful for ['a [cong] b] and ['z > 0], or /a/ \u226a 1, and ['z < 0]. 299 300Then we have Tricomi's approximation (given in simplified form in A&S 13.3.7) [link math_toolkit.hypergeometric.hypergeometric_refs (7)] 301 302[:[/\large $$_1F_1(a, b; z) = \Gamma(b) e^{ \frac{1}{2} z} \sum_{n=0}^{\infty} 2^{-n}z^n A_n(a,b) E_{b-1+n}(z(\frac{b}{2}-a)) $$] 303[$../equations/hypergeometric_1f1_4.svg]] 304 305with 306 307[:[/\large $$A_0(a,b) = 1, A_1(a,b) = 0, A_2(a,b) = \frac{b}{2}, A_3(a,b)= - \frac{1}{3}(b-2a) $$] 308[$../equations/hypergeometric_1f1_5.svg]] 309 310and 311 312[:[/\large $$(n+1)A_{n+1} = (n+b-1)A_{n-1} - (b-2a)A_{n-2} \quad;\quad n \geq 2 $$] 313[$../equations/hypergeometric_1f1_6.svg]] 314 315With ['E[sub v]] defined as: 316 317[:[/ 318\begin{equation*} 319\begin{split} 320E_v(z) & = z^{-\frac{1}{2}v} J_v(2 \sqrt{z})\\ 321E_v(-z) & = z^{-\frac{1}{2}v} I_v(2 \sqrt{z})\\ 322E_v(0) & = \frac{1}{\Gamma(v + 1)} 323\end{split} 324\end{equation*}] 325[$../equations/hypergeometric_1f1_7.svg]] 326 327This approximation is especially effective when ['a < 0]. 328 329For /a/ and /b/ both large and positive and somewhat similar in size then an expansion in terms of gamma functions can be used [link math_toolkit.hypergeometric.hypergeometric_refs (6)]: 330 331[:[/\large $$_1F_1(a, b; x) = \frac{1}{B(a, b-a)} e^x \sum_{k=0}^{\infty} \frac{(1-a)_k}{k!} \frac{\gamma(b-a+k, x)}{x^{b-a+k}} $$] 332[$../equations/hypergeometric_1f1_8.svg]] 333 334For /z/ large we have the asymptotic expansion: 335 336[:[/\large $$_1F_1(a, b; x) \approx \frac{e^x x^{a-b}}{\Gamma(a)} \sum_{k=0}^{\infty} \frac{(1-a)_k(b-a)_k}{k! x^k} $$] 337[$../equations/hypergeometric_1f1_9.svg]] 338 339which is a special case of the gamma function expansion above. 340 341In the situation where `ab/z` is reasonably small then Luke's rational approximation is used - this is too complex to document 342here, refer either to the code or the original paper [link math_toolkit.hypergeometric.hypergeometric_refs (3)]. 343 344The special case [sub 1]/F/[sub 1](1, /b/, -/z/) is treated via Luke's Pade approximation [link math_toolkit.hypergeometric.hypergeometric_refs (3)]. 345 346That effectively completes the "direct" methods used, the remaining areas are treated indirectly via recurrence relations. 347These require extreme care in their use, as they often change the direction of stability, or else are not stable at all. 348Sometimes this is a well defined and characterized change in direction: for example with /a,b/ and /z/ all positive recurrence on /a/ via 349 350[:[/\large $$(b-a) _1F_1(a-1, b; z) + (2a-b+z) _1F_1(a, b; z) -a _1F_1(a+1, b; z) = 0 $$] 351[$../equations/hypergeometric_1f1_10.svg]] 352 353abruptly changes from stable forwards recursion to stable backwards if /2a-b+z/ changes sign. 354Thus recurrence on /a/, even when [sub 1]/F/[sub 1](/a/+/N/, /b/, /z/) is strictly increasing, needs careful handling as /a/ \u2192 0. 355 356The transitory nature of the stable recurrence relations is well documented in the literature, 357for example [link math_toolkit.hypergeometric.hypergeometric_refs (10)] 358gives many examples, including the anomalous behaviour of recurrence on /a/ and /b/ for large /z/ as first documented by 359Gauchi [link math_toolkit.hypergeometric.hypergeometric_refs (12)]. 360Gauchi also provides the definitive reference on 3-term recurrence relations 361in general in [link math_toolkit.hypergeometric.hypergeometric_refs (11)]. 362 363Unfortunately, not all recurrence stabilities are so well defined. 364For example, when considering [sub 1]F[sub 1](/a/, -/b/, /z/) it would be convenient to use 365the continued fractions associated with the recurrence relations to calculate [sub 1]F[sub 1](/a/, -/b/, /z/) / [sub 1]F[sub 1](/a/, 1-/b/, /z/) 366and then normalize 367either by iterating forwards on /b/ until /b > 0/ and comparison with a reference value, or by using the Wronskian 368 369[:[/\large $$_1F_1(a, b; z) \frac{d}{dz}z^{1-b}_1F_1(1+a-b, 2-b; z) - z^{1-b}_1F_1(1+a-b, 2-b; z)\frac{d}{dz}{_1F_1}(a, b; z) = (1-b)z^{-b}e^z,$$] 370[$../equations/hypergeometric_1f1_11.svg]] 371 372which is of course the well known Miller's method [link math_toolkit.hypergeometric.hypergeometric_refs (12)]. 373 374Unfortunately, stable forwards recursion on /b/ is stable only for /|b| << |z|/, as /|b|/ increases in magnitude it passes through a region 375where recursion is unstable in both directions before eventually becoming stable in the backwards direction (at least for a while!). 376This transition is moderated not by a change of sign in the recurrence coefficients themselves, but by a change in the behaviour of the ['values] of [sub 1]F[sub 1] - 377from alternating in sign when ['|b|] is small to having the same sign when /b/ is larger. During the actual transition, [sub 1]F[sub 1] will either 378pass through a zero or a minima depending on whether b is even or odd (if there is a minima at [sub 1]F[sub 1](a, b, z) then there is necessarily a zero 379at [sub 1]F[sub 1](a+1, b+1, z) as per the [@https://dlmf.nist.gov/13.3#E15 derivative of the function]). 380As a result the behaviour of the recurrence relations 381is almost impossible to reason about in this area, and we are left to using heuristic based approaches to "guess" which region we are in. 382 383In this implementation we use recurrence relations as follows: 384 385For /a,b,z > 0/ and large we can either: 386 387* Make /0 < a < 1/ and /b > z/ and evaluate the defining series directly, or 388* The gamma function approximation has decent convergence and accuracy for /2b - 1 < a < 2b < z/, so we can move /a/ and /b/ into this region, or 389* We can recurse on /b/ alone so that /b-1 < a < b/ and use A&S 13.3.6 as long as /z/ is not too large. 390 391For ['b < 0] and ['a] tiny we would normally use A&S 13.3.6, but when that is non-convergent for some inputs we can use the forward recurrence relation on ['b] to 392calculate the ratio ['[sub 1]F[sub 1](a, -b, z)/[sub 1]F[sub 1](a, 1-b, z)] and then iterate forwards until ['b > 0] and compute a reference value 393and normalize (Millers Method). 394In the same domain when ['b] is near -1 we can use a single backwards recursion on /b/ to compute ['[sub 1]F[sub 1](a, -b, z)] 395from ['[sub 1]F[sub 1](a, 2-b, z)] and ['[sub 1]F[sub 1](/a/, 1-/b/, /z/)] even though technically we are recursing in the unstable direction. 396 397For ['[sub 1]F[sub 1](-N, b, z)] and integer /N/ then backwards recursion from ['[sub 1]F[sub 1](0, b, z)] and ['[sub 1]F[sub 1](-1, b, z)] works well. 398 399For /a/ or /b/ < 0 and if all the direct methods have failed, then we use various fallbacks: 400 401For ['[sub 1]F[sub 1](-a, b, z)] we can use backwards recursion on /a/ as long as ['b > z], otherwise a more complex scheme is required 402which starts from ['[sub 1]F[sub 1](-a + N, b + M, z)], and recurses backwards in up to 3 stages: first on /a/, then on /a/ and /b/ together, 403and finally on /b/ alone. 404 405For /b < 0/ we have no good methods in some domains (see the unsolved issues above). 406However in some circumstances we can either use: 407 408* 3-stage backwards recursion on both /a/, /a/ and /b/ and then /b/ as above. 409* Calculate the ratio ['[sub 1]F[sub 1](a, b, z) / ['[sub 1]F[sub 1](a-1, b-1, z)]] via backwards recurrence when z is small, and then normalize via the Wronskian above (Miller's method). 410* Calculate the ratio ['[sub 1]F[sub 1](a, b, z) / ['[sub 1]F[sub 1](a+1, b+1, z)]] via forwards recurrence when z is large, and then normalize by iterating until b > 1 and comparing to a reference value. 411 412The latter two methods use a lookup table to determine whether inputs are in either of the domains or neither. Unfortunately the methods are basically 413limited to double precision: calculation of the ratios require iteration ['towards] the no-mans-land between the two methods where iteration is unstable in 414both directions. As a result, only low-precision results which require few iterations to calculate the ratio are available. 415 416If all else fails, then we use a checked series summation which will raise an __evaluation_error if more than half the digits 417are destroyed by cancellation. 418 419[endsect] [/section:hyper_geometric_1f1 Hypergeometric [sub 1]/F/[sub 1]] 420 421[section:hypergeometric_pfq Hypergeometric [sub p]F[sub q]] 422 423 #include <boost/math/special_functions/hypergeometric_pFq.hpp> 424 425 namespace boost { namespace math { 426 427 template <class Seq, class Real> 428 ``__sf_result`` hypergeometric_pFq(const Seq& aj, const Seq& bj, const Real& z, Real* p_abs_error, const Policy& pol); 429 430 template <class Seq, class Real> 431 ``__sf_result`` hypergeometric_pFq(const Seq& aj, const Seq& bj, const Real& z, Real* p_abs_error = 0); 432 433 template <class R, class Real> 434 ``__sf_result`` hypergeometric_pFq(const std::initializer_list<R>& aj, const std::initializer_list<R>& bj, const Real& z, Real* p_abs_error, const Policy& pol); 435 436 template <class R, class Real> 437 ``__sf_result`` hypergeometric_pFq(const std::initializer_list<R>& aj, const std::initializer_list<R>& bj, const Real& z, Real* p_abs_error = 0); 438 439 template <class Seq, class Real, class Policy> 440 Real hypergeometric_pFq_precision(const Seq& aj, const Seq& bj, Real z, unsigned digits10, double timeout, const Policy& pol); 441 442 template <class Seq, class Real> 443 Real hypergeometric_pFq_precision(const Seq& aj, const Seq& bj, const Real& z, unsigned digits10, double timeout = 0.5); 444 445 template <class Real, class Policy> 446 Real hypergeometric_pFq_precision(const std::initializer_list<Real>& aj, const std::initializer_list<Real>& bj, const Real& z, unsigned digits10, double timeout, const Policy& pol); 447 448 template <class Real> 449 Real hypergeometric_pFq_precision(const std::initializer_list<Real>& aj, const std::initializer_list<Real>& bj, const Real& z, unsigned digits10, double timeout = 0.5); 450 451 }} // namespaces 452 453[h4 Description] 454 455The function `hypergeometric_pFq` returns the result of: 456 457[:[/\Large $$ {}_pF_q\left(\left\{a_1, \ldots, a_p\right\}, \left\{b_1, \ldots, b_q\right\}; z\right)=\sum_{n=0}^{\infty} \frac{\prod_{j=1}^{p}\left(a_j\right)_n}{\prod_{j=1}^{q}\left(b_{j}\right)_n} \frac{z^n}{n!} $$] 458[$../equations/hypergeometric_pfq_1.svg]] 459 460It is most naturally used via initializer lists as in: 461 462 double h = hypergeometric_pFq({2,3,4}, {5,6,7,8}, z); // Calculate 3F4 463 464The optional ['p_abs_error] argument calculates an estimate of the absolute error in the result from the 465L1 norm of the sum, plus some other factors (see implementation below). 466 467You should divide this value by the result to get an estimate of ['relative error]. 468 469It should be noted that the error estimates are rather crude - the error can be significantly underestimated 470in some circumstances, and over-estimated in others. 471 472The `hypergeometric_pFq_precision` functions will calculate `pFq` to a specified number of decimal places, and if `timeout` 473is reached then they raise an __evaluation_error. Note that while these functions are defined as templates, they 474require type `Real` to be a *variable-precision* floating-point type: in practice the only type currently supported 475(as of July 2019) is `boost::multiprecision::mpfr_float`. Typical usage would be: 476 477 478 typedef boost::multiprecision::mpfr_float mp_type; 479 // 480 // Calculate 2F2 to 20 decimal places using a 10 second timeout: 481 // 482 mp_type result = boost::math::hypergeometric_pFq_precision( 483 { mp_type(a1), mp_type(a2) }, { mp_type(b1), mp_type(b2) }, mp_type(z), 20, 10.0 484 ); 485 // 486 // Convert the result back to a double: 487 // 488 double d_result = static_cast<double>(result); 489 490[h4 Implementation] 491 492This function is implemented by direct summation of the series; summation normally starts with the zeroth term, 493but will skip forward and sum outward (ie in both forward and backward directions) from some term /N/ when 494required. This is particularly important when some of the /b/ arguments are negative, as in this situation 495the sum undergoes "false-convergence", and then diverges again as each b[sub j] passes the origin. Consequently, 496were you to plot the magnitude of the terms in the sum, you would see them pass through a series of 497minima and maxima before finally converging to zero at infinity. For some values of /p/ and /q/ we 498can compute where the maxima occur, and therefore sum only the terms that will have an impact on the 499result. For other /p/ and /q/ values, predicting the exact locations of the maxima is not so easy, and we 500simply restart summation at the point where each b[sub j] passes the origin: this will eventually reach 501all the significant terms of the sum, but the key word may well be "eventually". 502 503The error term /p_abs_error/ is normally simply the sum of the absolute values of each term multiplied 504by machine epsilon for type `Real`. However, 505if it was necessary for the summation to skip forward, then /p_abs_error/ is adjusted to account for the 506error inherent in calculating the N'th term via logarithms. 507 508[endsect] [/section:pFq Hypergeometric [sub p]F[sub q]] 509 510[section:hypergeometric_refs Hypergeometric References] 511 512# Beals, Richard, and Roderick Wong. ['Special functions: a graduate text.] Vol. 126. Cambridge University Press, 2010. 513# Pearson, John W., Sheehan Olver, and Mason A. Porter. ['Numerical methods for the computation of the confluent and Gauss hypergeometric functions.] Numerical Algorithms 74.3 (2017): 821-866. 514# Luke, Yudell L. ['Algorithms for Rational Approximations for a Confluent Hypergeometric Function II.] MISSOURI UNIV KANSAS CITY DEPT OF MATHEMATICS, 1976. 515# Derezinski, Jan. ['Hypergeometric type functions and their symmetries.] Annales Henri Poincar[eacute]. Vol. 15. No. 8. Springer Basel, 2014. 516# Keith E. Muller ['Computing the confluent hypergeometric function, M(a, b, x)]. Numer. Math. 90: 179-196 (2001). 517# Carlo Morosi, Livio Pizzocchero. ['On the expansion of the Kummer function in terms of incomplete Gamma functions.] Arch. Inequal. Appl. 2 (2004), 49-72. 518# Jose Luis Lopez, Nico M. Temme. ['Asymptotics and numerics of polynomials used in Tricomi and Buchholz expansions of Kummer functions]. Numerische Mathematik, August 2010. 519# Javier Sesma. ['The Temme's sum rule for confluent hypergeometric functions revisited]. Journal of Computational and Applied Mathematics 163 (2004) 429-431. 520# Javier Segura, Nico M. Temme. ['Numerically satisfactory solutions of Kummer recurrence relations]. Numer. Math. (2008) 111:109-119. 521# Alfredo Deano, Javier Segura. ['Transitory Minimal Solutions Of Hypergeometric Recursions And Pseudoconvergence of Associated Continued Fractions]. Mathematics of Computation, Volume 76, Number 258, April 2007. 522# W. Gautschi. ['Computational aspects of three-term recurrence relations]. SIAM Review 9, no.1 (1967) 24-82. 523# W. Gautschi. ['Anomalous convergence of a continued fraction for ratios of Kummer functions]. Math. Comput., 31, no.140 (1977) 994-999. 524# British Association for the Advancement of Science: ['Bessel functions, Part II, Mathematical Tables vol. X]. Cambridge (1952). 525 526 527[endsect] [/section:hypergeometric_refs Hypergeometric References] 528 529[endsect] [/section:hypergeometric Hypergeometric Functions] 530 531[/ Copyright 2017 John Maddock. 532 Distributed under the Boost Software License, Version 1.0. 533 (See accompanying file LICENSE_1_0.txt or copy at 534 http://www.boost.org/LICENSE_1_0.txt). 535]
