# Evaluating Lauricella functions in Mathematica?

Consider the Lauricella functions. The Mathematica documentation does not contain anything regarding Lauricella. Is there any known way (e.g. an unofficial package) to evaluate Lauricella functions in Mathematica? Thanks for any suggestion!

• Those are difficult to evaluate numerically in my experience. I have an experimental routine for $F_D$, but it's very finicky. Aug 24, 2017 at 15:03
• @J.M. I believe knowledge of your routine might make the world a better place! Could you perhaps post it, in case if it is not a secret? Aug 24, 2017 at 15:10
• I meant "finicky" in the sense of "difficult to use except by me and actual special functions experts"; so, it might take a while before I can post something I won't be embarrassed to show publicly. Hopefully you can wait. Aug 24, 2017 at 15:14
• @J.M. Yes, certainly! Having such a routine sounds very exciting and promising for the future. I'll be looking forward to it when you decide that it is ready to be posted. Aug 24, 2017 at 15:22
• I'm not an expert; that's why I deliberately separated "me" in my last comment. :) I am unfortunately not current on research on multivariate hypergeometrics; that sounds like a good question to ask on MO. Aug 24, 2017 at 17:05

I offer the following tentative answer with the same caveat I gave in the comments and this answer: the method is a bit finicky and not yet in a stage where it can be fully automatized, without needing manual intervention. Nevertheless, I am hoping someone more knowledgeable than me can build on this.

In fact, the method in my answer that I linked to (i.e. Cuyt's transformation + Wynn's $\varepsilon$ algorithm) is also perfectly admissible for computing the Lauricella functions; one merely needs to change the function describing the summand in the infinite sum definition, and then sum over the results of FrobeniusSolve[ConstantArray[1, n], k], where n is the number of variables.

However, I will instead present a special method for Lauricella $F_D$. The method is rooted in this paper by van Laarhoven and Kalker, who derives a compact closed form for the transformation of the multiple sum defining $F_D$ into a single one. Here is a Mathematica routine that expresses their notation compactly:

lauricellaTerm[m_Integer?Positive, a_, b_?VectorQ, c_, z_?VectorQ] :=
(Pochhammer[a, m]/Pochhammer[c, m]) *
Sum[BellY[m, k, Table[(j - 1)! b.z^j, {j, m}]], {k, m}]/m!


This can now be used within NSum[]; e.g.

With[{a = 2, b = {1, 2/3}, c = 3, z = {1, 2}/5},
1 + NSum[lauricellaTerm[m, a, b, c, z], {m, 1, ∞}, Compiled -> False,
Method -> {"WynnEpsilon", "ExtraTerms" -> 12},
NSumTerms -> 12, WorkingPrecision -> 30]]
1.4418435572012925151


Recalling that the bivariate case of $F_D$ is Appell's $F_1$:

N[AppellF1[2, 1, 2/3, 3, 1/5, 2/5], 20]
1.4418435572012925152


which is in excellent agreement. However, note that much higher precision was needed for the NSum[] evaluation; this is due to the unavoidable numerical cancellation during the application of the $\varepsilon$ algorithm. (The finicky part I was talking about is that I haven't yet found a good way to automatically adjust the precision so that the numerical evaluation is accurate, yet still reasonably efficient. So, one still needs to manually adjust the number of terms taken, and the evaluation precision, and check if the results being returned are consistent.)

As another demonstration, here is a trivariate special case of $F_D$:

With[{a = 1/2, b = {1/2, 1/2, 1}, c = 3/2, z = {1, 4/7, 1/3} Sin[π/5]^2},
1 + NSum[lauricellaTerm[m, a, b, c, z], {m, 1, ∞}, Compiled -> False,
Method -> {"WynnEpsilon", "ExtraTerms" -> 12},
NSumTerms -> 12, WorkingPrecision -> 30]]
1.159839217443075387343


and compare with the known closed form:

N[EllipticPi[1/3, π/5, 4/7]/Sin[π/5], 20]
1.1598392174430753873