Computing the volume of compact semi-algebraic sets

# Computing the volume of compact semi-algebraic sets

Pierre Lairez InriaFrance Marc Mezzarobba Sorbonne Université, CNRS,
Laboratoire d’Informatique de Paris 6, LIP6, Équipe PeQuaN
F-75252, Paris Cedex 05France
and  Mohab Safey El Din Sorbonne Université, CNRS, Inria,
Laboratoire d’Informatique de Paris 6, LIP6, Équipe PolSys
F-75252, Paris Cedex 05France
July 31, 20192019
###### Abstract.

Let be a compact basic semi-algebraic set defined as the real solution set of multivariate polynomial inequalities with rational coefficients. We design an algorithm which takes as input a polynomial system defining and an integer  and returns the -dimensional volume of at absolute precision .

Our algorithm relies on the relationship between volumes of semi-algebraic sets and periods of rational integrals. It makes use of algorithms computing the Picard-Fuchs differential equation of appropriate periods, properties of critical points, and high-precision numerical integration of differential equations.

The algorithm runs in essentially linear time with respect to . This improves upon the previous exponential bounds obtained by Monte-Carlo or moment-based methods. Assuming a conjecture of Dimca, the arithmetic cost of the algebraic subroutines for computing Picard-Fuchs equations and critical points is singly exponential in and polynomial in the maximum degree of the input.

Semi-algebraic sets; Picard-Fuchs equations; Symbolic-numeric algorithms
Marc Mezzarobba is supported in part by ANR grant ANR-14-CE25-0018-01 FastRelax. Mohab Safey El Din is supported by the ANR grants ANR-17-CE40-0009 Galop, ANR-18-CE33-0011 Sesame, the PGMO grant Gamma and the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement N° 813211 (POEMA).
copyright: rightsretainedjournalyear: 2019copyright: licensedothergovconference: International Symposium on Symbolic and Algebraic Computation; July 15–18, 2019; Beijing, Chinabooktitle: International Symposium on Symbolic and Algebraic Computation (ISSAC ’19), July 15–18, 2019, Beijing, Chinaprice: 15.00doi: 10.1145/3326229.3326262isbn: 978-1-4503-6084-5/19/07

## 1. Introduction

Semi-algebraic sets are the subsets of which are finite unions of real solution sets to polynomial systems of equations and inequalities with coefficients in . Starting from Tarski’s algorithm for quantifier elimination (Tarski, 1998) improved by Collins through the Cylindrical Algebraic Decomposition algorithm (Collins, 1975), effective real algebraic geometry yields numerous algorithmic innovations and asymptotically faster routines for problems like deciding the emptiness of semi-algebraic sets, answering connectivity queries or computing Betti numbers (e.g., Basu et al., 2006; Safey El Din and Schost, 2003; Canny, 1988; Safey El Din and Schost, 2017; Bürgisser et al., 2019). The output of all these algorithms is algebraic in nature. In this paper, we study the problem of computing the volume of a (basic) compact semi-algebraic set defined over . The output may be transcendental: for instance, the area of the unit circle in is .

Volumes of semi-algebraic sets actually lie in a special class of real numbers, for they are closely related to Kontsevich-Zagier periods introduced in (Kontsevich and Zagier, 2001). A (real) period is the value of an absolutely convergent integral of a rational function with rational coefficients over a semi-algebraic set defined by polynomials with rational coefficients. For example, algebraic numbers are periods, as are , , . Since , volumes of semi-algebraic sets defined over  are periods. Conversely, interpreting an integral as a “volume under a graph” shows that periods are differences of volumes of semi-algebraic sets defined over . In (Viu-Sos, 9 03), it is further shown that periods are differences of volumes of compact semi-algebraic sets defined over .

The problem we consider in this paper is thus a basic instance of the more general problem of integrating an algebraic function over a semi-algebraic set; it finds applications in numerous areas of engineering sciences. Performing these computations at high precision (hundreds to thousands of digits) is also relevant in experimental mathematics, especially for discovering formulas, as explained, for example, in (Bailey and Borwein, 2011). Most of the examples featured in this reference are periods, sometimes in disguise.

#### Prior work.

