I'm looking for a method to efficiently compute a numerical approximation of $$F^n_D(x_1,\ldots,x_n) = \sum_{m=0}^{\infty} \sum_{i_1 +\ldots+i_n=m}\frac{(a)_{m}(b_1)_{i_1}\ldots (b_n)_{i_n}}{(c)_{m}i_1!\ldots i_n!}x_1^{i_1}\ldots x_n^{i_n}$$
I'm restricted to the cases where $$\left\{\begin{array}{c} (x_1,\ldots,x_n) \in (-1,1)^{n}\\\\ (b_1,\ldots,b_n) \in \mathbb{R}^{+}\\\\ c=\sum_{i=0}^n b_i \\\\ a=\sum_{i=0}^{p<n} b_i \end{array}\right. $$ which may help speed up the computation. Also, typically $n \sim 6$.
Since $c > a > 0$, I'm aware of the integral representation $$F^n_D(x_1,\ldots,x_n)=\frac{\Gamma(c)}{\Gamma(a)\Gamma(c-a)}\int_0^1 \frac{t^{a-1}t^{a-c-1}}{(1-x_1t)^{b_1}\ldots (1-x_nt)^{b_n}}\mathrm{dt}$$
Currently, that is the fastest way I know of computing the series is to evaluate the integral directly by quadrature. The series under its current form converges very slowly due to the combinatorial explosion of terms. Writing the terms of the series as a graph, I tried to explore the graph with a heap, adding the largest terms first, but even that is much slower than the integral.
This article introduces an alternate expression of the series, but my understanding is that it hides a sum over the partition number, which grows exponentially with $m$ rather than polynomially as in the original expression.
Is there any fast method to approximate the series when all parameters and variables are real?
