Efficient Algorithms for Mixed Creative Telescoping
aa r X i v : . [ c s . S C ] M a y Efficient Algorithms for Mixed Creative Telescoping
Alin Bostan
InriaFrance
[email protected] Louis Dumont
InriaFrance
[email protected] Bruno Salvy
Inria, Laboratoire LIP(U. Lyon, CNRS, ENS Lyon, UCBL)France
ABSTRACT
Creative telescoping is a powerful computer algebra paradigm–initiated by Doron Zeilberger in the 90’s– for dealing withdefinite integrals and sums with parameters. We addressthe mixed continuous–discrete case, and focus on the inte-gration of bivariate hypergeometric-hyperexponential terms.We design a new creative telescoping algorithm operatingon this class of inputs, based on a Hermite-like reductionprocedure. The new algorithm has two nice features: it isefficient and it delivers, for a suitable representation of theinput, a minimal-order telescoper. Its analysis reveals tightbounds on the sizes of the telescoper it produces.
CCS Concepts • Computing methodologies → Algebraic algorithms;
Keywords symbolic integration; creative telescoping; hypergeometric-hyperexponential term; Hermite reduction
1. INTRODUCTION
Context.
Creative telescoping is an algorithmic approachintroduced in computer algebra by Zeilberger [23, 24, 21] toaddress definite summation and integration for a large classof functions and sequences involving parameters.In this article, we focus on the mixed continuous–discretecase. Given a term F n ( x ) that is both hypergeometric (i.e., F n +1 ( x ) /F n ( x ) is a rational function) and hyperexponential(i.e., F ′ n ( x ) /F n ( x ) is a rational function), the question is tofind a linear recurrence relation satisfied by the sequence ofintegrals I n = R γ F n ( x ) dx over a domain γ where F n ( x ) isintegrable. To do this, the method of creative telescopinglooks for polynomials c ( n ) , . . . , c r ( n ), not all zero, and fora rational function Q ( n, x ) such that G n ( x ) = Q ( n, x ) F n ( x ) ACM ISBN 978-1-4503-2138-9.DOI: http://dx.doi.org/10.1145/2930889.2930907 satisfies the telescoping relation L ( F n ( x )) def = r X i =0 c i ( n ) F n + i ( x ) = G ′ n ( x ) . (1)The recurrence operator L = P ri =0 c i ( n ) S in in the shift op-erator S n is called a telescoper for F n ( x ), and the rationalfunction Q ( n, x ) is called a certificate for the telescoper L .The integer r is the order of L and max j deg c j is its degree .Throughout the article, the ground field, denoted by k , isassumed to be of characteristic zero. Under suitable addi-tional assumptions, L ( I n ) = 0 is a recurrence relation satis-fied by the sequence of integrals I n = R γ F n ( x ) dx .A subtle point is that not all hypergeometric-hyperexponentialterms admit telescopers. A criterion for deciding when thisis the case has been given only recently [7, Section 6]: ahypergeometric-hyperexponential term is telescopable if andonly if it can be written as the sum of a derivative and of a proper term . In our context, proper terms are of the form P ( n, x ) · H ( x ) n · exp (cid:18)Z S ( x ) T ( x ) (cid:19) · Υ( n ) , (2)for P ∈ k [ n, x ], H ∈ k ( x ), S, T ∈ k [ x ] and Υ hypergeomet-ric.Several methods are known for computing a telescoper L and the corresponding certificate Q , even in much greatergenerality [1, 21, 11, 20, 2, 15]. Despite a very rich researchactivity during the last decade, not much is known about thecomplexity of creative telescoping methods when applied toan input of the type (2). In particular, few estimates areavailable in the literature on the order and degree of theminimal-order telescoper.If Υ( n ) = 1, one can compute a telescoper of the form (1)for the rest and then multiply the coefficient c i ( n ) by therational function Υ( n ) / Υ( n + i ); this process does not affectthe minimality of the telescoper. Thus we focus on the case F n ( x ) = P ( n, x ) · H ( x ) n · exp (cid:18)Z S ( x ) T ( x ) (cid:19) . (3) Previous work.
Among the various classes of creativetelescoping methods, three different approaches allow to ob-tain bounds on the sizes of the telescoper and the certifi-cate, and to control the algorithmic complexity. The firstapproach is based on elimination techniques [23, 21, 22]; itgenerally yields pessimistic bounds, not very efficient algo-rithms and telescopers of non-minimal order. The secondapproach, initiated by Apagodu and Zeilberger [2], providessharp bounds on the order and degree of telescopers, more ef-cient algorithms, but does not provide information on tele-scopers of minimal order and works under a restrictive gener-icity assumption. It was successfully applied by Kauers andco-authors [10, 9, 14] to (bivariate) hyperexponential termsand hypergeometric terms. The third approach, based onHermite-like reduction, is the only one that computes tele-scopers of minimal order, while guaranteeing a good controlon sizes and complexity. It has been introduced by Bostan et alii for the integration of bivariate rational functions [3],then extended to bivariate hyperexponential functions [4],multivariate rational functions [5, 17] and recently adaptedto summation [8] for bivariate hypergeometric terms. Thepresent work is part of the on-going effort in this direction.
Contributions.
We present the first Hermite-style algo-rithm in the mixed (continuous–discrete) setting. Our ap-proach is inspired by the proof of Manivel’s lemma [12, § polyno-mial rigidity conjecture .Our algorithm works on terms of the form (3). Its inputis F n ( x ) = P ( n, x )Φ( n, x ), given by P ( n, x ) , Φ ′ Φ ( n, x ) = A ( n, x ) B ( x ) = n H ′ ( x ) H ( x ) + S ( x ) T ( x ) , ( ⋆ )i.e., a polynomial in k [ n, x ] and a rational function in k ( n, x )of a special form given either in lowest terms ( A/B withgcd(
A, B ) = 1 and deg n B = 0) or decomposed as the sumof the logarithmic derivative of a rational function H ∈ k ( x )multiplied by n and another rational function S/T ∈ k ( x ).Note that, for a given term, several choices are available for P and Φ. If we further assume that Φ ′ / Φ has no positiveinteger residue, they become unique and we call them the minimal decomposition of F n ( x ).Our main results consist in bounds on the order and de-gree of telescopers for F n ( x ), that are summarized in Theo-rem 1 below (see Theorems 6 and 15 for more precise state-ments). A complexity analysis of the algorithms leading tothese bounds is conducted in section 4 (see in particularTheorem 21).In the statement of the following theorem, deg x H and deg n P denote the maximum of the degrees of their numerator anddenominator, and d H is the degree of H at infinity. Theorem 1
Given a decomposition as in ( ⋆ ) , F n ( x ) admitsa telescoper of order r bounded by δ , where δ := max(deg x A, deg x B − , (4) and degree bounded by r (deg x H + deg n P ) + deg x P. Moreover, if the decomposition is minimal then these boundsapply to a minimal-order telescoper for F n ( x ) .Algorithm MixedCT ( § d denotes an upper bound on the degrees of thenumerator and denominator of H , and of all the polynomialsin ( ⋆ ) , and if all these polynomials are square-free, then thetelescoper has arithmetic size O ( d ) and MixedCT computesit using at most ( ˜ O ( d ω +1 ) if d H < , ˜ O ( d ) if d H > arithmetic operations in k , where ω denotes a feasible matrixmultiplication exponent for k . Note that one can always choose the decomposition of F n ( x )in ( ⋆ ) to be minimal. Indeed, we may write Φ = Q ˜Φ where Q is a polynomial and ˜Φ ′ / ˜Φ has no positive integer residue.We then have a minimal decomposition F n ( x ) = ( P Q ) ˜Φ.The proof of Theorem 1 is achieved through two mainingredients: confinement and
Hermite reduction . Confine-ment is the property that given Φ, any polynomial P canbe reduced modulo derivatives to a polynomial R of degreeat most δ − P Φ = R Φ + Γ ′ for some Γ. It gives a finitedimensional vector space over k ( n ), where the computationwill be confined. Hermite reduction is more classical. In thiscontext, it consists in an algorithm that performs a reduc-tion P Φ( n + 1 , x ) = R Φ( n, x ) + Γ ′ for some Γ. By iteratedapplication of both these operations, all P ( n + i, x )Φ( n + i, x )for i = 0 , . . . , δ can be rewritten in a vector space of dimen-sion δ over k ( n ). Thus by (polynomial) linear algebra thereexists a non-trivial linear relation between them, i.e., a tele-scoper. A more careful study of the increase of the degreesin n during these rewritings gives the degree bound. Notation.
In all that follows, k and K will denote fieldsof characteristic 0. In our applications we will often set K = k ( n ). We denote by K [ x ] d the set of polynomials in K [ x ] of degree less than d .Rational functions are always written in reduced form,with monic denominator. Thus the numerator and denomi-nator are defined without ambiguity. If k is the degree of thenumerator and ℓ the degree of the denominator of a rationalfunction F , we say that F has rational degree ( k, ℓ ), that wedenote Rdeg( F ) = [ k ] / [ ℓ ]. We also define the regular degree deg( F ) of F as deg( F ) = max( k, ℓ ). Finally, the degree atinfinity of F is defined as deg ∞ ( F ) = k − ℓ . The variablewith respect to which the degree is taken will be indicatedas a subscript when there is an ambiguity.A polynomial is called square-free when its gcd with itsderivative is trivial. The square-free decomposition of amonic polynomial Q ∈ K [ x ] is a factorization Q = Q · · · Q mm ,with Q i ∈ K [ x ] monic and square-free, the Q i ’s pairwise co-prime and deg x ( Q m ) >
0. The square-free part of Q is thepolynomial Q ⋆ = Q · · · Q m . Structure of the article.
Section 2 gives the main prop-erties of the confinement and of the mixed Hermite-like re-duction, leading to the bound on the order of the telescoper.Section 3 gives a bound on the degree of the telescoper andanalyzes the evolution of the degrees during the reductions,preparing the complexity analyses in Section 4. We con-clude the article in Section 5 with a few applications of ourimplementation and experiments on the actual growth of theminimal-order telescopers.
2. ALGORITHMS AND ORDER BOUND
In this section, we introduce the algorithms with justenough information to prove their correctness and to obtaina bound on the order of the telescoper they compute. Amore thorough analysis of the degrees is in the next section.
In terms of integrals, the operation of confinement writes Z γ P ( n, x )Φ( n, x ) dx = Z γ R ( n, x )Φ( n, x ) dx, with R a polynomial of degree smaller than δ from Eq (4).This transformation is based on the following lemma. emma 2 Let Φ , A and B be as in Eq. ( ⋆ ) , and δ as inEq. (4) . Then for any polynomial P ∈ K [ x ] , there existunique polynomials R and Q in K [ x ] with deg R δ − and deg Q deg P − δ such that P Φ = R Φ + ( QB Φ) ′ . (5)Thus all polynomial multiples of Φ can be written moduloderivatives on a vector space of dimension δ . Proof.
Equation (5) rewrites QA + ( QB ) ′ = P − R. (6)Denote d = deg( P ) − δ , and consider the linear map f : K [ x ] d +1 → K [ x ] d +1 Q ( QA + ( QB ) ′ ) div x δ , where u div v denotes the quotient in the Euclidean divisionof u by v . For m d , f ( x m ) has degree at most m . Itscoefficient of degree m equals either lc A = 0 if δ = deg A > deg B −
1, or m + δ + 1 = 0 if δ = deg B − > deg A , orlc A + m + δ +1 if deg A = deg B −
1. For this case to occur, itis necessary that deg
S < deg T . In that case, write H ′ /H = U/V with
U, V polynomials such that deg U = deg V − A/B = ( nUT + V S ) / ( V T ), anddeg(
V S ) deg V + deg T − UT ) . From this we see that lc( A ) depends on n and we deducethat f ( x m ) has degree exactly m in all cases.Thus f is an isomorphism. It follows that Equation (6)is equivalent to Q being the unique polynomial such that f ( Q ) = P div x δ , and then R is P − QA − ( QB ) ′ . The degreesof Q and R directly follow from the construction.Algorithm Confinement implements the proof of Lemma 2.Its correctness follows from the fact that the relation usedin the loop is obtained by extracting the coefficient of x i + δ in Eq. (6).The key to the minimality in Theorem 1 is the followingproperty of the confinement. Proposition 3
With the notation of ( ⋆ ) , further assumethat Φ ′ / Φ has no positive integer residue. Then, for anypolynomial R ∈ K [ x ] such that deg x R < δ ∃ K ∈ K ( x ) R Φ = ( K Φ) ′ ⇔ R = 0 . Proof.
Only the direct implication is not obvious. Ifsuch a K exists, the equation rewrites R = K ′ + K AB . If K is a polynomial, this equality can only be satisfied if B divides K , in which case the result is a direct consequence ofthe uniqueness in Lemma 2. Now assume that K has a pole x of order v >
0. Then the equation shows that
A/B musthave a simple pole at x with residue v , which contradictsthe assumption on the residues of Φ ′ / Φ. In terms of integrals, our Algorithm
HermiteReduction letsone change H n +1 into H n in the integral, writing Z γ P ( n, x ) H ( x )Φ( n, x ) dx = Z γ ˜ P ( n, x )Φ( n, x ) dx, for some polynomial ˜ P . It relies on a sequence of elementarysteps ( BasicReduction ) based on the following lemma. Algorithm
Confinement ( P , F ) Input
A polynomial P ∈ K [ x ], a rational function F = A/B with gcd(
A, B ) = 1.
Output
A polynomial R ∈ K [ x ] of degree less thanmax(deg( A ) , deg( B ) −
1) such that P = R + ( QB ) ′ + QA for some Q ∈ K [ x ]. δ ← max(deg( A ) , deg( B ) − d ← deg( P ) − δ ;Write A = P i a i x i , B = P i b i x i , P = P i p i x i ; for i ← d to 0 by − do c ← a δ + ( δ + i + 1) b δ +1 ; q i ← c (cid:16) p δ + i − P δj =1 q i + j a δ − j − ( δ + i + 1) P δ +1 j =1 q i + j b δ +1 − j (cid:17) ; Q ← P di =0 q i x i ; return P − ( QB ) ′ − QA . Lemma 4
Let Φ , A and B be as in Eq. ( ⋆ ) . Then for any G dividing B and satisfying gcd( G, A + B ′ ) = 1 , there existpolynomials Q and R in K [ x ] such that P Φ = RG Φ + ( QB Φ) ′ . (7) Proof.
The hypothesis gcd(
G, A + B ′ ) = 1 implies theexistence of polynomials U, Q ∈ K [ x ] such that P = UG + Q ( A + B ′ ). Then the derivative ( QB Φ) ′ expands as( QB Φ) ′ Φ = Q ′ B + Q ( A + B ′ ) = ( Q ′ B/G − U ) G + P, which has exactly the form of Equation (7).The crucial condition gcd( G, A + B ′ ) = 1 required to applythis lemma does not hold for arbitrary A and B and divisor G of B , but when G is square-free, it is a consequence of thepresence of n in Eq. ( ⋆ ), as shown in the following. Lemma 5
Let Φ , A and B be as in Eq. ( ⋆ ) . Then for anysquare-free polynomial G in k [ x ] dividing the denominatorof H , gcd( G, A + B ′ ) = 1 . This is also true if A/B is replacedby the reduced form of
A/B + iG ′ /G for some i ∈ Z . Proof. If H ′ = 0 then the denominator of H is 1 andthen G = 1 and the property holds.Otherwise, let first G be an irreducible factor of the de-nominator of H , so that there exist an integer k and poly-Algorithm BasicReduction ( P , F , G , k ) Input
A polynomial P ∈ K [ x ], a rational function F = A/B ∈ K ( x ), a square-free factor G of B s.t.gcd( G, A + B ′ + iBG ′ /G ) = 1 for all i ∈ Z , a positiveinteger k . Output
A polynomial R ∈ K [ x ] such that P = G k ( R + q ′ + F q ) for some q ∈ K [ x ][ G − ]. R ← P ; for i ← k do C ← A + ( i − k − BG ′ /G + B ′ ;Write R = QC + V G with deg
Q < deg G ; R ← ( R − Q ′ B − QC ) /G ; ⊲ PG k Φ = RG k − i Φ + ( q Φ) ′ return R .lgorithm HermiteReduction ( P , H , S/T ) Input
A polynomial P ∈ k ( n )[ x ],two rational functions H and S/T in k ( x ). Output
A polynomial R ∈ k ( n )[ x ] such that P = R/H + nQH ′ /H + QS/T + Q ′ for some Q ∈ k ( n, x ).Compute the square-free decomposition of thedenominator of H : g = g g . . . g mm ;Compute the corresponding partial fractiondecomposition: H = U + P mk =1 U k /g kk ; R ← P U ; K ← S/T + ( n − H ′ /H ; for k ← m do r k ← BasicReduction ( P U k , K, g k , k ); return R + r + · · · + r m .nomials H , H ∈ k [ x ] such that AB = nk G ′ G + n H H + ST , with gcd(
G, H ) = 1 and B = lcm( G, H , T ). Write B = G ν ˜ B with gcd( G, ˜ B ) = 1. Then A = nkG ′ G ν − ˜ B + nH G ν ( ˜ B/H ) + S ( G ν ˜ B/T ) . Reducing A + B ′ modulo G then yields A + B ′ ≡ ( nk + ν ) G ′ G ν − ˜ B + S ( G ν ˜ B/T ) mod G. (8)Since G does not depend on n , G | A + B ′ would imply thatboth G | G ′ G ν − ˜ B and G | SG ν ˜ B/T . The first condition im-plies ν >
1, which forces that G ν | T , making G | S necessarytoo, a contradiction. This reasoning also proves the resultwhen adding integer multiples of G ′ /G to the fraction A/B ,which adds an integer to nk + ν in Eq. (8).If G is only assumed to be a square-free divisor of thedenominator of H , then the property holds for each of itsirreducible factors, and thus A + B ′ is invertible modulo theirproduct G by the Chinese remainder theorem. Correctness of
BasicReduction and
HermiteReduction.
Al-gorithm
HermiteReduction treats each square-free factor ofthe denominator of H separately, while the second part ofthe lemma is used in Algorithm BasicReduction to deal withmultiplicities. If G is a square-free factor of multiplicity k , then Lemma 4 is used successively with A/B the reducedform of the logarithmic derivatives of Φ /G k , Φ /G k − , . . . , Φ /G ,thus rewriting P Φ /G k as R Φ up to a derivative.
Combining confinement and Hermite reduction gives thefinal result.
Theorem 6
Let P be a polynomial in k [ n, x ] and Φ , A and B as in Eq. ( ⋆ ) and δ as in Eq. (4) . Then, F n = P ( n, x )Φ admits a telescoper of order bounded by δ . Proof.
By Lemma 4, Algorithm
HermiteReduction canbe used to rewrite all the shifts F n , F n +1 , F n +2 , . . . underthe form R Φ modulo derivatives, with R ∈ k ( n )[ x ]. The nec-essary condition to apply the algorithm is satisfied at eachstep thanks to Lemma 5. By Lemma 2, F n , F n +1 , . . . , F n + δ Algorithm
MixedCT ( P, H, S/T ) Input
A polynomial P ∈ k [ n, x ],two rational functions H and S/T in k ( x ). Output A r -tuple ( c , . . . , c r − ) such that P ( n + r, x ) H r − P r − i =0 c i P ( n + i, x ) H i = nQH ′ /H + QS/T + Q ′ for some Q ∈ k ( n, x ). F ← S/T + ( n − H ′ /H ; R ← Confinement ( P | n n − , F ); for k ← , . . . doif rank k ( n ) ( R , R , . . . , R k ) < k + 1 then Solve P k − i =0 c i R i = R k for c , . . . , c k − in k ( n ); return ( c , . . . , c k − ). P ← R k | n n +1 ; P ← HermiteReduction ( P, H, S/T ); R k +1 ← Confinement ( P, F );are linearly dependent modulo derivatives. A linear relationbetween them provides a telescoper of order at most δ .Algorithm MixedCT implements that proof. In practice,for efficiency purposes, the reduction of F n + i is obtainedby applying HermiteReduction to the shift of the confinedreduction of F n + i − . Theorem 7
With the notation of Theorem 6, further as-sume that Φ ′ / Φ does not have any positive integer residue.Then Algorithm MixedCT ( P , H , S/T ) computes a minimal-order telescoper for P Φ . Proof.
Consider the minimal-order monic telescoper L ( n, S n ) = S rn − r − X i =0 c i ( n ) S in of P Φ and its certificate C ∈ k ( n, x ) such that L ( n, S n ) · ( P Φ) = ( C Φ) ′ . By Lemma 2, there exist R and K satisfying L ( n, S n ) · ( P Φ) = R Φ + ( K Φ) ′ , where K ∈ k ( n, x ) and R = R r − r − X i =0 c i ( n ) R i is a linear combination of the R i ’s of Algorithm MixedCT . Itfollows that R Φ = (( C − K )Φ) ′ and Proposition 3 then im-plies that R = 0. Thus, this linear combination is detectedby the algorithm, producing the output ( c , . . . , c r − ). Certificates.
Algorithm
MixedCT as given here computesthe certificate, although not in a normalized form. We choseto only output the telescoper, but it would be possible toreturn the certificate as well (and normalize it or not).
3. DEGREE BOUNDS
We now review more precisely the algorithms and obtainbounds on the degrees at each step. The first part of thissection consists of technical results that are needed for thecomplexity analysis in the next section. Then, at the end ofsection 3.3, we give a bound on the degree of the telescoperproduced by Algorithm
MixedCT . .1 Confinement Lemma 8
Let A and B be as in Eq. ( ⋆ ) and let P be apolynomial in k ( n )[ x ] . Let R be the polynomial returned byAlgorithm Confinement . Then
Rdeg n R − Rdeg n P is at most ( (deg x P − deg x A + 1) [1][1] , if deg x B deg x A + 1 , hj deg x P − deg x B +1deg x B − deg x A − k + 1 i / [0] otherwise. Proof.
The proof is a case by case analysis of Algo-rithm
Confinement ; we use its notation.When δ = deg x A , the recurrence for q i has a summand a δ − q i +1 except when i = d , while c has a δ for summand.Thus by induction, the degree in n of the numerator anddenominator of q d − i increase by 1 at each step. Since thereare d +1 steps, Rdeg n Q − Rdeg n P is bounded by [ d ] / [ d +1].The result for R then follows from Equation (6).When δ = deg x B −
1, the coefficients a δ − j are zero for j <δ − deg x A and c has degree 0 in n . By induction on i , q d − i has degree that changes (by increases of 1) only when i ≡ δ − deg x A . More explicitly, Rdeg n q d − i − Rdeg n P isbounded by [ ⌊ i/ ( δ − deg x A ) ⌋ ] / [0]. Again, the conclusionfor R follows from Equation (6). In order to track the degrees in n of the polynomials in-volved in Algorithms BasicReduction and
HermiteReduction ,we need to look deeper into the modular inversions involved.This is done in the next lemma, using the same notation asin the discussion preceding Lemma 5.
Lemma 9
Let A , B and H be as in Eq. ( ⋆ ) . Let G be anirreducible factor of the denominator of H , and ν the G -adic valuation of B . Let P be in k ( n )[ x ] and let Q be thepolynomial such that deg x Q < deg x ( G ) and P ≡ Q ( A + B ′ ) mod G . Then Rdeg n Q − Rdeg n P is bounded by [0] / [0] , if ν > / [1] , if ν = 1 and G T ;[deg x ( G ) − / [deg x ( G )] , if ν = 1 and G | T . Proof.
The existence of Q follows from Lemma 5. Noticefirst that Rdeg n Q − Rdeg n P = Rdeg n ( A + B ′ ) − (where( A + B ) − denotes an inverse mod G ). Indeed, writing p ( n ) P = P ( x ) + · · · + P d ( x ) n d with p ( n ) the denomina-tor of P , we see that p ( n ) Q = ( A + B ′ ) − ( P mod C ) + · · · +( A + B ′ ) − n d ( P d mod C ) has rational degree in n at most[ d ] / [0] + Rdeg n ( A + B ′ ) − .Thus we just need to bound Rdeg n ( A + B ′ ) − . To do so,we take a closer look at Equation (8). If ν >
1, the equationbecomes A + B ′ ≡ S ˜ BT /G ν mod G, so that ( A + B ′ ) − does not depend on n in this case.If ν = 1 and G T , the equation becomes A + B ′ ≡ ( nk + 1) G ′ ˜ B. Then, ( A + B ′ ) − ≡ ( G ′ ˜ B ) − / ( nk + ν ) has rational degree[0] / [1]. Finally, if ν = 1 and G | T , the result follows fromLemma 10 below. Lemma 10
Let
A, B, C be polynomials in k [ x ] such that C is relatively prime with at least one of A or B . Let U, V ∈ k ( n )[ x ] be such that U ( An + B ) + V C, with deg x U < deg C. (9) Then
Rdeg n U is bounded by [deg C − / [deg C ] . Proof.
Let R ( n ) = 0 denote the resultant with respectto x of C and An + B . Then R = S ( An + B ) + T C for some
S, T in k [ n, x ] with deg x S < deg C . Moreover, by Cramer’srule applied to the Sylvester matrix of C and An + B , wehave deg R ≤ deg C and deg n S < deg C .Now denote by q ( n ) the monic denominator of U , and by˜ U ∈ k [ n, x ] its numerator. We need to prove that deg q ≤ deg C and deg n ˜ U < deg C . Equalities1 = U ( An + B ) + V C,
S/R ( An + B ) + ( T /R ) C imply by subtraction that C divides ( An + B )( U − S/R )in k ( n )[ x ]. Since C is coprime with An + B , this impliesthat C divides U − S/R . As deg
C < deg x ( U − S/R ), thisshows that U = S/R . In particular, q divides R . It followsthat deg q ≤ deg R ≤ deg C and deg n ˜ U = deg y ( qS/R ) ≤ deg n S < deg C , which concludes the proof. Lemma 11
With the same notation as in Lemma 9, andassuming that G has multiplicity k in the denominator of H ,the output R of BasicReduction ( P, A/B, G, k ) satisfies deg x R max(deg x P − k deg x G, deg x A − , deg x B − , and Rdeg n R − Rdeg n P is bounded by k [1] / [0] , if ν > k [1] / [1] , if ν = 1 and G T ; k deg x ( G )[1] / [1] , if ν = 1 and G | T . Proof.
Both bounds follow from Lemma 9. For the de-gree in n , the polynomial Q of the algorithm is obtainedfrom a modular inverse as above, that is multiplied by R ,and then by C that has rational degree [1] / [0] in n . Thebound directly follows in the first and third cases since thecondition on ν and G is preserved at each step. In the sec-ond case, the condition G ∤ T is not preserved, but writing J = iG ′ /G + S/T at each step shows that the degree in n still increases by [1] / [1] only. For the degree in x , the boundis obtained by bounding each term in the expression of R that is used in the algorithm. Lemma 12
With the same notation, write the denominator g of H as g = efh , where e = gcd( g, T, T ′ ) , f = gcd( g/e, T ) . Also denote m the highest multiplicity of the roots of g , andlet e = e e · · · e mm and h = h h · · · h mm be the square-freedecompositions of e and h .Then, the result of Algorithm HermiteReduction ( P, H, S/T ) is a polynomial in x of degree at most max(deg x P + deg ∞ x H, deg x P − , deg x A − , deg x B − . Seen as a rational function in n , it has degree at most Rdeg n P + [max e k =1 k + P h k =1 k + deg x f ][ P h k =1 k + deg x f ] . Moreover, the certificate Q in the algorithm satisfies Q = qB/ ( gH ) for some polynomial q such that deg x q < deg x g . roof. Following the notation of Algorithm
HermiteRe-duction , the partial fraction decomposition of H produces U with deg x U deg ∞ x H and U k with deg x U k < deg x g k . P U obviously satisfies the bounds. Now write the square-freedecomposition f = f f · · · f mm . By Lemma 11, Rdeg n r k Rdeg n P +[ k ( e k =1 + h k =1 +deg x f k )] / [ k ( h k =1 +deg x f k )].Normalizing R + r + r + · · · + r m then yields the result. Thebound for the degree in x is obtained by bounding separatelyeach term using Lemma 11. The form of Q follows from thefact that the certificate of the k -th call to BasicReduction hasthe form q k B/ ( g kk H ) with deg x q k < deg x ( g kk ). Lemma 13
With the notation of Eq. ( ⋆ ) , Eq. (4) , Algo-rithm MixedCT , and d H = deg ∞ x H , for all i in { , , . . . , δ } we have Rdeg n R i deg n P · [1][0] + α + i ( β + γ ) where α = ( max(deg x P − δ + 1 , · [1][1] if δ = deg x A , h max (cid:16)j deg x P − δδ − deg x A k + 1 , (cid:17)i / [0] otherwise. β = [max e k =1 k + P h k =1 k + deg x f ][ P h k =1 k + deg x f ] γ = ( d H + 1) · [1][1] if d H > and δ = deg x A , hj d H δ − deg x A k + 1 i / [0] if d H > and δ = deg x B − , otherwise . Proof.
By Lemma 8, the initial confinement increasesRdeg n P by α . HermiteReduction is then always used with aninput polynomial of degree less than δ . By Lemma 12 eachcall to HermiteReduction increases Rdeg n P by at most β ,and produces an output of degree at most δ if d H < δ + d H if d H >
0. Thus the confinement is only necessary in thelatter case. Plugging this bound into Lemma 8 shows thateach call to
Confinement increases Rdeg n P by at most γ . Lemma 14
With the notation of ( ⋆ ) , let ( c , . . . , c r − ) bethe output of MixedCT ( P , H , S/T ) and write H = f/g with gcd( f, g ) = 1 .Then there exists a polynomial Q ∈ k ( n )[ x ] such that S rn − r − X i =0 c i S in ! ( P Φ) = (cid:18) Qg r B Φ (cid:19) ′ , with deg x Q r deg x H + max(deg x P − δ, − . Proof.
Tracking the certificates of the various rewritingsin
MixedCT , it suffices to show that for all i ∈ { , . . . , r } P ( n + i − , x ) H i Φ H = R i Φ H + (cid:18) Q i g i B Φ H (cid:19) ′ for some polynomial Q i such thatdeg x Q i max(deg x P − δ,
0) + i deg x H − . This is obvious for i = 0 (initial confinement). Assumethis is true for i −
1, then the next Hermite reduction writes R i − ( n + 1 , x )Φ = ˜ R i Φ H + (cid:18) ˜ Q i g B Φ H (cid:19) ′ with deg x ˜ Q i < deg x g and deg x ˜ R i δ − δ H , where δ H =max(deg ∞ x H,
0) (see Lemma 12). As for the confinement,˜ R i Φ H = R i Φ H + (cid:18) Q i B Φ H (cid:19) ′ , with deg x Q i δ H − i with Q i = fQ i − ( n + 1 , x ) + ˜ Q i g i − + Q i g i , from which follows the bound on the degree of Q i . Theorem 15
With the notation of ( ⋆ ) , the telescoper L produced by Algorithm MixedCT ( P , H , S/T ) satisfies deg n L r (deg n P + deg x H ) + max(deg x P − δ, , where r is the order of L . Proof.
Write H = f/g with gcd( f, g ) = 1. RewritingLemma 14 in terms of polynomials yields P ( n + r ) f r − r − X i =0 c i P ( n + i ) f i g r − i = ( QB ) ′ − rQ Bg ′ g + QA for some Q = P si =0 q i x i with s = max(deg x P − δ,
0) + r deg x H −
1. This equation is a linear system with twoblocks of unknowns: c , . . . , c r − with coefficients of degreebounded by deg n P and q , . . . , q s with coefficients of degree1, which by Hadamard’s bound yields deg n c i r deg n P + s + 1, whence the theorem.
4. COMPLEXITY
We will rely on some classical complexity results for the ba-sic operations on polynomials and rational functions. Stan-dard references for these questions are the books [13] and [6].We will also use the fact that linear systems with polynomialcoefficients can be solved efficiently using Storjohann andVillard’s algorithm [19]. The needed results are summarizedin the following lemma.
Lemma 16
Addition, product and differentiation of ratio-nal functions in K ( x ) of regular degree less than d, as wellas extended gcd and square-free decomposition in K [ x ] d canbe performed using ˜ O ( d ) operations in K .The kernel of a s × ( s + 1) matrix with polynomial entriesin k [ x ] d can be solved using ˜ O ( s ω d ) operations in k . Lemma 17
With P ∈ K [ x ] , A/B ∈ K ( x ) as input, and δ = max(deg x A, deg x B − , Confinement performs at most O ( δ deg x P ) operations in K . Proof.
Each iteration of the loop performs O ( δ ) opera-tions in K and the loop is executed deg x P − δ +1 times. Thecomputation of R then needs ˜ O (deg x P ) operations in K . .2 Hermite Reduction Lemma 18
With the same notation as in the preceding lemma, G a square-free factor of B and k a positive integer, BasicReduction ( P , A/B , G , k ) performs at most ˜ O ( k (deg x P + δ )) operations in K . Proof.
The costly steps are the gcd computations, whichaccording to Lemma 16 can be performed using ˜ O (deg x P + δ ) operations in K . The result then follows since there are k gcd computations. Lemma 19
With the same notation as in Lemma 12, set ǫ = P g k =1 k . Then HermiteReduction ( P , H , S/T ) performsat most ˜ O (max( ǫ deg x P, d H ) + deg x g + ǫδ ) operations in K = k ( n ) . Proof.
By Lemma 16, the square-free decomposition of g can be computed in ˜ O (deg x g ) operations and the product P U is computed at a cost ˜ O (max(deg x P, deg ∞ x H )). Bythe preceding lemma, the k -th call to BasicReduction uses˜ O ( k (deg x P +deg x g k + δ )) operations. The announced boundis then obtained by summation. For the sake of simplicity, Algorithm
MixedCT searchesfor telescopers for all the possible orders, starting from 0.In practice, a more efficient variant consists in carrying adichotomic search of the order between 0 and δ . This waythe complexity is that of the last step up to a logarithmicfactor. Here, we analyze this variant. Lemma 20
With the same notation as in Section 3.3, and ǫ = P g k =1 k , the number of operations in K = k ( n ) per-formed by MixedCT ( P , H , S/T ) is1. ˜ O ( δ max( ǫ deg x P, d H ) + δ deg x g + ǫδ + δ ω ) if d H < ;2. ˜ O ( δ max( ǫ deg x P, δd H )+ δ deg x g + ǫδ + δ ) if d H > . Proof.
When d H <
0, by Lemma 12 there are no con-finement steps and the construction of the system to solveamounts to δ Hermite reductions. The algorithm then com-putes a vector in the kernel of a δ × ( δ + 1) matrix, whichby Lemma 16 can be performed in ˜ O ( δ ω ) operations in K ,hence the complexity.When d H >
0, we have to add the cost of the confinementsteps, which by Lemma 12 are performed on polynomials ofdegree at most δ + d H . There are at most δ + 1 calls toconfinement, so the result follows from Lemma 17. Theorem 21
With the same notation, set µ = max( a, b ) where deg n P · [1] / [0] + α + δ ( β + γ ) = [ a ] / [ b ] , the number ofoperations in k performed by Algorithm MixedCT is ˜ O ( µ ( ǫδ deg x P + δ deg x g + ǫδ + δ ω )) , if d H < , or ˜ O ( µ ( ǫδ deg x P + δ deg x g + ǫδ + δ + δ d H )) if d H > . Proof.
The result follows directly from the precedinglemma and the fact that all the elements of K = k ( n ) appear-ing in the construction of the linear system have numeratorand denominator of degree in n bounded by µ . The cost ofsolving is then ˜ O ( δ ω µ ) operations in k by Lemma 16.The complexity result from Theorem 1 follows directly.
5. EXPERIMENTS AND APPLICATIONS5.1 Various Integrals
Example 1.
The Jacobi polynomials have the following in-tegral representation, up to a factor that does not dependon n : I (cid:18) z − z − x ) (cid:19) n (1 − z ) α (1 + z ) β dzz − x , with a contour enclosing z = x once in the positive sense [18,18.10.8]. It is well-known that the Jacobi polynomials satisfya recurrence of order 2. Theorem 6 is tight in that case:it predicts a bound 2 on the order of the telescoper, sincethe logarithmic derivative of the integrand has numerator ofdegree 2 and denominator of degree 3.Note however that in such an example, a more direct andefficient way to obtain a recurrence is to change the vari-able z into x + u making the integral that of the extractionof the n -th coefficient in a hyperexponential term. That isachieved easily by translating the first order linear differen-tial equation it satisfies into a linear recurrence. The onlydifference is that our method guarantees the minimality ofthe telescoper. Example 2.
A typical example where several of the diffi-culties are met at once is the term (cid:16) xn +1 (cid:17) (cid:16) ( x +1) ( x − x − ( x − (cid:17) n p x − e x x ( x − x − . Our code finds a telescoper of order 9 and degree 90 in 1.4sec. The only other code we are aware of that can performthis computation is Koutschan’s HolonomicFunctions pack-age [16], which takes more than 3 min. This is by no meansa criticism of this excellent and very versatile package butrather to indicate the advantage of implementing specializedalgorithms like ours for this class.
Here is a generalization of Manivel’s lemma [12], whichled us to this work. It gives an efficient way to compute therecurrence satisfied by the coefficients of the compositionalinverse of a rational function. The starting point is Lagrangeinversion. Let f ∈ Q ( x ) be a rational function such that f (0) = 0, and denote by f ( − its compositional inverse. ByCauchy’s formula, the n -th coefficient u n of f ( − is givenby u n = 12 πi I f ( − ( x ) d xx n +1 , where the contour is a small circle around the origin. In-tegrating by parts and then using the change of variables The code and Maple worksheet for these examples can befound at http://mixedct.gforge.inria.fr. = f ( u ) yields u n = 12 πin I f ( − ′ ( x )d xx n = 12 πin I d uf ( u ) n . Thus a recurrence for u n can be computed with Algorithm MixedCT . Theorem 6 and Theorem 15 then provide boundsfor the order and degree of this recurrence.
Theorem 22
Let f ∈ Q ( x ) be a rational function such that f (0) = 0 . Write f = P/Q and denote P = P P · · · P mm thesquare-free decomposition of P . Also denote p, p ⋆ , q, q ⋆ thedegrees of P, P ⋆ , Q, Q ⋆ respectively. Then the Taylor coeffi-cients of f ( − satisfy a recurrence of order at most q ⋆ + p ⋆ − and of rational degree in n at most ( q ⋆ + p ⋆ )( q ⋆ + p ⋆ + 1)2 (max( q − p,
0) + X P k =1 k ) · [1][1] .k order degree coeffs time5 10 61 1759 2.496 12 88 2440 4.647 14 120 3778 13.368 16 157 4666 33.899 18 199 6192 88.3410 20 246 8364 260.5911 22 298 10146 628.2112 24 355 11802 1451.54 Table 1.
Order, degree, bit size and timings for Example 3
Example 3.
Experimental results on the family of rationalfunctions f k = xP k ( x ) /Q k ( x ) with P k and Q k two densepolynomials of degree k and integer coefficients of absolutevalue bounded by 100 are presented in Table 1. The firstcolumn gives the index k . The second one is the order ofthe minimal-order telescopers, which is as predicted by The-orem 22. The next one gives the degree of the telescoper;it displays a quadratic growth, as predicted by Theorem 22.The column “coeffs” gives the bit size of the largest coeffi-cient of the telescoper, whose growth seems slightly morethan quadratic. Finally, the time (in seconds) taken by ourimplementation is given in the last column. Acknowledgements.
We are grateful to the referees fortheir thorough work and helpful comments. This work hasbeen supported in part by FastRelax ANR-14-CE25-0018-01.
6. REFERENCES [1] G. Almkvist and D. Zeilberger. The method ofdifferentiating under the integral sign.
J. SymbolicComput. , 10(6):571–591, 1990.[2] M. Apagodu and D. Zeilberger. Multi-variable Zeilbergerand Almkvist-Zeilberger algorithms and the sharpening ofWilf-Zeilberger theory.
Adv. in Appl. Math. , 37(2):139–152,2006. [3] A. Bostan, S. Chen, F. Chyzak, and Z. Li. Complexity ofcreative telescoping for bivariate rational functions. In
ISSAC’10 , pages 203–210. ACM, New York, 2010.[4] A. Bostan, S. Chen, F. Chyzak, Z. Li, and G. Xin. Hermitereduction and creative telescoping for hyperexponentialfunctions. In
ISSAC’13 , pages 77–84. ACM, New York,2013.[5] A. Bostan, P. Lairez, and B. Salvy. Creative telescoping forrational functions using the Griffiths-Dwork method. In
ISSAC’13 , pages 93–100. ACM, New York, 2013.[6] P. B¨urgisser, M. Clausen, and M. A. Shokrollahi.
Algebraiccomplexity theory , volume 315 of
Grundlehren derMathematischen Wissenschaften . Springer, 1997.[7] S. Chen, F. Chyzak, R. Feng, G. Fu, and Z. Li. On theexistence of telescopers for mixed hypergeometric terms.
J.Symbolic Comput. , 68:1–26, 2015.[8] S. Chen, H. Huang, M. Kauers, and Z. Li. A modifiedAbramov-Petkovˇsek reduction and creative telescoping forhypergeometric terms. In
ISSAC’15 , pages 117–124. ACM,New York, 2015.[9] S. Chen and M. Kauers. Order-degree curves forhypergeometric creative telescoping. In
ISSAC’12 , pages122–129. ACM, New York, 2012.[10] S. Chen and M. Kauers. Trading order for degree in creativetelescoping.
J. Symbolic Comput. , 47(8):968–995, 2012.[11] F. Chyzak. An extension of Zeilberger’s fast algorithm togeneral holonomic functions.
Discrete Math. ,217(1-3):115–134, 2000.[12] J.-P. Furter. Polynomial composition rigidity and planepolynomial automorphisms.
J. Lond. Math. Soc. (2) ,91(1):180–202, 2015.[13] J. von zur Gathen and J. Gerhard.
Modern ComputerAlgebra . Cambridge Univ. Press, second edition, 2003.[14] M. Kauers and L. Yen. On the length of integers intelescopers for proper hypergeometric terms.
J. SymbolicComput. , 66:21–33, 2015.[15] C. Koutschan. Creative telescoping for holonomic functions.In C. Schneider and J. Bl¨umlein, editors,
Computer Algebrain Quantum Field Theory: Integration, Summation andSpecial Functions , Texts & Monographs in SymbolicComputation, pages 171–194. Springer, Wien, 2013.[16] C. Koutschan. Holonomic Functions in Mathematica.
ACMCommun. Comput. Algebra , 47(3/4):179–182, Jan. 2014.[17] P. Lairez. Computing periods of rational integrals.
Math.Comp. , 85(300):1719–1752.[18] F. W. J. Olver, D. W. Lozier, R. F. Boisvert, and C. W.Clark, editors.
NIST Handbook of Mathematical Functions .Cambridge University Press, 2010.[19] A. Storjohann and G. Villard. Computing the rank and asmall nullspace basis of a polynomial matrix. In
ISSAC’05 ,pages 309–316. ACM, New York, 2005.[20] A. Tefera. MultInt, a MAPLE package for multipleintegration by the WZ method.
J. Symbolic Comput. ,34(5):329–353, 2002.[21] H. S. Wilf and D. Zeilberger. An algorithmic proof theoryfor hypergeometric (ordinary and “ q ”) multisum/integralidentities. Invent. Math. , 108(3):575–633, 1992.[22] L. Yen.
Contributions to the proof theory of hypergeometricidentities . ProQuest LLC, Ann Arbor, MI, 1993. Thesis(Ph.D.)–University of Pennsylvania.[23] D. Zeilberger. A holonomic systems approach to specialfunctions identities.
J. Comput. Appl. Math. ,32(3):321–368, 1990.[24] D. Zeilberger. The method of creative telescoping.