# Rigorous Multiple-Precision Evaluation of D-Finite Functions in SageMath

Marc Mezzarobba

► To cite this version:

Marc Mezzarobba. Rigorous Multiple-Precision Evaluation of D-Finite Functions in SageMath. 5th International Congress on Mathematical Software (ICMS 2016), Jul 2016, Berlin, Germany. http://hal.science/hal-01342769

HAL Id: hal-01342769

https://hal.science/hal-01342769

Submitted on 6 Jul 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# Rigorous Multiple-Precision Evaluation of D-Finite Functions in SageMath

Marc Mezzarobba*

CNRS, LIP6, Université Pierre et Marie Curie**, Paris, France

marc@mezzarobba.net,

http://marc.mezzarobba.net/

**Abstract.** We present a new open source implementation in the Sage-Math computer algebra system of algorithms for the numerical solution of linear ODEs with polynomial coefficients. Our code supports regular singular connection problems and provides rigorous error bounds.

## 1 Introduction

Many special functions satisfy differential equations

$$p_r(x)f^{(r)}(x) + \dots + p_1(x)f'(x) + p_0(x)f(x) = 0 \quad (1)$$

whose coefficients $p_i$ depend polynomially on the variable $x$. In virtually all cases, such special functions can be defined by complementing (1) either with simple initial values $f(0), f'(0), \dots$ or with constraints on the asymptotic behavior of $f(x)$ as $x$ approaches a singular point. For example, the error function satisfies

$$\operatorname{erf}''(x) + 2x \operatorname{erf}'(x) = 0, \quad \operatorname{erf}(0) = 0, \quad \operatorname{erf}'(0) = \frac{2}{\sqrt{\pi}}, \quad (2)$$

while the modified Bessel function $K_0$ is the solution of

$$xK_0''(x) + K_0'(x) - xy(x) = 0 \quad \text{s.t.} \quad K_0(x) = -\log(x/2) - \gamma + O_{x \to 0}(x). \quad (3)$$

This observation has led to the idea of developing algorithms that deal with these functions in a uniform way, using the ODE (1) as a data structure [11,2].

In this context, solutions of linear ODEs with polynomial coefficients are called *D-finite* (or *holonomic*) functions. These names originate from combinatorics, where D-finite power series arise naturally as generating functions [17,4]. While classical special functions typically satisfy ODEs of order 2 to 4 with

<footnote>* Supported in part by ANR grant ANR-14-CE25-0018-01 (FastRelax).</footnote>

<footnote>** Sorbonne Universités, UPMC Univ Paris 06, CNRS, LIP6 UMR 7606, 4 place Jussieu 75005 Paris.</footnote>

<footnote>This article is in the public domain. In jurisdictions where this is not possible, any entity is granted the perpetual right to use this work for any purpose, without any conditions other than those required by law.</footnote>

<footnote>Extended abstract for a talk given at the 5th International Congress on Mathematical Software (ICMS 2016). Accepted for publication in the proceedings, but withdrawn due to a disagreement with Springer about the above public domain notice.</footnote>
simple polynomial coefficients, D-finite generating functions are usually given
by more general and sometimes very large equations. Other notable sources of
general D-finite functions include mathematical physics and number theory. Ex-
plicit examples in theoretical works on meromorphic ODEs also tend to have
polynomial coefficients.

Packages for manipulating D-finite power series symbolically exist for major
computer algebra systems, and include in particular gfun [16] and the rest of the
algolib bundle for Maple, GeneratingFunctions [12] and the RISCErgoSum suite
for Mathematica, and ore_algebra [9] for SageMath [18].

We extended ore_algebra with features dedicated to D-finite analytic func-
tions such as special functions, starting with numerical evaluation algorithms.
The need for such features is most evident in the context of special functions,
but they turn out to be useful in a wide range of applications, as we will see from
the examples in the sequel. We refer to our implementation as oa.analytic, af-
ter the name of the subpackage of ore_algebra where most of the code resides.
The present paper introduces the main features of oa.analytic, with example
applications. The reader is invited to consult the package's documentation for
more information. We further refer to Hille [5] for more background on complex
analytic ODEs, and to Kauers [8] for an overview of the field of D-finite func-
tions in computer algebra. Also note that oa.analytic is intended to supersede
NumGfun [14], a Maple package with similar features previously developed by
the author<fnref n="1" />.