The simplest semi-algebraic sets one can consider are polytopes. The computation of their volume has been extensively studied, with a focus on the complexity with respect to the dimension. It is known that even approximating the volume of a polytope deterministically is #P-hard (Dyer and Frieze, 1988; Khachiyan, 1989). The celebrated probabilistic approximation algorithm in (Dyer et al., 1991), which applies to more general convex sets, computes an -approximation in time polynomial in the dimension of the set and . A key ingredient for this algorithm is a Monte Carlo method for efficiently sampling points from a convex set. Since then, Monte Carlo schemes have been adopted as the framework of several volume estimation algorithms.

In contrast, we deal here with compact semi-algebraic sets which can be non-convex and even non-connected. Additionally, while volumes of polytopes are rational, the arithmetic nature of volumes of semi-algebraic sets is much less clear, as unclear as the nature of periods. This raises the question of the computational complexity of a volume, even taken as a single real number.

A simple Monte Carlo technique applies in our setting as well: one samples points uniformly in a box containing and estimates the probability that they lie in . This method is of practical interest at low precision but requires  samples to achieve an error bounded by with high probability. We refer to (Koiran, 1995) which deals with definable sets, a class which encompasses semi-algebraic sets.

In a different direction, numerical approximation schemes based on the moment problem and semi-definite programming have been designed in (Henrion et al., 2009). They are also of practical interest at low precision, and can provide rigorous error bounds, but the convergence is worse than exponential with respect to  (Korda and Henrion, 2018).

Another line of research, going back to the nineteenth century, is concerned with the computation of periods of algebraic varieties. In particular, we build on work (Chudnovsky and Chudnovsky, 1990) on the high-precision numerical solution of ODEs with polynomial coefficients which was motivated, among other things, by applications to periods of Abelian integrals (see Chudnovsky and Chudnovsky, 1990, p. 133).

#### Main result.

We describe a new strategy for computing volumes of semi-algebraic sets, at the crossroads of effective algebraic and real algebraic geometry, symbolic integration, and rigorous numerical computing. Our approach effectively reduces the volume computation to the setting of (Chudnovsky and Chudnovsky, 1990). It yields an algorithm that approximates the volume of a fixed, bounded basic semi-algebraic set in almost linear time with respect to the precision. More precisely, we prove the following bit complexity estimate.

###### Theorem 1 ().

Let be polynomials in , and let be the semi-algebraic set defined by . Assume that is compact. There exists an algorithm which computes, on input  and , an approximation  of the volume  of with . When are fixed, the algorithm runs in time (for any ) as .

The algorithm recursively computes integrals of volumes of sections of . Let  denote the -dimensional volume of , for some nonzero linear projection . In our setting,  is a piecewise analytic function and, except at finitely many , is solution of a linear differential equation with polynomial coefficients known as a Picard-Fuchs equation.

The problem points belong to the critical locus of the restriction of the projection to a certain hypersurface containing the boundary of and are found by solving appropriate polynomial systems. (Compare (Khachiyan, 1993) in the case of polytopes.) The Picard-Fuchs equation for is produced by algorithms from symbolic integration, in particular (Bostan et al., 2013; Lairez, 2016). To obtain the volume of , it then suffices to compute  with a rigorous numerical ODE solving algorithm, starting from values  at suitable points  obtained through recursive calls.

The complexity with respect to the dimension  of the ambient space and the number , maximum degree , and coefficient size of the polynomials is harder to analyze. We will see, though, that under reasonable assumptions, the “algebraic” steps (computing the critical loci and of the Picard-Fuchs equations) take at most arithmetic operations in .

#### Example.

The idea of the method is well illustrated by the example of a torus , here of major radius  and minor radius . Let

 S={(x,y,z)∈R3 ∣∣ (x2+y2+z2+3)2⩽16(x2+y2)}.

The area (2-dimensional volume) of a section  defines a function  (see Figure 1).

It is analytic, except maybe at the critical values  and  where the real locus of the curve is singular. On each interval on which  is analytic, it satisfies the Picard-Fuchs equation

 (1) (t−3)(t+3)(t−1)2(t+1)2t2v′′′(t)+(t2+9)(t−1)2(t+1)2tv′′(t)−(2t4+11t2−9)(t−1)(t+1)v′(t)+2(t2+3)t3v(t)=0,

