123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535 |
- [section:hypergeometric Hypergeometric Functions]
- [section:hypergeometric_1f0 Hypergeometric [sub 1]/F/[sub 0] ]
- #include <boost/math/special_functions/hypergeometric_1F0.hpp>
- namespace boost { namespace math {
- template <class T1, class T2>
- ``__sf_result`` hypergeometric_1F0(T1 a, T2 z);
- template <class T1, class T2, class ``__Policy``>
- ``__sf_result`` hypergeometric_1F0(T1 a, T2 z, const ``__Policy``&);
- }} // namespaces
- [h4 Description]
- The function `hypergeometric_1F0` returns the result of:
- [equation hyper_1f0]
- The return type of these functions is computed using the __arg_promotion_rules
- when `T1` and `T2` are different types.
- [optional_policy]
- The functions return the result of __domain_error whenever the result is
- undefined or complex. This occurs for `z == 1` or `1 - z < 0` and `a` not an integer.
- [h4 Implementation]
- The implementation is trivial:
- [expression ['[sub 1]F[sub 0](a, z) = (1-z)[super -a]]]
- [endsect] [/section:hyper_geometric_1f0 Hypergeometric [sub 1]/F/[sub 0] ]
- [section:hypergeometric_0f1 Hypergeometric [sub 0]/F/[sub 1] ]
- #include <boost/math/special_functions/hypergeometric_0F1.hpp>
- namespace boost { namespace math {
- template <class T1, class T2>
- ``__sf_result`` hypergeometric_0F1(T1 b, T2 z);
- template <class T1, class T2, class ``__Policy``>
- ``__sf_result`` hypergeometric_0F1(T1 b, T2 z, const ``__Policy``&);
- }} // namespaces
- [h4 Description]
- The function `hypergeometric_0F1` returns the result of
- [equation hyper_0f1]
- The return type of these functions is computed using the __arg_promotion_rules
- when `T1` and `T2` are different types.
- [optional_policy]
- The functions return the result of __domain_error whenever the result is
- undefined or complex.
- This occurs only when `b` is an integer <= 0.
- [h4 Implementation]
- The function is implemented via its defining series wherever convergent.
- For a divergent series we use the continued fraction as long as the result is not too small:
- [equation hyper_0f1_cf]
- and one of the Bessel function relations otherwise:
- [equation hyper_0f1_bessel_j]
- [equation hyper_0f1_bessel_i]
- [endsect] [/section:hypergeometric_0f1 Hypergeometric [sub 0]/F/[sub 1] ]
- [section:hypergeometric_2f0 Hypergeometric [sub 2]/F/[sub 0]]
- #include <boost/math/special_functions/hypergeometric_2F0.hpp>
- namespace boost { namespace math {
- template <class T1, class T2, class T3>
- ``__sf_result`` hypergeometric_2F0(T1 a1, T2 a2, T3 z);
- template <class T1, class T2, class T3, class ``__Policy``>
- ``__sf_result`` hypergeometric_2F0(T1 a1, T2 a2, T3 z, const ``__Policy``&);
- }}
- [h4 Description]
- The function `hypergeometric_2F0` returns the result of
- [equation hyper_2f0]
- The return type of these functions is computed using the __arg_promotion_rules
- when `T1` and `T2` are different types.
- [optional_policy]
- The functions return the result of __domain_error whenever the result is
- undefined or complex. The valid domain for this function occurs only when one of `a1` or
- `a2` is a negative integer: ie the polynomial case.
- [h4 Implementation]
- When `a1 == a2 - 0.5` then the function is implemented in terms of the Hermite polynomial:
- [equation hyper_2f0_hermite]
- When both `a1` and `a2` are integers then the function is implemented in terms of the associated-Laguerre polynomial:
- [equation hyper_2f0_laguerre]
- If the defining series is divergent, we use the continued fraction
- [equation hyper_2f0_cf]
- Otherwise we use the defining series.
- [endsect] [/section:hyper_geometric_1f0 Hypergeometric [sub 2]/F/[sub 0]]
- [section:hypergeometric_1f1 Hypergeometric [sub 1]/F/[sub 1]]
- #include <boost/math/special_functions/hypergeometric_1F1.hpp>
- namespace boost { namespace math {
- template <class T1, class T2, class T3>
- ``__sf_result`` hypergeometric_1F1(T1 a, T2 b, T3 z);
- template <class T1, class T2, class T3, class ``__Policy``>
- ``__sf_result`` hypergeometric_1F1(T1 a, T2 b, T3 z, const ``__Policy``&);
- template <class T1, class T2, class T3>
- ``__sf_result`` hypergeometric_1F1_regularized(T1 a, T2 b, T3 z);
- template <class T1, class T2, class T3, class ``__Policy``>
- ``__sf_result`` hypergeometric_1F1_regularized(T1 a, T2 b, T3 z, const ``__Policy``&);
- template <class T1, class T2, class T3>
- ``__sf_result`` log_hypergeometric_1F1(T1 a, T2 b, T3 z);
- template <class T1, class T2, class T3, class ``__Policy``>
- ``__sf_result`` log_hypergeometric_1F1(T1 a, T2 b, T3 z, const ``__Policy``&);
- template <class T1, class T2, class T3>
- ``__sf_result`` log_hypergeometric_1F1(T1 a, T2 b, T3 z, int* sign);
- template <class T1, class T2, class T3, class ``__Policy``>
- ``__sf_result`` log_hypergeometric_1F1(T1 a, T2 b, T3 z, int* sign, const ``__Policy``&);
- }} // namespaces
- [h4 Description]
- [? __build_html '''<?dbhtml-include href="'''__base_path__'''/graphs/hypergeometric_1f1/script_include.html"?>''']
- The function `hypergeometric_1F1(a, b, z)` returns the non-singular solution to
- [@https://en.wikipedia.org/wiki/Confluent_hypergeometric_function Kummer's equation]
- [:[/\large $$z \frac{d^2 w}{d z^2} + (b-z) \frac{dw}{dz} - aw = 0 $$]
- [$../equations/hypergeometric_1f1_2.svg]]
- which for |/z/| < 1 has the hypergeometric series expansion
- [:[$../equations/hypergeometric_1f1_1.svg]]
- where (/a/)[sub /n/] denotes rising factorial.
- This function has the same definition as Mathematica's `Hypergeometric1F1[a, b, z]` and Maple's `KummerM(a, b, z)`.
- The "regularized" versions of the function return:
- [:[/ \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!} $$]
- [$../equations/hypergeometric_1f1_17.svg]]
- The "log" versions of the function return:
- [:[/ \Large $$ \ln(|_1F_1(a, b, z)|) $$ ]
- [$../equations/hypergeometric_1f1_18.svg]]
- When 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/).
- Both the regularized and the logarithmic versions of these functions return results without the spurious
- under/overflow that plague naive implementations.
- [h4 Known Issues]
- This function is still very much the subject of active research,
- and a full range of methods capable of calculating the function
- over its entire domain are not yet available.
- We have worked extremely hard to ensure that problem domains result in an exception being thrown
- (an __evaluation_error) rather than return a garbage result.
- Domains that are still unsolved include:
- [table
- [[domain][comment][example]]
- [[ [/\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)]]
- [[ [/\large $$_1F_1(-a, -b; z)$$] [$../equations/hypergeometric_1f1_14.svg]][ /a < b/, /b > z/, this is really the same domain as above.][ ]]
- [[ [/\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)]]
- [[ [/\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
- are some values where either that series is non-convergent (most particularly for /b/ < 0)
- 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)]]
- ]
- The following graph illustrates the problem area for /b < 0/ and /a,z > 0/:
- [? __build_html '''<?dbhtml-include href="'''__base_path__'''/graphs/hypergeometric_1f1/negative_b_incalculable.html"?>''']
- [?! __build_html [$../graphs/hypergeometric_1f1/negative_b_incalculable.png]]
- [h4 Testing]
- There are 3 main groups of tests:
- * Spot tests for special inputs with known values.
- * Sanity checks which use integrals of hypergeometric functions of known value. These are particularly useful
- since for fixed ['a] and ['b] they evaluate ['[sub 1]F[sub 1](a,b,z)] for all /z/.
- * Testing against values precomputed using arbitrary precision arithmetic and the /pFq/ series.
- We 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
- is also used for the error rate plots below and has been extremely useful in fine-tuning the implementation.
- It should be noted however, that there are some domains for large negative /a/ and /b/ that have poor test coverage as we were
- simply unable to compute test values in reasonable time: it is not uncommon for the /pFq/ series to cancel many hundreds of digits
- and sometimes into the thousands of digits.
- [h4 Errors]
- We divide the domain of [sub 1]/F/[sub 1] up into several domains to explain the error rates.
- In the following graphs we ran 100000 random test cases over each domain, note that the scatter plots show only the very worst errors
- as otherwise the graphs are both incomprehensible and virtually unplottable (as in sudden browser death).
- For 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
- for performance reasons. Even so the errors rates are broadly acceptable, while the scatter graph shows where the worst errors are located:
- [$../graphs/hypergeometric_1f1/positive_abz_bins.svg]
- [? __build_html '''<?dbhtml-include href="'''__base_path__'''/graphs/hypergeometric_1f1/positive_abz.html"?>''']
- [?! __build_html [$../graphs/hypergeometric_1f1/positive_abz.png]]
- For -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:
- [$../graphs/hypergeometric_1f1/negative_a_bins.svg]
- [? __build_html '''<?dbhtml-include href="'''__base_path__'''/graphs/hypergeometric_1f1/negative_a.html"?>''']
- [?! __build_html [$../graphs/hypergeometric_1f1/negative_a.png]]
- For -1000 < /b/ < 0 and /a/,/z/ > 0 the errors are mostly low at double precision except near the "unimplementable zone".
- Note however, that the some of the methods used here fail to converge to much better than double precision.
- [$../graphs/hypergeometric_1f1/negative_b_bins.svg]
- [? __build_html '''<?dbhtml-include href="'''__base_path__'''/graphs/hypergeometric_1f1/negative_b.html"?>''']
- [?! __build_html [$../graphs/hypergeometric_1f1/negative_b.png]]
- For 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/:
- [$../graphs/hypergeometric_1f1/negative_ab_bins.svg]
- [? __build_html '''<?dbhtml-include href="'''__base_path__'''/graphs/hypergeometric_1f1/negative_ab.html"?>''']
- [?! __build_html [$../graphs/hypergeometric_1f1/negative_ab.png]]
- [h4 Implementation]
- The implementation for this function is remarkably complex;
- readers will have to refer to the code for the transitions between methods, as we can only give an overview here.
- In almost all cases where /z < 0/ we use [@https://en.wikipedia.org/wiki/Confluent_hypergeometric_function Kummer's relation]
- to make /z/ positive:
- [:[/\large $$_1F_1(a, b; -z) = e^{-z} {_1F_1}(a, b; z)$$]
- [$../equations/hypergeometric_1f1_12.svg]]
- The main series representation
- [:[$../equations/hypergeometric_1f1_1.svg]]
- is used only when
- * we are sure that it is convergent and not subject to excessive cancellation, and
- * there is no other better method available.
- The 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
- to life again as /b/ crosses the origin.
- A&S 13.3.6 gives
- [:[/\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})$$]
- [$../equations/hypergeometric_1f1_3.svg]]
- which is particularly useful for ['a [cong] b] and ['z > 0], or /a/ \u226a 1, and ['z < 0].
- Then we have Tricomi's approximation (given in simplified form in A&S 13.3.7) [link math_toolkit.hypergeometric.hypergeometric_refs (7)]
- [:[/\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)) $$]
- [$../equations/hypergeometric_1f1_4.svg]]
- with
- [:[/\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) $$]
- [$../equations/hypergeometric_1f1_5.svg]]
- and
- [:[/\large $$(n+1)A_{n+1} = (n+b-1)A_{n-1} - (b-2a)A_{n-2} \quad;\quad n \geq 2 $$]
- [$../equations/hypergeometric_1f1_6.svg]]
- With ['E[sub v]] defined as:
- [:[/
- \begin{equation*}
- \begin{split}
- E_v(z) & = z^{-\frac{1}{2}v} J_v(2 \sqrt{z})\\
- E_v(-z) & = z^{-\frac{1}{2}v} I_v(2 \sqrt{z})\\
- E_v(0) & = \frac{1}{\Gamma(v + 1)}
- \end{split}
- \end{equation*}]
- [$../equations/hypergeometric_1f1_7.svg]]
- This approximation is especially effective when ['a < 0].
- For /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)]:
- [:[/\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}} $$]
- [$../equations/hypergeometric_1f1_8.svg]]
- For /z/ large we have the asymptotic expansion:
- [:[/\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} $$]
- [$../equations/hypergeometric_1f1_9.svg]]
- which is a special case of the gamma function expansion above.
- In the situation where `ab/z` is reasonably small then Luke's rational approximation is used - this is too complex to document
- here, refer either to the code or the original paper [link math_toolkit.hypergeometric.hypergeometric_refs (3)].
- The special case [sub 1]/F/[sub 1](1, /b/, -/z/) is treated via Luke's Pade approximation [link math_toolkit.hypergeometric.hypergeometric_refs (3)].
- That effectively completes the "direct" methods used, the remaining areas are treated indirectly via recurrence relations.
- These require extreme care in their use, as they often change the direction of stability, or else are not stable at all.
- Sometimes this is a well defined and characterized change in direction: for example with /a,b/ and /z/ all positive recurrence on /a/ via
- [:[/\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 $$]
- [$../equations/hypergeometric_1f1_10.svg]]
- abruptly changes from stable forwards recursion to stable backwards if /2a-b+z/ changes sign.
- Thus recurrence on /a/, even when [sub 1]/F/[sub 1](/a/+/N/, /b/, /z/) is strictly increasing, needs careful handling as /a/ \u2192 0.
- The transitory nature of the stable recurrence relations is well documented in the literature,
- for example [link math_toolkit.hypergeometric.hypergeometric_refs (10)]
- gives many examples, including the anomalous behaviour of recurrence on /a/ and /b/ for large /z/ as first documented by
- Gauchi [link math_toolkit.hypergeometric.hypergeometric_refs (12)].
- Gauchi also provides the definitive reference on 3-term recurrence relations
- in general in [link math_toolkit.hypergeometric.hypergeometric_refs (11)].
- Unfortunately, not all recurrence stabilities are so well defined.
- For example, when considering [sub 1]F[sub 1](/a/, -/b/, /z/) it would be convenient to use
- the 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/)
- and then normalize
- either by iterating forwards on /b/ until /b > 0/ and comparison with a reference value, or by using the Wronskian
- [:[/\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,$$]
- [$../equations/hypergeometric_1f1_11.svg]]
- which is of course the well known Miller's method [link math_toolkit.hypergeometric.hypergeometric_refs (12)].
- Unfortunately, stable forwards recursion on /b/ is stable only for /|b| << |z|/, as /|b|/ increases in magnitude it passes through a region
- where recursion is unstable in both directions before eventually becoming stable in the backwards direction (at least for a while!).
- This 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] -
- from 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
- pass 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
- at [sub 1]F[sub 1](a+1, b+1, z) as per the [@https://dlmf.nist.gov/13.3#E15 derivative of the function]).
- As a result the behaviour of the recurrence relations
- is almost impossible to reason about in this area, and we are left to using heuristic based approaches to "guess" which region we are in.
- In this implementation we use recurrence relations as follows:
- For /a,b,z > 0/ and large we can either:
- * Make /0 < a < 1/ and /b > z/ and evaluate the defining series directly, or
- * 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
- * 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.
- For ['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
- calculate 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
- and normalize (Millers Method).
- In 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)]
- from ['[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.
- For ['[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.
- For /a/ or /b/ < 0 and if all the direct methods have failed, then we use various fallbacks:
- For ['[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
- which 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,
- and finally on /b/ alone.
- For /b < 0/ we have no good methods in some domains (see the unsolved issues above).
- However in some circumstances we can either use:
- * 3-stage backwards recursion on both /a/, /a/ and /b/ and then /b/ as above.
- * 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).
- * 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.
- The latter two methods use a lookup table to determine whether inputs are in either of the domains or neither. Unfortunately the methods are basically
- limited to double precision: calculation of the ratios require iteration ['towards] the no-mans-land between the two methods where iteration is unstable in
- both directions. As a result, only low-precision results which require few iterations to calculate the ratio are available.
- If all else fails, then we use a checked series summation which will raise an __evaluation_error if more than half the digits
- are destroyed by cancellation.
- [endsect] [/section:hyper_geometric_1f1 Hypergeometric [sub 1]/F/[sub 1]]
- [section:hypergeometric_pfq Hypergeometric [sub p]F[sub q]]
- #include <boost/math/special_functions/hypergeometric_pFq.hpp>
- namespace boost { namespace math {
- template <class Seq, class Real>
- ``__sf_result`` hypergeometric_pFq(const Seq& aj, const Seq& bj, const Real& z, Real* p_abs_error, const Policy& pol);
- template <class Seq, class Real>
- ``__sf_result`` hypergeometric_pFq(const Seq& aj, const Seq& bj, const Real& z, Real* p_abs_error = 0);
- template <class R, class Real>
- ``__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);
- template <class R, class Real>
- ``__sf_result`` hypergeometric_pFq(const std::initializer_list<R>& aj, const std::initializer_list<R>& bj, const Real& z, Real* p_abs_error = 0);
- template <class Seq, class Real, class Policy>
- Real hypergeometric_pFq_precision(const Seq& aj, const Seq& bj, Real z, unsigned digits10, double timeout, const Policy& pol);
- template <class Seq, class Real>
- Real hypergeometric_pFq_precision(const Seq& aj, const Seq& bj, const Real& z, unsigned digits10, double timeout = 0.5);
- template <class Real, class Policy>
- 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);
- template <class Real>
- 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);
- }} // namespaces
- [h4 Description]
- The function `hypergeometric_pFq` returns the result of:
- [:[/\Large $$ _pF_q(\{a_1...a_n\}, \{b_1...b_n\}; z) = \sum_{k=0}^{\infty} \frac{\Pi_{j=1}^n(a_j)_n z^k}{\Pi_{j=1}^n(b_j)_n k!} $$]
- [$../equations/hypergeometric_pfq_1.svg]]
- It is most naturally used via initializer lists as in:
- double h = hypergeometric_pFq({2,3,4}, {5,6,7,8}, z); // Calculate 3F4
- The optional ['p_abs_error] argument calculates an estimate of the absolute error in the result from the
- L1 norm of the sum, plus some other factors (see implementation below).
- You should divide this value by the result to get an estimate of ['relative error].
- It should be noted that the error estimates are rather crude - the error can be significantly underestimated
- in some circumstances, and over-estimated in others.
- The `hypergeometric_pFq_precision` functions will calculate `pFq` to a specified number of decimal places, and if `timeout`
- is reached then they raise an __evaluation_error. Note that while these functions are defined as templates, they
- require type `Real` to be a *variable-precision* floating-point type: in practice the only type currently supported
- (as of July 2019) is `boost::multiprecision::mpfr_float`. Typical usage would be:
- typedef boost::multiprecision::mpfr_float mp_type;
- //
- // Calculate 2F2 to 20 decimal places using a 10 second timeout:
- //
- mp_type result = boost::math::hypergeometric_pFq_precision(
- { mp_type(a1), mp_type(a2) }, { mp_type(b1), m_type(b2) }, mp_type(z), 20, 10.0
- );
- //
- // Convert the result back to a double:
- //
- double d_result = static_cast<double>(result);
- [h4 Implementation]
- This function is implemented by direct summation of the series; summation normally starts with the zeroth term,
- but will skip forward and sum outward (ie in both forward and backward directions) from some term /N/ when
- required. This is particularly important when some of the /b/ arguments are negative, as in this situation
- the sum undergoes "false-convergence", and then diverges again as each b[sub j] passes the origin. Consequently,
- were you to plot the magnitude of the terms in the sum, you would see them pass through a series of
- minima and maxima before finally converging to zero at infinity. For some values of /p/ and /q/ we
- can compute where the maxima occur, and therefore sum only the terms that will have an impact on the
- result. For other /p/ and /q/ values, predicting the exact locations of the maxima is not so easy, and we
- simply restart summation at the point where each b[sub j] passes the origin: this will eventually reach
- all the significant terms of the sum, but the key word may well be "eventually".
- The error term /p_abs_error/ is normally simply the sum of the absolute values of each term multiplied
- by machine epsilon for type `Real`. However,
- if it was necessary for the summation to skip forward, then /p_abs_error/ is adjusted to account for the
- error inherent in calculating the N'th term via logarithms.
- [endsect] [/section:pFq Hypergeometric [sub p]F[sub q]]
- [section:hypergeometric_refs Hypergeometric References]
- # Beals, Richard, and Roderick Wong. ['Special functions: a graduate text.] Vol. 126. Cambridge University Press, 2010.
- # 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.
- # Luke, Yudell L. ['Algorithms for Rational Approximations for a Confluent Hypergeometric Function II.] MISSOURI UNIV KANSAS CITY DEPT OF MATHEMATICS, 1976.
- # Derezinski, Jan. ['Hypergeometric type functions and their symmetries.] Annales Henri Poincar[eacute]. Vol. 15. No. 8. Springer Basel, 2014.
- # Keith E. Muller ['Computing the confluent hypergeometric function, M(a, b, x)]. Numer. Math. 90: 179-196 (2001).
- # Carlo Morosi, Livio Pizzocchero. ['On the expansion of the Kummer function in terms of incomplete Gamma functions.] Arch. Inequal. Appl. 2 (2004), 49-72.
- # Jose Luis Lopez, Nico M. Temme. ['Asymptotics and numerics of polynomials used in Tricomi and Buchholz expansions of Kummer functions]. Numerische Mathematik, August 2010.
- # Javier Sesma. ['The Temme's sum rule for confluent hypergeometric functions revisited]. Journal of Computational and Applied Mathematics 163 (2004) 429-431.
- # Javier Segura, Nico M. Temme. ['Numerically satisfactory solutions of Kummer recurrence relations]. Numer. Math. (2008) 111:109-119.
- # 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.
- # W. Gautschi. ['Computational aspects of three-term recurrence relations]. SIAM Review 9, no.1 (1967) 24-82.
- # W. Gautschi. ['Anomalous convergence of a continued fraction for ratios of Kummer functions]. Math. Comput., 31, no.140 (1977) 994-999.
- # British Association for the Advancement of Science: ['Bessel functions, Part II, Mathematical Tables vol. X]. Cambridge (1952).
- [endsect] [/section:hypergeometric_refs Hypergeometric References]
- [endsect] [/section:hypergeometric Hypergeometric Functions]
- [/ Copyright 2017 John Maddock.
- Distributed under the Boost Software License, Version 1.0.
- (See accompanying file LICENSE_1_0.txt or copy at
- http://www.boost.org/LICENSE_1_0.txt).
- ]
|