Since the main feature of oa.analytic is to compute numerical values of
solutions of differential equations, it is, at its core, a numerical ODE solver.
However, it is limited to linear ODEs with polynomial coefficients, and differs
from typical ODE solvers in a number of other respects. First of all, the results
it computes come with *rigorous error bounds*, making it a “validated” ODE
solver. It works in *arbitrary precision*, and implements specific algorithms to
make evaluations at very high precisions feasible. Last but not least, it offers
extensive support for connection problems between *regular singular points of*
differential operators. These distinctive features all rely on *symbolic-numeric
algorithms* that make use of specific properties of D-finite functions.

The development of oa.analytic is still very much a work in progress. At
the time of writing, there is no formal release of ore_algebra containing the
analytic subpackage yet. The latest development version of oa.analytic can
be found at http://marc.mezzarobba.net/code/ore_algebra-analytic. It is
distributed under the GNU General Public Licence, version 2 or later. Comments,
bug reports and feature requests are welcome.

## 2 Differential Operators in ore_algebra

In ore_algebra, differential operators are represented algebraically as *Ore polynomials*. More specifically, the operators we are interested in are elements of the

<footnote n="1">¹ Many of the examples in this article are adapted from earlier work on NumGfun [15].</footnote>
ring, denoted $\mathbb{K}[x] \langle D_x \rangle$, of polynomials over $\mathbb{K}[x]$ for some number field $\mathbb{K} \subset \mathbb{C}$ in a variable $D_x$ subject to the commutation rule $D_x x = x D_x + 1$ (which expresses the equality $(x f(x))' = x f'(x) + f(x)$). Elements of $\mathbb{K}[x] \langle D_x \rangle$ act on differentiable functions in a natural way, namely by $x \cdot f = (\xi \mapsto \xi f(\xi))$ and $D_x \cdot f = f'$, and can thus be identified with differential operators.

The first thing to do in order to use `ore_algebra` in a Sage session is to install it in the Python search path and `import` it. The rings `K[x]` and `K[x]⟨D_x⟩` and their respective generators `x` and `D_x` can then be created using standard Sage syntax as follows:

```python
sage: from ore_algebra import OreAlgebra
sage: Pols.<x> = PolynomialRing(QQ); Dops.<Dx> = OreAlgebra(Pols)
```

Note that `ore_algebra` supports general Ore polynomials, but basic arithmetic with differential operators in one variable is all we are going to need here.

After these commands, the Python variable `Dx` contains an object of type "Ore polynomial" that satisfies

```python
sage: Dx*x
x*Dx + 1
```

Most features of `ore_algebra`, including those presented in this article, are available as methods of Ore polynomial objects. Some functions of the `analytic` subpackage can also be called directly in advanced usage.

## 3 Numerical Analytic Continuation

Let $L = p_r(x)D_x^r + \dots + p_0(x) \in \mathbb{K}[x] \langle D_x \rangle$ be a differential operator, and assume that the origin is an ordinary point of $L$, i.e. that $p_r(0) \neq 0$. By the classical Cauchy existence theorem [5], the equation $L \cdot f = 0$ then admits an $r$-dimensional space of complex analytic solutions defined in a neighborhood of 0 and characterized by $r$ initial values at 0. The `numericalSolution()` method of differential operators computes values of D-finite functions defined in this way. Perhaps the simplest example is the computation of $e = \exp(1)$ as the value at 1 of the solution of $f' - f = 0$ such that $f(0) = 1$:

```python
sage: (Dx-1).numericalSolution(ini=[1], path=[0,1], eps=1e-40)
[2.71828182845904523536028747135266249775725 +/- 7.53e-42]
```

The last argument, `eps`, provides an indication of the accuracy the user is interested in. The method outputs a "mid-rad interval", or *ball*, consisting of a multiple-precision center and a low-precision radius. In normal circumstances, the diameter of this *ball* will be of the order of magnitude of `eps`, but this is not guaranteed to be the case. In contrast, barring any bug, the mathematical result is guaranteed to lie within the bounds of the *ball*.

Evaluations on the complex plane are also supported. Based on (2), we may for instance evaluate the error function as follows:

```python
sage: (Dx^2 + 2*x*Dx).numericalSolution([0, 2/sqrt(pi)], [0, i])
[+/- 2.91e-25] + [1.6504257587975429 +/- 2.70e-17]*I
```
Allowed evaluation points and initial values include rational and floating-point numbers, algebraic numbers, and constant symbolic expressions.