which we compute in 2 seconds on a laptop using the algorithm of (Lairez, 2016) and Theorem 8.

We know some special values of , namely , and . Additionally, we have  as . These properties characterize the analytic function  in the -dimensional space of analytic solutions of the differential equation (1) on , and similarly for . (Our algorithm actually uses recursive calls at generic points instead of these ad hoc conditions.) The rigorous ODE solver part of the Sage package ore_algebra (Mezzarobba, 2016) determines in less than a second that

 ∫3−3v(t)dt=39.4784176043[...]25056533975±10−60.

And indeed, it is not hard to see in this case that . We can obtain 1000 digits in less than a minute.

#### Outline.

The remainder of this article is organized as follows. In Section 2, we give a high-level description of the main algorithm. As sketched above, the algorithm relies on the computation of critical points, Picard-Fuchs equations, and numerical solutions of these equations. In Section 3, we discuss the computation of Picard-Fuchs equations and critical points, relating these objects with analyticity properties of the “section volume” function. Then, in Section 4, we describe the numerical solution process and study its complexity with respect to the precision. Finally, in Section 5, we conclude the proof of Theorem 1 and state partial results on the complexity of the algorithm with respect to , , and .

#### Acknowledgements.

We would like to thank the anonymous reviewers for their careful reading and valuable comments.

## 2. Volumes of semi-algebraic sets

We start by designing an algorithm which deals with the case of a union of connected components of a semi-algebraic set defined by a single inequality. Next, we will use a deformation technique to handle semi-algebraic sets defined by several inequalities.

### 2.1. Sets defined by a single inequality

Let  and be the semi-algebraic set

 A≜{(ρ,x)∈R×Rn ∣∣ f(ρ,x)⩾0}.

Let  be the projection on the -coordinate. We want to compute the volume of a union  of connected components of  starting from the volumes of suitable fibers . For technical reasons, we first consider the slightly more general situation where  is a union of connected components of for some open interval . From a computational point of view, we assume that  is described by a semi-algebraic formula , that is,

 U={(ρ,x)∈A ∣∣ ΘU(ρ,x)},

where is a finite disjunction of conjunctions of polynomial inequalities with (in our setting) rational coefficients.

For , let  and . Let  (we will often omit the subscript ) be the set of exceptional values

 (2) Σf≜{ρ∈R∣∃x∈Rn,f(ρ,x)=0∧∀i,∂∂xif(ρ,x)=0}.

Thus, when is square-free, exceptional values are either critical values of the restriction to the hypersurface of the projection , or images of singular points of . By definition of , for any , the zero set of is a smooth submanifold of .

Further, we say that assumption (R) holds for if

 (R) {z∈Rn+1 ∣∣ f(z)=0∧∂∂tf(z)=0∧∀i,∂∂xif(z)=0}=∅.

Observe that by Sard’s theorem (e.g. Basu et al., 2006, Theorem 5.56), when (R) holds, the exceptional set  is finite.

The mainstay of the method is the next result, to be proved in §3. Let denote the set of Fuchsian linear differential operators with coefficients in  whose local exponents at singular points are rational (see §4 for reminders on Fuchsian operators and their exponents).

###### Theorem 2 ().

If is bounded and , then the function  is solution of a computable differential equation of the form , where depends only on .

We will also use the following proposition, which summarizes the results of Proposition 12 and Lemma 13 in §4. The complete definition of “good initial conditions” is given there as well. Up to technical details, this simply means a system  of linear equations of the form that suffices to characterize a particular solution  among the solutions of . An -approximation of  is made of the same equations with each right-hand side  replaced by an enclosure  of diameter .

###### Proposition 2 ().

Let have order , and let be a real interval with algebraic endpoints. Let be a solution of with a finite limit at  and be a system of good initial conditions for  on  defining .

1. Given , a precision  and a -approximation of , one can compute an interval of width  (as for fixed , , and ) containing .

2. Given , one can compute such that the form a system of good initial conditions for  on .

