Computing the volume of compact semi-algebraic sets
CComputing the volume of compact semi-algebraic sets
Pierre Lairez
Marc Mezzarobba
Sorbonne Université, CNRS,Laboratoire d’Informatique de Paris 6,LIP6, Équipe PeQuaNF-75252, Paris Cedex 05, [email protected]
Mohab Safey El Din
Sorbonne Université, CNRS, Inria,Laboratoire d’Informatique de Paris 6,LIP6, Équipe PolSysF-75252, Paris Cedex 05, [email protected]
ABSTRACT
Let S ⊂ R n be a compact basic semi-algebraic set defined as thereal solution set of multivariate polynomial inequalities with ra-tional coefficients. We design an algorithm which takes as input apolynomial system defining S and an integer p ⩾ n -dimensional volume of S at absolute precision 2 − p .Our algorithm relies on the relationship between volumes ofsemi-algebraic sets and periods of rational integrals. It makes useof algorithms computing the Picard-Fuchs differential equation ofappropriate periods, properties of critical points, and high-precisionnumerical integration of differential equations.The algorithm runs in essentially linear time with respect to p .This improves upon the previous exponential bounds obtained byMonte-Carlo or moment-based methods. Assuming a conjecture ofDimca, the arithmetic cost of the algebraic subroutines for comput-ing Picard-Fuchs equations and critical points is singly exponentialin n and polynomial in the maximum degree of the input. KEYWORDS
Semi-algebraic sets; Picard-Fuchs equations; Symbolic-numeric al-gorithms
ACM Reference Format:
Pierre Lairez, Marc Mezzarobba, and Mohab Safey El Din. 2019. Computingthe volume of compact semi-algebraic sets. In
International Symposium onSymbolic and Algebraic Computation (ISSAC ’19), July 15–18, 2019, Beijing,China.
ACM, New York, NY, USA, 8 pages. https://doi.org/10.1145/3326229.3326262
Semi-algebraic sets are the subsets of R n which are finite unions ofreal solution sets to polynomial systems of equations and inequa-lities with coefficients in R . Starting from Tarski’s algorithm forquantifier elimination [39] improved by Collins through the Cylin-drical Algebraic Decomposition algorithm [11], effective real al-gebraic geometry yields numerous algorithmic innovations and 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 Horizon2020 research and innovation programme under the Marie Skłodowska-Curie grantagreement N° 813211 (POEMA)..Publication rights licensed to ACM. ACM acknowledges that this contribution wasauthored or co-authored by an employee, contractor or affiliate of a national govern-ment. As such, the Government retains a nonexclusive, royalty-free right to publish orreproduce this article, or to allow others to do so, for Government purposes only.
ISSAC ’19, July 15–18, 2019, Beijing, China © 2019 Copyright held by the owner/author(s). Publication rights licensed to theAssociation for Computing Machinery.ACM ISBN 978-1-4503-6084-5/19/07...$15.00https://doi.org/10.1145/3326229.3326262 asymptotically faster routines for problems like deciding the empti-ness of semi-algebraic sets, answering connectivity queries or com-puting Betti numbers [e.g., 3, 6, 7, 36, 37]. The output of all thesealgorithms is algebraic in nature. In this paper, we study the prob-lem of computing the volume of a (basic) compact semi-algebraicset S ⊂ R n defined over Q . The output may be transcendental: forinstance, the area of the unit circle in R is π .Volumes of semi-algebraic sets actually lie in a special class of realnumbers, for they are closely related to Kontsevich-Zagier periods introduced in [23]. A (real) period is the value of an absolutelyconvergent integral of a rational function with rational coefficientsover a semi-algebraic set defined by polynomials with rationalcoefficients. For example, algebraic numbers are periods, as are π , log 2, ζ ( ) . Since vol S = ∫ S x , volumes of semi-algebraic setsdefined over Q are periods. Conversely, interpreting an integralas a “volume under a graph” shows that periods are differences ofvolumes of semi-algebraic sets defined over Q . In [41], it is furthershown that periods are differences of volumes of compact semi-algebraic sets defined over Q .The problem we consider in this paper is thus a basic instance ofthe more general problem of integrating an algebraic function overa semi-algebraic set; it finds applications in numerous areas of engi-neering sciences. Performing these computations at high precision(hundreds to thousands of digits) is also relevant in experimentalmathematics, especially for discovering formulas, as explained, forexample, in [1]. Most of the examples featured in this reference areperiods, sometimes in disguise. Prior work.
The simplest semi-algebraic sets one can consider arepolytopes. The computation of their volume has been extensivelystudied, with a focus on the complexity with respect to the dimen-sion. It is known that even approximating the volume of a polytopedeterministically is ε -approximation in time polynomial inthe dimension of the set and 1 / ε . A key ingredient for this algo-rithm is a Monte Carlo method for efficiently sampling points froma convex set. Since then, Monte Carlo schemes have been adoptedas the framework of several volume estimation algorithms.In contrast, we deal here with compact semi-algebraic sets whichcan be non-convex and even non-connected. Additionally, whilevolumes of polytopes are rational, the arithmetic nature of volumesof semi-algebraic sets is much less clear, as unclear as the nature ofperiods. This raises the question of the computational complexityof 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 S and estimates a r X i v : . [ c s . S C ] A p r SSAC ’19, July 15–18, 2019, Beijing, China Pierre Lairez, Marc Mezzarobba, and Mohab Safey El Din the probability that they lie in S . This method is of practical interestat low precision but requires 2 Ω ( p ) samples to achieve an errorbounded by 2 − p with high probability. We refer to [22] which dealswith definable sets, a class which encompasses semi-algebraic sets.In a different direction, numerical approximation schemes basedon the moment problem and semi-definite programming have beendesigned in [16]. They are also of practical interest at low precision,and can provide rigorous error bounds, but the convergence isworse than exponential with respect to p [24].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 [9] on the high-precision numericalsolution of ODEs with polynomial coefficients which was moti-vated, among other things, by applications to periods of Abelianintegrals [see 9, p. 133]. Main result.
We describe a new strategy for computing volumesof semi-algebraic sets, at the crossroads of effective algebraic andreal algebraic geometry, symbolic integration, and rigorous nu-merical computing. Our approach effectively reduces the volumecomputation to the setting of [9]. It yields an algorithm that ap-proximates the volume of a fixed , bounded basic semi-algebraic setin almost linear time with respect to the precision. More precisely,we prove the following bit complexity estimate.Theorem 1.
Let f , . . . , f r be polynomials in Q [ x , . . . , x n ] , andlet S ⊂ R n be the semi-algebraic set defined by f ⩾ , . . . , f r ⩾ .Assume that S is compact. There exists an algorithm which computes,on input p ⩾ and ( f , . . . , f r ) , an approximation ˜ V of the volume V of S with | ˜ V − V | ⩽ − p . When f , . . . , f r are fixed, the algorithmruns in time O ( p log ( p ) + ε ) (for any ε > ) as p → ∞ . The algorithm recursively computes integrals of volumes ofsections of S . Let v ( t ) denote the ( n − ) -dimensional volume of S ∩ pr − ( t ) , for some nonzero linear projection pr : R n → R . In oursetting, v is a piecewise analytic function and, except at finitelymany t , is solution of a linear differential equation with polynomialcoefficients known as a Picard-Fuchs equation.The problem points belong to the critical locus of the restric-tion of the projection pr to a certain hypersurface containing theboundary of S and are found by solving appropriate polynomialsystems. (Compare [20] in the case of polytopes.) The Picard-Fuchsequation for v is produced by algorithms from symbolic integration,in particular [5, 26]. To obtain the volume of S , it then suffices tocompute ∫ R v with a rigorous numerical ODE solving algorithm,starting from values v ( ρ i ) at suitable points ρ i obtained throughrecursive calls.The complexity with respect to the dimension n of the ambientspace and the number r , maximum degree D , and coefficient size ofthe polynomials f i is harder to analyze. We will see, though, thatunder reasonable assumptions, the “algebraic” steps (computingthe critical loci and of the Picard-Fuchs equations) take at most ( rD ) O ( n ) arithmetic operations in Q . Example.
The idea of the method is well illustrated by the exam-ple of a torus S , here of major radius 2 and minor radius 1. Let S = (cid:8) ( x , y , z ) ∈ R (cid:12)(cid:12) ( x + y + z + ) ⩽ ( x + y ) (cid:9) . Figure 1: Volume of the sections of the torus S as a functionof the parameter t . In red, a singular section. The area (2-dimensional volume) of a section S ∩ { x = t } definesa function v : R → R (see Figure 1). It is analytic, except maybeat the critical values t = ± t = ± ( t + y + z + ) = ( t + y ) is singular. On each intervalon which v is analytic, it satisfies the Picard-Fuchs equation ( t − )( t + )( t − ) ( t + ) t v ′′′ ( t ) + ( t + )( t − ) ( t + ) tv ′′ ( t )− ( t + t − )( t − )( t + ) v ′ ( t ) + ( t + ) t v ( t ) = , (1)which we compute in 2 seconds on a laptop using the algorithmof [26] and Theorem 9.We know some special values of v , namely v ( ) = π , v (± ) = v (± ) =
0. Additionally, we have v ( ± t ) = O ( t ) as t →∓
0. These properties characterize the analytic function v |(− , ) inthe 2-dimensional space of analytic solutions of the differentialequation (1) on (− , ) , and similarly for v |( , ) . (Our algorithmactually uses recursive calls at generic points instead of these adhoc conditions.) The rigorous ODE solver part of the Sage packageore_algebra [30] determines in less than a second that ∫ − v ( t ) d t = . [ ... ] ± − . And indeed, it is not hard to see in this case that vol S = π . Wecan obtain 1000 digits in less than a minute. Outline.
The remainder of this article is organized as follows. InSection 2, we give a high-level description of the main algorithm. Assketched above, the algorithm relies on the computation of criticalpoints, Picard-Fuchs equations, and numerical solutions of theseequations. In Section 3, we discuss the computation of Picard-Fuchsequations and critical points, relating these objects with analyticityproperties of the “section volume” function. Then, in Section 4, wedescribe the numerical solution process and study its complexitywith respect to the precision. Finally, in Section 5, we conclude theproof of Theorem 1 and state partial results on the complexity ofthe algorithm with respect to n , r , and D . Acknowledgements.
We would like to thank the anonymous re-viewers for their careful reading and valuable comments. omputing the volume of compact semi-algebraic sets ISSAC ’19, July 15–18, 2019, Beijing, China
We start by designing an algorithm which deals with the case of aunion of connected components of a semi-algebraic set defined bya single inequality. Next, we will use a deformation technique tohandle semi-algebraic sets defined by several inequalities.
Let f ∈ Q [ t , x , . . . , x n ] and A be the semi-algebraic set A ≜ (cid:8) ( ρ , x ) ∈ R × R n (cid:12)(cid:12) f ( ρ , x ) ⩾ (cid:9) . Let pr : R n + → R be the projection on the t -coordinate. We wantto compute the volume of a union U of connected components of A starting from the volumes of suitable fibers U ∩ pr − ( ρ ) . For tech-nical reasons, we first consider the slightly more general situationwhere U is a union of connected components of A ∩ pr − ( I ) forsome open interval I ⊆ R . From a computational point of view, weassume that U is described by a semi-algebraic formula Θ U , that is, U = {( ρ , x ) ∈ A | Θ U ( ρ , x )} , where Θ U is a finite disjunction of conjunctions of polynomialinequalities with (in our setting) rational coefficients.For ρ ∈ I , let U ρ ≜ U ∩ pr − ( ρ ) and v ( ρ ) ≜ vol n U ρ . Let Σ f ⊆ R (we will often omit the subscript f ) be the set of exceptional values Σ f ≜ (cid:110) ρ ∈ R | ∃ x ∈ R n , f ( ρ , x ) = ∧ ∀ i , ∂∂ x i f ( ρ , x ) = (cid:111) . (2)Thus, when f is square-free, exceptional values are either criticalvalues of the restriction to the hypersurface { f = } of the pro-jection pr, or images of singular points of { f = } . By definitionof Σ , for any ρ ∈ R \ Σ , the zero set of f ρ = f ( ρ , −) is a smoothsubmanifold of R n .Further, we say that assumption (R) holds for f if (cid:110) z ∈ R n + (cid:12)(cid:12)(cid:12) f ( z ) = ∧ ∂∂ t f ( z ) = ∧ ∀ i , ∂∂ x i f ( z ) = (cid:111) = ∅ . (R)Observe that by Sard’s theorem [e.g. 3, 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 D ⊂ Q [ t ][ dd t ] denote the set of Fuchsian linear differential op-erators with coefficients in Q [ t ] whose local exponents at singularpoints are rational (see §4 for reminders on Fuchsian operators andtheir exponents).Theorem 2. If U is bounded and I ∩ Σ = ∅ , then the function v | I is solution of a computable differential equation of the form P ( v ) = ,where P ∈ D depends only on f . We will also use the following proposition, which summarizesthe results of Proposition 14 and Lemma 15 in §4. The completedefinition of “good initial conditions” is given there as well. Up totechnical details, this simply means a system I of linear equationsof the form y ( k ) ( u ) = s that suffices to characterize a particularsolution y among the solutions of P ( y ) =
0. An ε -approximation of I is made of the same equations with each right-hand side s replaced by an enclosure ˜ s ∋ s of diameter ⩽ ε .Proposition 3. Let P ∈ D have order m , and let J = ( α , β ) be areal interval with algebraic endpoints. Let y : J → R be a solution of P ( y ) = with a finite limit at α and I be a system of good initialconditions for P on J defining y . Algorithm 1
Volume of U at precision O ( − p ) procedure Volume1( f , Θ U , ( t , x , . . . , x n ) , p ) if n = then return UnivariateVolume ( f , Θ U , p ) . ( α , . . . , α ℓ ) ← CriticalValues ( f , t ) P ← PicardFuchs ( f , t ) for ⩽ i ⩽ ℓ − do ▷ ˜ s j , ˜ S i are intervals ( ρ , . . . , ρ m ) ← PickGoodPoints ( P , α i , α i + ) for ⩽ j ⩽ m do ˜ s j ← Volume1 ( f | t = ρ j , Θ U | t = ρ j , ( x , . . . , x n ) , p ) ˜ I ← [ y ′ ( ρ ) = ˜ s , . . . , y ′ ( ρ m ) = ˜ s m , y ( α i + ) = ] ˜ S i ← − DSolve ( P dd t , ˜ I , α i , p ) return ˜ S + · · · + ˜ S ℓ (1) Given P , α , a precision p ∈ N and a − p -approximation ˜ I of I , one can compute an interval of width O ( − p ) (as p → ∞ for fixed P , α , and I ) containing lim t → α y ( t ) .(2) Given P , α , β , one can compute ρ , . . . , ρ m ∈ J ∩ Q such thatthe y ( ρ j ) form a system of good initial conditions for P on J . Assume now that U is a bounded union of connected componentsof A (i.e., that we can take I = R above), and that (R) holds for f .The algorithm is recursive. Starting with input f , Θ U , and p , it firstcomputes the set Σ = { α ⩽ · · · ⩽ α ℓ } of exceptional values soas to decompose R − Σ into intervals over which the function v satisfies the differential equation P ( y ) = U is bounded, one hasvol n + U = ℓ − (cid:213) i = vol n + (cid:16) U ∩ pr − ( α i , α i + ) (cid:17) = ℓ − (cid:213) i = ∫ α i + α i v ( t ) d t . Fix i and consider the interval J = ( α i , α i + ) . Since v | J is anni-hilated by P , its anti-derivative w : J → R vanishing at α i + isannihilated by the operator P dd t , which belongs to D since P does.Additionally, if [ v ( ρ j ) = s j ] j is a system of good initial conditionsfor P that defines v | J , then [ w ′ ( ρ j ) = s j ] j ∪[ w ( α i + ) = ] is a systemof good initial conditions for P dd t defining w (see Lemma 13 in §4).Thus, by Proposition 3, to compute w ( α i ) to absolute precision p , itsuffices to compute v ( ρ j ) , 1 ⩽ j ⩽ m , to precision p + O ( ) .By definition of Σ , since ρ j (cid:60) Σ , there is no solution to the system f ( ρ j , −) = ∂∂ x f ( ρ j , −) = · · · = ∂∂ x n f ( ρ j , −) = f ( ρ j , −) . Additionally, U ∩ pr − ( ρ j ) is a bounded union of connected components of A ∩ pr − ( ρ j ) . Hence,the values v ( ρ j ) can be obtained by recursive calls to the algorithmwith t instantiated to ρ j .The process terminates since each recursive call handles one lessvariable. In the base case, we are left with the problem of computingthe length of a union of real intervals encoded by a semi-algebraicformula. This is classically done using basic univariate polynomialarithmetic and real root isolation [3, Chap. 10].The complete procedure is formalized in Algorithm 1. The quan-tities denoted with a tilde in the pseudo-code are understood to berepresented by intervals, and the operations involving them followthe semantics of interval arithmetic. Additionally, we assume thatwe have at our disposal the following subroutines: SSAC ’19, July 15–18, 2019, Beijing, China Pierre Lairez, Marc Mezzarobba, and Mohab Safey El Din • PicardFuchs ( f , t ) , DSolve ( P , ˜ I , α , p ) , and PickGoodPoints ( P , α , β ) , which implement the algorithms implied, re-spectively, by Theorem 2 and Proposition 3 (1) and (2); • CriticalValues ( f , t ) , which returns an encoding for a finiteset of real algebraic numbers containing the exceptionalvalues associated to f , sorted in increasing order; • UnivariateVolume ( д , Θ U , p ) where д ∈ Q [ t ] and Θ U is asemi-algebraic formula describing a union U of connectedcomponents of { д ⩾ } , which returns an interval of width ⩽ − p containing vol U .The following result summarizes the above discussion.Theorem 4. Assume that U is a bounded union of connectedcomponents of A and that (R) holds. Then, on input f , Θ U , p and ( t , x , . . . , x n ) , Algorithm 1 (Volume1) returns a real interval of width O ( − p ) (for fixed f ) containing vol n + U . Now, we show how to compute the volume of a basic semi-algebraicset S ⊂ R n defined by f ⩾ , . . . , f r ⩾ , f i ∈ Q [ x , . . . , x n ] , assuming that S is compact.We set f = f · · · f r − t ∈ Q [ t , x , . . . , x n ] , and consider thesemi-algebraic set A ⊂ R n + defined by f ⩾
0. Observe that thepolynomial f satisfies (R) because ∂ f ∂ t = −
1. We can hence choosean interval I = ( , α ) with α ∈ Q that contains no element of Σ f .Let U ≜ A ∩ ( I × S ) and pr be the projection on the t -coordinate.For fixed ρ ∈ I , the set U ∩ pr − ( ρ ) can be viewed as a boundedsubset of S , whose volume v ( ρ ) = vol n ( U ∩ pr − ( ρ )) tends to vol n S as ρ → U itself is bounded and the formula Θ U = f ⩾ ∧ · · · ∧ f r ⩾ ∧ < t < α defines U in A . In addition, U is a union of connected componentsof A ∩ pr − ( I ) . Indeed, for any point ( ρ , x ) ∈ A with ρ ∈ I , it holdsthat f ( x ) · · · f r ( x ) >
0. This implies that U = A ∩ ( I × ˚ S ) where ˚ S isthe interior of S . Therefore, U is both relatively closed (as the traceof R × S ) and open (as that of R × ˚ S ) in A ∩ pr − ( I ) .We are hence in the setting of the previous subsection. Since I ∩ Σ f = ∅ by definition of I , Theorem 2 applies, and the func-tion v : I → R is annihilated by an operator P ∈ D which iscomputed using the routine PicardFuchs introduced earlier. ByProposition 3, one can choose rational points ρ j ∈ I such that thevalues of v at these points characterize it among the solutions of P ,and, given sufficiently precise approximations of v ( ρ j ) , one cancompute vol n S = lim t → v ( t ) to any desired accuracy.The “initial conditions” v ( ρ j ) are computed by calls to Algo-rithm 1 with f and Θ U specialized to t = ρ j . In the notation of§2.1, this corresponds to taking A = A ( ρ j ) = { f · · · f r ⩾ ρ j } and U = U ( ρ j ) = A ( ρ j ) ∩ S . Thus, U ( ρ j ) is compact, and, since no f i can change sign on a connected component of A ( ρ ) for ρ > A ( ρ j ) where f , . . . , f r ⩾
0. Additionally, (R) holds for f ( ρ j , −) since ρ j (cid:60) Σ f .Therefore, the assumptions of Theorem 4 are satisfied.We obtain Algorithm 2 (which uses the same subroutines andconventions as Algorithm 1) and the following correctness theorem. Algorithm 2
Volume of S procedure Volume( ( f , . . . , f r ) , p ) f ← f · · · f r − t ( α , . . . , α ℓ ) ← CriticalValues ( f , t ) α ← a rational s.t. 0 < α < min ({ α i | α i > } ∪ { }) Θ U ← f ⩾ ∧ · · · ∧ f r ⩾ ∧ < t < α P ← PicardFuchs ( f , t ) ( ρ , . . . , ρ m ) ← PickGoodPoints ( P , , α ) for ⩽ j ⩽ m do ▷ ˜ s j are intervals ˜ s j ← Volume1 ( f | t = ρ j , ( Θ U ) | t = ρ j , ( x , . . . , x n ) , p ) return DSolve ( P , [ y ( ρ j ) = ˜ s j ] mj = , p ) Theorem 5.
Let f , . . . , f r ∈ Q [ x , . . . , x n ] . Let S be the semi-algebraic set defined by f ⩾ , . . . , f r ⩾ . Assume that S is bounded.Then, given ( f , . . . , f r ) and a working precision p ∈ N , Algorithm 2(Volume) computes an interval containing vol n ( S ) of width O ( − p ) as p → ∞ for fixed f , . . . , f r Remark 6.
In case S has empty interior, Algorithm 2 returns zero.When S is contained in a linear subspace of dimension k < n , onecould in principle obtain the k -volume of S by computing linearequations defining the subspace (using quantifier elimination asin [21, 38]) and eliminating n − k variables. The new system wouldin general have algebraic instead of rational coefficients, though.Lastly, we note that a more direct symbolic computation of inte-grals on general semi-algebraic sets depending on a parameter ispossible with Oaku’s algorithm [32], based on the effective theoryof D -modules. Let us now discuss in more detail the main black boxes used by thevolume computation algorithm. In this section, we study how thevolume of a section U ∩ pr − ( ρ ) varies with the parameter t = ρ . Let R ( t , x , . . . , x n ) be a rational function. A period of the parameter-dependent rational integral ∮ R ( t , x , . . . , x n ) d x · · · d x n is an ana-lytic function ϕ : Ω → C , for some open subset Ω of R or C suchthat for any s ∈ Ω there is an n -cycle γ ⊂ C n and a neighbor-hood Ω ′ ⊂ Ω of s such that for any t ∈ Ω ′ , γ is disjoint from thepoles of R ( t , −) and ϕ ( t ) = ∫ γ R ( t , x , . . . , x n ) d x · · · d x n . (3)Recall that an n -cycle is a compact n -dimensional real submanifoldof C n and that such an integral is invariant under a continuousdeformation of the integration domain γ as long as it stays awayfrom the poles of R ( t , −) , as a consequence of Stokes’ theorem. It isalso well known that such a function ϕ depends analytically on t ,by Morera’s theorem for example.For instance, algebraic functions are periods: if ϕ : Ω → C satisfies a nontrivial relation P ( t , ϕ ( t )) =
0, with square-free P ∈ C [ t , x ] , then ϕ ( t ) is a period by the residue theorem applied to ϕ ( t ) = πi ∮ γ xP ( t , x ) ∂ P ∂ x ( t , x ) d x omputing the volume of compact semi-algebraic sets ISSAC ’19, July 15–18, 2019, Beijing, China where γ ⊂ C encloses ϕ ( t ) and no other root of P . Indeed, the inte-grand decomposes as (cid:205) deg x Pi = x /( x − ψ i ( t )) , where the functions ψ i parametrize the roots of P ( t , −) , and, w.l.o.g., ϕ = ψ .Periods of rational functions are solutions of Fuchsian linear dif-ferential equations with polynomial coefficients known as Picard-Fuchs equations. This was proved in [34] in the case of three vari-ables at most and a parameter and generalized later, using either thefiniteness of the algebraic De Rham cohomology [e.g. 8, 15, 31] orthe theory of D-finite functions [28]. The regularity of Picard-Fuchsequations is due to Griffiths [see 18].Theorem 7. If ϕ : Ω → C is the period of a rational integralthen ϕ is solution of a nontrivial linear differential equation withpolynomial coefficients P ( ϕ ) = , where the operator P belongs to theclass D introduced in §2.1. Several algorithms are known and implemented to compute suchPicard-Fuchs equations [10, 25, 26].Theorem 8 ([5]).
A period of the form (3) is solution of a differen-tial equation of order at most D n where D is the degree of R ; and onecan compute such an equation in D O ( n ) operations in Q . Note however that the algorithm underlying this result might notreturn the equation of minimal order, but rather a left multiple of the
Picard-Fuchs equation. So there is no guarantee that the computedoperator belongs to D . On the other hand, Lairez’s algorithm [26]can compute a sequence of operators with non-increasing orderwhich eventually stabilizes to the minimal order operator. In partic-ular, as long as the computed operator is not in D , we can computethe next one, with the guarantee that this procedure terminates. Aconjecture of Dimca [12] ensures that it terminates after at most n steps, leading to a D O ( n ) complexity bound as in Theorem 8. We prove Theorem 2 as a consequence of Theorem 7 and the fol-lowing result. It is probably well known to experts but it is stillworth an explicit proof. We use the notation of §2.Theorem 9. If I ∩ Σ = ∅ and if U is bounded then the function ρ ∈ I (cid:55)→ vol n U ρ is a period of the rational integral iπ ∮ x f ρ ∂ f ρ ∂ x d x · · · d x n . Proof. Let ρ ∈ I . By Stokes’ formula,vol n U ρ = ∫ U ρ d x · · · d x n = ∮ ∂ U ρ x d x · · · d x n , where ∂ U ρ is the boundary of U ρ . Due to the regularity assumption ρ (cid:60) Σ , the gradient ∇ p f ρ does not vanish on the real zero locus of f ρ ,denoted V ( f ρ ) . Because U ρ is a union of connected componentsof A ∩ pr − ( ρ ) , it follows that ∂ U ρ is a compact ( n − ) -dimensionalsubmanifold of R n + contained in V ( f ρ ) .For ε >
0, let τ ( ρ ) be the Leray tube defined by τ ( ρ ) ≜ (cid:8) p + u ∇ p f ρ (cid:12)(cid:12) p ∈ ∂ U ρ , u ∈ C and | u | = ε (cid:9) . This is an n -dimensional submanifold of C n . We choose ε smallenough that τ ( ρ ) ∩ V ( f ρ ) = ∅ : this is possible because ∇ p f doesnot vanish on ∂ U ρ which is compact. Let R ( ρ , x , . . . , x n ) = x f − ρ ∂ f ρ / ∂ x ; observe that τ ( ρ ) doesnot cancel the denominator of R ( ρ , x , . . . , x n ) . Leray’s residue the-orem [27] shows that2 πi ∮ ∂ U t x d x · · · d x n = ∮ τ ( ρ ) d f ρ f ρ ∧ ( x d x · · · d x n ) = ∮ τ ( ρ ) R ( ρ , x , . . . , x n ) d x · · · d x n . (In Pham’s [33, Thm. III.2.4] notation, we have γ = ∂ U ρ , δγ = τ ( ρ ) , φ = f − ρ d f ρ ∧ ( x d x · · · d x n ) , and res [ φ ] = x d x · · · d x n .)To match the definition of a period and conclude the proof, it isenough to prove that, locally, the integration domain τ ( ρ ) can bemade independent of ρ . And indeed, since U is a union of connectedcomponents of A ∩ pr − ( I ) , we have ∂ U ⊆ f =
0. Therefore, since I is connected and I ∩ Σ = ∅ , the restriction of the projection prdefines a submersive map from ∂ U ∩ pr − ( I ) onto I . Additionally, ∂ U is compact, hence this map is proper. Ehresmann’s theorem thenimplies that there exists a continuous map h : I × ∂ U ρ → R n suchthat h ( σ , −) induces a homeomorphism ∂ U ρ ≃ ∂ U σ for any σ ∈ I .In particular, we have τ ( σ ) = (cid:110) h ( σ , p ) + u ∇ h ( σ , p ) f σ (cid:12)(cid:12)(cid:12) p ∈ ∂ U σ , u ∈ C and | u | = ε (cid:111) . This formulation makes it clear that τ ( σ ) deforms continuouslyinto τ ( ρ ) as σ varies. Since τ ( σ ) does not intersect the polar locus V ( f σ ) of R ( σ , −) , neither does τ ( ρ ) when σ and ρ are close enough,by compactness of τ ( σ ) and continuity of the deformation. There-fore, given any ρ ∈ I , we have ∮ τ ( σ ) R ( s , −) = ∮ τ ( ρ ) R ( σ , −) for σ close enough to ρ . □ The choice of x d x . . . d x n as a primitive of d x . . . d x n in The-orem 9 is arbitrary, but of little consequence, since the Picard-Fuchsequation only depends on the cohomology class of the integrand. Theorem 2 does not guarantee that v satisfies the Picard-Fuchsequation on the whole domain where the equation is nonsingular.It could happen that the solutions extend analytically across anexceptional point, or that some of them have singularities betweentwo consecutive exceptional points. As a consequence, we need toexplicitly compute Σ .Lemma 10. There exists an algorithm which, given on input apolynomial f ∈ Q [ t , x , . . . , x n ] of degree D satisfying (R) , computesa polynomial д ∈ Q [ t ] − { } of degree D O ( n ) whose set of real rootscontains Σ , using D O ( n ) operations in Q . Proof. Recall that, when (R) holds, the set Σ is finite. Our goalis to write Σ as the root set of a univariate polynomial д . Considerthe polynomial h = f + ( ∂ f / ∂ x ) + · · · + ( ∂ f / ∂ x n ) . We start bycomputing at least one point in each connected component of thereal algebraic set defined by h = D O ( n ) operations. It returns arational parametrization: polynomials P , F , G , . . . , G n in Q [ y ] ofdegree ⩽ D O ( n ) such that P is square-free and the set of points (cid:8) P ′ ( ξ ) − (cid:0) F ( ξ ) , G ( ξ ) , . . . , G n ( ξ ) (cid:1) ∈ R n + (cid:12)(cid:12) ξ ∈ R , P ( ξ ) = (cid:9) meets every connected component of the zero set of h . In particular, Σ = { F ( ξ )/ P ′ ( ξ ) | ξ ∈ R , P ( ξ ) = } . As a polynomial д , we take theresultant with respect to y of P ( y ) and F ( y ) − tP ′ ( y ) : its set of roots SSAC ’19, July 15–18, 2019, Beijing, China Pierre Lairez, Marc Mezzarobba, and Mohab Safey El Din contains Σ . Since P and F have degree D O ( n ) , this last step also uses D O ( n ) operations in Q [42]. □ Let us turn to the numerical part of the main algorithm. It isknown [9, 40] that Fuchsian differential equations with coefficientsin Q [ t ] can be solved numerically in quasi-linear time w.r.t. the pre-cision. Yet, some minor technical points must be addressed to applythe results of the literature to our setting. We start with reminderson the theory of linear ODEs in the complex domain [e.g. 17, 35].Consider a linear differential operator P = p m ( t ) d m d t m + · · · + p ( t ) dd t + p ( t ) (4)of order m with coefficients in Q [ t ] .Recall that u ∈ C is a singular point of P when the leadingcoefficient p m of P vanishes at u . A point that is not a singularpoint is called ordinary . Singular points are traditionally classifiedin two categories: a singular point u ∈ C is a regular singular pointof P if, for 0 ⩽ i < m , its multiplicity as a pole of p i / p m is atmost m − i , and an irregular singular point otherwise. The point atinfinity in P ( C ) is said to be ordinary, singular, etc., depending onthe nature of 0 after the change of variable t (cid:55)→ t − . An operatorwith no irregular singular point in P ( C ) is called Fuchsian .Fix a simply connected domain Ω ⊆ C containing only ordinarypoints of P , and let W be the space of analytic solutions y : Ω → C of the differential equation P ( y ) =
0. According to the Cauchyexistence theorem for linear analytic ODEs, W is a complex vectorspace of dimension m . A particular solution y ∈ W is determinedby the initial values y ( u ) , y ′ ( u ) , . . . , y ( m − ) ( u ) at any point u ∈ Ω .At a singular point, there may not be any nonzero analytic solu-tion. Yet, if u is a regular singular point, the differential equationstill admits m linearly independent solutions defined in the slit disk { u + ζ | | ζ | < η , ζ (cid:60) R − } for small enough η and each of the form y ( u + ζ ) = ζ γ ℓ (cid:213) k = y k ( ζ ) log ( ζ ) k = ℓ (cid:213) k = ∞ (cid:213) ν ∈ γ + N y k , ν ζ ν log ( ζ ) k (5)where γ ∈ ¯ Q , ℓ ∈ N , and y k , γ = y k ( ) (cid:44) k [35,§16]. The functions y k are analytic for | ζ | < η (including at 0). Thealgebraic numbers γ are called the exponents of P at u .Suppose now that u is either an ordinary point of P lying inthe topological closure ¯ Ω of Ω , or a regular singular point of P situated on the boundary of Ω . As a result of the previous discus-sion, we can choose a distinguished basis B u = ( ϕ u , , . . . , ϕ u , m ) of W in which each ϕ u , i is characterized by the leading monomial ( t − u ) γ log ( t − u ) k of its local expansion (5) at u . At an ordinarypoint u for instance, the coefficients of the decomposition of asolution y on B u are y ( i ) ( u )/ i !, that is, essentially the classical ini-tial values. Observe that when no two exponents γ have the sameimaginary part, the elements of B α all have distinct asymptoticbehaviours as t → u . In particular, at most one of them tends to anonzero finite limit. As Picard-Fuchs operators have real exponentsaccording to Theorem 7, this observation applies to them. More precisely, denoting λ k , ν ( y ) = y k , ν in (5), there are m computable pairs ( γ i , k i ) such that, for all i , we have λ ki , γi ( ϕ u , i ) = , λ kj , γj ( ϕ u , i ) = for j (cid:44) i ,and λ k , ν ( ϕ u , i ) = whenever ν − γ i (cid:60) N . Let u ′ ∈ ¯ Ω be a second point subject to the same restrictions as u .Let ∆ ( u , u ′ ) ∈ C m × m be the transformation matrix from B u to B u ′ .The key to the quasi-linear complexity of our algorithm is thatthe entries of this matrix can be computed efficiently, by solvingthe ODE with a Taylor method in which sums of Taylor series arecomputed by binary splitting [4, item 178], [9]. The exact resultwe require is due to van der Hoeven [40, Theorems 2.4 and 4.1];see also [29] for a detailed algorithm and some further refinements.Denote by M ( n ) the complexity of n -bit integer multiplication.Theorem 11 ([40]). For a fixed operator P and fixed algebraicnumbers u , u ′ as above, one can compute the matrix ∆ ( u , u ′ ) with anentry-wise error bounded by − p in O ( M ( p ( log p ) )) operations. Since P is linear, this result suffices to implement the procedureDSolve required by the main algorithm. More precisely, supposethat TransitionMatrix ( P , u , u ′ , p ) returns a matrix of complexintervals of width O ( − p ) that encloses ∆ ( u , u ′ ) entry-wise. Definition 12.
A system of good initial conditions for P on Ω ,denoted [ λ j ( y ) = s j ] m ′ i = , is a finite family of pairs ( λ j , s j ) where s j ∈ C and λ j is a linear form that belongs to the dual basis of B u for some algebraic point u ∈ ¯ Ω (which may depend on j ), with theproperty that λ , . . . , λ m ′ span the dual space of W .A system of good initial conditions on ( α , β ) ⊂ R is a system ofgood initial conditions on ( α , β ) + i ( , ε ) for some ε > P whose valuesdetermine at most one solution, and of prescribed values for thesecoefficients. When the system is compatible, we say that it defines the unique solution of P that satisfies all the constraints. Let us notein passing the following fact, which was used in §2.2.Lemma 13. Let u , . . . , u m ′ be ordinary points of P such that I = [ y ( u i ) = s i ] i is a system of good initial conditions for P on Ω , and let u ∈ ¯ Ω . Then I ′ = [ y ( u ) = s , y ′ ( u ) = s , . . . , y ′ ( u m ′ ) = s m ′ ] isa system of good initial conditions for P dd t on Ω . Proof. The derivative y (cid:55)→ y ′ maps the solution space of P dd t to that of P , and its kernel consists exactly of the constant functions.By assumption, a solution of P is completely defined by its values at u , . . . , u m , hence a solution of P dd t is characterized by the values ofits derivative at the same points, along with its limit at u . Because P dd t has order at least 2 (otherwise, I would not be a system of goodinitial conditions), the conditions y ′ ( u i ) are of the form λ ( y ) = s with λ belonging to the dual basis of some B u , as required. So is thecondition y ( u ) = s since P dd t has solutions with a nonzero finitelimit at u . □ Algorithm 3 evaluates the solution of an operator P given by asystem of good initial conditions. Note that the algorithm is allowedto fail. It fails if the intervals ˜ Λ i are not accurate enough for thelinear algebra step on line 10 to succeed, or if the linear system,which is in general over-determined, has no solution. The followingproposition assumes a large enough working precision p to ensurethat this does not happen. Additionally, we only require that theoutput be accurate to within O ( − p ) , so as to absorb any loss ofprecision resulting from numerical stability issues or from the useof interval arithmetic. omputing the volume of compact semi-algebraic sets ISSAC ’19, July 15–18, 2019, Beijing, China Algorithm 3
Solution of P ( y ) = procedure DSolve( P , [ ϕ ∗ u i , j i ( y ) = ˜ s i ] mi = , α , p ) ▷ ϕ ∗ u , j is the linear form dual to the element ϕ u , j of B u if B α has an element of leading monomial 1 then j ← its index else return u ← α ; ˜ ∆ ← id for ⩽ i ⩽ m do ▷ using interval arithmetic ˜ ∆ i ← TransitionMatrix ( P , u i − , u i , p ) · ˜ ∆ i − ˜ Λ i ← j i th row of ∆ i solve the linear system ˜ Λ i · ˜ c = ˜ s i , 1 ⩽ i ⩽ m (or fail) return the real part of ˜ c j Proposition 14.
Suppose that the operator P is Fuchsian withreal exponents. Let α < β be real algebraic numbers, and let y be areal analytic solution of P ( y ) = on the interval ( α , β ) such that y ( t ) tends to a finite limit as t → α . Let I = [ λ i ( y ) = s i ] be a system ofgood initial conditions for P on ( α , β ) that defines y .Given the operator P , the point α , a large enough working preci-sion p ∈ N , and an approximation ˜ I = [ λ i ( y ) = ˜ s i ] of I where ˜ s i is an interval of width at most − p containing s i , DSolve ( P , ˜ I , α , p ) (Algorithm 3) computes a real interval of width O ( − p ) containing lim t → α y ( t ) in time O ( M ( p ( log p ) )) . Proof. At the end of the loop, we have ∆ ( α , u i ) ∈ ˜ ∆ i , 1 ⩽ i ⩽ m ,and the entries of ˜ ∆ i are intervals of width O ( − p ) . The coefficients c = ( c i ) i of the decomposition of y in the basis B α satisfy Λ i · c = s i for all i , where Λ i is the j i th row of ∆ ( α , u i ) . As I is a system ofgood initial conditions, the linear system ( Λ i · x = s i ) i has no othersolution. Step 10 hence succeeds in solving the interval version assoon as the ˜ s i and the entries of the ˜ ∆ i are thin enough intervals.It then returns intervals of width O ( − p ) .We assumed that y tends to a finite limit at α . It follows that thedecomposition of y on B α only involves the basis elements with afinite limit at α . Either B α contains an element ϕ α , j that tends to 1,in which case lim α y = c j , or every solution that converges tendsto zero, and then the limit is zero. Since, by assumption, lim α y isreal, we can ignore the imaginary part of the computed value. Inboth cases, the algorithm, when it succeeds, returns a real intervalof width O ( − p ) containing lim α y .As for the complexity analysis, all u i including α are algebraic,hence Theorem 11 applies and shows that each call to TransitionMatrix runs in time O ( M ( p ( log p ) )) . The matrix multiplicationsat step 8 take O ( M ( p )) operations. The cost of solving the linearsystem (which is of bounded size) is O ( M ( p )) as well. The cost ofthe remaining steps is independent of p . □ It remains to show how to implement PickGoodPoints. Choos-ing the points at random works with probability one. The proceduredescribed below has the advantage of being deterministic and im-plying (at least in principle) bounds on the bit size of the u i .Lemma 15. Given P and two real numbers α < β , one can deter-ministically select m points u , . . . , u m ∈ ( α , β ) ∩ Q such that theevaluations y (cid:55)→ y ( u i ) are good initial conditions for P on ( α , β ) . Proof. A sufficient condition for y (cid:55)→ y ( u i ) to be good initialconditions is that the matrix M = ( ψ j ( u i )) i , j , for some basis ( ψ j ) of W , be invertible. Let K ⊂ ( α , β ) be a closed interval with rationalendpoints containing only ordinary points. Let u = min K . Assumewithout loss of generality u =
0, and take ( ψ j ) = B u . The matrix M is then of the form ( u j − i + η j ( u i )) mi , j = where, for all j , η j ( u ) = O ( u m ) as u →
0. In fact, there exists a computable [e.g., 40] constant C such that | η j ( u )| ⩽ C | u | m for all u ∈ K . Therefore, one can computea value ε > M is invertible for any distinct u , . . . , u m in ( , ε ) . The result follows. □ In practice, one can reduce the number of recursive calls in themain algorithm by replacing, when possible, some of the conditions y ( u i ) = s i by conditions that result from the continuity of v ( t ) atexceptional points, or from its analyticity at singular points of thePicard-Fuchs operator lying in R \ Σ . For instance, a solution thatis analytic at u must lie in the subspace spanned by the elementsof B u of leading term ( z − u ) γ with γ ∈ N and no logarithmic part. Let us finally study the complexity of Algorithm 2 to conclude theproof of Theorem 1. For fixed ( f , . . . , f r ) , all intermediate data(Picard-Fuchs equations, critical values and specialization pointschosen for the recursive calls) are fixed thanks to the deterministicbehaviour of PickGoodPoints (Lemma 15). Thus, the number ofrecursive calls does not depend on p .Now, the main point is to observe that, by Proposition 3, perform-ing recursive calls with precision p + O ( ) is enough. One can makethe width of the output interval smaller than 2 − p by doubling p andre-running the algorithm (if necessary with a more accurate ap-proximation of I ) a bounded number of times. By Proposition 14,the total cost of the calls to DSolve is O ( M ( p log ( p ) )) . The onlyother step whose complexity depends on p is the computation ofreal roots of fixed univariate polynomials in the base case, whichtakes O ( M ( p )) operations using Newton’s method. Using the bound M ( p ) = O ( p log ( p ) + ε ) , Theorem 1 follows.This theorem ignores the dependency of the cost on the dimen-sion n of the ambient space or the maximum degree D of the inputpolynomials. Under some assumptions, one can bound the num-ber of recursive calls arithmetic cost of computing Picard-Fuchsequations and critical values as follows. First consider Algorithm 1,and let δ be the degree of f . By Lemma 10, the number of criticalvalues and the cost of computing them are bounded by δ O ( n ) ; inthe notation of the algorithm, this shows that ℓ ⩽ δ O ( n ) .Under Dimca’s conjecture [12], the cost of computing the Picard-Fuchs equation is δ O ( n ) and it has order m ⩽ δ n according to thediscussion following Theorem 8. One can likely obtain the samebounds without this conjecture by replacing the deformed equationof Section 2.2 by f − t (cid:205) i x δ + ni , which permits using the “regularcase” of [5]. Solving the recurrence C ( n + , δ ) = δ O ( n ) C ( n , δ ) showsthat the algebraic steps of Algorithm 1 take δ O ( n ) operations in Q .Turning to Algorithm 2, Lemma 10 and Theorem 8 show thatthe cost of the calls to CriticalValues and PicardFuchs are domi-nated by that of the calls to Algorithm 1 (with an input polynomialof degree δ ⩽ rD ). Therefore, the algebraic steps use ( rD ) O ( n ) operations in Q in total, as announced in §1. SSAC ’19, July 15–18, 2019, Beijing, China Pierre Lairez, Marc Mezzarobba, and Mohab Safey El Din
We leave for future research the question of analyzing the booleancost of the full algorithm with respect to n , D , and the bit size ofthe input coefficients. This requires significantly more work, asone first needs to control the bit size of the points picked by Pick-GoodPoints in the recursive calls. Additionally, to the best of ourknowledge, no analogue of Theorem 11 fully taking into accountthe order, degree, and coefficient size of the operator P is availablein the literature. Our algorithm generalizes to non-basic bounded semi-algebraicsets since their volume can be written as a linear combination with ± basic semi-algebraic sets.An important question that we leave for future work is that ofthe practicality of our approach. While the worst-case complexitybound is exponential in n , there are a number of opportunitiesto exploit special features of the input that could help handlingnontrivial examples in practice. In particular: (1) the number ofrecursive calls only depends on the number of real critical points;(2) as already noted, it can be reduced by exploiting some knowledgeof the continuity of the slice volume function or its analyticity atexceptional points; (3) it turns out that, in our case, the integralappearing in Theorem 9 always is singular at infinity, and, as aconsequence, the Picard-Fuchs equations we encounter do not reachthe worst-case degree bounds. Ideally, one may hope to refine thecomplexity analysis to reflect some of these observations.Another natural question is to extend the algorithm to unboundedsemi-algebraic sets of finite volume, or even real periods in general,using the ideas in [41]. Note also that, using quantifier elimination[e.g., 2], boundedness can be verified in boolean time q ( rD ) O ( n ) where q bounds the bit size of the input coefficients.Finally, it is plausible that an algorithm of a similar structure butusing numerical quadrature recursively instead of solving Picard-Fuchs equations would also have polynomial complexity in theprecision for fixed n and be faster at medium precision. REFERENCES [1] D. Bailey and J. M. Borwein. 2011. High-precision numerical integration: progressand challenges.
Journal of Symbolic Computation
46, 7 (2011), 741–754.[2] S. Basu, R. Pollack, and M.-F. Roy. 1996. On the combinatorial and algebraiccomplexity of quantifier elimination.
J. ACM
43, 6 (1996), 1002–1045.[3] S. Basu, R. Pollack, and M.-F. Roy. 2006.
Algorithms in real algebraic geometry (second ed.). Algorithms and Computation in Mathematics, Vol. 10. Springer.[4] M. Beeler, R. W. Gosper, and R. Schroeppel. 1972.
Hakmem . AI Memo 239. MITArtificial Intelligence Laboratory.[5] A. Bostan, P. Lairez, and B. Salvy. 2013. Creative telescoping for rational functionsusing the Griffiths–Dwork method. In
ISSAC 2013 . ACM, 93–100.[6] P. Bürgisser, F. Cucker, and P. Lairez. 2019. Computing the homology of basicsemialgebraic sets in weak exponential time.
J. ACM
66, 1 (2019), 5.[7] J. Canny. 1988.
The complexity of robot motion planning . MIT press.[8] G. Christol. 1985. Diagonales de fractions rationnelles et équations de Picard-Fuchs.
Groupe de travail d’analyse ultramétrique
12, 1 (1985).[9] D. V. Chudnovsky and G. V. Chudnovsky. 1990. Computer algebra in the service ofmathematical physics and number theory. In
Computers in mathematics (Stanford,CA, 1986) . Lecture Notes in Pure and Appl. Math., Vol. 125. Dekker, 109–232.[10] F. Chyzak. 2000. An extension of Zeilberger’s fast algorithm to general holonomicfunctions.
Discrete Mathematics
Automata Theory and Formal Languages 2nd GI ConferenceKaiserslautern, May 20–23, 1975 . Springer, 134–183.[12] A. Dimca. 1991. On the de Rham cohomology of a hypersurface complement.
American Journal of Mathematics
SIAM J. Comput.
17, 5 (1988), 967–974.[14] M. Dyer, A. Frieze, and R. Kannan. 1991. A random polynomial-time algorithmfor approximating the volume of convex bodies.
J. ACM
38, 1 (1991), 1–17.[15] A. Grothendieck. 1966. On the de Rham cohomology of algebraic varieties.
Institutdes Hautes Études Scientifiques. Publications Mathématiques
29 (1966), 95–103.[16] D. Henrion, J.-B. Lasserre, and C. Savorgnan. 2009. Approximate volume andintegration for basic semialgebraic sets.
SIAM review
51, 4 (2009), 722–743.[17] E. Hille. 1976.
Ordinary differential equations in the complex domain . Wiley.Dover reprint, 1997.[18] N. M. Katz. 1971. The regularity theorem in algebraic geometry. In
Actes duCongrès International des Mathématiciens (Nice, 1970), Tome 1 . Gauthier-Villars,437–443.[19] L. Khachiyan. 1989. The problem of computing the volume of polytopes isNP-hard.
Uspekhi Mat. Nauk
44, 3 (1989), 199–200.[20] L. Khachiyan. 1993. Complexity of polytope volume computation. In
New trendsin discrete and computational geometry . Springer, 91–101.[21] L. Khachiyan and L. Porkolab. 2000. Integer Optimization on Convex Semialge-braic Sets.
Discrete & Computational Geometry
23, 2 (Feb. 2000), 207–224.[22] P. Koiran. 1995. Approximating the volume of definable sets. In
Foundations ofComputer Science . IEEE, 134–141.[23] M. Kontsevich and D. Zagier. 2001. Periods. In
Mathematics unlimited . Springer,771–808.[24] M. Korda and D. Henrion. 2018. Convergence rates of moment-sum-of-squareshierarchies for volume approximation of semialgebraic sets.
Optimization Letters
12, 3 (2018), 435–442.[25] C. Koutschan. 2010. A fast approach to creative telescoping.
Mathematics inComputer Science
4, 2-3 (2010), 259–266.[26] P. Lairez. 2016. Computing periods of rational integrals.
Math. Comp.
85, 300(2016), 1719–1752.[27] J. Leray. 1959. Le calcul différentiel et intégral sur une variété analytique complexe(Problème de Cauchy, III).
Bulletin de la Société mathématique de France
87 (1959),81–180.[28] L. Lipshitz. 1988. The diagonal of a D-finite power series is D-finite.
Journal ofAlgebra
Autour de l’évaluation numérique des fonctions D-finies .Thèse de doctorat. École polytechnique.[30] M. Mezzarobba. 2016. Rigorous multiple-precision evaluation of D-finite func-tions in SageMath. (2016). arXiv:1607.01967[31] P. Monsky. 1972. Finiteness of de Rham cohomology.
American Journal ofMathematics
94 (1972), 237–245.[32] T. Oaku. 2013. Algorithms for integrals of holonomic functions over domainsdefined by polynomial inequalities.
Journal of Symbolic Computation
50 (2013),1–27.[33] F. Pham. 2011.
Singularities of integrals: homology, hyperfunctions and microlocalanalysis . Springer ; EDP Sciences.[34] E. Picard. 1902. Sur les périodes des intégrales doubles et sur une classed’équations différentielles linéaires. In
Comptes rendus hebdomadaires des séancesde l’Académie des sciences , Vol. 134. Gauthier-Villars, 69–71.[35] E. G. C. Poole. 1936.
Introduction to the theory of linear differential equations .Clarendon Press.[36] M. Safey El Din and É. Schost. 2003. Polar varieties and computation of one pointin each connected component of a smooth real algebraic set. In
ISSAC 2003 . ACM,224–231.[37] M. Safey El Din and É. Schost. 2017. A nearly optimal algorithm for decidingconnectivity queries in smooth and bounded real algebraic sets.
J. ACM
63, 6(2017), 48.[38] M. Safey El Din and L. Zhi. 2010. Computing Rational Points in Convex Semial-gebraic Sets and Sum of Squares Decompositions.
SIAM Journal on Optimization
20, 6 (Jan. 2010), 2876–2889.[39] A. Tarski. 1998. A decision method for elementary algebra and geometry. In
Quantifier elimination and cylindrical algebraic decomposition . Springer, 24–84.[40] J. van der Hoeven. 2001. Fast Evaluation of Holonomic Functions Near and inRegular Singularities.
Journal of Symbolic Computation
31, 6 (2001), 717–743.[41] J. Viu-Sos. 2015-09-03. A semi-canonical reduction for periods of Kontsevich-Zagier. (2015-09-03). arXiv:1509.01097[42] J. von zur Gathen and J. Gerhard. 1999.