Complexity of Creative Telescoping for Bivariate Rational Functions
aa r X i v : . [ c s . S C ] J a n Complexity of Creative Telescopingfor Bivariate Rational Functions ∗ Alin Bostan, Shaoshi Chen, Fr ´ed ´eric Chyzak
Algorithms Project-Team, INRIAParis-Rocquencourt78153 Le Chesnay (France) { alin.bostan,shaoshi.chen,frederic.chyzak } @inria.fr Ziming Li
Key Laboratory of Mathematics Mechanization,Academy of Mathematics and System Sciences100190 Beijing (China) [email protected]
ABSTRACT
The long-term goal initiated in this work is to obtain fastalgorithms and implementations for definite integration inAlmkvist and Zeilberger’s framework of (differential) cre-ative telescoping. Our complexity-driven approach is to ob-tain tight degree bounds on the various expressions involvedin the method. To make the problem more tractable, we re-strict to bivariate rational functions. By considering thisconstrained class of inputs, we are able to blend the generalmethod of creative telescoping with the well-known Hermitereduction. We then use our new method to compute diago-nals of rational power series arising from combinatorics.
Categories and Subject Descriptors:
I.1.2 [
Computing Methodologies ]: Symbolic and Alge-braic Manipulations —
Algebraic Algorithms
General Terms:
Algorithms, Theory.
Keywords:
Hermite reduction, creative telescoping.
1. INTRODUCTION
The long-term goal of the research initiated in the presentwork is to obtain fast algorithms and implementations forthe definite integration of general special functions, in acomplexity-driven perspective.As most special-function integrals cannot be expressed inclosed form, their evaluation cannot be based on table look-ups only, and even when closed forms are available, theymay prove to be intractable in further manipulations. Inboth cases, the difficulty can be mitigated by representingfunctions by annihilating differential operators. This moti-vated Zeilberger to introduce a method now known as cre-ative telescoping [20], which applies to a large class of specialfunctions: the D-finite functions [14] defined by sets of lineardifferential equations of any order, with polynomial coeffi-cients. Zeilberger’s method applies in general to multipleintegrals and sums. ∗ We warmly thank the referees for their very helpful comments.— AB and FC were supported in part by the Microsoft Research –Inria Joint Centre, and SC and ZL by a grant of the NationalNatural Science Foundation of China (No. 60821002).
Permission to make digital or hard copies of all or part of this work forpersonal or classroom use is granted without fee provided that copies arenot made or distributed for profit or commercial advantage and that copiesbear this notice and the full citation on the first page. To copy otherwise, torepublish, to post on servers or to redistribute to lists, requires prior specificpermission and/or a fee.
ISSAC XXXX,
Copyright 2010 ACM ... $ A sketch of Zeilberger’s method is as follows. Given a D-finite function f of the variables x and y , the definite integral F ( x ) = R βα f ( x, y ) dy is D-finite, and a linear differentialequation satisfied by F can be constructed [20]. To explainthis, let k be a field of characteristic zero, D x and D y bethe usual derivations on the rational-function field k ( x, y ),both restricting to zero on k , and let k ( x, y ) h D x , D y i be thering of linear differential operators over k ( x, y ). The heart ofthe method is to solve the differential telescoping equation (1) below for L ∈ k [ x ] h D x i \ { } and g = R ( f ) for some R ∈ k ( x, y ) h D x , D y i . The operator L is called a telescoper for f , and g a certificate of L for f . Under the assumptionlim y → α g ( x, y ) = lim y → β g ( x, y ) for x in some domain ,L ( x, D x ) is then proved to be an annihilator of F .The main emphasis in works since the 1990’s has beenon finding telescopers of order minimal over all telescopersfor f , which are called minimal telescopers . (Two minimaltelescopers differ by a multiplicative factor in k ( x ).) In viewof the computational difficulty of solving (1), there has beenspecial attention to subclasses of inputs. Of particular im-portance is the case of hyperexponential functions, definedby first-order differential equations, studied by Almkvist andZeilberger in [1]. Their method is a direct differential ana-logue of Zeilberger’s algorithm for the recurrence case [21].On the other hand, very little is known about the com-plexity of creative telescoping: the only related result seemsto be an analysis in [9] of an algorithm for hyperexponen-tial indefinite integration. In order to get complexity esti-mates, we simplify the problem by restricting to a smallerclass of inputs, namely that of bivariate rational functions.Although restricted, this class already has many applica-tions, for instance in combinatorics, where many nontrivialproblems are encoded as diagonals of rational formal powerseries, themselves expressible as integrals. Our goal thusreads as follows. Problem
Given f = P/Q ∈ k ( x, y ) \{ } , find a pair ( L, g ) with L = P ρi =0 η i ( x ) D ix in k [ x ] h D x i \ { } and g in k ( x, y ) such that L ( x, D x )( f ) = D y ( g ) . (1)By considering this more constrained class of inputs, weare indeed able to blend the general method of creative tele-scoping with the well-known Hermite reduction [10].Essentially two algorithms for minimal telescopers can befound in the literature: The classical way [1] is to apply aethod deg D x ( L ) deg x ( L ) deg x ( g ) deg y ( g ) ComplexityMinimal Hermite reduction (new) ≤ d y O ( d x d y ) O ( d x d y ) O ( d y ) ˜ O ( d x d ω +3 y ) Las VegasTelescoper Almkvist and Zeilberger ≤ d y O ( d x d y ) O ( d x d y ) O ( d y ) ˜ O ( d x d ω +2 y ) Las VegasNonminimal Lipshitz elimination ≤ d x + 1)( d y + 1) O ( d x d y ) O ( d x d y ) O ( d x d y ) O ( d ωx d ωy ) deterministicTelescoper Cubic size ≤ d y O ( d x d y ) O ( d x d y ) O ( d y ) O ( d ωx d ωy ) deterministic Figure 1: Complexity of creative telescoping methods (under Hyp. (H’)), together with bounds on output differential analogue of Gosper’s indefinite summation algo-rithm, which reduces the problem to solving an auxiliarylinear differential equation for polynomial solutions. An al-gorithm developed later in [7] (see also [12]) performs Her-mite reduction on f to get an additive decomposition of theform f = D y ( a ) + P mi =1 u i /v i , where the u i and v i are in k ( x )[ y ] and the v i are squarefree. Then, the algorithm in [1]is applied to each u i /v i to get a telescoper L i minimal for it.The least common left multiple of the L i ’s is then proved tobe a minimal telescoper for f . This algorithm performs wellonly for specific inputs (both in practice and from the com-plexity viewpoint), but it inspired our Lemma 22 via [12].As a first contribution in this article, we present a new,provably faster algorithm for computing minimal telescopersfor bivariate rational functions. Instead of a single use ofHermite reduction as in [12], we apply Hermite reduction tothe D ix ( f )’s, iteratively for i = 0 , , . . . , which yields D ix ( f ) = D y ( g i ) + w i w (2)for some factor w of the squarefree part of the denominatorof f . If η , . . . , η ρ ∈ k ( x ) are not all zero and such that P ρi =0 η i w i = 0, then the operator P ρi =0 η i D ix is a telescoperfor f , and more specifically, the first nontrivial linear relationobtained in this way yields a minimal telescoper for f .As a second contribution, we give the first proof of a poly-nomial complexity for creative telescoping on a specific classof inputs, namely on bivariate rational functions. For min-imal telescopers, only a polynomial bound on d x (but noneon d y ) was given for special inputs in [7]; more specifically,we derive complexity estimates for all mentioned methods(see Fig. 1), showing that our approach is faster. Further-more, we analyse the bidegrees of non minimal telescop-ers generated by other approaches: Lipshitz’ work [13] canbe rephrased into an existence theorem for telescopers withpolynomial size; the approach followed in the recent work onalgebraic functions [3] leads to telescopers of smaller degreesizes. These are new instances of the philosophy, promotedin [3], that relaxing minimality can produce smaller outputs.A third contribution is a fast Maple implementation [22],incorporating a careful implementation of the original Her-mite reduction algorithm, making use of the special formof w i /w in (2) and of usual modular techniques (probabilis-tic rank estimate) to determine when to invoke the solverfor linear algebraic equations. Experimental results indicatethat our implementation outperforms Maple’s core routine.Note that for the fastest method we propose, denotedby H1 in Tables 1–3, we chose to output the certificate asa mere sum of (small) rational functions, without any formof normalisation. This choice seems to be uncommon forcreative-telescoping algorithms, but a motivation is how thecertificate is used in practice: Very often, like for appli-cations to diagonals in §
5, the certificate is actually notneeded. In other applications, the next step of the method of creative telescoping is to integrate (1) between α and β ,leading to L ( F )( x ) = g ( x, α ) − g ( x, β ). Therefore, only eval-uations of the certificate are really needed, and normalisa-tion can be postponed to after specialising at α and β .The end of this section, § k ( x ) in § § §
4: cubic for nonminimal order insteadof quartic for minimal order. See the summary in Figure 1,where the low complexity of algorithms for minimal tele-scopers relies on Storjohann and Villard’s algorithms [19],thus inducing a certified probabilistic feature. We apply ourresults to the calculation of diagonals in §
5, and describe ourimplementation and comment on execution timings in § We recall basic notation and complexity facts for later use.Let k be again a field of characteristic zero. Unless other-wise specified, all complexity estimates are given in terms ofarithmetical operations in k , which we denote by “ops”. Let k [ x ] m × n ≤ d be the set of m × n matrices with coefficients in k [ x ]of degree at most d . Let ω ∈ [2 ,
3] be a feasible exponent ofmatrix multiplication, so that two matrices from k n × n canbe multiplied using O ( n ω ) ops. Facts 1 and 2 below show thecomplexity of multipoint evaluation, rational interpolation,and algebraic operations on polynomial matrices using fastarithmetic, where the notation ˜ O ( · ) indicates cost estimateswith hidden logarithmic factors [6, Def. 25.8]. Fact 1
For p ∈ k [ x ] of degree less than n , pairwise distinct u , . . . , u n − in k , and v , . . . , v n − ∈ k , we have:(i) Evaluating p at the u i ’s takes ˜ O ( n ) ops.(ii) For m ∈ { , . . . , n } , constructing f = s/t ∈ k ( x ) with deg x ( s ) < m and deg x ( t ) ≤ n − m such that t ( u i ) = 0 and f ( u i ) = v i for ≤ i ≤ n − takes ˜ O ( n ) ops. Fact 2
For M in k [ x ] m × n ≤ d , d > , we have:(i) If M = (cid:0) M M (cid:1) is an invertible n × n matrix with M i ∈ k [ x ] n × n i ≤ d i , where i = 1 , and n + n = n , thenthe degree of det( M ) is at most n d + n d .(ii) If M = (cid:0) M M (cid:1) is not of full rank and with M i ∈ k [ x ] m × n i ≤ d i , where i = 1 , and n + n = n , then thereexists a nonzero u ∈ k [ x ] n with coefficients of degreeat most n d + n d such that Mu = 0 .(iii) The rank r and a basis of the null space of M can becomputed using ˜ O ( nmr ω − d ) ops. For proofs, see [6, Cor. 10.8, 5.18, 11.6] and [19, Th. 7.3].)We call squarefree factorisation of Q ∈ k [ x, y ] \ k [ x ] w.r.t. y the unique product qQ Q · · · Q mm equal to Q for q ∈ k [ x ] and Q i ∈ k [ x, y ] satisfying deg y ( Q m ) > Q i ’sare primitive, squarefree, and pairwise coprime. The square-free part Q ∗ of Q w.r.t. y is the product Q Q · · · Q m . Let Q − denote the polynomial Q/Q ∗ , and lc y ( Q ) the leadingcoefficient of Q w.r.t. y . The following two formulas about Q , Q ∗ , and Q − can be proved by mere calculations. Fact 3
Let ˆ Q i denote Q ∗ /Q i . Then we have(i) Q ∗ D y ( Q − ) /Q − = P mi =1 ( i −
1) ˆ Q i D y ( Q i ) ∈ k [ x, y ] ;(ii) D y ( Q ) /Q − = P mi =1 i ˆ Q i D y ( Q i ) ∈ k [ x, y ] . Let f = P/Q be a nonzero element in k ( x, y ), where P, Q are two coprime polynomials in k [ x, y ]. The degree of f in x is defined to be max { deg x ( P ) , deg x ( Q ) } , and denotedby deg x ( f ). The degree of f in y is defined similarly. The bidegree of f is the pair (deg x ( f ) , deg y ( f )), which is denotedby bideg( f ). The bidegree of f is said to be bounded (above)by ( α, β ), written bideg( f ) ≤ ( α, β ), when deg x ( f ) ≤ α and deg y ( f ) ≤ β .We say that f = P/Q is proper if the degree of P in y is less than that of Q . For creative telescoping, we mayalways assume w.l.o.g. that f = P/Q is proper. If not,rewrite f = D y ( p ) + ¯ f with p ∈ k ( x )[ y ] and ¯ f proper. Atelescoper L for ¯ f with certificate ¯ g is a telescoper for f withcertificate L ( p ) + ¯ g . Hypothesis (H)
From now on, P and Q are assumed to benonzero polynomials in k [ x, y ] such that deg y ( P ) < deg y ( Q ) , gcd( P, Q ) = 1 , and Q is primitive w.r.t. y . Notation
From now on, we write ( d x , d y ) , ( d ∗ x , d ∗ y ) , and ( d − x , d − y ) for the bidegrees of Q , Q ∗ , and Q − , respectively. The following hypothesis makes our estimates concise.
Hypothesis (H’)
Occasionally, we shall require the ex-tended hypothesis: Hypothesis (H) and deg x ( P ) ≤ d x .
2. HERMITE REDUCTION
Let K be a field of characteristic zero, either k or k ( x )in what follows. Let K ( y ) be the field of rational functionsin y over K , and D y be the usual derivation on it. For a ra-tional function f ∈ K ( y ), Hermite reduction [10] computesrational functions g and r = a/b in K ( y ) satisfying f = D y ( g ) + r, deg y ( a ) < deg y ( b ) , b is squarefree. (3)Horowitz and Ostrogradsky’s method [15, 11] computes thesame decomposition as in (3) by solving a linear system. Forthe details of those methods, see [4, Chapter 2]. Lemma 4 If f is proper, a pair ( g, r ) satisfying (3) forproper g, r is unique. Proof.
This is a consequence of [11, Theorem 2.10] afterwriting r as a sum P mi =1 α i / ( x − b i ) and integrating. Lemma 5
Let f be a nonzero rational function in K ( y ) ofdegree at most n in y , then Hermite reduction on f can beperformed using ˜ O ( n ) operations in K . Proof.
See [6, Theorem 22.7]. In contrast, the method of Horowitz and Ostrogradsky takes O ( n ω ) operations in K [6, § K = k ( x ) and analyse the complexityof Hermite reduction over k ( x ) in terms of operations in k .To this end, we use an evaluation-interpolation approach. We derive an upper bound on the bidegrees of g and r satisfying (3) by studying the linear system in [11].Analysing Hermite reduction (under (H)) shows the exis-tence of A, a ∈ k ( x )[ y ] with deg y ( A ) < d − y , deg y ( a ) < d ∗ y and PQ = D y (cid:18) AQ − (cid:19) + aQ ∗ . (4)In order to bound the bidegrees of A and a , we reformu-late (4) into the equivalent form P = Q ∗ D y ( A ) − (cid:18) Q ∗ D y ( Q − ) Q − (cid:19) A + Q − a, (5)where Q ∗ D y ( Q − ) /Q − is a polynomial in k [ x, y ] of bidegreeat most ( d ∗ x , d ∗ y −
1) by Fact 3. Viewing A and a as polyno-mials in k ( x )[ y ] with undetermined coefficients, we form thefollowing linear system, equivalent to (5), (cid:0) H H (cid:1) (cid:18) ˆ A ˆ a (cid:19) = ˆ P , (6)where H ∈ k [ x ] d y × d − y ≤ d ∗ x , H ∈ k [ x ] d y × d ∗ y ≤ d − x , and ˆ A , ˆ a , and ˆ P are the coefficient vectors of A , a , and P with sizes d − y , d ∗ y ,and d y , respectively. Under the constraint of properness of A/Q − and a/Q ∗ , ( A, a ) is unique by Lemma 4. Then (6)has a unique solution, which leads to the following lemma.
Lemma 6
The matrix (cid:0) H H (cid:1) is invertible over k ( x ) . As the matrix (cid:0) H H (cid:1) is uniquely defined by Q , we callit the matrix associated with Q , denoted by H ( Q ). Let δ be its determinant, so that deg x ( δ ) ≤ µ := d ∗ x d − y + d − x d ∗ y byFact 2 (i) . For later use, we also define δ ′ as the determinantof H ( Q ∗ ), so that deg x ( δ ′ ) ≤ µ ′ := 2 d ∗ x d ∗ y by Fact 2 (i) andsince ( Q ∗ ) − = Q ∗ . Lemma 7
There exist
B, b ∈ k [ x, y ] with deg y ( B ) < d − y and deg y ( b ) < d ∗ y , and such that:(i) PQ = D y (cid:16) BδQ − (cid:17) + bδQ ∗ ;(ii) deg x ( B ) ≤ µ − d ∗ x + deg x ( P ) and deg x ( b ) ≤ µ − d − x +deg x ( P ) . Proof.
Applying Cramer’s rule to (6) leads to (i) . As-sertion (ii) next follows by determinant expansions.In what follows, we shall encounter proper rational func-tions with denominator Q satisfying Q = Q ∗ . The followinglemma is an easy corollary of Lemma 7 for such functions. Corollary 8
Assuming Q = Q ∗ in addition to Hypothe-sis (H), there exist B, b ∈ k [ x, y ] with deg y ( B ) and deg y ( b ) less than d ∗ y , and such that(i) PQ ∗ = D y (cid:16) Bδ ′ Q ∗ (cid:17) + bδ ′ Q ∗ ;(ii) deg x ( B ) and deg x ( b ) are bounded by µ ′ − d ∗ x +deg x ( P ) . .2 Algorithm by evaluation and interpolation We observe that an asymptotically optimal complexitycan be achieved by evaluation and interpolation at each stepof Hermite reduction over k ( x ). This inspires us to adaptGerhard’s modular method [8, 9] to k ( x, y ). Recall that, byHyp. (H), Q ∈ k [ x, y ] is nonzero and primitive over k [ x ]. Definition
An element x ∈ k is lucky if lc y ( Q )( x ) = 0 and deg y (gcd( Q ( x , y ) , D y ( Q ( x , y )))) = d − y . Lemma 9
There are at most d x (2 d ∗ y − unlucky points. Proof.
Let σ ∈ k [ x ] be the d − y th subresultant w.r.t. y of Q and D y ( Q ). By [9, Corollary 5.5], all unlucky pointsare in the set U = { x ∈ k | σ ( x ) = 0 } . By [9, Corollary3.2 (ii) ], deg x ( σ ) ≤ d x (2 d ∗ y − Lemma 10
Let B , b , and δ be the same as in Lemma 7, andlet x ∈ k be lucky. Then δ ( x ) = 0 and ( B ( x , y ) , b ( x , y )) is the unique pair such that P ( x , y ) Q ( x , y ) = D y (cid:18) B ( x , y ) δ ( x ) Q − ( x , y ) (cid:19) + b ( x , y ) δ ( x ) Q ∗ ( x , y ) . (7) Proof.
By the luckiness of x , deg y ( Q ( x , y )) = d y and Q ( x , y ) − = Q − ( x , y ), so Q ( x , y ) ∗ = Q ∗ ( x , y ). This im-plies H ( Q )( x , y ) = H ( Q ( x , y )), which, by Lemma 6, isinvertible over k ( x ). Hence δ ( x ) = 0, and the evaluationat x = x of the equality in Lemma 7 (i) is well-defined.Thus, ( B ( x , y ) , b ( x , y )) is a solution of (7). Uniquenessfollows from Lemma 4. Theorem 11
Algorithm
HermiteEvalInterp in Figure 2 iscorrect and takes ˜ O ( d x d y + deg x ( P ) d y ) ops. Proof.
Set ν to d x (2 d ∗ y − λ + 1 lucky points found in Step 3 are all less than λ + ν + 1. By Lemmas 4 and 7 (i) , A = B/δ and a = b/δ . ByLemma 10, A = B ( x , y ) /δ ( x ) and a = b ( x , y ) /δ ( x ).By Lemma 7 (ii) and since deg x ( δ ) ≤ µ , it suffices to ratio-nally interpolate A and a from values at λ + 1 lucky points.This shows the correctness. The dominant computation inStep 1 is the gcd, which takes ˜ O ( d x d y ) ops by [6, Cor. 11.9].For each integer i ≤ λ + ν , testing luckiness amounts to eval-uations at x and computing gcd( Q ( x , y ) , D y ( Q ( x , y ))),which takes ˜ O ( d y ) ops by Fact 1 (i) and [6, Cor. 11.6]. Then,generating S in Step 3 costs ˜ O (( λ + ν + 1) d y ) ops. ByFact 1 (i) , evaluations in Step 4 take ˜ O (( λ + 1) d y ) ops. Foreach x ∈ S , the cost of the Hermite reduction in Step 4 is˜ O ( d y ) ops by Lemma 5. Thus, the total cost of Step 4 is˜ O (( λ + 1) d y ) ops. By Fact 1 (ii) , Step 5 takes ˜ O (( λ + 1) d y )ops. Since λ ≤ d x d y + deg x ( P ) and ν ≤ d x d y , the totalcost is as announced.As the generic output size of Hermite reduction is propor-tional to λd y , which is O (( d x d y + deg x ( P )) d y ), Algorithm HermiteEvalInterp has quasi-optimal complexity.
3. MINIMAL TELESCOPERS
We analyse two algorithms for constructing minimal tele-scopers for bivariate rational functions and their certificates. Algorithm
HermiteEvalInterp ( P, Q ) Input : P, Q ∈ k [ x, y ] satisfying Hypothesis (H). Output : (
A, a ) ∈ k ( x )[ y ] solving (4).1. Compute Q − := gcd( Q, D y ( Q )) and Q ∗ := Q/Q − ;2. Set λ := 2( d ∗ x d − y + d ∗ y d − x ) + deg x ( P ) − min { d − x , d ∗ x } ;3. Set S to the set of λ +1 smallest nonnegative integersthat are lucky for Q ;4. For each x ∈ S , compute ( A , a ) ∈ k [ y ] such that P ( x , y ) Q ( x , y ) = D y (cid:18) A Q − ( x , y ) (cid:19) + a Q ∗ ( x , y )using Hermite reduction over k ;5. Compute ( A, a ) ∈ k ( x )[ y ] by rational interpolationand return this pair. Figure 2: Hermite reduction over k ( x ) via evaluationand interpolation. We design a new algorithm, presented in Figure 3, to com-pute minimal telescopers for rational functions by basing onHermite reduction. For f = P/Q ∈ k ( x, y ) and i ∈ N , Her-mite reduction decomposes D ix ( f ) into D ix ( f ) = D y ( g i ) + r i , (8)where g i , r i ∈ k ( x, y ) are proper. Since the squarefree partof the denominator of D ix ( f ) divides Q ∗ , so does the de-nominator of r i . The following lemma shows that (8) re-combines into telescopers and certificates; next, Lemma 13implies that the first pair obtained in this way by Algorithm HermiteTelescoping in Figure 3 yields a minimal telescoper.
Lemma 12
The rational functions r , . . . , r d ∗ y are linearlydependent over k ( x ) . Proof.
The constraints on r i imply deg y ( r i Q ∗ ) < d ∗ y forall i ∈ N , from which follows the existence of a nontriviallinear dependence among the r i ’s over k ( x ). Lemma 13
An integer ρ is minimal such that P ρi =0 η i r i =0 for η , . . . , η ρ ∈ k ( x ) not all zero if and only if P ρi =0 η i D ix is a minimal telescoper for f with certificate P ρi =0 η i g i . Proof.
Multiplying (8) by η i before summing yields L ( f ) = D y (cid:18) ρ X i =0 η i g i (cid:19) + ρ X i =0 η i r i for L := ρ X i =0 η i D ix , where the first two sums are proper. Thus, by Lemma 4, L isa telescoper of order ρ for f with certificate P ρi =0 η i g i if andonly if P ρi =0 η i r i = 0 with η ρ = 0. The lemma follows. Lemmas 12 and 13 combine into an upper bound on theorder of minimal telescopers for f . Corollary 14
Minimal telescopers have order at most d ∗ y . The bound 6 d y is shown in [3] for rational functions ofthe form yD y ( Q ) /Q with Q ∈ k [ x, y ]. Apagodu and Zeil-berger [2] obtain a similar bound for a class of nonrationalyperexponential functions, but their proof does not seemto apply to rational functions, as it heavily relies on thepresence of a nontrivial exponential part.We also derive a lower bound on the order of the minimaltelescoper, to be used as an optimisation at the end of § x ∈ k , next applying Hermite reductionin k ( y ) to D ix ( f )( x , y ), yields D ix ( f )( x , y ) = D y ( g ,i ) + r ,i , (9)where g ,i , r ,i ∈ k ( y ) are proper and the denominator of r ,i divides Q ∗ ( x , y ). Let ρ be the smallest integer such that r , , . . . , r ,ρ are linearly dependent over k . Lemma 15
A minimal telescoper has order at least ρ . Proof.
We first claim that r ,i = r i ( x , y ), for r i asin (8). Note that the squarefree part w.r.t. y of the denomi-nator of D ix ( f ) divides Q ∗ for all i ∈ N . By [9, Cor. 5.5], x islucky for the denominator of D ix ( f ) for all i ∈ N . Then, theclaim on r ,i follows from Lemma 10 applied to D ix ( f ). Let ρ be the minimal order of a telescoper, then r , . . . , r ρ are lin-early dependent over k ( x ) by Lemma 13. Thus r , , . . . , r ,ρ are linearly dependent over k , which implies ρ ≤ ρ . To derive degree bounds for g i and r i in (8), let δ , δ ′ , µ ,and µ ′ be defined as before Lemma 7, and set µ ′′ = µ + µ ′ − Lemma 16
Let W be in k [ x, y ] with deg y ( W ) < d ∗ y . Then,for all i ∈ N , there exist B, b ∈ k [ x, y ] with both bideg( B ) and bideg( b ) bounded by (deg x ( W ) + µ ′′ , d ∗ y − , such that D x (cid:18) Wδ i +1 δ ′ i Q ∗ (cid:19) = D y (cid:18) Bδ i +2 δ ′ i +1 Q ∗ (cid:19) + bδ i +2 δ ′ i +1 Q ∗ . Proof.
A straightforward calculation leads to D x (cid:18) Wδ i +1 δ ′ i Q ∗ (cid:19) = ˜ Wδ i +2 δ ′ i +1 Q ∗ − δ i +1 δ ′ i W D x ( Q ∗ ) Q ∗ , where bideg( ˜ W ) ≤ (deg x ( W ) + µ ′′ , d ∗ y − B, ˜ b ∈ k [ x, y ] such that1 δ i +1 δ ′ i W D x ( Q ∗ ) Q ∗ = 1 δ i +2 δ ′ i +1 (cid:18) D y (cid:18) δ ˜ BQ ∗ (cid:19) + δ ˜ bQ ∗ (cid:19) , with bideg( ˜ B ) and bideg(˜ b ) bounded by (deg x ( W ) + µ ′ − ,d ∗ y − B, b ) = ( − δ ˜ B, ˜ W − δ ˜ b ) ends the proof. Lemma 17
For i ∈ N , there exist B i , b i ∈ k [ x, y ] such that D ix ( f ) = D y (cid:18) B i δ i +1 δ ′ i Q ∗ i Q − (cid:19) + b i δ i +1 δ ′ i Q ∗ . (10) Moreover, bideg( B i ) ≤ (deg x ( P ) + µ + iµ ′′ + ( i − d ∗ x , id ∗ y + d − y − and bideg( b i ) ≤ (deg x ( P ) + µ + iµ ′′ − d − x , d ∗ y − . Proof.
We proceed by induction on i . For i = 0, theclaim follows from Lemma 7. Assume that i > i . For brevity, weset γ = deg x ( P ) + µ , F i − = B i − / ( δ i δ ′ i − Q ∗ i − Q − ), and G i − = b i − / ( δ i δ ′ i − Q ∗ ). The induction hypothesis implies D ix ( f ) = D y D x ( F i − ) + D x ( G i − ) , with bidegree bounds on B i − and b i − . Fact 3 (i) impliesthat ˜ Q := Q ∗ D x ( Q − ) /Q − is in k [ x, y ], with bideg( ˜ Q ) ≤ Algorithm
HermiteTelescoping ( f ) Input : f = P/Q ∈ k ( x, y ) satisfying Hypothesis (H). Output : A minimal telescoper L ∈ k [ x ] h D x i with cer-tificate g ∈ k ( x, y ).1. Apply HermiteEvalInterp to f to get ( g , a ) suchthat f = D y ( g ) + a /Q ∗ . If a = 0, return (1 , g ).2. For i from 1 to deg y ( Q ∗ ) do(a) Apply HermiteEvalInterp to − a i − D x ( Q ∗ ) /Q ∗ to express it as D y (˜ g i ) + ˜ a i /Q ∗ .(b) Set g i = D x ( g i − ) + ˜ g i and a i = D x ( a i − ) + ˜ a i .(c) Solve P ij =0 η j a j = 0 for η j ∈ k ( x ) using [19].If there exists a nontrivial solution, then set( L, g ) := (cid:0)P ij =0 η j D jx , P ij =0 η j g j (cid:1) , and break.3. Compute the content c of L and return ( c − L, c − g ). Figure 3: Creative telescoping by Hermite reduction ( d ∗ x − , d ∗ y ). Hence D x (1 /Q − ) = − ˜ Q/Q . This observationand an easy calculation imply that D x ( F i − ) = ˜ B i − δ i +1 δ ′ i Q ∗ i Q − , where ˜ B i − ∈ k [ x, y ] and deg x ( ˜ B i − ) ≤ deg x ( B i − ) + µ ′′ + d ∗ x . Furthermore, by Lemma 16 there are ¯ B i , ¯ b i ∈ k [ x, y ]with bidegrees at most (deg x ( b i − ) + µ ′′ , d ∗ y − D x ( G i − ) = D y (cid:18) ¯ B i δ i +1 δ ′ i Q ∗ (cid:19) + ¯ b i δ i +1 δ ′ i Q ∗ . Setting B i = ˜ B i − + ¯ B i Q ∗ i − Q − and b i = ¯ b i , we arriveat (10). It remains to verify the degree bounds. The induc-tion hypothesis implies that both deg x ( ¯ B i ) and deg x ( b i ) arebounded by γ + iµ ′′ − d − x . It follows that deg x ( ¯ B i Q ∗ i − Q − )is bounded by γ + iµ ′′ + ( i − d ∗ x . Similarly, deg x ( ˜ B i − ) isbounded by γ + iµ ′′ + ( i − d ∗ x , and so is deg x ( B i ). Thebounds on degrees in y are obvious.We next derive degree bounds for the minimal telescop-ers obtained at an intermediate stage of HermiteTelescoping ;refined bounds on the output will be given by Theorem 25.
Lemma 18
Under (H’), Step 2(c) of Algorithm
Hermite-Telescoping computes a minimal telescoper L ∈ k [ x ] h D x i with order ρ and a certificate g ∈ k ( x, y ) for P/Q with deg x ( L ) ∈ O ( d x d y ρ ) and bideg( g ) ∈ O ( d x d y ρ ) × O ( d y ρ ) . Proof.
By Lemma 13, we exhibit a minimal telescoperby considering the first nontrivial linear dependence amongthe a i ’s in (10). Let M be the coefficient matrix of thesystem in ( η i ) obtained from P ρi =0 η i a i = 0. By Lemma 17, M is of size at most ( ρ + 1) × d ∗ y and with coefficients ofdegree at most σ := d x + µ + ρµ ′′ − d − x in x . Hence, thereexists a solution ( η , . . . , η ρ ) ∈ k [ x ] ρ +1 of degree at most σρ in x by Fact 2 (ii) . Since µ, µ ′′ ∈ O ( d x d y ) and d ∗ y ≤ d y , thedegree estimates of L and g are as announced. We proceed to analyse the complexity of the algorithm inFigure 3 and of an optimisation. heorem 19
Under Hyp. (H’), Algorithm
HermiteTelescop-ing in Figure 3 is correct and takes ˜ O ( ρ ω +1 d x d y ) ops, where ρ is the order of the minimal telescoper. Proof.
The formulas in Step 2(a) create the loop invari-ant D ix ( f ) = D y ( g i ) + a i /Q ∗ . Correctness then follows fromLemmas 12 and 20. Step 1 takes ˜ O ( d x d y ) ops by Theo-rem 11 under (H’). By Lemma 17, deg x ( − a i − D x ( Q ∗ )) ∈O ( id x d y ). So the cost for performing Hermite reduction on − a i − D x ( Q ∗ ) /Q ∗ in Step 2(a) is ˜ O ( id x d y ) ops by Theo-rem 11. The bidegrees of g i and a i in Step 2(b) are in O ( id x d y ) × O ( id y ) by Lemma 17. Since adding and differ-entiating have linear complexity, Step 2(b) takes ˜ O ( i d x d y )ops. For each i , the coefficient matrix of P ij =0 η j a j = 0in Step 2(c) is of size at most ( i + 1) × d ∗ y and with coef-ficients of degree at most deg x ( a i ) ∈ O ( id x d y ). Moreover,the rank of this matrix is either i or i + 1. Then, Step 2(c)takes ˜ O ( i ω d x d y ) ops by Fact 2 (iii) . Computing the contentand divisions in Step 3 has complexity ˜ O ( d x d y ρ ). If thealgorithm returns when i = ρ , then the total cost is in ρ X i =0 ˜ O ( i d x d y ) + ρ X i =1 ˜ O ( i ω d x d y ) ⊂ ˜ O ( ρ ω +1 d x d y ) ops , (11)which is as announced.An optimisation, based on Lemma 15, consists in guessingthe order ρ so as to perform Step 2(c) a few times only: Asa preprocessing step, choose x ∈ k lucky for Q , then detectlinear dependence of { r , , . . . , r ,j } in (9). The minimal j for dependence is a lower bound ρ on ρ . So Step 2(c) isthen performed only when i ≥ ρ . In practice, the lowerbound ρ computed in this way almost always coincides withthe actual order ρ . So normalising the g i ’s becomes thedominant step, as observed in experiments. We analyse thisoptimisation by first estimating the cost for computing ρ . Lemma 20
Under Hypothesis (H’), computing a lower or-der bound ρ for minimal telescopers takes ˜ O ( d x d y ρ ) ops. Proof.
Since differentiating has linear complexity, thederivative D ix ( f ) takes ˜ O ( i d x d y ) ops. By Fact 1 (i) , theevaluation D ix ( f )( x , y ) takes as much. The cost of Hermitereduction on D ix ( f )( x , y ) is ˜ O ( id y ) ops by Lemma 5. ByFact 2 (iii) with d = 1, computing the rank of the coefficientmatrix of P ij =0 η j r ,j , with r ,j as in (9), takes ˜ O ( d y i ω − )ops. Thus, the total cost for computing a lower bound on ρ is P ρ i =0 ˜ O ( i d x d y ) ∈ ˜ O ( d x d y ρ ) ops. Corollary 21
For runs such that ρ = ρ − O (1) , the previ-ous optimisation of HermiteTelescoping takes ˜ O ( ρ d x d y ) ops. Proof.
In view of Lemma 20, the estimate (11) becomes˜ O ( d x d y ρ ) + P ρi =0 ˜ O ( i d x d y ) + P ρi = ρ ˜ O ( i ω d x d y ), which is˜ O ( ρ d x d y ) + ˜ O (( ρ − ρ ) ρ ω d x d y ) ops, whence the result. We analyse the complexity of Almkvist and Zeilberger’salgorithm [1] when restricted to bivariate rational functions.In order to get a telescoper whose order ρ is minimal, the re-sulting algorithm, denoted RatAZ , solves (1) for increasing,prescribed values of ρ until it gets a solution ( η , . . . , η ρ , g ) ∈ k ( x ) ρ +1 × k ( x, y ) with the η i ’s not all zero. For the analysis,we start by studying the parameterisation of the differentialGosper algorithm of [1] under the same restriction to k ( x, y ). Definition ([9])
Let K be a field and a, b ∈ K [ y ] be non-zero polynomials. A triple ( p, q, r ) ∈ K [ y ] is said to be a differential Gosper form of the rational function a/b if ab = D y ( p ) p + qr and gcd( r, q − τ D y ( r )) = 1 for all τ ∈ N . For hyperexponential f , a key step in [1] is to compute adifferential Gosper form of the logarithmic derivative of F = P ρi =0 η i D ix ( f ), where the η i ’s are undetermined from k ( x ).In the analogue RatAZ , this form is predicted by Lemma 22below, which is a technical generalisation of a result byLe [12] on F when f has a squarefree denominator.Write Q = t ( y ) T ( x, y ), splitting content and primitivepart w.r.t. x . By an easy induction, D ix ( f ) = N i / ( QT ∗ i )for N i ∈ k [ x, y ]. For this section, set F = P ρi =0 η i D ix ( f ), N = P ρi =0 η i N i T ∗ ρ − i , and H = − D y ( Q ) /Q − − ρt ∗ D y ( T ∗ ). Lemma 22 If F is nonzero, the triple ( N, H, Q ∗ ) is a dif-ferential Gosper form of D y ( F ) /F . Proof.
First, observe F = N/ ( QT ∗ ρ ) and Q ∗ = t ∗ T ∗ .Next, D y ( F ) /F = D y ( N ) /N − D y ( Q ) /Q − ρD y ( T ∗ ) /T ∗ is D y ( N ) /N + H/Q ∗ . There remains to prove gcd( Q ∗ , H − τ D y ( Q ∗ )) = 1, for any τ ∈ N . Recall that the squarefreepart Q ∗ of Q is the product Q Q · · · Q m and that ˆ Q i de-notes Q ∗ /Q i . By Fact 3 (ii) , Z := H − τ D y ( Q ∗ ) = − ρt ∗ D y ( T ∗ ) − m X i =1 ( i + τ ) ˆ Q i D y ( Q i ) . If Q j divides t ∗ , Z reduces to − ( j + τ ) ˆ Q j D y ( Q j ) modulo Q j .If not, it reduces to − ( j + τ ) ˆ Q j D y ( Q j ) − ρt ∗ ( D y ( Q j ) T ∗ /Q j ),which rewrites to − ( j + τ + ρ ) ˆ Q j D y ( Q j ) modulo Q j . In bothcases, Z is coprime with Q ∗ , as j > τ ≥
0, and ρ ≥ N i ) ≤ (deg x ( P )+ i deg x ( T ∗ ) − i, d y + i deg y ( T ∗ ) − N ) ≤ (deg x ( P ) + ρ deg x ( T ∗ ) − ρ, d y + ρ deg y ( T ∗ ) − RatAZ is, for fixed ρ , to reduce (1) bythe change of unknown g = z/ ( Q − T ∗ ρ ), so as to determineall ( η i ) ∈ k ( x ) ρ +1 for which the differential equation in z ρ X i =0 η i N i T ∗ ρ − i = Q ∗ D y ( z ) + ( D y ( Q ∗ ) + H ) z (12)has a polynomial solution in k ( x )[ y ]. For later use, we recallthe following consequence of [9, Corollary 9.6]. Lemma 23
Let a, b ∈ K [ y ] be such that β = − lc y ( b ) / lc y ( a ) is a nonnegative integer and deg y ( b ) = deg y ( a ) − . Let c ∈ K [ y ] be such that β ≥ deg y ( c ) − deg y ( a ) + 1 . If u is apolynomial solution of aD y ( z ) + bz = c , then deg y ( u ) ≤ β . The following lemma generalises [12, Lemma 2] to present adegree bound for z . Lemma 24 If u ∈ k ( x )[ y ] is a solution of (12) for ( η i ) ∈ k ( x ) ρ +1 , then deg y ( u ) is bounded by β = d − y + ρ deg y ( T ∗ ) . Proof.
Let a = Q ∗ and b = D y ( Q ∗ ) + H . By the defi-nition of H , b = − Q ∗ D y ( Q − ) /Q − − ρt ∗ D y ( T ∗ ). Fact 3 (i) implies that lc y ( b ) = − ( d − y + ρ deg y ( T ∗ )) lc y ( a ). Therefore, β = − lc y ( b ) / lc y ( a ) = d − y + ρ deg y ( T ∗ ). As deg y ( N ) P, P, d − y , − Q ∗ D y ( Q ) /Q );3. For ℓ = 0 , , . . . do(a) Set z to P βj =0 z j y j , extract the linear system M (cid:0) η i z j (cid:1) T = 0 from (12) (for ρ = ℓ ) and com-pute a basis S of the null space of M by [19].(b) If S contains a solution ( η , . . . , η ℓ , s ) such that η , . . . , η ℓ are not all nonzero, then set ( L, g ) := (cid:0)P ℓi =0 η i D ix , s/ ( Q − T ∗ ℓ ) (cid:1) , and go to Step 4;(c) Update ˜ N := D x ( ˜ N ) T ∗ − ˜ N (cid:0) T ∗ D x ( T ) /T + iD x ( T ∗ ) (cid:1) , N := NT ∗ + η ℓ +1 ˜ N , β := β +deg y ( T ∗ ), and H := H − t ∗ D y ( T ∗ ).4. Compute the content c of L and return ( c − L, c − g ). Figure 4: Improved Almkvist–Zeilberger algorithm We end the present section using the approach of Almkvistand Zeilberger to provide tight degree bounds on the outputsfrom Algorithms HermiteTelescoping and RatAZ . Theorem 25 Under Hypothesis (H’), there exists a mini-mal telescoper L ∈ k [ x ] h D x i with certificate g ∈ k ( x, y ) with deg x ( L ) ∈ O ( d x d y d ∗ y ) and bideg( g ) ∈ O ( d x d y d ∗ y ) ×O ( d y d ∗ y ) . Proof. By Corollary 14, there exists a smallest ρ ∈ N at most d ∗ y , for which (1) has a solution with the η i ’s notall zero. For this ρ , we estimate the size of the polynomialmatrix M derived from (12) by undetermined coefficients.By the remark on N after Lemma 22, we have bideg( N ) ≤ ( n x , n y ) where n x := d x + ρ deg x ( T ∗ ) − ρ ∈ O ( ρd x ) and n y := d y + ρ deg y ( T ∗ ) − ∈ O ( ρd y ). The matrix M contains twoblocks M ∈ k [ x ] ( n y +1) × ( ρ +1) ≤ n x and M ∈ k [ x ] ( n y +1) × ( β +1) ≤ d x ,where β ∈ O ( ρd y ) is the same as in Lemma 24. By theminimality of ρ , the dimension of the null space of M is 1.So there exists u ∈ k [ x ] n y +1 with coefficients of degree atmost n x ( ρ + 1) + d x ( β + 1) ∈ O ( d x d y d ∗ y ) in x such that M (cid:0) η z (cid:1) T = 0, which implies degree bounds in x for L and g .The degree bound in y for g is obvious.We now analyse the complexity of the algorithm in Fig. 4. Theorem 26 Under Hypothesis (H’), Algorithm RatAZ inFigure 4 is correct and takes ˜ O ( d x d ωy ρ ω +2 ) ops, where ρ isthe order of the minimal telescoper. Proof. By the existence of a telescoper, Corollary 14,and Lemma 24, the algorithm always terminates and returnsa minimal telescoper L , of order ρ at most d ∗ y . Gcd computa-tions dominate the cost of Steps 1 and 2, which take ˜ O ( d x d y )ops. For each ℓ ∈ N , the dominating cost in Step 3 is com-puting the null space of M . Let n y = d y + ℓ deg y ( T ∗ ) − ∈O ( ℓd y ) and n x = d x + ℓ deg x ( T ∗ ) ∈ O ( ℓd x ). By the sameargument as in the proof of Theorem 25, the matrix M is of size at most ( n y + 1) × ( ℓ + β + 2) and with coeffi-cients of degree at most n x . Let r be the rank of M , whichis either ℓ + β + 2 or ℓ + β + 1 by construction. Thus,a basis of the null space of M can be computed within˜ O ( n x ( n y + 1)( ℓ + β + 2) r ω − ) ops by Fact 2 (iii) . Since β ∈ O ( ℓd y ), ˜ O ( n x ( n y + 1)( ℓ + β + 2) r ω − ) is included in˜ O ( d x d ωy ℓ ω +1 ). Since Step 3 terminates at ℓ = ρ , the totalcost of the algorithm is P ρℓ =0 d x d ωy ℓ ω +1 ops. This is withinthe announced complexity, ˜ O ( d x d ωy ρ ω +2 ) ops. Corollary 27 Algorithms HermiteTelescoping and RatAZ inFig. 3 and 4 both output the primitive minimal telescoper L together with its certificate g , which satisfy deg D x ( L ) ≤ d ∗ y , deg x ( L ) , deg x ( g ) ∈ O ( d x d y d ∗ y ) , and deg y ( g ) ∈ O ( d y d ∗ y ) . Proof. Both algorithms output the primitive minimaltelescoper, as they compute a minimal telescoper at an inter-mediate step, and owing to their last step of content removal.Bounds follow from Corollary 14 and Theorem 25. 4. NONMINIMAL TELESCOPERS Here, we discard Hypothesis (H) and trade the minimalityof telescopers for smaller total output sizes. To this end, weadapt and slightly extend the arguments in [13] and [3, § f = P/Q ∈ k ( x, y ) of bidegree ( d x , d y ), our goalis to find a (possibly nonminimal) telescoper for f . It issufficient to find a nonzero differential operator A ( x, D x , D y )that annihilates f . Indeed, any A ∈ k [ x ] h D x , D y i \ { } suchthat A ( f ) = 0 can be written A = D ry ( L + D y R ), where L is nonzero in k [ x ] h D x i and R ∈ k [ x ] h D x , D y i . If r = 0, thenclearly L is a telescoper for f ; otherwise, A ( f ) = 0 yields L ( f ) = D y ( − R ( f ) − P r − i =0 a i i +1 y i +1 ) for some a i ∈ k ( x ),which implies that L is again a telescoper for f . Moreover, inboth cases, deg x ( L ) ≤ deg x ( A ) and deg D x ( L ) ≤ deg D x ( A ).Furthermore, for any ( i, j, ℓ ) ∈ N , a direct calculation yields x i D jx D ℓy ( f ) = H i,j,ℓ Q j + ℓ +1 , (13)where H i,j,ℓ ∈ k [ x, y ] and deg x ( H i,j,ℓ ) ≤ ( j + ℓ + 1) d x + i − j and deg y ( H i,j,ℓ ) ≤ ( j + ℓ + 1) d y − ℓ . From these inequali-ties, we derive the size and complexity estimates in Figure 1(bottom half), using two different filtrations of k [ x ] h D x , D y i . Lipshitz’s filtration ([13]). Let F ν be the k -vector spaceof dimension f ν := (cid:0) ν +33 (cid:1) spanned by { x i D jx D ℓy | i + j + ℓ ≤ ν } . By (13), F ν ( f ) is contained in the vector space of di-mension g ν := (( ν + 1) d x + ν + 1) (( ν + 1) d y + 1) spannedby (cid:8) x i y j Q ν +1 | i ≤ ( ν + 1) d x + ν, j ≤ ( ν + 1) d y (cid:9) . Choos-ing ν = 6( d x + 1)( d y + 1) yields f ν > g ν ; therefore, thereexists A in k h x, D x , D y i \ { } with total degree at most6( d x +1)( d y +1) in x , D x , and D y that annihilates f . More-over, A is found by linear algebra in dimension O (( d x d y ) ). A better filtration ([3]). Instead of taking total de-gree, set F κ,ν to the k -vector space of dimension f κ,ν :=( κ + 1) (cid:0) ν +22 (cid:1) generated by { x i D jx D ℓy | i ≤ κ, j + ℓ ≤ ν } .By (13), F κ,ν ( f ) is contained in the vector space of dimen-sion g κ,ν := (( ν + 1) d x + κ + 1)(( ν + 1) d y + 1) spanned by (cid:8) x i y j Q ν +1 | i ≤ ( ν + 1) d x + κ, j ≤ ( ν + 1) d y (cid:9) . Choosing κ = 3 d x d y and ν = 6 d y results in f κ,ν > g κ,ν . This impliesthe existence of A in k h x, D x , D y i \ { } with total degree atmost 6 d y in D x and D y and degree at most 3 d x d y in x thatannihilates f . Again, A is found by linear algebra over k ,but in smaller dimension O ( d x d y ). o. AZ Abr RAZ H1 H2 HO EI MG 29 44 72 32 28 36 20 608 52843 52 76 36 20 24 32 652 58446 4268 1436 784 492 1288 752 343413 1894549 474269 34694 20977 10336 36254 22417 ∞ Table 1: Creative telescoping on random instances Timings in ms for algorithms in Table 3 (stopped after 30 min). 5. IMPLEMENTATION AND TIMINGS We implemented in Maple 13 all the algorithms described;as we used Maple’s generic solver SolveTools:-Linear , allof our implementations are deterministic.The evaluation-interpolation algorithm HermiteEvalInterp for Hermite reduction (Fig. 2) does not perform well, mainlybecause Maple’s rational interpolation routines are far tooslow. We thus implemented Algorithm HermiteReduce (orig-inal version) in [4, § HermiteTe-lescoping in Figure 3, using HermiteReduce in place of Her-miteEvalInterp , and including the optimisation at the end of § HermiteTelescoping re-turns the minimal telescoper L and the certificate g . Thealgorithm separates the computation for L from that for g .Indeed, g is formed by the coefficients of L , g , the ˜ g i andtheir derivatives given in Figure 3. This feature enables usto either return the certificate g as a sum of unnormalisedrational functions, or a normalised rational function.A selection of timings by this implementation and oth-ers are given in Table 1; our code, the full table, as well asthe random inputs are given in [22]. For our experiments,we exhaustively considered all 49 bidegree patterns in fac-torisations of denominators Q · · · Q mm ( m ≤ 5) that add upto bidegree (5,5), and generated corresponding random de-nominators, imposing the integers of the expanded forms tohave around 26 digits. Numerators were generated as ran-dom bidegree-(5,5) polynomials with coefficients of 26 digits. Application to diagonals. The diagonal of a formal powerseries f = P i,j ≥ f i,j x i y j in k [[ x, y ]] is defined to be thepower series ∆( f ) := P ∞ i =0 f i,i x i . For a D-finite power se-ries f , it is known to be D-finite [13], and it is even algebraicfor a bivariate rational function f ∈ k ( x, y ) ∩ k [[ x, y ]] [18, § L ∈ k ( x ) h D x i that an-nihilates ∆( f ) can then be computed via rational-functiontelescoping, owing to the following classical lemma from [13]. Lemma 28 Any telescoper for f ( y, xy ) /y annihilates ∆( f ) . By this lemma, it suffices to compute a telescoper without itscertificate to get an annihilator. Algorithm HermiteTelescop-ing is suitable for this task, since it separates computation oftelescopers and certificates. Alternatively, for f = P/Q , wecan compute an annihilator of ∆( f ) either as the differen-tial resolvent of the resultant Res y ( Q, P − τ D y Q ), or simply guess it from the first terms of the series expansion of ∆( f ).We compare the various algorithms on an example bor-rowed from [5] (timings of execution are given in Table 2): f = 11 − x − y − xy (1 − x d ) , where d ∈ N . (14)All computer calculations have been performed on a Quad-Core Intel Xeon X5482 processor at 3.20GHz, with 3GB ofRAM, using up to 6.5GB of memory allocated by Maple. d AZ Abr RAZ H1 H2 HO RR GHP ∞ Table 2: Computation of the diagonals of (14) Timings in ms by creative telescoping of f ( y, x/y ) /y (upper half)or f ( y/x, x ) /x (second half). Algorithms listed in Table 3. AZ DETools[Zeilberger]Abr AZ with Abramov’s denominator bound by option gosper_freeRAZ Algorithm RatAZ of Fig. 4, with lower-bound prediction H1 our Hermite-based approach, without certificate normalisation H2 H1 , but with normalised certificate HO RAZ , solving (1) by Horowitz–Ostrogradsky EI H1 with evaluation and interpolation for calculations over k ( x ) MG Mgfun ’s creative telescoping for general D-finite functions RR telescoper computation by resultant and differential resolvent GHP telescoper guessing by diagonal expansion and Hermite–Pad´e Table 3: List of the algorithms for the experiments 6. REFERENCES [1] G. Almkvist and D. Zeilberger. The method of differentiatingunder the integral sign. J. Symb. Comput. , 10:571–591, 1990.[2] M. Apagodu and D. Zeilberger. Multi-variable Zeilberger andAlmkvist-Zeilberger algorithms and the sharpening of Wilf-Zeilberger theory. Adv. in Appl. Math. , 37(2):139–152, 2006.[3] A. Bostan, F. Chyzak, B. Salvy, G. Lecerf, and ´E. Schost.Differential equations for algebraic functions. In ISSAC’07 ,pages 25–32. ACM, New York, 2007.[4] M. Bronstein. Symbolic Integration I: Transcendentalfunctions , volume 1 of Algorithms and Computation inMathematics . Springer-Verlag, Berlin, second edition, 2005.[5] A. Flaxman, A. W. Harrow, and G. B. Sorkin. Strings withmaximally many distinct subsequences and substrings. Electron. J. Combin. , 11(1):R8, 10 pp., 2004.[6] J. von zur Gathen and J. Gerhard. Modern Computer Algebra .Cambridge University Press, Cambridge, second edition, 2003.[7] K. O. Geddes and H. Q. Le. An algorithm to compute theminimal telescopers for rational functions (differential-integralcase). In Mathematical Software , pages 453–463. WSP, 2002.[8] J. Gerhard. Fast modular algorithms for squarefreefactorization and Hermite integration. Appl. Algebra Engrg.Comm. Comput. , 11(3):203–226, 2001.[9] J. Gerhard. Modular Algorithms in Symbolic Summation andSymbolic Integration (LNCS) . SpringerVerlag, 2004.[10] C. Hermite. Sur l’int´egration des fractions rationnelles. Ann.Sci. ´Ecole Norm. Sup. (2) , 1:215–218, 1872.[11] E. Horowitz. Algorithms for partial fraction decomposition andrational function integration. In SYMSAC’71 , pages 441–457,New York, USA, 1971. ACM.[12] H. Q. Le. On the differential-integral analogue of Zeilberger’salgorithm to rational functions. In Proc. of the 2000 AsianSymposium on Computer Mathematics , pages 204–213, 2000.[13] L. Lipshitz. The diagonal of a D-finite power series is D-finite. J. Algebra , 113(2):373–378, 1988.[14] L. Lipshitz. D-finite power series. J. Algebra , 122(2):353–373,1989.[15] M. Ostrogradsky. De l’int´egration des fractions rationnelles. Bull. de la classe physico-math´ematique de l’Acad. Imp´erialedes Sciences de Saint-P´etersbourg , 4:145–167, 286–300, 1845.[16] R. H. Risch. The problem of integration in finite terms. Trans.Amer. Math. Soc. , 139:167–189, 1969.[17] R. H. Risch. The solution of the problem of integration in finiteterms. Bull. Amer. Math. Soc. , 76:605–608, 1970.[18] R. P. Stanley. Enumerative Combinatorics. Vol. 2 , volume 62of Cambridge Studies in Advanced Mathematics . CUP, 1999.[19] A. Storjohann and G. Villard. Computing the rank and a smallnullspace basis of a polynomial matrix. In ISSAC’05 , pages309–316. ACM, New York, 2005.[20] D. Zeilberger. A holonomic systems approach to specialfunctions identities. J. Comput. Appl. Math. , 32:321–368, 1990.[21] D. Zeilberger. The method of creative telescoping. J. SymbolicComput. , 11(3):195–204, 1991.[22] http://algo.inria.fr/chen/BivRatCT/http://algo.inria.fr/chen/BivRatCT/