Assume now that is a bounded union of connected components of  (i.e., that we can take above), and that (R) holds for . The algorithm is recursive. Starting with input , and , it first computes the set of exceptional values so as to decompose into intervals over which the function  satisfies the differential equation given by Theorem 2. Since  is bounded, one has

 voln+1U=ℓ−1∑i=1voln+1(U∩pr−1(αi,αi+1))=ℓ−1∑i=1∫αi+1αiv(t)dt.

Fix  and consider the interval . Since is annihilated by , its anti-derivative vanishing at  is annihilated by the operator , which belongs to  since  does. Additionally, if is a system of good initial conditions for  that defines , then is a system of good initial conditions for defining  (see Lemma 12 in §4). Thus, by Proposition 2, to compute to absolute precision , it suffices to compute , , to precision .

By definition of , since , there is no solution to the system

 f(ρj,−)=∂∂x1f(ρj,−)=⋯=∂∂xnf(ρj,−)=0

which means that (R) holds for . Additionally, is a bounded union of connected components of . Hence, the values can be obtained by recursive calls to the algorithm with instantiated to .

The process terminates since each recursive call handles one less variable. In the base case, we are left with the problem of computing the length of a union of real intervals encoded by a semi-algebraic formula. This is classically done using basic univariate polynomial arithmetic and real root isolation (Basu et al., 2006, Chap. 10).

The complete procedure is formalized in Algorithm 1. The quantities denoted with a tilde in the pseudo-code are understood to be represented by intervals, and the operations involving them follow the semantics of interval arithmetic. Additionally, we assume that we have at our disposal the following subroutines:

• , , and , which implement the algorithms implied, respectively, by Theorem 2 and Proposition 2 (1) and (2);

• , which returns an encoding for a finite set of real algebraic numbers containing the exceptional values associated to , sorted in increasing order;

• where and is a semi-algebraic formula describing a union  of connected components of , which returns an interval of width containing .

The following result summarizes the above discussion.

###### Theorem 3 ().

Assume that  is a bounded union of connected components of  and that (R) holds. Then, on input and , Algorithm 1 (Volume1) returns a real interval of width (for fixed ) containing .

### 2.2. Sets defined by several inequalities

Now, we show how to compute the volume of a basic semi-algebraic set defined by

 f1⩾0,…,fr⩾0,fi∈Q[x1,…,xn],

assuming that is compact.

We set , and consider the semi-algebraic set defined by . Observe that the polynomial  satisfies (R) because . We can hence choose an interval with that contains no element of . Let and be the projection on the -coordinate. For fixed , the set can be viewed as a bounded subset of , whose volume tends to as .

The set itself is bounded and the formula

 ΘU=f1⩾0∧⋯∧fr⩾0∧0

defines in . In addition, is a union of connected components of . Indeed, for any point  with , it holds that . This implies that where is the interior of . Therefore, is both relatively closed (as the trace of ) and open (as that of ) in .

We are hence in the setting of the previous subsection. Since by definition of , Theorem 2 applies, and the function is annihilated by an operator which is computed using the routine PicardFuchs introduced earlier. By Proposition 2, one can choose rational points  such that the values of  at these points characterize it among the solutions of , and, given sufficiently precise approximations of , one can compute to any desired accuracy.

The “initial conditions” are computed by calls to Algorithm 1 with  and  specialized to . In the notation of §2.1, this corresponds to taking and . Thus, is compact, and, since no  can change sign on a connected component of for , it is the union of those connected components of  where . Additionally, (R) holds for since . Therefore, the assumptions of Theorem 3 are satisfied.

We obtain Algorithm 2 (which uses the same subroutines and conventions as Algorithm 1) and the following correctness theorem.

###### Theorem 4 ().

Let . Let be the semi-algebraic set defined by . Assume that  is bounded. Then, given and a working precision , Algorithm 2 (Volume) computes an interval containing of width as for fixed

###### Remark 5 ().

In case  has empty interior, Algorithm 2 returns zero. When is contained in a linear subspace of dimension , one could in principle obtain the -volume of  by computing linear equations defining the subspace (using quantifier elimination as in (Khachiyan and Porkolab, 2000; Safey El Din and Zhi, 2010)) and eliminating variables. The new system would in general have algebraic instead of rational coefficients, though.

