PrecisionLauricella: package for numerical computation of Lauricella functions depending on a parameter

Abstract

We introduce the PrecisionLauricella package, a computational tool developed in Wolfram Mathematica for high-precision numerical evaluations of Lauricella functions with indices linearly dependent on a parameter, $\varepsilon$. The package leverages a method based on analytical continuation via Frobenius generalized power series, providing an efficient and accurate alternative to conventional approaches relying on multi-dimensional series expansions or Mellin–Barnes representations. This one-dimensional approach is particularly advantageous for high-precision calculations and facilitates further optimization through $\varepsilon$-dependent reconstruction from evaluations at specific numerical values, enabling efficient parallelization. The underlying mathematical framework for this method has been detailed in our previous work, while the current paper focuses on the design, implementation, and practical applications of the PrecisionLauricella package.

PROGRAM SUMMARY

Program Title: PrecisionLauricella
Developer’s repository link: https://bitbucket.org/BezuglovMaxim/precisionlauricella-package/src/main/
Licensing provisions(please choose one): GPLv3
Programming language: Wolfram Mathematica
Supplementary material: PrecisionLauricella_Examples.nb
Nature of problem(approx. 50-250 words): Lauricella functions, generalizations of hypergeometric functions, appearing in physics and mathematics, including Feynman integrals and string theory. When their indices depend linearly on a parameter $\varepsilon$, their numerical evaluation becomes challenging due to the complexity of high-dimensional series and singularities. Traditional methods, like hypergeometric re-expansion or Mellin–Barnes integrals, often lack efficiency and precision.

Managing multi-dimensional sums exacerbates computational costs, especially for high-precision requirements, making these approaches unsuitable for many practical applications. Thus, there is a pressing need for efficient, scalable methods capable of maintaining numerical accuracy and effectively handling parameter dependencies.
Solution method(approx. 50-250 words): Our method uses the Frobenius approach to achieve analytical continuations of Lauricella functions through generalized power series. Representing the functions as one-dimensional series simplifies high-precision numerical evaluation compared to traditional methods relying on multi-dimensional expansions or Mellin–Barnes integrals.

We further optimize calculations by reconstructing $\varepsilon$ dependencies from evaluations at specific values, enabling efficient parallelization and reducing computational costs.

A comprehensive mathematical exposition of the method is provided in our previous work [1].

1 Introduction

Hypergeometric functions, both in one and several variables, play a fundamental role in physics and mathematics, particularly in quantum field theory and the computation of Feynman integrals. Numerous approaches have been devised to express Feynman integrals through hypergeometric functions, including the Mellin–Barnes method, the DRA method [4, 5, 6], based on dimensional recurrence relation and analytical properties, the method of functional equations [7], the exact Frobenius method [8, 9, 10], and the Gelfand–Kapranov–Zelevinsky (GKZ) approach to Feynman integrals [13, 14, 15, 16, 17, 18, 19, 20, 21, 22]. These hypergeometric functions often depend on parameters linearly related to the dimensional regularization parameter $\varepsilon$, and their Laurent expansions in $\varepsilon$ are essential for practical applications. Several computational tools have been developed for these expansions, both numerically and analytically [32, 33, 34, 35, 36, 37, 38, 39].

This paper introduces PrecisionLauricella, a software package written in Wolfram Mathematica, designed for the high-precision numerical evaluation of Lauricella functions with indices linearly dependent on the parameter $\varepsilon$. The mathematical foundations of the algorithm implemented in this package have been detailed in our previous work [40]. The package enables computations of these functions’ Laurent expansions about $\varepsilon=0$ for arbitrary argument values. The primary challenge lies in the analytical continuation of the defining series representations of Lauricella functions from their convergence domains to the entire $\mathbb{C}^{n}$ space, where $n$ is the number of arguments. While classical integral representations such as those of Euler [41] and Mellin–Barnes [41, 42] can, in principle, address this continuation, they are less practical for high-precision numerical evaluations. Instead, series representations that solve the corresponding partial differential equations and are valid in subdomains covering the entire $\mathbb{C}^{n}$ space are more effective. These representations not only facilitate accurate numerical computations, but also align with the classical approaches to analytic continuation developed by Kummer and Riemann [43, 44], which emphasize the monodromy group.