The reader may have noticed that the argument used to specify the evaluation point is called *path*. The main reason why `numericalSolution()` does not simply take an evaluation point is that in general, solutions of differential equations with singular points are not single-valued analytic functions. What `diffop.numericalSolution(ini, path)` really computes is the analytic continuation along *path* of the solution of `diffop` defined in the neighborhood of the starting point by the initial conditions *ini*. Consider for example the operator $L = xD_x + D_x$, obtained by differentiating the inhomogeneous equation $xy'(x) = 1$. A basis of solutions of $L$ is $\{1, \log(x)\}$, corresponding respectively to the initial values $(y(1), y'(1)) = (1, 0)$ and $(y(1), y'(1)) = (0, 1)$. Analytic continuation of the latter to $x = -1$ along a path that passes above the singular point $x = 0$ corresponds to the standard determination of the complex logarithm, while a path that passes below the origin yields another branch:

```
sage: (x*Dx^2 + Dx).numericalSolution([0, 1], [1, i, -1])
[+/- 5.11e-18] + [3.1415926535897932 +/- 4.38e-17]*I
```

```
sage: (x*Dx^2 + Dx).numericalSolution([0, 1], [1, -i, -1])
[+/- 5.11e-18] + [-3.1415926535897932 +/- 4.38e-17]*I
```

More generally, it is possible to compute the matrix that maps a vector

$$ [f(x_0), f'(x_0), \frac{1}{2}f''(x_0), \dots, \frac{1}{(r-1)!}f^{(r-1)}(x_0)] \quad (4) $$

of initial conditions<fnref n="2" /> at $x_0$ to the vector $[f(x_1), f'(x_1), \dots]$ of “initial” conditions at some other point $x_1$ that define the analytic continuation of the same solution $f$ along a given path. Following van der Hoeven [19], we call this matrix the *transition matrix* along that path.

```
sage: (x*Dx^2 + Dx).numerical_transition_matrix([1, 2], 1e-10)
[ 1.000000000000000 [0.69314718056 +/- 6.99e-14]]
[ 0 [0.500000000000000 +/- 9.1e-15]]
```

**Example 1.** Transition matrices along loops that turn around a singular point exactly once are called monodromy matrices. The first step of van Enckevort and van Straten's experimental study of differential equations of Calabi-Yau type [21] consists in computing high-precision approximations of all monodromy matrices starting from a common base point. The `numerical_transition_matrix()` method makes this computation straightforward. In the case of the singularity at 1 of the operator studied in Example 1 from their article, we obtain<fnref n="3" />:

```
sage: diffop = ((x*Dx)^4 - x*(65*(x*Dx)^4+130*(x*Dx)^3 + 105*(x*Dx)^2 + 40*x*Dx + 6) + 4*x^2*(4*x*Dx+3)*(x*Dx+1))^2*(4*x*Dx+5))
```

<footnote n="2">² Note the inverse factorials in our definition of initial conditions. This convention also applies to the *ini* argument of `numericalSolution()`. It is convenient in view of the generalization to regular singular points discussed in the next section.</footnote>

<footnote n="3">³ In situations such as this one where the singular point is regular singular, it is also possible, and typically more efficient, to express the monodromy matrix using the primitives described in the next section, as suggested by van der Hoeven [20].</footnote>
```
sage: mat = diffop.numerical_transition_matrix([1/2,1-i/2,3/2,1+i/2,1/2],
1e-1000)
```

```
sage: mat[0,0] \# 990 digits omitted
[1.000...000 +/- 6.45e-998] + [0.733...6314 +/- 8.99e-998]*I
```

The analytic continuation algorithm behind these functions is based on work by Chudnovsky and Chudnovsky [3] with some later improvements [19,15]. In the language of numerical methods for ODEs, it is a high-order Taylor series method with a step size of the order of the local radius of convergence of solutions. The computation of the transition matrix for each step boils down to that of sums of power series whose coefficients are generated by linear recurrence relations. Two algorithms for this task are currently implemented: a naïve iterative summation and a simple binary splitting method. In both cases, the entries of a given column are actually computed all at once, by working with truncated power series in a manner similar to direct-mode automatic differentiation. Once the terms of the sum have decreased below the target accuracy, the algorithm of [13] is used to get a rigorous bound on the tail of the series. Both the computation of the numerical result and that of error bounds rely on Arb [7] for controlling the round-off errors in basic arithmetic operations.

# 4 Regular Singular Connection Problems

Singular connection problems are similar to the analytic continuation problems of the previous section, except that one or both ends of the analytic continuation path are singular points of the operator (points where the leading coefficient vanishes). At a singular point, the Cauchy theorem does not apply, and solutions of the operator may diverge or display a branch point.

Our implementation supports the situation where the singularities involved are regular singular points. A singular point $x_0$ of $L \in \mathbb{K}[x] \langle D_x \rangle$ is regular when $L$ admits a full basis of solutions each of the form

$$ (x - x_0)^{\alpha} \left( g_0(x) + g_1(x) \log(x - x_0) + \dots + g_p(x) \frac{(\log(x - x_0))^p}{p!} \right) \quad (5) $$

where $g_0, \dots, g_p$ are analytic functions and $\alpha$ is an algebraic number.

A basic example is that of the modified Bessel equation (3), which possesses a regular singular point at $x = 0$ and whose standard solutions are defined by their asymptotic behavior there. Starting from (3), we can compute $K_0(1)$ by:

```
sage: ini = [-1, log(2)-euler_gamma]
sage: (x*Dx^2 + Dx - x).numericalSolution(ini, [0, 1], 1e-10)
[0.42102443824 +/- 2.11e-12]
```

This command is to be understood as follows. It is classical that, when $x_0$ is a regular singular point of $L$, a solution of $L$ is characterized by the coefficients of the monomials $(x - x_0)^\nu \log(x - x_0)^k/k!$ in its expansion (5) about $x_0$, where $\nu$ ranges over the roots of the so-called indicial polynomial of $L$ at $x_0$.
The generating series $a(x) = \sum_{n=0}^{\infty} a_n x^n$ and $b(x) = \sum_{n=0}^{\infty} b_n x^n$ are solutions of the following operator:

```
sage: diffop = (x^2*(x^2-34*x+1)*Dx^4 + 5*x*(2*x^2-51*x+1)*Dx^3 
             + (25*x^2-418*x+4)*Dx^2 + (15*x-117)*Dx + 1)
```

This operator has three regular singular points: $0, \alpha = (\sqrt{2} + 1)^4 \approx 33.971$, and $\alpha^{-1} = (\sqrt{2} - 1)^4 \approx 0.029$.

```
sage: s = diffop.leadingCoefficient().roots(AA, multiplicities=False); s
[0, 0.02943725152285942?, 33.97056274847714?]
```

This implies that the radius of convergence of $a(x)$ is in $\{\alpha, \alpha^{-1}, \infty\}$, so that $a_n$ must be one of the forms $\alpha^{-n+o(n)}$, $\alpha^{n+o(n)}$ and $e^{o(n)}$ as $n \to \infty$. The same holds for $b_n$. A key step of Apéry's proof is to show that $a_n = \alpha^{n+o(n)}$ and $b_n = \alpha^{n+o(n)}$, but $b_n - \zeta(3)a_n = \alpha^{-n+o(n)}$.

Analytic solutions at $x=0$ are determined by the first two coefficients of their Taylor expansions. At $x = \alpha^{-1}$, the canonical local basis takes the form $(f_0, f_1, f_2, f_3)$ where $f_0, f_2, f_3$ are analytic and $f_1(x) \sim \sqrt{x} - \alpha^{-1}$.

```
sage: diffop.local_basis_monomials(0)
[1/2*log(x)^2, log(x), 1, x]
sage: diffop.local_basis_monomials(s[1])
[1, sqrt(x - 0.02943725152285942?), x - 0.02943725152285942?,
(x - 0.02943725152285942?)^2]
```

To prove that $a_n = \alpha^{n+o(n)}$, it is enough to check that the coefficient $c_1$ in the decomposition $a(x) = \sum c_i f_i(x)$ is nonzero.

```
sage: mat = diffop.numerical_transition_matrix([0, sing[1]], 1e-40)
```

From the initial values $a_0 = 1, a_1 = 5$, we see that $c_1$ is equal to

```
sage: mat[1,2] + 5*mat[1,3]
[4.546376247522844600239593024915161553303 +/- 9.85e-41]*I
```

With a little more work, using the method of singularity analysis [4], this is actually enough to prove that $a_n = n^{-3/2}\alpha^n(\lambda + O(1/n))$ where $\lambda$ is given by

```
sage: CBF(-i/(2*sqrt(pi)))*sqrt(sing[1])*\(mat[1,2] + 5*mat[1,3])
[0.220043767112643 +/- 2.06e-16]
```

This result agrees with the closed form $a_n \sim \sqrt{\alpha}(2\pi\sqrt{2})^{-3/2}n^{-3/2}\alpha^n$ from [6].

In the same way, using the initial terms $b_0 = 0, b_1 = 6$ of $(b_n)$, we can compute the coefficient of $f_1$ in the decomposition of $b(x)$ and check that the ratio of the two values we obtained is $\zeta(3)$.

```
sage: 6*mat[1,3]/\(mat[1,2] + 5*mat[1,3])
[1.2020569031595942853997381615114499907650 +/- 6.08e-41]
```

**Acknowledgements.** This work was prompted by discussions with Fredrik Johansson. I am also indebted to the Sage developers, in particular Vincent Delecroix, Jeroen Demeyer and Clemens Heuberger, who reviewed or helped writing a number of Sage patches vital to this project, and to Christoph Lauter for comments on a version of this text.
## References

1. R. Apéry. Irrationalité de $ζ(2)$ et $ζ(3)$. *Astérisque*, 61, page 11–13, 1979.

2. A. Benoit, F. Chyzak, A. Darrasse, S. Gerhold, M. Mezzarobba, and B. Salvy. The dynamic dictionary of mathematical functions (DDMF). In K. Fukuda et al., editors, *Mathematical Software — ICMS 2010*, volume 6327 of *Lecture Notes in Computer Science*, page 35–41. Springer, 2010.

3. D. V. Chudnovsky and G. V. Chudnovsky. Computer algebra in the service of mathematical physics and number theory. In D. V. Chudnovsky and R. D. Jenks, editors, *Computers in Mathematics*, volume 125 of *Lecture Notes in Pure and Applied Mathematics*, page 109–232, Stanford University, 1986. Dekker.

4. P. Flajolet and R. Sedgewick. *Analytic Combinatorics*. Cambridge University Press, 2009.

5. E. Hille. *Ordinary differential equations in the complex domain*. Wiley, 1976. Dover reprint, 1997.

6. M. D. Hirschhorn. Estimating the Apéry numbers. *The Fibonacci Quarterly*, 50(2):129–131, 2012.

7. F. Johansson. Arb: a C library for ball arithmetic. *ACM Communications in Computer Algebra*, 47(4):166–169, 2013. http://fredrikj.net/arb/

8. M. Kauers. The holonomic toolkit. *Computer Algebra in Quantum Field Theory: Integration, Summation and Special Functions, Texts and Monographs in Symbolic Computation*, 2013.

9. M. Kauers, M. Jaroschek, and F. Johansson. Ore polynomials in Sage. In J. Gutierrez, J. Schicho, and M. Weimann, editors, *Computer Algebra and Polynomials*, page 105–125. Springer, 2015.

10. C. Koutschan. Lattice Green's functions of the higher-dimensional face-centered cubic lattices. *Journal of Physics A: Mathematical and Theoretical*, 46(12):125005, 2013.

11. C. Lanczos. *Applied analysis*. Prentice-Hall, 1956.

12. C. Mallinger. Algorithmic manipulation and transformations of univariate holonomic functions and sequences. Diplomarbeit, RISC-Linz, 1996.

13. M. Mezzarobba. Truncation bounds for D-finite series. In preparation.

14. M. Mezzarobba. NumGfun: a package for numerical and analytic computation with D-finite functions. In S. M. Watt, editor, *ISSAC '10*, page 139–146. ACM, 2010.

15. M. Mezzarobba. *Autour de l'évaluation numérique des fonctions D-finies*. Thèse de doctorat, École polytechnique, 2011.

16. B. Salvy and P. Zimmermann. Gfun: A Maple package for the manipulation of generating and holonomic functions in one variable. *ACM Transactions on Mathematical Software*, 20(2):163–177, 1994.

17. R. P. Stanley. *Enumerative combinatorics*, volume 2. Cambridge University Press, 1999.

18. The SageMath Developers. SageMath mathematics software, 2005–. http://www.sagemath.org/

19. J. van der Hoeven. Fast evaluation of holonomic functions. *Theoretical Computer Science*, 210(1):199–216, 1999.

20. J. van der Hoeven. Fast evaluation of holonomic functions near and in regular singularities. *Journal of Symbolic Computation*, 31(6):717–743, 2001.

21. C. van Enckevort and D. van Straten. Monodromy calculations of fourth order equations of Calabi-Yau type. In J. D. Lewis, S.-T. Yau, and N. Yui, editors, *Mirror Symmetry V*, volume 38 of *AMS/IP Studies in Advanced Mathematics*, page 539–559. International Press, 2006.