...one of the most highly
regarded and expertly designed C++ library projects in the
world.

— Herb Sutter and Andrei
Alexandrescu, C++
Coding Standards

#include <boost/math/special_functions/hypergeometric_pFq.hpp> namespace boost { namespace math { template <class Seq, class Real>calculated-result-typehypergeometric_pFq(const Seq& aj, const Seq& bj, const Real& z, Real* p_abs_error, const Policy& pol); template <class Seq, class Real>calculated-result-typehypergeometric_pFq(const Seq& aj, const Seq& bj, const Real& z, Real* p_abs_error = 0); template <class R, class Real>calculated-result-typehypergeometric_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>calculated-result-typehypergeometric_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

The function `hypergeometric_pFq`

returns the result of:

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), mp_type(b2) }, mp_type(z), 20, 10.0 ); // // Convert the result back to a double: // double d_result = static_cast<double>(result);

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_{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_{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.