Existing methods for the analytic continuation of Lauricella functions can be categorized into two main approaches: re-expansion of hypergeometric series using known continuations of simpler functions [42, 45, 46, 41, 47, 48] and Mellin–Barnes integral representations [49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59].

In this work, we address the analytic continuation problem by leveraging the Frobenius method to construct generalized power series solutions of the Pfaffian systems satisfied by Lauricella functions along chosen analytic-continuation paths. This approach produces one-dimensional series, which are computationally more efficient than the multi-dimensional series arising in other methods, such as Mellin–Barnes representations. PrecisionLauricella automates the entire process, offering an accessible and efficient tool for researchers dealing with Lauricella functions and their applications.

The remainder of the paper is organized as follows. In Section 2, we briefly recall the definition of Lauricella functions and associated systems of differential equations. Section 3, we provide a brief description of all the steps of the algorithm and present the program’s flowchart. Section 4 contains details on the usage of the PrecisionLauricella package and its performance. Finally, in Section 5, we present our conclusion.

2 Lauricella Functions and associated differential equations

In this section, we discuss the class of functions that are implemented in the presented package. Specifically, we focus on three Lauricella functions: $F_{A}^{(n)}$, $F_{B}^{(n)}$, and $F_{D}^{(n)}$ for $n\leq 3$. Generally, these functions are defined as hypergeometric series of several variables [60, 61]:

$$F_{A}^{(n)}(\alpha;\beta_{1},\dots,\beta_{n};\gamma_{1},\dots,\gamma_{n};x_{1},\dots,x_{n})$$

$$= \sum\limits_{m_{1},\dots,m_{n}=0}^{\infty}\frac{(\alpha)_{m_{1}+\dots+m_{n}}(\beta_{1})_{m_{1}}\dots(\beta_{n})_{m_{n}}}{(\gamma_{1})_{m_{1}}\dots(\gamma_{n})_{m_{n}}\,m_{1}!\dots m_{n}!}\,x_{1}^{m_{1}}\dots x_{n}^{m_{n}}\,,$$

$$F_{B}^{(n)}(\alpha_{1},\dots,\alpha_{n};\beta_{1},\dots,\beta_{n};\gamma;x_{1},\dots,x_{n})$$

$$= \sum\limits_{m_{1},\dots,m_{n}=0}^{\infty}\frac{(\alpha_{1})_{m_{1}}\dots(\alpha_{n})_{m_{n}}(\beta_{1})_{m_{1}}\dots(\beta_{n})_{m_{n}}}{(\gamma)_{m_{1}+\dots+m_{n}}\,m_{1}!\dots m_{n}!}\,x_{1}^{m_{1}}\dots x_{n}^{m_{n}}\,,$$

$$F_{D}^{(n)}(\alpha;\beta_{1},\dots,\beta_{n};\gamma;x_{1},\dots,x_{n})$$ (1)

$$= \sum\limits_{m_{1},\dots,m_{n}=0}^{\infty}\frac{(\alpha)_{m_{1}+\dots+m_{n}}(\beta_{1})_{m_{1}}\dots(\beta_{n})_{m_{n}}}{(\gamma)_{m_{1}+\dots+m_{n}}\,m_{1}!\dots m_{n}!}\,x_{1}^{m_{1}}\dots x_{n}^{m_{n}}\,.$$

All these series converge absolutely within specific regions. For the functions $F_{B}^{(n)}$ and $F_{D}^{(n)}$, the region of convergence is given by $|x_{1}|<1$, $|x_{2}|<1$, …, $|x_{n}|<1$. For the function $F_{A}^{(n)}$, the region is given by $|x_{1}|+|x_{2}|+\dots+|x_{n}|<1$.

When reduced to a single variable, these series correspond to the well-known hypergeometric function:

$$F_{A}^{(1)}=F_{B}^{(1)}=F_{D}^{(1)}={}_{2}F_{1}\,.$$ (2)

Calculating these functions within their convergence regions is straightforward and does not pose significant challenges. However, our primary interest lies in their analytic continuation. To achieve this, we utilize partial differential equations associated with these functions. These partial differential equations can be derived directly from the series definitions.