Lastly, we note that a more direct symbolic computation of integrals on general semi-algebraic sets depending on a parameter is possible with Oaku’s algorithm (Oaku, 2013), based on the effective theory of -modules.

## 3. Periods depending on a parameter

Let us now discuss in more detail the main black boxes used by the volume computation algorithm. In this section, we study how the volume of a section varies with the parameter .

### 3.1. Picard-Fuchs equations

Let  be a rational function. A period of the parameter-dependent rational integral is an analytic function , for some open subset of  or  such that for any  there is an -cycle and a neighborhood  of  such that for any , is disjoint from the poles of  and

 (3) ϕ(t)=∫γR(t,x1,…,xn)dx1⋯dxn.

Recall that an -cycle is a compact -dimensional real submanifold of and that such an integral is invariant under a continuous deformation of the integration domain  as long as it stays away from the poles of , as a consequence of Stokes’ theorem. It is also well known that such a function  depends analytically on , by Morera’s theorem for example.

For instance, algebraic functions are periods: if  satisfies a nontrivial relation , with square-free , then is a period by the residue theorem applied to

 ϕ(t)=12πi∮γxP(t,x)∂P∂x(t,x)dx

where encloses  and no other root of . Indeed, the integrand decomposes as , where the functions  parametrize the roots of , and, w.l.o.g., .

Periods of rational functions are solutions of Fuchsian linear differential equations with polynomial coefficients known as Picard-Fuchs equations. This was proved in (Picard, 1902) in the case of three variables at most and a parameter and generalized later, using either the finiteness of the algebraic De Rham cohomology (e.g. Grothendieck, 1966; Monsky, 1972; Christol, 1985) or the theory of D-finite functions (Lipshitz, 1988). The regularity of Picard-Fuchs equations is due to Griffiths (see Katz, 1971).

###### Theorem 6 ().

If  is the period of a rational integral then  is solution of a nontrivial linear differential equation with polynomial coefficients , where the operator  belongs to the class  introduced in §2.1.

Several algorithms are known and implemented to compute such Picard-Fuchs equations (Lairez, 2016; Koutschan, 2010; Chyzak, 2000).

###### Theorem 7 ((Bostan et al., 2013)).

A period of the form (3) is solution of a differential equation of order at most  where  is the degree of ; and one can compute such an equation in operations in .

Note however that the algorithm underlying this result might not return the equation of minimal order, but rather a left multiple of the Picard-Fuchs equation. So there is no guarantee that the computed operator belongs to . On the other hand, Lairez’s algorithm (Lairez, 2016) can compute a sequence of operators with non-increasing order which eventually stabilizes to the minimal order operator. In particular, as long as the computed operator is not in , we can compute the next one, with the guarantee that this procedure terminates. A conjecture of Dimca (Dimca, 1991) ensures that it terminates after at most  steps, leading to a  complexity bound as in Theorem 7.

### 3.2. Volume of a section and proof of Theorem 2

We prove Theorem 2 as a consequence of Theorem 6 and the following result. It is probably well known to experts but it is still worth an explicit proof. We use the notation of §2.

###### Theorem 8 ().

If and if is bounded then the function  is a period of the rational integral

 12iπ∮x1fρ∂fρ∂x1dx1⋯dxn.
###### Proof.

Let . By Stokes’ formula,

 volnUρ=∫Uρdx1⋯dxn=∮∂Uρx1dx2⋯dxn,

where  is the boundary of . Due to the regularity assumption , the gradient does not vanish on the real zero locus of , denoted . Because is a union of connected components of , it follows that is a compact -dimensional submanifold of  contained in .

For , let be the Leray tube defined by

 τ(ρ)≜{p+u∇pfρ ∣∣ p∈∂Uρ,u∈C and |u|=ε}.

This is an -dimensional submanifold of . We choose small enough that : this is possible because does not vanish on  which is compact.

Let ; observe that does not cancel the denominator of . Leray’s residue theorem (Leray, 1959) shows that

 2πi∮∂Utx1dx2⋯dxn= ∮τ(ρ)dfρfρ∧(x1dx2⋯dxn) = ∮τ(ρ)R(ρ,x1,…,xn)dx1⋯dxn.

(In Pham’s (Pham, 2011, Thm. III.2.4) notation, we have , , , and .)

