#include <boost/math/special_functions/hypergeometric_pFq.hpp>

namespace boost { namespace math {

template <class Seq, class Real>
calculated-result-type hypergeometric_pFq(const Seq& aj, const Seq& bj, const Real& z, Real* p_abs_error, const Policy& pol);

template <class Seq, class Real>
calculated-result-type hypergeometric_pFq(const Seq& aj, const Seq& bj, const Real& z, Real* p_abs_error = 0);

template <class R, class Real>
calculated-result-type 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>
calculated-result-type 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

Description

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);

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