Consider a general series of the form

$$f=\sum A_{i_{1},i_{2},\dots,i_{n}}x_{1}^{i_{1}}x_{2}^{i_{2}}\dots x_{n}^{i_{n}}\,,$$ (3)

where

$$\frac{A_{i_{1},\dots,i_{k}+1,\dots,i_{n}}}{A_{i_{1},\dots,i_{k},\dots,i_{n}}}=\frac{P_{k}(i_{1},i_{2},\dots,i_{n})}{P^{\prime}_{k}(i_{1},i_{2},\dots,i_{n})}\,,$$ (4)

and $P_{k}$, $P^{\prime}_{k}$ are polynomials in their variables. The function $f$ then satisfies the following system of partial differential equations:

$$\left[P^{\prime}_{k}\left(\theta_{x_{1}},\theta_{x_{2}},\dots,\theta_{x_{n}}\right)\frac{1}{x_{k}}-P_{k}\left(\theta_{x_{1}},\theta_{x_{2}},\dots,\theta_{x_{n}}\right)\right]f=0\,,\qquad k=1,\dots,n\,,$$ (5)

where $\theta_{a}=a\partial/\partial a$.

For example, for the function $F_{B}^{(2)}(\alpha_{1},\alpha_{2};\beta_{1},\beta_{2};\gamma;x,y)=F_{3}$, the partial differential equations are:

$$\small\left[x(1-x)\frac{\partial^{2}}{\partial x^{2}}+y\frac{\partial^{2}}{\partial x\partial y}+\left(\gamma-(\alpha_{1}+\beta_{1}+1)x\right)\frac{\partial}{\partial x}-\alpha_{1}\beta_{1}y\frac{\partial}{\partial y}-\alpha\beta_{1}\right]F_{3} = 0\,,$$

$$\left[y(1-y)\frac{\partial^{2}}{\partial y^{2}}+x\frac{\partial^{2}}{\partial x\partial y}+\left(\gamma-(\alpha_{2}+\beta_{2}+1)y\right)\frac{\partial}{\partial y}-\alpha_{2}\beta_{2}x\frac{\partial}{\partial x}-\alpha\beta_{2}\right]F_{3} = 0\,.$$

While obtaining such partial differential systems is straightforward, applying the Frobenius method to them can be challenging. To facilitate this, we require a complete system of equations in Pfaffian form:

$$dJ=\left(M_{x}\,dx+M_{y}\,dy\right)J\,.$$ (7)

If we choose a basis as

$$J=\left\{F_{3},~{}\theta_{x}F_{3},~{}\theta_{y}F_{3},~{}\theta_{x}\theta_{y}F_{3}\right\}\,,$$ (8)

then the matrices will have the form

$$\scriptsize M_{x} = \left(\begin{array}[]{cccc}0&\frac{1}{x}&0&0\\ -\frac{\alpha_{1}\beta_{1}}{x-1}&-\frac{-\gamma+\alpha_{1}x+\beta_{1}x+1}{(x-1)x}&0&\frac{1}{(x-1)x}\\ 0&0&0&\frac{1}{x}\\ 0&-\frac{\alpha_{2}\beta_{2}y}{x(xy-x-y)}&-\frac{\alpha_{1}\beta_{1}(y-1)}{xy-x-y}&-\frac{-\alpha_{1}x-\beta_{1}x+\alpha_{1}xy+\beta_{1}xy+\alpha_{2}y+\beta_{2}y-\gamma y+y}{x(xy-x-y)}\\ \end{array}\right)\,,$$ (13)

$$M_{y} = \left(\begin{array}[]{cccc}0&0&\frac{1}{y}&0\\ 0&0&0&\frac{1}{y}\\ -\frac{\alpha_{2}\beta_{2}}{y-1}&0&-\frac{-\gamma+\alpha_{2}y+\beta_{2}y+1}{(y-1)y}&\frac{1}{(y-1)y}\\ 0&-\frac{\alpha_{2}\beta_{2}(x-1)}{xy-x-y}&-\frac{\alpha_{1}\beta_{1}x}{y(xy-x-y)}&-\frac{\alpha_{1}x+\beta_{1}x-\gamma x+\alpha_{2}xy+\beta_{2}xy+x-\alpha_{2}y-\beta_{2}y}{y(xy-x-y)}\\ \end{array}\right)\,.$$ (18)