To match the definition of a period and conclude the proof, it is enough to prove that, locally, the integration domain  can be made independent of . And indeed, since  is a union of connected components of , we have . Therefore, since  is connected and , the restriction of the projection  defines a submersive map from onto . Additionally, is compact, hence this map is proper. Ehresmann’s theorem then implies that there exists a continuous map such that  induces a homeomorphism  for any . In particular, we have

 τ(σ)={h(σ,p)+u∇h(σ,p)fσ ∣∣ p∈∂Uσ,u∈C and |u|=ε}.

This formulation makes it clear that  deforms continuously into  as  varies. Since does not intersect the polar locus of , neither does when  and  are close enough, by compactness of  and continuity of the deformation. Therefore, given any , we have for close enough to . ∎

The choice of as a primitive of in Theorem 8 is arbitrary, but of little consequence, since the Picard-Fuchs equation only depends on the cohomology class of the integrand.

### 3.3. Critical values

Theorem 2 does not guarantee that  satisfies the Picard-Fuchs equation on the whole domain where the equation is nonsingular. It could happen that the solutions extend analytically across an exceptional point, or that some of them have singularities between two consecutive exceptional points. As a consequence, we need to explicitly compute .

###### Lemma 9 ().

There exists an algorithm which, given on input a polynomial of degree satisfying (R), computes a polynomial of degree whose set of real roots contains , using operations in .

###### Proof.

Recall that, when (R) holds, the set  is finite. Our goal is to write as the root set of a univariate polynomial . Consider the polynomial We start by computing at least one point in each connected component of the real algebraic set defined by using (Basu et al., 2006, Algorithm 13.3). By (Basu et al., 2006, Theorem 13.22), this algorithm uses operations. It returns a rational parametrization: polynomials , , in of degree such that is square-free and the set of points

meets every connected component of the zero set of . In particular, . As a polynomial , we take the resultant with respect to  of  and : its set of roots contains . Since and have degree , this last step also uses operations in (von zur Gathen and Gerhard, 1999). ∎

## 4. Numerics

Let us turn to the numerical part of the main algorithm. It is known (Chudnovsky and Chudnovsky, 1990; van der Hoeven, 2001) that Fuchsian differential equations with coefficients in  can be solved numerically in quasi-linear time w.r.t. the precision. Yet, some minor technical points must be addressed to apply the results of the literature to our setting. We start with reminders on the theory of linear ODEs in the complex domain (e.g. Poole, 1936; Hille, 1976). Consider a linear differential operator

 (4) P=pm(t)dmdtm+⋯+p1(t)ddt+p0(t)

of order  with coefficients in .

Recall that is a singular point of  when the leading coefficient  of  vanishes at . A point that is not a singular point is called ordinary. Singular points are traditionally classified in two categories: a singular point is a regular singular point of  if, for , its multiplicity as a pole of is at most , and an irregular singular point otherwise. The point at infinity in is said to be ordinary, singular, etc., depending on the nature of  after the change of variable . An operator with no irregular singular point in is called Fuchsian.

Fix a simply connected domain containing only ordinary points of , and let be the space of analytic solutions of the differential equation . According to the Cauchy existence theorem for linear analytic ODEs, is a complex vector space of dimension . A particular solution is determined by the initial values at any point .

At a singular point, there may not be any nonzero analytic solution. Yet, if  is a regular singular point, the differential equation still admits  linearly independent solutions defined in the slit disk for small enough  and each of the form

 (5) y(u+ζ)=ζγℓ∑k=0yk(ζ)log(ζ)k=ℓ∑k=0∞∑ν∈γ+Nyk,νζνlog(ζ)k

where , , and for exactly one  (Poole, 1936, §16). The functions  are analytic for  (including at ). The algebraic numbers  are called the exponents of  at .

Suppose now that  is either an ordinary point of  lying in the topological closure  of , or a regular singular point of  situated on the boundary of . As a result of the previous discussion, we can choose a distinguished basis of  in which each  is characterized by the leading monomial111 More precisely, denoting in (5), there are  computable pairs such that, for all , we have , for , and whenever . of its local expansion (5) at . At an ordinary point  for instance, the coefficients of the decomposition of a solution  on  are , that is, essentially the classical initial values. Observe that when no two exponents  have the same imaginary part, the elements of  all have distinct asymptotic behaviours as . In particular, at most one of them tends to a nonzero finite limit. As Picard-Fuchs operators have real exponents according to Theorem 6, this observation applies to them.