In our work, for the functions under consideration, we utilize pre-derived Pfaffian systems from Ref. [39] with the following function bases:

$$\left\{\theta_{x_{j_{1}}}\dots\theta_{x_{j_{k}}}F_{N}\,\Big{|}\,0\leq k\leq n,\,j_{1}<j_{2}<\dots<j_{k}\right\}\,,\qquad N=A,B\,,$$

$$\left\{F_{D},\,\theta_{x_{j}}F_{D}\,\Big{|}\,j=1,\dots,n\right\}\,.$$ (20)

Our aim is to evaluate the value of a chosen Lauricella function at the point $\vec{x}$ given its initial value at the point $\vec{x}_{0}$. This can be most easily done by considering the above differential equation system along some path connecting these two points. Suppose we want to get a solution along the path $f$ parameterized by the parameter $t$, such that

$$x_{1}=f_{1}(t)\,,\quad x_{2}=f_{2}(t)\,,\quad\dots\,,\quad x_{n}=f_{n}(t)\,.$$ (21)

In this way, the differential equation system in Pfaffian form restricted to this path takes the form

$$\frac{dJ}{dt}=M_{t}J\,,\qquad M_{t}=\frac{\partial x_{1}}{\partial t}M_{1}+\frac{\partial x_{2}}{\partial t}M_{2}+\dots+\frac{\partial x_{n}}{\partial t}M_{n}\,.$$ (22)

In our work, we choose the contour $f$ as a line passing through the center of coordinates,

$$x_{1}=\kappa_{1}t\,,\quad x_{2}=\kappa_{2}t\,,\quad\dots\,,\quad x_{n}=\kappa_{n}t\,.$$ (23)

3 Algorithm overview

Figure 1: Flowchart for the PrecisionLauricella package.

In this section, we provide a general overview of the algorithm and its main stages. The flowchart of the algorithm is presented in Fig. 1. The algorithm takes as input the target function, the required precision, and the desired number of terms in the $\varepsilon$ expansion. One of the key features of the algorithm is its treatment of $\varepsilon$ dependencies. Instead of directly expanding the system of equations as a series in $\varepsilon$, we evaluate the system for several numerical values of $\varepsilon$ at specific grid points. The grid is constructed based on the required accuracy and the desired order of expansion in $\varepsilon$. All computations are then performed independently for each grid point, simplifying the calculations and enabling efficient parallelization. In the initial step, we determine the parameters of the grid along with other essential solution parameters, such as boundary conditions.

Next, we select a path in the multivariable space, as described in the previous section. During this step, we impose the condition that all $\kappa_{i}$ values must be non-negative real numbers. As a result, this step may need to be repeated several times, albeit never more than four times. The rationale behind this condition is discussed in detail in Ref. [40]. The system of equations in Pfaffian form is preconfigured in the program for the considered Lauricella functions, which simplifies this process. This step ultimately reduces the problem to a system of differential equations in a single variable.

We then determine whether analytic continuation is required. If so, the analytic continuation is performed according to the following procedure. First, we generate circular regions within a rectangular area containing the starting and endpoint. The radii of these regions are determined by their distance from the nearest singularities. Next, we compute the intersection graph: each region is treated as a node, with special nodes representing the starting and end points. Nodes are connected if the minimum distance between regions is less than three quarters of the radius of the region. Branch cuts are also taken into account, with regions on opposite sides of a cut considered non-intersecting. The intersection graph is then used to reduce the continuation problem to a graph traversal between the starting and end nodes, which is solved using a built-in Wolfram Mathematica function. If no path is found, the grid is refined, and the process is repeated until a path is identified. Examples of this procedure are shown in Figs. 2 and 3.

Following this, we compute solutions to the systems of differential equations in each region along the analytic continuation path. These systems have the general form

$$\frac{dJ}{dt}=M(t)J\,.$$ (24)

Due to our choice of basis for the hypergeometric functions, we can always write the matrix of the differential equation in the form

$$M(t)=\frac{A_{0}}{t}+R(t)\,,$$ (25)

where $R(t)$ is a rational function which has no poles at the point $t=0$. We are looking for the fundamental solution matrix for each of these equations in the form

$$U=\sum\limits_{\lambda\in S}t^{\lambda}\sum\limits_{n=0}^{\infty}\sum_{k=0}^{m_{\lambda}}c_{n}(\lambda,k)t^{n}\log^{k}(t)\,,$$ (26)

where the set $S$ consists of the non-resonant $A_{0}$ matrix eigenvalues and $m_{\lambda}$ is the number of eigenvalues in resonance with eigenvalue $\lambda$, that is their differences are integers. If several eigenvalues are in resonance, then the set $S$ will include only the smallest of them. Accordingly, the next two steps consist of deriving difference equations for the coefficients of the Frobenius series and solving these equations. From these solutions, one can obtain a numerical expression for the fundamental solution matrix for each of the regions of analytic continuation and each $\varepsilon$ value from the lattice.

Finally, we check whether we have reached the end point in the multi-variable plane. If not, we shift the origin of the system to the current point and repeat the previous steps again. If we have reached the end point, we collect all the fundamental solution matrices together to obtain the answer for the end point. Finally, we reconstruct the Laurent series coefficients for the $\varepsilon$ expansion using Lagrange interpolation polynomials.

We can estimate the numerical error of our method as

$$\text{max}\left(h^{2n-2\lfloor k/2\rfloor},\Delta_{\text{Frob}}\right)\,,$$ (27)

where $h$ is the $\varepsilon$ lattice step size, $2n$ is the number of lattice elements, $k$ is the number of terms in the $\varepsilon$ expansion and $\Delta_{\text{Frob}}$ is the error arising from the truncation of the Frobenius generalized power series.

Figure 2: On the left, one of the possible paths with complex singular points is shown. The singular points are shown in red, the end point in green, the staring point is at the origin, and shaded circles represent convergence regions. Horizontal dotted lines indicate branch cuts of singular points. On the right, the intersection graph of the expansion regions is displayed. The nodes in the intersection graph include expansion regions at regular points. The nodes and edges highlighted in red show the found path of analytic continuation. Note also, that the intersection graph may have several connected components.

Figure 3: On the left, one of the possible paths with complex singular points is shown. The singular points are shown in red, the end point in green, the staring point is at the origin, and shaded circles represent convergence regions. Horizontal dotted lines indicate branch cuts of singular points. On the right, the intersection graph of the expansion regions is displayed. The nodes in the intersection graph include expansion regions at regular points. The nodes and edges highlighted in red show the found path of analytic continuation. Note also, that the intersection graph may have several connected components.

4 PrecisionLauricella package

The PrecisionLauricella package can be freely downloaded from the bitbucket repository https://bitbucket.org/BezuglovMaxim/precisionlauricella-package/src/main/. The entire package consists of one file PrecisionLauricella.wl and, provided the path is set correctly, is loaded with the command

⬇

<< PrecisionLauricella‘

The universal function for the numerical expansion of Lauricella functions in Laurent series in parameter $\varepsilon$ is NExpandHypergeometry. This function takes three arguments. The first is a hypergeometric function that needs to be expanded. The second is the list $\{\varepsilon,k\}$, where the first argument $\varepsilon$ is the variable with respect to which the expansion is performed, and the second is the required number of terms in the $\varepsilon$ expansion. The last argument specifies the accuracy of the result with the desired number of decimal places. For those functions that are not included by default in Wolfram Mathematica, we introduce our own notations: AppellF2, AppellF3, LauricellaFA, LauricellaFB and LauricellaFD. The arguments of these functions are the same as in their definitions. In the case of Lauricella functions, numbered indices and variables are collected into lists. More details about the package functionality can be found in the Mathematica notebook with examples.

Figure 4: Analytical continuation with option $\texttt{DeltaPrescription}\rightarrow-I$ on the left and option $\texttt{DeltaPrescription}\rightarrow+I$ on the right. The singular points are shown in red and the end point in green, the staring point is at the origin and shaded circles represent convergence regions. Horizontal dotted lines indicate branch cuts of singular points.