Let be a second point subject to the same restrictions as . Let be the transformation matrix from  to . The key to the quasi-linear complexity of our algorithm is that the entries of this matrix can be computed efficiently, by solving the ODE with a Taylor method in which sums of Taylor series are computed by binary splitting (Beeler et al., 1972, item 178), (Chudnovsky and Chudnovsky, 1990). The exact result we require is due to van der Hoeven (van der Hoeven, 2001, Theorems 2.4 and 4.1); see also (Mezzarobba, 2011) for a detailed algorithm and some further refinements. Denote by the complexity of -bit integer multiplication.

###### Theorem 10 ((van der Hoeven, 2001)).

For a fixed operator  and fixed algebraic numbers  as above, one can compute the matrix with an entry-wise error bounded by  in operations.

Since is linear, this result suffices to implement the procedure DSolve required by the main algorithm. More precisely, suppose that returns a matrix of complex intervals of width  that encloses entry-wise.

###### Definition 11 ().

A system of good initial conditions for  on , denoted , is a finite family of pairs where and is a linear form that belongs to the dual basis of  for some algebraic point  (which may depend on ), with the property that span the dual space of .

A system of good initial conditions on is a system of good initial conditions on for some .

In other words, a system of good initial conditions is a choice of coefficients of local decompositions of a solution of  whose values determine at most one solution, and of prescribed values for these coefficients. When the system is compatible, we say that it defines the unique solution of  that satisfies all the constraints. Let us note in passing the following fact, which was used in §2.2.

###### Lemma 12 ().

Let be ordinary points of  such that is a system of good initial conditions for on , and let . Then is a system of good initial conditions for  on .

###### Proof.

The derivative maps the solution space of to that of , and its kernel consists exactly of the constant functions. By assumption, a solution of  is completely defined by its values at , hence a solution of  is characterized by the values of its derivative at the same points, along with its limit at . Because has order at least  (otherwise, would not be a system of good initial conditions), the conditions  are of the form with belonging to the dual basis of some , as required. So is the condition since has solutions with a nonzero finite limit at . ∎

Algorithm 3 evaluates the solution of an operator  given by a system of good initial conditions. Note that the algorithm is allowed to fail. It fails if the intervals are not accurate enough for the linear algebra step on line 10 to succeed, or if the linear system, which is in general over-determined, has no solution. The following proposition assumes a large enough working precision  to ensure that this does not happen. Additionally, we only require that the output be accurate to within , so as to absorb any loss of precision resulting from numerical stability issues or from the use of interval arithmetic.

###### Proposition 12 ().

Suppose that the operator  is Fuchsian with real exponents. Let be real algebraic numbers, and let be a real analytic solution of on the interval such that tends to a finite limit as . Let  be a system of good initial conditions for  on that defines .

Given the operator , the point , a large enough working precision , and an approximation of  where is an interval of width at most containing , (Algorithm 3) computes a real interval of width containing in time .

###### Proof.

At the end of the loop, we have , , and the entries of are intervals of width . The coefficients of the decomposition of  in the basis  satisfy for all , where is the th row of . As  is a system of good initial conditions, the linear system has no other solution. Step 10 hence succeeds in solving the interval version as soon as the  and the entries of the  are thin enough intervals. It then returns intervals of width .

We assumed that tends to a finite limit at . It follows that the decomposition of  on  only involves the basis elements with a finite limit at . Either  contains an element that tends to , in which case , or every solution that converges tends to zero, and then the limit is zero. Since, by assumption, is real, we can ignore the imaginary part of the computed value. In both cases, the algorithm, when it succeeds, returns a real interval of width containing .

As for the complexity analysis, all including  are algebraic, hence Theorem 10 applies and shows that each call to TransitionMatrix runs in time . The matrix multiplications at step 8 take operations. The cost of solving the linear system (which is of bounded size) is as well. The cost of the remaining steps is independent of