The NExpandHypergeometry function has the following set of additional options:

• SimpleAnalyticContinuation - can take the values True or False. The default value is False. If False, uses the more complex four-step analytic continuation scheme described in [40]. If True, uses a simpler two-step scheme in which we only require that $\kappa_{i}$ in Eq. (23) be real numbers, rather than real non-negative numbers as in the four-step scheme. This scheme may be more useful in some cases and for comparing results with an answer from polylogarithms if the latter is possible.

• UseParallelComputing - can take the values True or False. The default value is True. If True, uses all available kernels for calculations. Parallelization is implemented using the built-in function ParallelTable. Otherwise, it will simply execute sequentially.

• FrobeniusNumberTerms - can take the values "Auto" or non-negative integer. The default value is "Auto". Specifies the number of elements of the Frobenius series that the program needs to calculate. In the case of "Auto", the program determines this number automatically based on the desired accuracy.

• InternalPrecision - can take the values "Auto" or non-negative integer. The default value is "Auto". Specifies the precision with which the individual coefficients in the Frobenius expansion must be calculated. This number must be significantly greater than the required precision in order to obtain a reliable final answer. In the case of "Auto", the program determines this number automatically based on the desired accuracy.

• DeltaPrescription - can take the values $\pm I$. The default value is $-I$. Determines on which side of the cut the end point is taken if its value pertains to the cut. The analytical continuation procedure in these two cases is different. The corresponding example is shown in Figure 4.

Figure 5: Average time in seconds required to expand $F_{1}\left(\frac{1}{2};1,\varepsilon;\frac{3}{2};\frac{4}{3},\frac{7}{4}\right)$ using 16 parallel kernels on the left and 8 parallel kernels on the right. The horizontal axis shows the number of terms in $\varepsilon$, and the vertical axis shows the evaluation time in seconds. Solid lines represent the best linear fit.

Figure 6: Average time in seconds required to expand $F_{2}\left(1,\frac{2}{3}\varepsilon,1,1+\frac{3}{2}\varepsilon,1-\frac{15}{7}\varepsilon;\frac{3}{2},4\right)$ on the left and $F_{D}^{(3)}(\frac{1}{2}-\varepsilon;1,\varepsilon,\varepsilon;1+2\varepsilon;\frac{4}{3},\frac{3}{4},\frac{8}{5})$ on the right using 8 parallel kernels. The horizontal axis shows the number of terms in $\varepsilon$, and the vertical axis shows the evaluation time in seconds. Solid lines represent the best linear fit.

The computation time for certain functions is also of particular interest. Examples of these test calculations are provided in Figs. 5 and 6. The test calculations were performed on a standard laptop with a 12th Gen Intel® Core™ i7-12700H × 20 processor. The most obvious consequence is that, as the number of terms in the $\varepsilon$ expansion increases, the program’s execution time grows linearly. This is a direct result of our approach to $\varepsilon$ dependencies.

5 Conclusion

In this work, we presented PrecisionLauricella, a software package written in Wolfram Mathematica, for the high-precision numerical evaluation of Lauricella functions with indices linearly dependent on a parameter $\varepsilon$. The package automates the computation of Laurent series expansions of these functions about $\varepsilon=0$, leveraging their analytical continuation in terms of Frobenius generalized power series. Unlike multi-dimensional series arising in other methods, these one-dimensional series offer significant advantages in terms of computational efficiency and precision, making them especially suitable for large-scale calculations.

Our approach includes a dedicated treatment of the $\varepsilon$ dependences in these expansions, which not only accelerates computations but also facilitates parallelization, making the package scalable for high-performance applications. We validated the accuracy and efficiency of PrecisionLauricella by comparing its results, where applicable, with existing tools for hypergeometric functions, such as those presented in Refs. [62, 48] and Refs. [38, 39]. These comparisons demonstrated the robustness and reliability of the presented approach.

We hope that PrecisionLauricella will serve as a reliable and efficient tool for researchers working with Lauricella and similar hypergeometric functions. We also look forward to exploring further extensions of the package to other classes of hypergeometric functions and additional applications in both mathematics and physics.