THÈSE

présentée à l’École Polytechnique

Pour obtenir le titre de Docteur en Sciences

Spécialités : Informatique

Shaoshi Chen

Quelques applications de l’algèbre différentielle et

aux différences pour le télescopage créatif

Soutenue le 16 février 2011 devant le jury composé de :

William Y. C.Chen,  Rapporteur,
Frédéric Chyzak,  Directeur de thèse,
Joris van der Hoeven,  Examinateur,
Manuel Kauers,  Rapporteur,
Ziming Li,  Directeur de thèse,
François Morain,  Examinateur,
Marko Petkovšek,  Rapporteur,
Felix Ulmer,  Examinateur.

2

Doctoral Thesis

presented to l’École Polytechnique

For the French degree of Doctor in Sciences

Specialty: Computer Science

Shaoshi Chen

Some applications of differential-difference algebra

to creative telescoping

Defended on February 16, 2011 in front of the committee consisting of:

William Y. C.Chen,  Reviewer,
Frédéric Chyzak,  Thesis advisor,
Joris van der Hoeven,  Examiner,
Manuel Kauers,  Reviewer,
Ziming Li,  Thesis advisor,
François Morain,  Examiner,
Marko Petkovšek,  Reviewer,
Felix Ulmer,  Examiner.

0

Abstract

Since the 1990’s, Zeilberger’s method of creative telescoping has played an important role in the automatic verification of special-function identities. The long-term goal initiated in this work is to obtain fast algorithms and implementations for definite integration and summation in the framework of this method. Our contributions include new practical algorithms, complexity analyses of algorithms, and theoretical criteria for the termination of existing algorithms.

On the practical side, we present a new algorithm for computing minimal telescopers for bivariate rational functions. This algorithm is based on Hermite reduction. We also improve the classical Almkvist and Zeilberger’s algorithm for rational-function inputs. The Hermitereduction based algorithm and improved Almkvist and Zeilberger’s algorithm are analyzed in terms of field operations. Both complexity analysis and experimental results show that our algorithms are superior to other known ones for rational-function inputs.

On the theoretic side, we present a structure theorem for multivariate hyperexponentialhypergeometric functions. This theorem is based on (multivariate) Christopher’s theorem for hyperexponential functions, the Ore-Sato theorem for hypergeometric terms, and our generalization of a recent result by Feng, Singer, and Wu on compatible bivariate continuous-discrete rational functions. The structure theorem allows us to decompose a hyperexponential-hypergeometric function as a product of a rational function, several exponential and power functions, and factorial terms. Furthermore, we derive two criteria for the existence of telescopers for bivariate hyperexponential-hypergeometric functions: one is with respect to the continuous variable, and the other with respect to the discrete one. The two criteria solve the termination problem of the continuous-discrete analogue of Zeilberger’s algorithms.

2

Acknowledgements

I would like to express my gratitude to my two co-supervisors: Frédéric Chyzak and Ziming Li, for their academic guidance and constant supports. Thanks go to Frédéric Chyzak for his motivated lessons on mathematics, English, and French with great patience. Thanks are due to Ziming Li for his long-term seminar in which I could report my understanding and misunderstanding of each sentence of Bronstein’s book Symbolic Integration I.

Special thanks go to Alin Bostan from whom I learnt the precise knowledge on complexity analysis. He spent a lot of time in reading the draft of my thesis and offered me precise and constructive comments. Without him, I could not set up the “teagroup" in the Algorithms team. I am deeply grateful to Manuel Kauers for his constant encouragement and for vivid talks. He also impressed me with his wide interest and strong ability in communicating with others.

I also want to thank Philippe Flajolet and Bruno Salvy for running the Algorithms team in an atmosphere of freedom. I thank Philippe Flajolet for all the croissants. Bruno Salvy is always willing to listen to my questions and present me great answers, which has motivated my further work. I really appreciate all of my colleagues of the Algorithms team for their assistances, particularly, our secretary Virginie Collette for helping me with dealing with all kinds of problems in my daily life, Marc Mezzarobba for answering my questions on Linux, and Alexandre Benoit for moving my books from Roquencourt to Orsay and back.

I appreciate many members of the Key Laboratory of Mathematics Mechanization for their help, particularly, Wentsün Wu for his great influence on mathematics-mechanization, Lihong Zhi for making the meeting between my wife and me possible, and Xiao-shan Gao for recommending me to join the PhD Research Program between Algorithms and KLMM. I thank Ruyong Feng, Guofeng Fu, Chunming Yuan and many other “lab-mates" for their assistances.

I thank Peter Paule for inviting me to visit RISC-Linz where I presented my first seminar talk and made a lot of new friends. My thanks also go to Moulay Barkatou for inviting my wife and me to visit his university where I gave my first conference talk. I am very grateful to Bill Chen and Qinghu Hou for offering me many opportunities to present my work in Nankai University. My thanks also goes to Zhengyi Lu for many encouragements during my PhD study.

I thank George Andrews for sending me Payne’s thesis as well as Garth Payne for the updated version. Thanks go to Doron Zeilberger for a motivated discussion about his method and for his kind recommendation letter for my postdoctoral application. I am also very grateful to Michael Singer for offering me a postdoctoral position in the North Carolina State University.

I would like to express my gratitude to Bill Chen, Manuel Kauers, and Marko Petkovšek for agreeing to be the reviewers of this thesis and their valuable and constructive comments. My thanks also go to Joris van der Hoeven, François Morain, and Felix Ulmer for agreeing to serve as the examiners of this thesis.

Last but not least, I would like to thank my parents for their education and unconditional support, my uncles and their families in Paris for the hospitality, and my wife for her love, understanding and appreciation of my work.

This work was partially supported by INRIA Paris-Rocquencourt and KLMM. Two grants of the National Natural Science Foundation of China (no. 60821002/F02 and no. 76596100) financially supported my graduate study in Beijing and several international trips.

Contents

1 Introduction 9 1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 1.2 Main results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 1.3 Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 1.4 Notation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

2 Preliminaries 15 2.1 Differential and difference rings . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.2 Ore polynomials . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.3 Integration of rational functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

2.3.1 Hermite reduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.3.2 Ostrogradsky and Horowitz’s method . . . . . . . . . . . . . . . . . . . . . 21 2.3.3 Residues and Rothstein–Trager resultants . . . . . . . . . . . . . . . . . . 22

2.4 Background on complexity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

3 Hermite Reduction for Rational-Function Telescoping 27 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 3.2 Notation and hypotheses . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 3.3 Hermite reduction for bivariate rational functions . . . . . . . . . . . . . . . . . . 30

3.3.1 Output size estimates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 3.3.2 Algorithm by evaluation and interpolation . . . . . . . . . . . . . . . . . . 32

3.4 Minimal telescopers for bivariate rational functions . . . . . . . . . . . . . . . . . 34 3.4.1 Hermite-reduction based method . . . . . . . . . . . . . . . . . . . . . . . 34

3.4.2 Improved Almkvist and Zeilberger’s method . . . . . . . . . . . . . . . . . 39 3.5 Non-minimal telescopers for bivariate rational functions . . . . . . . . . . . . . . 43 3.6 Implementation and experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . 45

3.6.1 Implementation and examples . . . . . . . . . . . . . . . . . . . . . . . . . 45 3.6.2 Experimental results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47

4 Structure of Multivariate Hyperexponential-Hypergeometric Functions 51 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 4.2 Algebraic setting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52

4.2.1 Hyperexponential-hypergeometric functions . . . . . . . . . . . . . . . . . 52 4.2.2 First-order fully integrable systems . . . . . . . . . . . . . . . . . . . . . . 54

4.3 Rational normal forms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 4.3.1 Differential and shift rational normal forms . . . . . . . . . . . . . . . . . 56 4.3.2 𝕐-rational normal forms . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58

4.4 Structure of compatible rational functions . . . . . . . . . . . . . . . . . . . . . . 63 4.4.1 The Ore-Sato theorem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 4.4.2 Alternative proof of multivariate Christopher’s theorem . . . . . . . . . . 65 4.4.3 Multivariate extension of Feng-Singer-Wu’s lemma . . . . . . . . . . . . . 68

4.5 Multiplicative structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72

5 Existence of Telescopers for Hyperexponential-Hypergeometric Functions 75 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 5.2 Preliminaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77

5.2.1 Ring of sequences . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 5.2.2 Ring of differential-recurrence operators . . . . . . . . . . . . . . . . . . . 81 5.2.3 Split polynomials . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81

5.3 Existence problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 5.4 Standard representations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 5.5 Two additive decompositions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87

5.5.1 Additive decomposition with respect to x . . . . . . . . . . . . . . . . . . 87 5.5.2 Additive decomposition with respect to n . . . . . . . . . . . . . . . . . . 94

5.5.3 Applying operators to additive decompositions . . . . . . . . . . . . . . . 98 5.6 Two existence criteria . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100

5.6.1 Existence of telescopers implies properness . . . . . . . . . . . . . . . . . . 100 5.6.2 Properness implies the existence of telescopers . . . . . . . . . . . . . . . . 103

5.7 Algorithms and examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108 5.7.1 Algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108 5.7.2 Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109

6 Conclusion and Perspectives 111 6.1 Extensions of the Hermite-reduction based method . . . . . . . . . . . . . . . . . 111 6.2 Existence criteria for telescopers in the q-setting . . . . . . . . . . . . . . . . . . . 113 6.3 Wilf and Zeilberger’s conjecture . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 8

Chapter 1

Introduction

As most special-function integrals and sums cannot be expressed in closed form, their evaluation cannot be based on table look-ups only. Even when closed forms are available, they may prove to be intractable in further manipulations. In both cases, the difficulty can be mitigated by representing functions by annihilating differential and difference operators. This motivated Zeilberger to introduce a method now known as creative telescoping [95], which applies to a large class of special functions. Zeilberger’s method has been used extensively in the automatic proofs of special-function identities.

Let us illustrate the basic idea of Zeilberger’s method by proving the identity

∑_{n=0}^{+∞} f(x, n) = \frac{1}{\sqrt{1 − 4x}}, where f(x, n) = \binom{2n}{n}x^{n} and x ∈ [0, 1/4).  (1.1)

Let D_{x} and Δ_{n} denote the derivation with respect to x, and the difference operator with respect to n, respectively. The first step is to find a differential equation satisfied by the sum on the left-hand side of (1.1). To this end, we try to find a nonzero operator L ∈ ℚ(x)⟨D_{x}⟩ such that

L(x, D_{x})(f) = Δ_{n}(g),  (1.2)

where g is a function in x and n. Moreover, the method requires that the ratio between g and f is a rational function. For this specific example, we have

L = 2 − (1 − 4x)D_{x} and g = \frac{n}{x} · f = nx^{n−1}\binom{2n}{n}.

Note that taking sums with respect to n commutes with the application of L(x, D_{x}) and, additionally, that

∑_{n=0}^{+∞} Δ_{n}(g) = 0.

This implies that the sum on the left-hand side of (1.1) satisfies

2y(x) − (1 − 4x)\frac{dy(x)}{dx} = 0.

It is easy to verify that the function on the right-hand side of (1.1) also satisfies the above differential equation. Moreover, the identity holds when x = 0. By Cauchy’s theorem, the identity (1.1) holds for all well-defined values x ∈ [0, 1/4). The operator L in (1.2) is called a telescoper for f with respect to n, and g the certificate of L for f.

This thesis focuses on the bivariate case, in which there are four kinds of telescoping problems related to different integration and summation problems.

∫ / ∑ problems Telescoping equations
(P1) ∫_{a}^{b} f(x, y) dy L(x, D_{x})(f) = D_{y}(g)
(P2) ∫_{a}^{b} f(n, y) dy L(n, S_{n})(f) = D_{y}(g)
(P3) ∑_{m=a}^{b} f(m, x) L(x, D_{x})(f) = Δ_{m}(g)
(P4) ∑_{m=a}^{b} f(n, m) L(n, S_{n})(f) = Δ_{m}(g)

Table 1.1: Four telescoping problems

1.1 Motivation

Since the 1990’s, Zeilberger’s method has been extensively studied in the literature [95, 96, 15, 92, 63, 6, 7, 65, 28, 16]. The main focus of those studies is on the efficiency and termination of creative-telescoping algorithms.

Related to efficiency, progress and improvements have since been made by Abramov [7], Apagodu [16], Chyzak [35, 33], Le [63, 46, 65], Takayama [85], etc. However, very little is known about the complexity of creative-telescoping algorithms. We believe that complexity analysis helps us understand algorithms, and indicates good ways of improving and implementing them.

This is one of the motivations for Chapter 3. The existing algorithms bind the construction of telescopers with certificates. In certain applications, only telescopers are needed. For bivariate rational-function inputs, our complexity analysis shows that the arithmetic size of certificates is asymptotically larger than that of telescopers. So the complexity could be lower if one can avoid the calculations for certificates in certain applications. This motivates us to find a way in which one could choose to compute or not to compute the certificates according to the user requirements.

By the functionality of Zeilberger-style algorithms, their termination is equivalent to the existence of telescopers. In his pioneering work [95], Zeilberger has proved that telescopers exist for holonomic functions. In particular, Wilf and Zeilberger [92] have presented an elementary and constructive proof of the existence of telescopers for proper hyperexponential-hypergeometric functions by basing on ideas of Fasenmyer [40] and Verbaeten [89]. By a proper hyperexponentialhypergeometric function, we mean a function that can be written as a product of a polynomial, exponential functions, power functions, and factorial terms. Properness can be detected from the certificates of a hyperexponential-hypergeometric function. Moreover, Wilf and Zeilberger conjectured in [92, page 585] that a hyperexponential-hypergeometric function is holonomic if and only if it is proper. If this conjecture were verified, then one could algorithmically detect the holonomicity of hyperexponential-hypergeometric functions. In the case of several discrete variables, a slightly modified version of the conjecture has been proved independently by Payne in his Ph.D. thesis [73] and by Abramov and Petkovšek [14]. In particular, the case of two variables has also been shown by Hou [54, 55] and by Abramov and Petkovšek [12]. However, holonomicity is only a sufficient condition for the existence of telescopers. In fact, Chyzak, Kauers, and Salvy [34] have listed certain classes of functions that are not holonomic but still have telescopers. Therefore, the challenge is to find theoretical criteria that enable us to algorithmically detect the existence of telescopers.

In view of the theoretical difficulty, special attention has been focused mainly on the subclass of hyperexponential-hypergeometric functions. In the continuous case, the works by Bernstein [17] and Lipshitz [67] show that every hyperexponential function has a telescoper. This implies that Zeilberger’s algorithm always succeeds on those inputs. However, the situation in other cases turns out to be more involved. In the discrete case and its q-analogue, the first complete solution to the termination problem is Abramov and Le’s criterion [64, 9], which decides whether telescopers exist for a given bivariate rational function in the discrete variables m and n. According to their criterion, the rational function

f = \frac{1}{m^{2} + n^{2}}

has no telescoper. Soon, the criterion was extended to the general case of bivariate hypergeometric terms by Abramov [5, 6]. Basically, Abramov proves that a hypergeometric term can be written as a sum of a hypergeometric-summable term and a proper one if it has a telescoper [6, Theorem 10]. Similar results have been obtained in the general q-shift case by Chen, Hou and Mu [28]. These results are fundamental for predicting the termination of Zeilberger’s algorithm.

A continuous-discrete analogue of Zeilberger’s algorithm was presented by Almkvist and Zeilberger [15]. This analogue has been shown to be very useful in the study of orthogonal polynomials [60, Chapters 10–13]. In this setting, not all hyperexponential-hypergeometric functions have telescopers. For example, we will show in this thesis that the rather simple rational function in the continuous variable x and the discrete variable n

f = \frac{1}{x + n}

has no telescoper with respect to either the continuous variable x or the discrete variable n. Therefore, an existence criterion is also needed in this mixed setting.

1.2 Main results

For a bivariate rational function, we present a new algorithm to compute its minimal telescoper, which is based on Hermite reduction. We also obtain some improvements over the classical algorithm by Almkvist and Zeilberger by extending the idea of Geddes and Le to general rationalfunction inputs. Moreover, we give the first proof of a polynomial complexity (in terms of field operations) for creative telescoping on this specific class of inputs. Our algorithms are proved to be faster concerning both theoretical complexity and actual performance.

Motivated by a conjecture of Wilf and Zeilberger [92, page 585], we study the possible form of a multivariate hyperexponential-hypergeometric function. For such a function, we prove a structure theorem for its certificates, which generalizes a recent result by Feng, Singer, and Wu [42, Proposition 5]. Combining our result with the result on multivariate hyperexponential functions in [29, 97] and the Ore-Sato theorem on multivariate hypergeometric terms, we obtain a structure theorem for multivariate hyperexponential-hypergeometric functions, which says that a multivariate hyperexponential-hypergeometric function can be written as a product of a rational function, an exponential function, power functions and factorial terms.

With the help of the result by Feng, Singer, and Wu [42, Proposition 5], we derive two criteria for the existence of telescopers for bivariate hyperexponential-hypergeometric functions. We show that a hyperexponential-hypergeometric function can be written as a sum of a hypergeometricsummable, resp. hyperexponential-integrable, function and a proper one if it has a telescoper with respect to the discrete, resp. continuous, variable. Our criteria are based on standard representations and the two adapted additive decompositions. With them, we can decide in advance the termination of the continuous-discrete analogue of Zeilberger’s algorithm for bivariate hyperexponential-hypergeometric functions.

1.3 Outline

In this section, we provide the reader with an outline of this thesis.

Chapter 2. We recall basic notation and facts on differential rings, difference rings and Ore polynomials. Our new algorithm for the construction of minimal telescopers will be based on the Hermite reduction for the integration of the rational function. So we review the classical algorithms for rational-function integration, including Hermite reduction, Ostrogradsky–Horowitz’s method, and Rothstein–Trager’s algorithm. We also summarize some complexity results for later use.

Chapter 3. We focus on deriving a fast algorithm for computing minimal telescopers for bivariate rational functions. First, we present an optimal algorithm for the Hermite reduction over k(x) by fast evaluation-interpolation strategy. Second, we present a new algorithm for computing the minimal telescoper of a bivariate rational function, based on a bivariate extension of Hermite reduction. Moreover, some improvements over the classical method by Almkvist and Zeilberger [15] are achieved in this chapter. We also analyze the arithmetic complexity of those algorithms. Third, we adapt and slightly extend the arguments by Lipshitz [67] and by Bostan and others [20] to derive smaller total degree sizes of telescopers. At last, we describe our implementation and show experimental results. Chapter 4. We first present an algebraic setting for hyperexponential-hypergeometric functions. After that, we review various normal forms of rational functions and introduce a new kind of rational normal forms, which enables us to extend the existing results to multivariate continuousdiscrete setting. Our main result is Theorem 4.4.6 on the structure of certificates of a multivariate hyperexponential-hypergeometric function. This result is a generalization of a result by Feng, Singer, and Wu [42, Proposition 5]. At last, we describe a multiplicative form of multivariate hyperexponential-hypergeometric functions. Chapter 5. We study the existence of telescopers for hyperexponential-hypergeometric functions of one continuous variable and one discrete variable. First, we review the construction of a ring of sequences from [41, 43], which allows us to study the existence problem in an algebraic manner. After that, we introduce standard representations for bivariate hyperexponential-hypergeometric functions and then adapt two additive decompositions [13, 47] to bivariate hyperexponentialhypergeometric functions described by their standard representations. At last, we describe our existence criteria and its algorithmic description with examples. Chapter 6. We present some conclusions and propose some topics for future work.

1.4 Notation

We shall use the letters ℂ, ℕ, ℚ, ℝ, and ℤ to denote the set of complex numbers, non-negative integers, rational numbers, real numbers, and integers, respectively. For a field F, we denote its algebraic closure by \overline{F}.

Preliminaries

In this chapter, we recall basic notation and facts on differential and difference rings and Ore polynomials from [56, 21, 70]. Moreover, we review classical techniques for rational-function integration and background on complexity analysis.

2.1 Differential and difference rings

Let R be a commutative ring. A map δ : R → R is said to be a derivation on R if

δ(a + b) = δ(a) + δ(b), δ(ab) = δ(a)b + aδ(b) for all a, b ∈ R.

The pair (R, δ) is called a differential ring. Moreover, if R is a field, R is called a differential field. An element c ∈ R is called a constant with respect to δ if δ(c) = 0. All constants with respect to δ form a subring of R, denoted by C_{δ,R}. If R is a field, C_{δ,R} is also a field. Some basic facts concerning derivations are collected in the next lemma whose proof can be found in any book about differential algebra, e.g. [21] Lemma 2.1.1. Let (R, δ) be a differential ring. For a, b ∈ R and n ∈ ℕ, we have

(i) δ(1) = 0;

(ii) δ(a^{n}) = na^{n−1}δ(a);

(iii) If b is invertible, then

δ (\frac{a}{b}) = \frac{δ(a)b − aδ(b)}{b^{2}}.

(iv) Logarithmic derivative identity: if a and b are invertible, then

\frac{δ(a^{m}b^{n})}{a^{m}b^{n}} = m\frac{δ(a)}{a} + n\frac{δ(b)}{b} for all m, n ∈ ℤ.

Let σ : R → R be a monomorphism on R. Then the pair (R, σ) is called a difference ring. Moreover, if R is a field, R is called a difference field. An element c ∈ R is called a constant with respect to σ if σ(c) = c. All constants with respect to σ form a subring of R, denoted by C_{σ,R}. If R is a field, C_{σ,R} is also a field. The triple (R, δ, σ) is called a differential-difference ring. (R, δ, σ) is said to be orthogonal [62] if the derivation δ commutes with the monomorphism σ. Unless otherwise specified, all fields in the thesis are of characteristic zero and all differentialdifference rings are orthogonal.

2.2 Ore polynomials

Ore polynomials are a common abstraction for linear differential operators, linear difference operators, shifts, and q-shifts [70, 23]. Let (F, δ, σ) be a differential-difference field. The ring of Ore polynomials F [x; δ, σ] over F is the univariate polynomial ring with indeterminate x, addition defined as usual, and multiplication being associative, distributive, and satisfying the commutation rule:

xf = σ(f)x + δ(f), for any f ∈ F.  (2.1)

By the commutation rule (2.1), the product of two polynomials in F [x; δ, σ] can be decided by association and distribution.

Example 2.2.1. Let F (t) be the univariate rational-function field in t over F.

(i) The usual polynomial ring F (t)[x] on which δ is the zero map and σ is the identity map.

(ii) The ring of differential operators F (t)[x; d/dt, 1], where 1 is the identical map.

(iii) The ring of recurrence operators F (t)[x; 0, σ], where σ(f(t)) = f(t + 1) for any f ∈ F (t).

(iv) The ring of difference operators F (t)[x; Δ, σ], where Δ = σ − 1.

In this thesis, we also denote the ring of differential operators over F (t) by F (t)⟨D_{t}⟩ and the ring of difference operators over F (t) by F (t)⟨S_{t}⟩.

If σ is an automorphism, we can perform both left and right Euclidean division on the ring F [x; δ, σ] [70], that is, for any f, g ∈ F [x; δ, σ], there exist q, r ∈ F [x; δ, σ] such that f = qg+r or f = gq + r, where deg(r) < deg(g). This implies that F [x; δ, σ] is a left and right principal ideal domain. A multivariate extension of Ore polynomials is studied in [35, 32, 93].

2.3 Integration of rational functions

Symbolic evaluation of integrals is an active research domain of computer algebra, which revives the study of the integration problems from the algorithmic point of view. The first landmark algorithm, designed by Risch [77, 78] in the late 1960’s, can decide whether indefinite integrals of elementary functions are elementary or not. State of the art about integrating transcendental functions has been presented in Bronstein’s book [21]. For later use, we review some basic methods for rational-function integration. For more intensive presentations, see the books [21, Chapter 2], [45, Chapter 11], and [90, Chapter 22].

In this section, let F be a field of characteristic zero and F (x) be the rational-function field in x over F. On the field F (x), the derivation δ is defined by setting δ(x) = 1 and δ(c) = 0 for all c ∈ F. In most undergraduate calculus textbooks, we can find that any rational function f ∈ ℂ(x) over the field ℂ of complex numbers has an elementary integral of the form

∫ f dx = g + ∑_{i=1}^{n} β_{i} log(x − α_{i})  (2.2)

where g ∈ ℂ(x) and α_{i}, β_{i} ∈ ℂ for i = 1, …, n. Here, we call g the rational part of the integral and the sum of logarithms the logarithmic part of the integral. In this analytic setting when the ground field is the algebraically closed field ℂ, the problem of rational function integration is well-understood. However, the algorithmic question becomes more involved when the ground field is not algebraically closed. Over an arbitrary field F, one may need to compute objects over an algebraic extension of F. In this case, the method in calculus books has many practical difficulties. For any rational function f ∈ F (x), it has been shown by Hermite [52] that the rational part of the integral ∫ f dx still lies in F (x) and he presented a method to compute this part without introducing any algebraic extension of F. His method is named Hermite reduction in [21, 45, 90]. For the logarithmic part, in general it is required to introduce algebraic extensions of F and then the question is how to compute the minimal algebraic extension of F that is sufficient to express the integral. This question has been independently solved by Trager [87] and Rothstein [80].

2.3.1 Hermite reduction

Hermite reduction is a practical method for computing the rational part, which uses the method of integration by parts in order to reduce the denominator of the integrand to a square-free polynomial. More explicitly, Hermite reduction decomposes f ∈ F (x) into

f = δ(g) + \frac{a}{b},  (2.3)

where g ∈ F (x) and a, b ∈ F [x] with deg(a) < deg(b) and b squarefree. Such a pair (g, a/b) in (2.3) is called an additive decomposition for f (with respect to x). The rational functions g and a/b are called the rational and logarithmic parts of f, respectively.

We summarize the idea of Hermite reduction following the treatment in [21, Chapter 2.2]. Write the integrand as f = A/D with A, D ∈ F [x] and gcd(A, D) = 1. Let D = D_{1}D_{2}^{2} … D_{n}^{n} be the squarefree factorization of D. By computing the partial fraction decomposition for f, we have

f = P + ∑_{i=1}^{n} \frac{A_{i}}{D_{i}^{i}}  (2.4)

where P and the A_{i}’s are in F [x] and deg(A_{i}) < deg(D_{i}^{i}) for each i. By the linearity of δ, it is sufficient to consider the same decomposition problem for any fraction of the form

\frac{P}{Q^{m}} ∈ F [x], with m > 1, deg(P) < m deg(Q) and Q squarefree.

The square-freeness of Q implies the existence of a Bézout relation

Sδ(Q) + T Q = 1,  (2.5)

where S, T ∈ F [x] can be obtained by the extended Euclidean algorithm. Furthermore, upon multiplication and division with remainder, we get

P = SP δ(Q) + T P Q = (GQ + B)δ(Q) + T P Q = Bδ(Q) + CQ  (2.6)

where C = Gδ(Q) + T P with deg(C) < (m − 1) deg(Q). Now, integration by parts yields

\frac{P}{Q^{m}} = \frac{Bδ(Q) + CQ}{Q^{m}} = δ (\frac{−(m − 1)^{−1}B}{Q^{m−1}} ) + \frac{C + (m − 1)^{−1}δ(B)}{Q^{m−1}}  (2.7)

where deg(C + (m − 1)^{−1}δ(B)) < (m − 1) deg(Q) since deg(B) < deg(Q) and deg(C) < (m − 1) deg(Q). This process is repeated until the denominator is square-free.

In 1975, Mack [69] introduced a variant of Hermite reduction that requires neither partial fraction decomposition nor squarefree factorization, but only extended GCD computation.

The following lemma due to Ostrogradsky [71] shows the uniqueness (up to an additive constant) of additive decompositions. For proofs, one can follow the argument in [45, Theorem 11.4].

Lemma 2.3.1. Let f = a/b be a nonzero rational function in F (x) such that a, b ∈ F [x], gcd(a, b) = 1, deg(a) < deg(b), and b is squarefree. Then there is no g ∈ F (x) such that f = δ(g).

Corollary 2.3.2. Let f be a rational function in F (x). Then the pair (g, a/b) satisfying (2.3) for f is unique up to adding a constant to g.

Proof. Assume that (g_{1}, a_{1}/b_{1}) and (g_{2}, a_{2}/b_{2}) are two additive decompositions for f with a_{1}/b_{1} ≠ a_{2}/b_{2}. The difference of logarithmic parts can be written as

0 ≠ a_{1}/b_{1} − a_{2}/b_{2} = \frac{A}{B} , where A, B ∈ F [x], gcd(A, B) = 1, and deg(A) < deg(B).

Since both b_{1} and b_{2} are squarefree and B divides lcm(b_{1}, b_{2}), B is also squarefree. However,

\frac{A}{B} = δ(g_{2} − g_{1}),

which contradicts Lemma 2.3.1. So a_{1}/b_{1} equals a_{2}/b_{2} and then g_{1} and g_{2} differ by an additive constant.

For later use, we recall a fact on the logarithmic derivatives of rational functions.

Lemma 2.3.3. Let f = a/b ∈ F (x) be such that a, b ∈ F [x] and gcd(a, b) = 1, then

\frac{δf}{f} = \frac{p}{a^{∗}b^{∗}} ,

where a^{∗} and b^{∗} are, respectively, the squarefree parts of a and b, and p ∈ F [x] with deg(p) < deg(a^{∗}b^{∗}) and gcd(p, a^{∗}b^{∗}) = 1.

Proof. Let a = a_{1}a_{2}^{2} ⋯ a_{m}^{m} and b = b_{1}b_{2}^{2} ⋯ b_{n}^{n} be the squarefree factorizations of a and b, respectively. Then a^{∗} = a_{1}a_{2} ⋯ a_{m} and b^{∗} = b_{1}b_{2} ⋯ b_{n}. By Lemma 2.1.1 (iv), we have

\frac{δf}{f} = \frac{p}{a^{∗}b^{∗}} , where p = b^{∗} ∑_{i=1}^{m} \frac{iδ(a_{i})a^{∗}}{a_{i}} − a^{∗} ∑_{j=1}^{n} \frac{jδ(b_{j})b^{∗}}{b_{i}} ∈ F [x].

It is easy to see that deg(p) < deg(a^{∗}b^{∗}). Since the a_{i}’s and b_{j}’s are squarefree and pairwise coprime, we have

gcd (a^{∗}, ∑_{i=1}^{m} \frac{iδ(a_{i})a^{∗}}{a_{i}} ) = gcd (b^{∗}, ∑_{j=1}^{n} \frac{jδ(b_{j})b^{∗}}{b_{i}} ) = 1,

which further implies gcd(p, a^{∗}b^{∗}) = 1. This completes the proof.

The following lemma will be useful to verify uniqueness properties, in particular, in the proofs of Theorem 4.4.4 and Lemma 4.3.10. Lemma 2.3.4. Assume that f is a rational function in F (x), p_{1}, …, p_{n} are pairwise coprime polynomials in F [x], and c_{1}, …, c_{n} are constants in F. If

δ(f) = ∑_{i=1}^{n} c_{i} \frac{δ(p_{i})}{p_{i}},

then f ∈ F and either c_{i} = 0 or p_{i} ∈ F for all i such that 1 ≤ i ≤ n. Proof. Let p_{i}^{∗} be the squarefree part of p_{i} for all i such that 1 ≤ i ≤ n. By Lemma 2.3.3, we have

∑_{i=1}^{n} c_{i} \frac{δ(p_{i})}{p_{i}} = ∑_{i=1}^{n} c_{i} \frac{q_{i}}{p_{i}^{∗}} =: \frac{a}{b},

where a, b ∈ F [x] with gcd(a, b) = 1 and the q_{i}’s are polynomials in F [x] with gcd(q_{i}, p_{i}^{∗}) = 1 and deg(q_{i}) < deg(p_{i}^{∗}). Since the p_{i}’s are pairwise coprime, so are the p_{i}^{∗}’s. Hence, the denominator b is squarefree and deg(a) < deg(b). By Lemma 2.3.1, a/b must be equal to zero and then f ∈ F. By the uniqueness of squarefree partial fraction decomposition, all the fractions c_{i}q_{i}/p_{i}^{∗} are equal to zero. This implies that either c_{i} = 0 or p_{i} ∈ F for each i such that 1 ≤ i ≤ n.

Corollary 2.3.5. Let c_{1}, …, c_{n} ∈ F be linearly independent over ℤ. If there exist rational functions f_{1}, …, f_{n} ∈ F (x) such that

∑_{i=1}^{n} c_{i} \frac{δ(f_{i})}{f_{i}} = 0,

then f_{1}, …, f_{n} belong to F.

Proof. Since every element of F is a constant with respect to δ, we may suppose that none of f_{1}, …, f_{n} belongs to F, and look for a contradiction. Under this assumption, the sum can be rewritten as

∑_{j=1}^{m} c̄_{j} \frac{δ(p_{j})}{p_{j}} = 0

where c̄_{j} is a ℤ-linear combination of c_{1}, …, c_{n}, and the p_{j}’s are nontrivial distinct irreducible factors of the denominators of the δ(f_{i})/f_{i}. Since every c̄_{j} is nonzero, Lemma 2.3.4 implies that every p_{j} is in F, which is a contradiction.

2.3.2 Ostrogradsky and Horowitz’s method

Ostrogradsky and Horowitz’s method [71, 53] computes the additive decomposition of a rational function by solving a linear system. Although this method has asymptotically higher complexity than that of Hermite reduction [90, Section 22.2], it is useful for our complexity analyses in the sequel.

Let f = P/Q ∈ F (x) be such that P, Q ∈ F [x] and gcd(P, Q) = 1. After reading out the polynomial part of f, we further assume that deg(P) < deg(Q). Let Q^{∗} be the squarefree part of Q and Q^{−} = Q/Q^{∗}. Denote d_{x}^{∗} = deg_{x}(Q^{∗}) and d_{x}^{−} = deg_{x}(Q^{−}). According to the functionality of Hermite reduction, one can find two unique polynomials A and a in F [x] with deg_{x} A < deg(Q^{−}) and deg_{x} a < deg(Q^{∗}) such that

\frac{P}{Q} = δ (\frac{A}{Q^{−}} ) + \frac{a}{Q^{∗}} .  (2.8)

Note that A and a satisfy (2.8) if and only if they satisfy the equation

P = δ(A)Q^{∗} − AQ̃ + aQ^{−},  (2.9)

where Q̃ = Q^{∗}δ(Q^{−})/Q^{−} is a polynomial in F [x] of degree deg(Q^{∗}) − 1. Now, equation (2.9) can be solved by the method of undetermined coefficients. Write P = ∑_{l=0}^{d_{x}^{∗}+d_{x}^{−}−1} P_{l}x^{l} and set A = ∑_{i=0}^{d_{x}^{−}−1} A_{i}x^{i}, a = ∑_{j=0}^{d_{x}^{∗}−1} a_{j}x^{j} with undetermined coefficients A_{i} and a_{j}. Then (2.9) holds if and only if

(A_{d_{x}^{−}−1}, …, A_{0}, a_{d_{x}^{∗}−1}, …, a_{0}) M = (P_{d_{x}^{∗}+d_{x}^{−}−1}, …, P_{0}) ,  (2.10)

where M is a deg(Q) × deg(Q) matrix over F obtained by equating the likewise powers of x in (2.9). By the uniqueness of A/Q^{−} and a/Q^{∗}, the system (2.10) has a unique solution, which leads to the following lemma.

Lemma 2.3.6. The matrix M in (2.10) is invertible over F.

We call the linear system (2.10) the Ostrogradsky–Horowitz system and M the Ostrogradsky– Horowitz matrix associated with Q. Note that M is uniquely determined by Q. So we denote this matrix by M(Q).

2.3.3 Residues and Rothstein–Trager resultants

After additive decomposition, the integration problem of rational functions is reduced to computing the integrals of the form

∫ \frac{a}{b} dx, where a, b ∈ F [x] with gcd(a, b) = 1, deg(a) < deg(b) and b squarefree. (2.11)

Over the algebraic closure \overline{F} of F, the rational function a/b above decomposes into

\frac{a}{b} = ∑_{i=1}^{n} \frac{β_{i}}{x − α_{i}} , where α_{i}, β_{i} ∈ \overline{F} and b(α_{i}) = 0 for 1 ≤ i ≤ n.

Consequently, the integral of a/b can be expressed by

∫ \frac{a}{b} dx = ∑_{i=1}^{n} β_{i} log(x − α_{i}).

By convention, the value β_{i} is called the residue of a/b at the point x = α_{i}. According to the Lagrange formula ([44, page 38] or [45, Exercise 11.8]), the residue of a/b at α_{i} is

β_{i} = \frac{a}{δ(b)}(α_{i}) ∈ F (α_{i}).  (2.12)

So we can always express the integrals in (2.11) over the splitting field of the denominator b. However, it is not necessary to compute splitting fields for obtaining the integrals. For instance, the integral below can be expressed without any algebraic extension:

∫ \frac{2x}{x^{2} − 2} dx = log(x + \sqrt{2}) + log(x − \sqrt{2}) = log(x^{2} − 2).

In fact, Rothstein [80] and Trager [87] have shown that the minimal algebraic extension of F for expressing the integrals is the splitting field of the following polynomial

R(z) = resultant_{x}(b, a − z · δ(b)) ∈ F [z].  (2.13)

In the literature, R(z) above is called the Rothstein–Trager resultant of a/b with respect to x, denoted henceforth by RT_{x}(a/b). By the formula (2.12), all the residues of a/b at its poles are the roots of R(z). Moreover, we have the following lemma, which appears implicitly in the literature [45, 21, 90]. Lemma 2.3.7. Let a, b ∈ F [x] be such that deg(a) < deg(b), gcd(a, b) = 1 and b is squarefree in F [x]. Let R(z) be the Rothstein–Trager resultant of a/b with respect to x. Then we have

(i) all roots of R(z) are nonzero;

(ii) if α_{1} and α_{2} are two distinct roots of R(z), then p_{1} and p_{2} are coprime over \overline{F}, where p_{i} = gcd(b, a − α_{i}δ(b)) ∈ F (α_{i})[x], for i = 1, 2.

Proof. The first assertion follows from the fact that gcd(a, b) = 1. For the second one, we suppose p_{1} and p_{2} are not coprime over \overline{F}. Then there exists β ∈ \overline{F} such that p_{1}(β) = p_{2}(β) = 0. By the definition of p_{1} and p_{2}, b(β) = 0 and (a − α_{i}δ(b))(β) = 0 for i = 1, 2. Since b is squarefree in F [x], δ(b)(β) ≠ 0, which implies

α_{1} = α_{2} = \frac{a}{δb}(β),

which is a contradiction with α_{1} ≠ α_{2}.

The next theorem is fundamental to design the algorithm for computing the logarithmic part of the integral of a rational function. For proofs, see [21, Theorem 2.4.1] or [90, Theorem 22.8]. Theorem 2.3.8 (Rothstein and Trager’s theorem). Let a, b ∈ F [x] be such that deg(a) < deg(b), gcd(a, b) = 1 and b is squarefree in F [x]. Let R(z) be the Rothstein–Trager resultant of a/b with respect to x. Then the integral of rational function a/b can be expressed as

∫ \frac{a}{b} dx = ∑_{i=1}^{n} c_{i} log(g_{i}),

where c_{1}, … c_{n} ∈ \overline{F} are the distinct roots of R(z) and g_{i} = gcd(b, a − c_{i}δ(b)) ∈ F (c_{i})[x] for all 1 ≤ i ≤ n. Moreover, if E is an algebraic extension of F such that we have

∫ \frac{a}{b} dx = ∑_{i=1}^{m} c̃_{i} log(g̃_{i}),

where c̃_{1}, …, c̃_{m} ∈ E\{0} and g̃_{1}, …, g̃_{m} ∈ E[x]\E are monic, squarefree, and pairwise coprime, then E contains all the roots of R(z).

We conclude this section by the following theorem on the integration of rational functions.

Theorem 2.3.9. Let f be a nonzero element in F (x). Then there exist g ∈ F (x), nonzero elements c_{j} ∈ \overline{F} and p_{j} ∈ F (c_{j})[x] \ F (c_{j}) for j = 1, …, n such that

f = δ(g) + ∑_{j=1}^{n} c_{j} \frac{δ(p_{j})}{p_{j}}.

Moreover, p_{j} and p_{j^{∗}} are coprime over \overline{F} whenever j ≠ j^{∗}.

Proof. By Hermite reduction, there exist f ∈ F (x), a, b ∈ F [x] with deg a < deg b, gcd(a, b) = 1 and b being squarefree such that

f = δ(g) + \frac{a}{b}.

If a is zero, then there is nothing to prove. Assume that a is nonzero and let Λ be the set of distinct roots of the Rothstein-Trager resultant R(z) of a/b. By Lemma 2.3.7 (i), all elements in Λ are nonzero. By Theorem 2.3.8,

f = δ(g) + ∑_{λ∈Λ} λ\frac{δ(p_{λ})}{p_{λ}},

where p_{λ} = gcd(b, a − λδb) ∈ F (λ)[x] \ F (λ). By Lemma 2.3.7 (ii), p_{λ} and p_{μ} are coprime if λ and μ are two distinct elements of Λ. This completes the proof.

2.4 Background on complexity

In this section, we recall basic notation and complexity results from [90]. These results are useful to analyze algorithms in Chapter 3. The complexity results are expressed using the “big Oh" notation [90, Section 25.7]. For instance, we say that the Karatsuba algorithm for multiplying two polynomials in ℚ[x] of degree at most n takes O(n^{log_{2} 3}) arithmetic operations in ℚ. Furthermore, we use the notation Õ(·) to indicate cost estimates with hidden logarithmic factors [90, Definition 25.8].

Let k be a field of characteristic zero. Unless otherwise specified, all complexity estimates are given in terms of arithmetic operations in k, which we denote by “ops”. Let k[x]_{≤d}^{m×n} be the set of m × n matrices with coefficients in k[x] of degree at most d. Let ω ∈ [2, 3] be a feasible exponent of matrix multiplication, so that two matrices from k^{n×n} can be multiplied using O(n^{ω}) ops.

Facts 2.4.1 and 2.4.2 below show the complexity of multi-point evaluation, rational interpolation, and algebraic operations on polynomial matrices using fast arithmetic. For proofs, see [90, Corollaries 10.8, 5.18, 11.6] and [84, Theorem 7.3]. Fact 2.4.1. For a polynomial p ∈ k[x] of degree less than n, pairwise distinct u_{0}, …, u_{n−1} in k, and v_{0}, …, v_{n−1} ∈ k, we have:

(i) Fast evaluation: evaluating p at the u_{i}’s takes Õ(n) ops. (ii) Fast interpolation: for m ∈ {1, …, 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 0 ≤ i ≤ n − 1 takes Õ(n) ops. Fact 2.4.2. For M in k[x]_{≤d}^{m×n} with d > 0, we have:

(i) If M = (M_{1} M_{2}) is an invertible n × n matrix with M_{i} ∈ k[x]_{≤d_{i}}^{n×n_{i}}, where i = 1, 2 and n_{1} + n_{2} = n, then the degree of det(M) is at most n_{1}d_{1} + n_{2}d_{2}.

(ii) If M = (M_{1} M_{2}) is not of full rank and with M_{i} ∈ k[x]_{≤d_{i}}^{m×n_{i}}, where i = 1, 2 and n_{1}+n_{2} = n, then there exists a nonzero u ∈ k[x]^{n} with coefficients of degree at most n_{1}d_{1} +n_{2}d_{2} such that Mu = 0.

(iii) The rank r and a basis of the null space of M can be computed using Õ(nmr^{ω−2}d) ops.

The complexity in terms of arithmetic operations in k of Hermite reduction has been analyzed by Yun in [94] and also [90, Theorem 22.7]. Lemma 2.4.3 (Yun, 1977). Let f be a nonzero rational function in k(x) of degree at most n in x, then Hermite reduction on f can be performed using Õ(n) operations in k.

26

Hermite Reduction for

Rational-Function Telescoping

3.1 Introduction

Although creative-telescoping algorithms have been now used extensively in modern computer algebra, very little is known about their complexity: some related results seem to be the complexity analyses of an algorithm for hyperexponential indefinite integration in [49], of algorithms for rational and hypergeometric summation in [50, 19], and of an algorithm by Takayama for finding a recurrence for hypergeometric series of a special form in [86]. In order to get complexity estimates, we simplify the problem by restricting to a smaller class of inputs, namely that of bivariate rational functions. Therefore, our goal reads as follows.

Problem 1. Given a rational function f ∈ k(x, y) \ {0}, find a nonzero operator L in k(x)⟨D_{x}⟩ and a rational function g in k(x, y) such that

L(x, D_{x})(f) = D_{y}(g).  (3.1)

Here, L is called a telescoper for f with respect to y and g the certificate of L for f. Since D_{y} commutes with any element in k(x)⟨D_{x}⟩, we can always choose telescopers with polynomial coefficients in k[x].

By considering this constrained class of inputs, we are indeed able to blend the general method of creative telescoping with the well-known Hermite reduction [52]. Although restricted, this class already has many applications. In 1827, Abel observed that an algebraic function satisfied a linear differential equation [1]. The annihilating differential equations are important in the study of algebraic functions [37, 30, 31]. In [20], it was shown that differential equations for algebraic functions can be computed via rational-function telescoping. In combinatorics, many nontrivial problems are encoded as diagonals of rational formal power series [83, 74]. Differential equations for diagonals of a bivariate rational function can also be constructed via rationalfunction telescoping.

Essentially two algorithms for computing telescopers of minimal order can be found in the literature: the classical way [15] is to apply a differential analogue of Gosper’s indefinite summation algorithm [51], which reduces the problem to solving an auxiliary linear differential equation for polynomial solutions. An algorithm developed later in [46] (see also [63]) performs Hermite reduction on f to get an additive decomposition of the form

f = D_{y}(g) + ∑_{i=1}^{m} \frac{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 [15] is applied to each u_{i}/v_{i} to get a minimal telescoper L_{i} for it. The least common left multiple of the L_{i}’s is proved to be a minimal telescoper for f. This algorithm performs well only for specific inputs (both in practice and from the complexity viewpoint), but it inspired our improvements over Almkvist and Zeilberger’s method.

Our contribution. For bivariate rational functions, we present a new and provably faster algorithm for computing minimal telescopers, which is based on Hermite reduction. Over the classical method by Almkvist and Zeilberger, we make some improvements avoiding unnecessary resultant calculations and integer-root finding. We derive complexity estimates of all algorithms described in Figure 3.1, showing that our approach is faster. These complexity results give the first proof of a polynomial complexity (in terms of field operations) for creative telescoping on bivariate rational-function inputs. In the computer algebra system Maple, we have implemented our algorithms for rational-function telescoping, which has been integrated into the existing Maple library Algolib 13.0 (for Maple 13). Our experimental results show that the Hermitereduction based algorithm outperforms all other known algorithms concerning both the worstcase complexity and the actual timings of implementations.

Algorithms  deg_{D_{x}}(L) deg_{x}(L) deg_{x}(g) deg_{y}(g) Complexity
HermiteTelescoping ≤ d_{y}  O(d_{x}d_{y}^{2}) O(d_{x}d_{y}^{2}) O(d_{y}^{2}) Õ(d_{x}d_{y}^{ω+3})
RationalAZ ≤ d_{y} O(d_{x}d_{y}^{2}) O(d_{x}d_{y}^{2}) O(d_{y}^{2}) Õ(d_{x}d_{y}^{2ω+2})

Table 3.1: Complexity for rational-function telescoping

The rest of this chapter is organized as follows. In Section 3.2, we recall notation about bivariate polynomials and set some hypotheses used in this chapter. We study Hermite reduction over k(x) in Section 3.3, proving tight output degree bounds and an optimal algorithm via fast evaluation and interpolation. We present a Hermite-reduction based algorithm for minimal telescopers and some improvements over the classical method in Section 3.4. In Section 3.5, we show the existence of telescopers and estimate their order bounds. At last, we describe in Section 3.6 an implementation and show some experimental results.

The work in this chapter is published in [18], which is a joint work with Alin Bostan, Frédéric Chyzak and Ziming Li.

3.2 Notation and hypotheses

In this section, we review some basic notation and introduce some hypotheses for later use.

Let Q be a bivariate polynomial in k[x, y]\k[x]. The squarefree factorisation of Q with respect to y is the unique product qQ_{1}Q_{2}^{2} ⋯ Q_{m}^{m} where q ∈ k[x] and Q_{i} ∈ k[x, y] satisfying deg_{y}(Q_{m}) > 0 and such that the Q_{i}’s are primitive, squarefree, and pairwise coprime. The squarefree part Q^{∗} of Q with respect to y is the product Q_{1}Q_{2} ⋯ Q_{m}. Let Q^{−} denote the polynomial Q/Q^{∗}, and lc_{y}(Q) the leading coefficient of Q with respect to y. The following two formulas about Q, Q^{∗}, and Q^{−} can be proved by mere calculations.

Fact 3.2.1. Let Q̂_{i} denote Q^{∗}/Q_{i}. Then we have

(i) Q^{∗}D_{y}(Q^{−})/Q^{−} = ∑_{i=1}^{m}(i − 1)Q̂_{i}D_{y}(Q_{i}) ∈ k[x, y];

(ii) D_{y}(Q)/Q^{−} = ∑_{i=1}^{m} iQ̂_{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 denoted by 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 denoted by 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 may always 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. A telescoper L for f̄ with certificate ḡ is a telescoper for f with certificate L(p) + ḡ. So we introduce the following hypothesis.

Hypothesis (H) From now on, P and Q are assumed to be nonzero polynomials in k[x, y] such that deg_{y}(P) < deg_{y}(Q), gcd(P, Q) = 1, and Q is primitive with respect to 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.

Sometimes, we use the following hypothesis in order to make estimates concise.

Hypothesis (H’) Hypothesis (H) and deg_{x}(P) ≤ d_{x}.

3.3 Hermite reduction for bivariate rational functions

In this section, we will apply Hermite reduction to a bivariate rational function of k(x, y) with respect to y. We present a quasi-optimal algorithm to perform the Hermite reduction by using a fast evaluation-interpolation approach.

3.3.1 Output size estimates

In order to use the evaluation-interpolation strategy, we first derive an upper bound on the degrees of the outputs of the Hermite reduction on a rational function. To this end, we will use the linear system introduced in Ostrogradsky and Horowitz’s method.

After reading out the polynomial part, we may assume that f = P/Q ∈ k(x, y) with P and Q in k[x, y] such that gcd(P, Q) = 1 and deg_{y}(P) < deg_{y}(Q). Let D_{x} and D_{y} denote the usual derivations ∂/∂x and ∂/∂y on k(x, y), respectively. Recall that Q^{∗} denotes the squarefree part of Q with respect to y and Q^{−} = Q/Q^{∗}. Set

d_{x}^{∗} = deg_{x}(Q^{∗}), d_{y}^{∗} = deg_{y}(Q^{∗}), d_{x}^{−} = deg_{x}(Q^{−}), and d_{y}^{−} = deg_{y}(Q^{−}).

According to Ostrogradsky and Horowitz’s method, we make the ansatz

\frac{P}{Q} = D_{y} (\frac{A}{Q^{−}} ) + \frac{a}{Q^{∗}} ,  (3.2)

where A, a ∈ k(x)[y] with deg_{y}(A) < d_{y}^{−} and deg_{y}(a) < d_{y}^{∗}. In order to bound the bidegrees of A and a, we reformulate (3.2) into the equivalent form

P = Q^{∗}D_{y}(A) − (\frac{Q^{∗}D_{y}(Q^{−})}{Q^{−}}) A + Q^{−}a,  (3.3)

where Q̃ = Q^{∗}D_{y}(Q^{−})/Q^{−} is a polynomial in k[x, y] of bidegree at most (d_{x}^{∗}, d_{y}^{∗} −1) by Fact 3.2.1. Viewing A and a as polynomials in k(x)[y] with undetermined coefficients, we form the following Ostrogradsky–Horowitz system,

M (Â ; â) = P̂,  (3.4)

where M is the Ostrogradsky–Horowitz matrix associated to Q and Â, â, and P̂ are the coefficient vectors of A, a, and P with sizes d_{y}^{−}, d_{y}^{∗}, and d_{y}, respectively. From equation (3.3), the matrix M is of the form

(M_{1} M_{2}) , where M_{1} ∈ k[x]_{≤d_{x}^{∗}}^{d_{y}×d_{y}^{−}} and M_{2} ∈ k[x]_{≤d_{x}^{−}}^{d_{y}×d_{y}^{∗}}.

Let Δ be the determinant of M, so that deg_{x}(Δ) ≤ μ := d_{x}^{∗}d_{y}^{−} + d_{x}^{−}d_{y}^{∗} by Fact 2.4.2(i). For later use, we also define Δ′ as the determinant of M(Q^{∗}^{2}), so that deg_{x}(Δ′) ≤ μ′ := 2d_{x}^{∗}d_{y}^{∗} by Fact 2.4.2 (i) and since (Q^{∗}^{2})^{−} = Q^{∗}.

Lemma 3.3.1. Let f = P/Q ∈ k(x, y) with P, Q ∈ k[x, y], gcd(P, Q) = 1 and deg_{y}(P) < deg_{y}(Q). Then there exist B, b ∈ k[x, y] with deg_{y}(B) < d_{y}^{−} and deg_{y}(b) < d_{y}^{∗}, and such that:

(i) f = D_{y} (\frac{B}{ΔQ^{−}} ) + \frac{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 (3.4) leads to the first assertion. The second assertion follows by determinant expansions.

Algorithm HermiteEvalInterp Input: P, Q ∈ k[x, y] satisfying Hypothesis (H). Output: (A, a) ∈ k(x)[y]^{2} solving (3.2).

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 integers that are lucky for Q;

4. For each x_{0} ∈ S, compute (A_{0}, a_{0}) ∈ k[y]^{2} such that

\frac{P (x_{0}, y)}{Q(x_{0}, y)} = D_{y} (\frac{A_{0}}{Q^{−}(x_{0}, y)}) + \frac{a_{0}}{Q^{∗}(x_{0}, y)}

using Hermite reduction over k;

5. Compute (A, a) ∈ k(x)[y] by rational interpolation and return this pair. Figure 3.1: Hermite reduction over k(x) via evaluation and interpolation.

In what follows, we shall encounter proper rational functions with denominator Q satisfying Q = Q^{∗}^{2}. The following lemma is an easy corollary of Lemma 3.3.1 for such functions.

Corollary 3.3.2. Assume that Q = Q^{∗}^{2} additionally in Lemma 3.3.1. Then there exist B, b ∈ k[x, y] with deg_{y}(B) and deg_{y}(b) less than d_{y}^{∗}, and such that

(i) \frac{P}{Q^{∗}^{2}} = D_{y} (\frac{B}{Δ′Q^{∗}} ) + \frac{b}{δ′Q^{∗}};

(ii) deg_{x}(B) and deg_{x}(b) are bounded by μ′ − d_{x}^{∗} + deg_{x}(P).

3.3.2 Algorithm by evaluation and interpolation

We observe that an asymptotically optimal complexity can be achieved by evaluation and interpolation at each step of Hermite reduction over k(x). This inspires us to adapt Gerhard’s modular method [48, 49] to k(x, y). For simplicity, we further assume that Q ∈ k[x, y] is nonzero and primitive over k[x].

Definition 3.3.3. An element x_{0} ∈ k is said to be lucky if lc_{y}(Q)(x_{0}) ≠ 0 and the degree in y of gcd(Q(x_{0}, y), D_{y}Q(x_{0}, y)) is equal to deg_{y}(Q^{−}).

Lemma 3.3.4. There are at most d_{x}(2d_{y}^{∗} − 1) unlucky points. Proof. Let σ ∈ k[x] be the d_{y}^{−}th subresultant with respect to y of Q and D_{y}(Q). By Corollary 5.5 in [49], all unlucky points are in the set U = { x_{0} ∈ k | σ(x_{0}) = 0 }. By Corollary 3.2 (ii) in [49], deg_{x}(σ) ≤ d_{x}(2d_{y}^{∗} − 1).

Lemma 3.3.5. Let B, b, and Δ be the same as in Lemma 3.3.1, and let x_{0} ∈ k be lucky. Then Δ(x_{0}) ≠ 0 and (B(x_{0}, y), b(x_{0}, y)) is the unique pair such that

\frac{P (x_{0}, y)}{Q(x_{0}, y)} = D_{y} (\frac{B(x_{0}, y)}{Δ(x_{0})Q^{−}(x_{0}, y)}) + \frac{b(x_{0}, y)}{Δ(x_{0})Q^{∗}(x_{0}, y)}.  (3.5)

Proof. By the luckiness of x_{0}, deg_{y}(Q(x_{0}, y)) = d_{y} and Q(x_{0}, y)^{−} = Q^{−}(x_{0}, y). Then we have Q(x_{0}, y)^{∗} = Q^{∗}(x_{0}, y). This implies M(Q)(x_{0}, y) = M(Q(x_{0}, y)), which, by Lemma 2.3.6, is invertible over k(x). Hence Δ(x_{0}) ≠ 0, and the evaluation at x = x_{0} of the equality in Lemma 3.3.1 (i) is well-defined. Thus, (B(x_{0}, y), b(x_{0}, y)) is a solution of (3.5). Uniqueness follows from Corollary 2.3.2.

Theorem 3.3.6. Algorithm HermiteEvalInterp in Figure 3.1 is correct and takes Õ(d_{x}d_{y}^{2} + deg_{x}(P)d_{y}) operations in k. Proof. Set ν to d_{x}(2d_{y}^{∗} − 1). Lemma 3.3.4 implies that the λ + 1 lucky points found in Step 3 are all less than λ+ν +1. By Corollary 2.3.2 and 3.3.1 (i), A = B/Δ and a = b/Δ. By Lemma 3.3.5, A_{0} = B(x_{0}, y)/Δ(x_{0}) and a_{0} = b(x_{0}, y)/Δ(x_{0}). By Lemma 3.3.1 (ii) and since deg_{x}(Δ) ≤ μ, it suffices to rationally interpolate A and a from values at λ + 1 lucky points. This shows the correctness. The dominant computation in Step 1 is the gcd, which takes Õ(d_{x}d_{y}) ops by Corollary 11.9 in [90]. For each integer i ≤ λ + ν, testing luckiness amounts to evaluations at x_{0} and computing gcd(Q(x_{0}, y), D_{y}(Q(x_{0}, y))), which takes Õ(d_{y}) ops by Fact 2.4.1 (i) and Corollary 11.6 in [90]. Then, generating S in Step 3 costs Õ((λ + ν + 1)d_{y}) ops. By Fact 2.4.1 (i), evaluations in Step 4 take Õ((λ + 1)d_{y}) ops. For each x_{0} ∈ S, the cost of the Hermite reduction in Step 4 is Õ(d_{y}) ops by Lemma 2.4.3. Thus, the total cost of Step 4 is Õ((λ + 1)d_{y}) ops. By Fact 2.4.1(ii), Step 5 takes Õ((λ + 1)d_{y}) ops. Since λ ≤ 2d_{x}d_{y} + deg_{x}(P) and ν ≤ 2d_{x}d_{y}, the total cost is as announced.

Remark 3.3.7. As the generic output size of Hermite reduction is proportional to λd_{y}, which is equal to O((d_{x}d_{y} + deg_{x}(P))d_{y}), Algorithm HermiteEvalInterp has quasi-optimal complexity.

3.4 Minimal telescopers for bivariate rational functions

Since the 1990’s, the main emphasis in existing works [95, 15, 46, 63] has been on finding telescopers of order minimal over all telescopers for f, which are called minimal telescopers. For a rational function f ∈ k(x, y), all its telescopers form a left ideal I_{f} in the left principal ideal domain k(x)⟨D_{x}⟩. So a minimal telescoper is a generator of I_{f} and two minimal telescopers differ by a multiplicative left factor in k(x). Therefore, a rational function has a unique monic minimal telescoper. In this section, we present a new method for computing minimal telescopers for rational functions in k(x, y), which is based on Hermite reduction. Also, we make some improvements over the classical method by Almkvist and Zeilberger. We will show that the arithmetic complexity of Hermite reduction approach is lower than that of the classical one.

3.4.1 Hermite-reduction based method

In some applications, we only need to compute telescopers without certificates. In the next section, we will show that the arithmetic size of certificates is asymptotically larger than that of telescopers. This motivates us to find a way in which one could choose to compute or not to compute the certificates according to the user requirements. To this end, we design a new algorithm, presented in Figure 3.2, to compute minimal telescopers for rational functions by basing on Hermite reduction.

For f = P/Q ∈ k(x, y) and i ∈ ℕ, Hermite reduction decomposes D_{x}^{i}(f) into

D_{x}^{i}(f) = D_{y}(g_{i}) + r_{i},  (3.6)

where g_{i}, r_{i} ∈ k(x, y) are proper. Since the squarefree part of the denominator of D_{x}^{i}(f) divides Q^{∗}, so does the denominator of r_{i}. The following lemma shows that (3.6) recombines into telescopers and certificates; next, Lemma 3.4.2 implies that the first pair obtained in this way by Algorithm HermiteTelescoping in Figure 3.2 yields a minimal telescoper.

Lemma 3.4.1. The rational functions r_{0}, …, r_{d_{y}^{∗}} are linearly dependent over k(x).

Proof. The constraints on r_{i} imply deg_{y}(r_{i}Q^{∗}) < d_{y}^{∗} for all i ∈ ℕ, from which follows the existence of a nontrivial linear dependence among the r_{i}’s over k(x).

Algorithm HermiteTelescoping Input: f = P/Q ∈ k(x, y) satisfying Hypothesis (H). Output: A minimal telescoper L ∈ k[x]⟨D_{x}⟩ with certificate g ∈ k(x, y).

1. Apply HermiteEvalInterp to f to get (g_{0}, a_{0}) such that f = D_{y}(g_{0}) + a_{0}/Q^{∗}. If a_{0} = 0, return (1, g_{0}).

2. For i from 1 to deg_{y}(Q^{∗}) do

(a) Apply HermiteEvalInterp to −a_{i−1}D_{x}(Q^{∗})/Q^{∗}^{2} to express it as D_{y}(g̃_{i}) + ã_{i}/Q^{∗}.

(b) Set g_{i} = D_{x}(g_{i−1}) + g̃_{i} and a_{i} = D_{x}(a_{i−1}) + ã_{i}.

(c) Solve ∑_{j=0}^{i} η_{j}a_{j} = 0 for η_{j} ∈ k(x) using [84, Algorithm Nullspace].

If there exists a nontrivial solution, then set (L, g) := (∑_{j=0}^{i} η_{j}D_{x}^{j}, ∑_{j=0}^{i} η_{j}g_{j}), and break.

3. Compute the content c of L and return (c^{−1}L, c^{−1}g).

Figure 3.2: Creative telescoping by Hermite reduction

Lemma 3.4.2. An integer ρ is minimal such that ∑_{i=0}^{ρ} η_{i}r_{i} = 0 for η_{0}, …, η_{ρ} ∈ k(x) not all zero if and only if ∑_{i=0}^{ρ} η_{i}D_{x}^{i} is a minimal telescoper for f with certificate ∑_{i=0}^{ρ} η_{i}g_{i}.

Proof. Multiplying (3.6) by η_{i} before summing yields

L(f) = D_{y}(∑_{i=0}^{ρ} η_{i}g_{i}) + ∑_{i=0}^{ρ} η_{i}r_{i} for L := ∑_{i=0}^{ρ} η_{i}D_{x}^{i},

where the sum ∑_{i=0}^{ρ} η_{i}r_{i} is a proper fraction in y with a squarefree denominator with respect to y. Thus, by Corollary 2.3.2, L is a telescoper of order ρ for f with certificate ∑_{i=0}^{ρ} η_{i}g_{i} if and only if ∑_{i=0}^{ρ} η_{i}r_{i} = 0 with η_{ρ} ≠ 0. The lemma follows.

Remark 3.4.3. In the algorithm HermiteTelescoping described in Figure 3.2, one can choose not to compute the certificate g, because the minimal telescoper and this corresponding certificate are computed separately.

Order bounds for minimal telescopers

Lemmas 3.4.1 and 3.4.2 combine into a tight upper bound on the order of minimal telescopers for f.

Corollary 3.4.4. Minimal telescopers have order at most d_{y}^{∗}.

We also derive a lower bound on the order of the minimal telescoper, to be used as an optimisation trick in our implementation. First, we choose a lucky element x_{0} ∈ k, and then apply Hermite reduction in k(y) to D_{x}^{i}(f)(x_{0}, y), we get

D_{x}^{i}(f)(x_{0}, y) = D_{y}(g_{0,i}) + r_{0,i},  (3.7)

where g_{0,i}, r_{0,i} ∈ k(y) are proper and the denominator of r_{0,i} divides Q^{∗}(x_{0}, y). Let ρ_{0} be the smallest integer such that r_{0,0}, …, r_{0,ρ_{0}} are linearly dependent over k.

Lemma 3.4.5. A minimal telescoper has order at least ρ_{0}.

Proof. We first claim that r_{0,i} = r_{i}(x_{0}, y), for r_{i} as in (3.6). Note that the squarefree part with respect to y of the denominator of D_{x}^{i}(f) divides Q^{∗} for all i ∈ ℕ. By Corollary 5.5 in [49], x_{0} is lucky for the denominator of D_{x}^{i}(f) for all i ∈ ℕ. Then the claim on r_{0,i} follows from Lemma 3.3.5 applied to D_{x}^{i}(f). Let ρ be the minimal order of a telescoper, then r_{0}, …, r_{ρ} are linearly dependent over k(x) by Lemma 3.4.2. Thus the rational functions r_{0,0}, …, r_{0,ρ} are linearly dependent over k, which implies ρ_{0} ≤ ρ.

Degree bounds for minimal telescopers

To derive degree bounds for g_{i} and r_{i} in (3.6), let Δ, Δ′, μ, and μ′ be defined as in the paragraph before Lemma 3.3.1, and set μ′′ = μ + μ′ − 1.

Lemma 3.4.6. Let W be in k[x, y] with deg_{y}(W) < d_{y}^{∗}. Then, for all i ∈ ℕ, there exist B, b ∈ k[x, y] with both bideg(B) and bideg(b) bounded by (deg_{x}(W) + μ′′, d_{y}^{∗} − 1), such that

D_{x} (\frac{W}{Δ^{i+1}Δ′^{i}Q^{∗}} ) = D_{y} (\frac{B}{Δ^{i+2}Δ′^{i+1}Q^{∗}} ) + \frac{b}{Δ^{i+2}Δ′^{i+1}Q^{∗}} .

Proof. A straightforward calculation leads to

D_{x} (\frac{W}{Δ^{i+1}Δ′^{i}Q^{∗}} ) = \frac{W̃}{Δ^{i+2}Δ′^{i+1}Q^{∗}} − \frac{1}{Δ^{i+1}Δ′^{i}} \frac{W D_{x}(Q^{∗})}{Q^{∗}^{2}},

where bideg(W̃) ≤ (deg_{x}(W) + μ′′, d_{y}^{∗} − 1). By Corollary 3.3.2, there exist B̃, b̃ ∈ k[x, y] such that

\frac{1}{Δ^{i+1}Δ′^{i}} \frac{W D_{x}(Q^{∗})}{Q^{∗}^{2}} = \frac{1}{Δ^{i+2}Δ′^{i+1}} (D_{y} (\frac{ΔB̃}{Q^{∗}} ) + \frac{Δb̃}{Q^{∗}}) ,

where bideg(B̃) and bideg(b̃) are bounded by (deg_{x}(W) + μ′ − 1, d_{y}^{∗} − 1). The proof is completed by setting (B, b) = (−ΔB̃, W̃ − Δb̃).

Lemma 3.4.7. For i ∈ ℕ, there exist B_{i}, b_{i} ∈ k[x, y] such that

D_{x}^{i}(f) = D_{y} (\frac{B_{i}}{Δ^{i+1}Δ′^{i}Q^{∗}^{i}Q^{−}} ) + \frac{b_{i}}{Δ^{i+1}Δ′^{i}Q^{∗}} .  (3.8)

Moreover, bideg(B_{i}) ≤ (deg_{x}(P) + μ + iμ′′ + (i − 1)d_{x}^{∗}, id_{y}^{∗} + d_{y}^{−} − 1) and bideg(b_{i}) ≤ (deg_{x}(P) + μ + iμ′′ − d_{x}^{−}, d_{y}^{∗} − 1).

Proof. We proceed by induction on i. For i = 0, the claim follows from Lemma 3.3.1. Assume that i > 0 and that the claim holds for the values less than i. For brevity, we set γ = deg_{x}(P)+μ, F_{i−1} = B_{i−1}/(Δ^{i}Δ′^{i−1}Q^{∗}^{i−1}Q^{−}), and G_{i−1} = b_{i−1}/(Δ^{i}Δ′^{i−1}Q^{∗}). The induction hypothesis implies

D_{x}^{i}(f) = D_{y}D_{x}(F_{i−1}) + D_{x}(G_{i−1}),

with bidegree bounds on B_{i−1} and b_{i−1}. Fact 3.2.1(i) implies that Q̃ := Q^{∗}D_{x}(Q^{−})/Q^{−} is in k[x, y], with bideg(Q̃) ≤ (d_{x}^{∗} − 1, d_{y}^{∗}). Hence D_{x}(1/Q^{−}) = −Q̃/Q. This observation and an easy calculation imply that

D_{x}(F_{i−1}) = \frac{B̃_{i−1}}{Δ^{i+1}Δ′^{i}Q^{∗}^{i}Q^{−}} ,

where B̃_{i−1} ∈ k[x, y] and deg_{x}(B̃_{i−1}) ≤ deg_{x}(B_{i−1}) + μ′′ + d_{x}^{∗}. Furthermore, by Lemma 3.4.6 there are B̄_{i}, b̄_{i} ∈ k[x, y] with bidegrees at most (deg_{x}(b_{i−1}) + μ′′, d_{y}^{∗} − 1), such that

D_{x}(G_{i−1}) = D_{y} (\frac{B̄_{i}}{Δ^{i+1}Δ′^{i}Q^{∗}} ) + \frac{b̄_{i}}{Δ^{i+1}Δ′^{i}Q^{∗}} .

Setting B_{i} = B̃_{i−1} + B̄_{i}Q^{∗}^{i−1}Q^{−} and b_{i} = b̄_{i}, we arrive at (3.8). It remains to verify the degree bounds. The induction hypothesis implies that both deg_{x}(B̄_{i}) and deg_{x}(b_{i}) are bounded by γ + iμ′′−d_{x}^{−}. It follows that deg_{x}(B̄_{i}Q^{∗}^{i−1}Q^{−}) is bounded by γ+iμ′′+(i−1)d_{x}^{∗}. Similarly, deg_{x}(B̃_{i−1}) is bounded by γ +iμ′′+(i−1)d_{x}^{∗}, and so is deg_{x}(B_{i}). The bounds on degrees in y are obvious.

We next derive degree bounds for the minimal telescopers obtained at an intermediate stage of HermiteTelescoping; refined bounds on the output will be given by Theorem 3.4.16.

Lemma 3.4.8. Under (H’), Step 2(c) of Algorithm HermiteTelescoping computes a minimal telescoper L ∈ k[x]⟨D_{x}⟩ with order ρ and a certificate g ∈ k(x, y) for P/Q with deg_{x}(L) ∈ O(d_{x}d_{y}ρ^{2}) and bideg(g) ∈ O(d_{x}d_{y}ρ^{2}) × O(d_{y}ρ).

Proof. By Lemma 3.4.2, we exhibit a minimal telescoper by considering the first nontrivial linear dependence among the a_{i}’s in (3.8). Let M be the coefficient matrix of the system in (η_{i}) obtained from ∑_{i=0}^{ρ} η_{i}a_{i} = 0. By Lemma 3.4.7, M is of size at most (ρ + 1) × d_{y}^{∗} and with coefficients of degree at most σ := d_{x} + μ + ρμ′′ − d_{x}^{−} in x. Hence, there exists a solution (η_{0}, …, η_{ρ}) ∈ k[x]^{ρ+1} of degree at most σρ in x by Fact 2.4.2(ii). Since μ, μ′′ ∈ O(d_{x}d_{y}) and d_{y}^{∗} ≤ d_{y}, the degree estimates of L and g are as announced.

Complexity estimates

We proceed to analyse the complexity of the algorithm in Figure 3.2.

Theorem 3.4.9. Under Hyp. (H’), Algorithm HermiteTelescoping in Figure 3.2 is correct and takes Õ(ρ^{ω+1}d_{x}d_{y}^{2}) ops, where ρ is the order of the minimal telescoper.

Proof. The formulas in Step 2(a) create the loop invariant D_{x}^{i}(f) = D_{y}(g_{i}) + a_{i}/Q^{∗}. Correctness then follows from Lemmas 3.4.1 and 3.4.2. Step 1 takes Õ(d_{x}d_{y}^{2}) ops by Theorem 3.3.6 under (H’). By Lemma 3.4.7, deg_{x}(−a_{i−1}D_{x}(Q^{∗})) ∈ O(id_{x}d_{y}). So the cost for performing Hermite reduction on −a_{i−1}D_{x}(Q^{∗})/Q^{∗}^{2} in Step 2(a) is Õ(id_{x}d_{y}^{2}) ops by Theorem 3.3.6. The bidegrees of g_{i} and a_{i} in Step 2(b) are in O(id_{x}d_{y}) × O(id_{y}) by Lemma 3.4.7. Since adding and differentiating have linear complexity, Step 2(b) takes Õ(i^{2}d_{x}d_{y}^{2}) ops. For each i, the coefficient matrix of ∑_{j=0}^{i} η_{j}a_{j} = 0 in Step 2(c) is of size at most (i + 1) × d_{y}^{∗} and with coefficients 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 Õ(i^{ω}d_{x}d_{y}^{2}) ops by Fact 2.4.2(iii). Computing the content and divisions in Step 3 has complexity Õ(d_{x}d_{y}ρ^{3}). If the algorithm returns when i = ρ, then the total cost is in

∑_{i=0}^{ρ} Õ(i^{2}d_{x}d_{y}^{2}) + ∑_{i=1}^{ρ} Õ(i^{ω}d_{x}d_{y}^{2}) ⊂ Õ(ρ^{ω+1}d_{x}d_{y}^{2}) ops,  (3.9)

which is as announced.

An optimisation, based on Lemma 3.4.5, consists in guessing the order ρ so as to perform Step 2(c) a few times only: As a preprocessing step, choose x_{0} ∈ k lucky for Q, then detect linear dependence of {r_{0,0}, …, r_{0,j}} in (3.7). The minimal j with dependence is a lower bound ρ_{0} on ρ. So Step 2(c) is then performed only when i ≥ ρ_{0}. In practice, the lower bound ρ_{0} computed in this way almost always coincides with the actual order ρ. So normalising the g_{i}’s becomes the dominant step, as observed in experiments. We analyse this optimisation by first estimating the cost for computing ρ_{0}.

Lemma 3.4.10. Under Hypothesis (H’), computing a lower order bound ρ_{0} for minimal telescopers takes Õ(d_{x}d_{y}ρ_{0}^{3}) ops.

Proof. Since differentiating has linear complexity, the derivative D_{x}^{i}(f) takes Õ(i^{2}d_{x}d_{y}) ops. By Fact 2.4.1(i), the evaluation D_{x}^{i}(f)(x_{0}, y) takes as much. The cost of Hermite reduction on D_{x}^{i}(f)(x_{0}, y) is Õ(id_{y}) ops by Lemma 2.4.3. By Fact 2.4.2(iii) with d = 1, computing the rank of the coefficient matrix of ∑_{j=0}^{i} η_{j}r_{0,j}, with r_{0,j} as in (3.7), takes Õ(d_{y}i^{ω−1}) ops. Thus, the cost for computing a lower bound on ρ_{0} is ∑_{i=0}^{ρ_{0}} Õ(i^{2}d_{x}d_{y}) ∈ Õ(d_{x}d_{y}ρ_{0}^{3}) ops.

Corollary 3.4.11. Assume that ρ_{0} = ρ − O(1). Then the previous optimisation of HermiteTelescoping takes Õ(ρ^{3}d_{x}d_{y}^{2}) ops.

Proof. In view of Lemma 3.4.10, the estimate (3.9) becomes

Õ(d_{x}d_{y}ρ_{0}^{3}) + ∑_{i=0}^{ρ} Õ(i^{2}d_{x}d_{y}^{2}) + ∑_{i=ρ_{0}}^{ρ} Õ(i^{ω}d_{x}d_{y}^{2}),

which is Õ(ρ^{3}d_{x}d_{y}^{2}) + Õ((ρ − ρ_{0})ρ^{ω}d_{x}d_{y}^{2}) ops, whence the result.

3.4.2 Improved Almkvist and Zeilberger’s method

We analyse the complexity of Almkvist and Zeilberger’s algorithm [15] when restricted to bivariate rational functions. In order to get a telescoper whose order ρ is minimal, the resulting algorithm, denoted RationalAZ, solves the telescoping equation

∑_{i=0}^{ρ} η_{i}D_{x}^{i}(f) = D_{y}(g)

for increasing, prescribed values of ρ until it gets a solution (η_{0}, …, η_{ρ}, g) ∈ k(x)^{ρ+1} × k(x, y) with the η_{i}’s not all zero. For the analysis, we start by studying the parametrisation of the differential Gosper algorithm of [15] under the same restriction to k(x, y).

Definition 3.4.12 ([49]). Let K be a field and a, b ∈ K[y] be nonzero polynomials. A triple (p, q, r) ∈ K[y]^{3} is said to be a differential Gosper form of the rational function a/b if

\frac{a}{b} = \frac{D_{y}(p)}{p} + \frac{q}{r} and gcd(r, q − τD_{y}(r)) = 1 for all τ ∈ ℕ.

For hyperexponential f, a key step in [15] is to compute a differential Gosper form of the logarithmic derivative of F = ∑_{i=0}^{ρ} η_{i}D_{x}^{i}(f), where the η_{i}’s are undetermined from k(x). In the analogue RationalAZ, this form is predicted by Lemma 3.4.13 below, which is a technical generalisation of a result by Le [63] on F when f has a squarefree denominator.

For f = P/Q, write Q = t(y)T (x, y), splitting content and primitive part with respect to x. By an easy induction, D_{x}^{i}(f) = N_{i}/(QT^{∗}^{i}) for N_{i} ∈ k[x, y]. For this section, set

F = ∑_{i=0}^{ρ} η_{i}D_{x}^{i}(f), N = ∑_{i=0}^{ρ} η_{i}N_{i}T^{∗}^{ρ−i}, and H = −D_{y}(Q)/Q^{−} − ρt^{∗}D_{y}(T^{∗}).

Lemma 3.4.13. If F is nonzero, the triple (N, H, Q^{∗}) is a differential 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^{∗}. It remains to prove gcd(Q^{∗}, H −τD_{y}(Q^{∗})) = 1, for any τ ∈ ℕ. Recall that the squarefree part Q^{∗} of Q is the product Q_{1}Q_{2} ⋯ Q_{m} and that Q̂_{i} denotes Q^{∗}/Q_{i}. By Fact 3.2.1(ii),

Z := H − τD_{y}(Q^{∗}) = −ρt^{∗}D_{y}(T^{∗}) − ∑_{i=1}^{m}(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 both cases, Z is coprime with Q^{∗}, as j > 0, τ ≥ 0, and ρ ≥ 0.

By induction, we have bideg(N_{i}) ≤ (deg_{x}(P) + i deg_{x}(T^{∗}) − i, d_{y} + i deg_{y}(T^{∗}) − 1), so that bideg(N) ≤ (deg_{x}(P) + ρ deg_{x}(T^{∗}) − ρ, d_{y} + ρ deg_{y}(T^{∗}) − 1).

Algorithm RationalAZ Input: f = P/Q ∈ k(x, y) satisfying Hypothesis (H). Output: A minimal telescoper L ∈ k[x]⟨D_{x}⟩ with certificate g ∈ k(x, y).

1. Compute Q^{−} = gcd(Q, D_{y}(Q)), Q^{∗} = Q/Q^{−}, and T, T^{∗} primitive parts of Q, Q^{∗} with respect to x, respectively;

2. Set (Ñ, N, β, H) to (P, P, d_{y}^{−}, −Q^{∗}D_{y}(Q)/Q);

3. For ℓ = 0, 1, … do

(a) Set z to ∑_{j=0}^{β} z_{j}y^{j}, extract the linear system M (η_{i} z_{j})^{T} = 0 from (3.10) (for ρ = ℓ) and compute a basis S of the null space of M by [84].

(b) If S contains a solution (η_{0}, …, η_{ℓ}, s) such that η_{0}, …, η_{ℓ} are not all nonzero, then set (L, g) := (∑_{i=0}^{ℓ} η_{i}D_{x}^{i}, s/(Q^{−}T^{∗}^{ℓ})), and go to Step 4;

(c) Update Ñ := D_{x}(Ñ)T^{∗} − Ñ(T^{∗}D_{x}(T)/T + iD_{x}(T^{∗})), N := NT^{∗} + η_{ℓ+1}Ñ,

β := β + deg_{y}(T^{∗}), and H := H − t^{∗}D_{y}(T^{∗}).

4. Compute the content c of L and return (c^{−1}L, c^{−1}g).

Figure 3.3: Improved Almkvist–Zeilberger algorithm

The next step in RationalAZ is, for fixed ρ, to reduce (3.1) by the change of unknown g = z/(Q^{−}T^{∗}^{ρ}), so as to determine all (η_{i}) ∈ k(x)^{ρ+1} for which the differential equation in z

∑_{i=0}^{ρ} η_{i}N_{i}T^{∗}^{ρ−i} = Q^{∗}D_{y}(z) + (D_{y}(Q^{∗}) + H) z  (3.10)

has a polynomial solution in k(x)[y]. For later use, we recall the following consequence of Corollary 9.6 in [49].

Lemma 3.4.14. 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) − 1. Let c ∈ K[y] be such that β ≥ deg_{y}(c) − deg_{y}(a) + 1. If u is a polynomial solution of aD_{y}(z) + bz = c, then deg_{y}(u) ≤ β.

The following lemma generalizes Lemma 2 in [63] to present a degree bound for z.

Lemma 3.4.15. If u ∈ k(x)[y] is a solution of (3.10) 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 definition of H, b = −Q^{∗}D_{y}(Q^{−})/Q^{−} − ρt^{∗}D_{y}(T^{∗}). Fact 3.2.1(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) < d_{y} + ρ deg_{y}(T^{∗}) and d_{y} = d_{y}^{∗} + d_{y}^{−}, β ≥ deg_{y}(N) − d_{y}^{∗} + 1. The lemma holds by Lemma 3.4.14.

We end the present section using the approach of Almkvist and Zeilberger to provide tight degree bounds on the outputs from Algorithms HermiteTelescoping and RationalAZ.

Theorem 3.4.16. Under Hypothesis (H’), there exists a minimal telescoper L ∈ k[x]⟨D_{x}⟩ 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 3.4.4, there exists a smallest ρ ∈ ℕ at most d_{y}^{∗}, for which (3.1) has a solution with the η_{i}’s not all zero. For this ρ, we estimate the size of the polynomial matrix M derived from (3.10) by undetermined coefficients. By the remark on N after Lemma 3.4.13, 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^{∗}) − 1 ∈ O(ρd_{y}). The matrix M contains two blocks M_{1} ∈ k[x]_{≤n_{x}}^{(n_{y}+1)×(ρ+1)} and M_{2} ∈ k[x]_{≤d_{x}}^{(n_{y}+1)×(β+1)}, where β ∈ O(ρd_{y}) is the same as in Lemma 3.4.15. By the minimality of ρ, the dimension of the null space of M is 1. So there exists u ∈ k[x]^{n_{y}+1} with coefficients of degree at most n_{x}(ρ + 1) + d_{x}(β + 1) ∈ O(d_{x}d_{y}d_{y}^{∗}) in x such that M (η z)^{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. 3.3.

Theorem 3.4.17. Under Hypothesis (H’), Algorithm RationalAZ in Figure 3.3 works correctly and takes Õ(d_{x}d_{y}^{ω}ρ^{ω+2}) ops, where ρ is the order of the minimal telescoper.

Proof. By the existence of a telescoper, Corollary 3.4.4, and Lemma 3.4.15, the algorithm always terminates and returns a minimal telescoper L, of order ρ at most d_{y}^{∗}. Gcd computations dominate the cost of Steps 1 and 2, which take Õ(d_{x}d_{y}^{2}) ops. For each ℓ ∈ ℕ, the dominating cost in Step 3 is computing the null space of M. Let n_{y} = d_{y} + ℓ deg_{y}(T^{∗}) − 1 ∈ O(ℓd_{y}) and n_{x} = d_{x} + ℓ deg_{x}(T^{∗}) ∈ O(ℓd_{x}). By the same argument as in the proof of Theorem 3.4.16, the matrix M is of size at most (n_{y} +1)×(ℓ+β +2) and with coefficients of degree at most n_{x}. Let r be the rank of M, which is either ℓ + β + 2 or ℓ + β + 1 by construction. Thus, a basis of the null space of M can be computed within Õ(n_{x}(n_{y} + 1)(ℓ + β + 2)r^{ω−2}) ops by Fact 2.4.2(iii). Since β ∈ O(ℓd_{y}), Õ(n_{x}(n_{y} + 1)(ℓ + β + 2)r^{ω−2}) is included in Õ(d_{x}d_{y}^{ω}ℓ^{ω+1}). Since Step 3 terminates at ℓ = ρ, the total cost of the algorithm is ∑_{ℓ=0}^{ρ} d_{x}d_{y}^{ω}ℓ^{ω+1} ops. This is within the announced complexity, Õ(d_{x}d_{y}^{ω}ρ^{ω+2}) ops.

Corollary 3.4.18. Algorithms HermiteTelescoping and RationalAZ in Figures 3.2 and 3.3 both output the primitive minimal telescoper L together with its certificate g, which satisfy deg_{D_{x}}(L) ≤ d_{y}^{∗}, deg_{x}(L) and 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 minimal telescoper, as they compute a minimal telescoper at an intermediate step, and owing to their last step of content removal. Bounds follow from Corollary 3.4.4 and Theorem 3.4.16.

3.5 Non-minimal telescopers for bivariate rational functions

In this section, we trade the minimality of telescopers for smaller total output sizes. To this end, we adapt and slightly extend the arguments in [67] and [20, Section 3].

Let f = P/Q be a nonzero element in k(x, y) such that P, Q ∈ k[x, y] and gcd(P, Q) = 1. Denote d_{x} = max{deg_{x}(P), deg_{x}(Q)} and d_{y} = max{deg_{y}(P), deg_{y}(Q)}. Our next goal is to find a nonzero operator A(x, D_{x}, D_{y}) in k[x]⟨D_{x}, D_{y}⟩ such that A(f) = 0. To this end, we consider the k-vector space 𝕎_{N} generated by the set { x^{i}D_{x}^{j}D_{y}^{ℓ} | i + j + ℓ ≤ N } over k. By an easy combinatorial counting, the dimension of 𝕎_{N} is \binom{N+3}{3} over k. Furthermore, for any (i, j, ℓ) ∈ ℕ^{3}, a direct calculation yields

x^{i}D_{x}^{j}D_{y}^{ℓ}(f) = \frac{P_{i,j,ℓ}}{Q^{i+j+ℓ+1}} ,  (3.11)

where P_{i,j,ℓ} ∈ k[x, y] with deg_{x}(P_{i,j,ℓ}) ≤ (i + j + ℓ + 1)d_{x} + i and deg_{y}(P_{i,j,ℓ}) ≤ (i + j + ℓ + 1)d_{y}. So the set 𝕎_{N}(f) is included in the set

𝕍_{N} = span_{k} \{\frac{x^{i}y^{j}}{Q^{N+1}} | i ≤ (N + 1)d_{x} + N, j ≤ (N + 1)d_{y}\} ,

where the dimension of 𝕍_{N} is (N + 1)(d_{x} + 1)((N + 1)d_{y} + 1) over k. Define a linear map ϕ : 𝕎_{N} → 𝕍_{N} by ϕ(L) = L(f) for any L ∈ 𝕎_{N}. Choosing N = 6(d_{x} + 1)(d_{y} + 1) yields the inequality

\binom{N + 3}{3} > (N + 1)(d_{x} + 1)((N + 1)d_{y} + 1),

which implies the kernel of ϕ is nontrivial whenever N ≥ 6(d_{x} + 1)(d_{y} + 1). Therefore, there exists a nonzero operator A in k[x]⟨D_{x}, D_{y}⟩ with total degree at most 6(d_{x} + 1)(d_{y} + 1) in x, D_{x}, and D_{y} that annihilates f. Moreover, A can be found by solving a linear system of size O((d_{x}d_{y})^{3}) over k. The following lemma shows that one can construct a telescoper for f from any y-free annihilator A(x, D_{x}, D_{y}) of f.

Lemma 3.5.1. Let f be a nonzero rational function in k(x, y) and A ∈ k[x]⟨D_{x}, D_{y}⟩ be a nonzero operator such that A(f) = 0. Then there exists a nonzero operator L(x, D_{x}) in k[x]⟨D_{x}⟩ such that L(f) = D_{y}(g) for some g ∈ k(x, y).

Proof. Since D_{y} commutes with x and D_{x}, we can write A = D_{y}^{m}(L(x, D_{x})+D_{y}M), where m ∈ ℕ, M ∈ k[x]⟨D_{x}, D_{y}⟩, and L is a nonzero operator in k[x]⟨D_{x}⟩. By a differential extension of a “non-commutative” trick, used by Wegschaider in [91, Theorem 3.2], there exist w ∈ k[y] and nonzero r ∈ k such that

wD_{y}^{m} = D_{y}Q + r,  (3.12)

where Q ∈ k[y]⟨D_{y}⟩. In particular, r = (−1)^{m}m! ≠ 0 if we take w = y^{m}. Using the fact rD_{y} = D_{y}r and (3.12), we find

\frac{y^{m}}{(−1)^{m}m!}A = L + D_{y}G, where G ∈ k[x, y]⟨D_{x}, D_{y}⟩.

Since A(f) = 0, L(f) = D_{y}(−G(f)). Note that g = −G(f) is still in k(x, y). This completes the proof.

From Lipshitz’s argument, we see that there exists a telescoper of order in O(d_{x}d_{y}). Now, we slightly extend the argument in [20, Section 3] to show that a telescoper of order linear in d_{y} exists. Instead of taking total degree, set 𝕎_{N_{x},N_{∂}} to the k-vector space

𝕎_{N_{x},N_{∂}} = span_{k} { x^{i}D_{x}^{j}D_{y}^{ℓ} | i ≤ N_{x}, j + ℓ ≤ N_{∂} } .

where the dimension of 𝕎_{N_{x},N_{∂}} is (N_{x} + 1)\binom{N_{∂}+2}{2} over k. By (3.11), the set 𝕎_{N_{x},N_{∂}}(f) is included in the set

𝕍_{N_{x},N_{∂}} = span_{k} \{\frac{x^{i}y^{j}}{Q^{N+1}} | i ≤ (N_{∂} + 1)d_{x} + N_{x}, j ≤ (N_{∂} + 1)d_{y}\} .

where the dimension of 𝕍_{N_{x},N_{∂}} is ((N_{∂} + 1)d_{x} + N_{x} + 1)((N_{∂} + 1)d_{y} + 1) over k. Choosing N_{x} = 3d_{x}d_{y} and N_{∂} = 6d_{y} yields the inequality

(N_{x} + 1)\binom{N_{∂} + 2}{2} > ((N_{∂} + 1)d_{x} + N_{x} + 1)((N_{∂} + 1)d_{y} + 1).

Therefore, there exists a nonzero operator A in k[x]⟨D_{x}, D_{y}⟩ with degree at most 3d_{x}d_{y} in x and total degree at most 6d_{y} in D_{x} and D_{y} that annihilates f. Again, A can be found by solving a linear system, but of smaller size O(d_{x}d_{y}^{3}) over k. By the construction in Lemma 3.5.1, the order of L is bounded by 6d_{y}, which is linear in d_{y}.

Remark 3.5.2. The bound 6d_{y} has been shown in [20] for rational functions yD_{y}(Q)/Q with Q ∈ k[x, y]. Apagodu and Zeilberger [16] obtain a similar bound for a class of nonrational hyperexponential functions, but their proof does not seem to apply to rational functions, as it heavily relies on the presence of a nontrivial exponential part.

3.6 Implementation and experiments

We describe in this section an implementation of algorithms HermiteTelecoping and RationalAZ, and compare ours with Maple’s routine.

3.6.1 Implementation and examples

The evaluation-interpolation algorithm HermiteEvalInterp for Hermite reduction in Figure 3.1 does not perform well, mainly because Maple’s rational interpolation routines are far too slow. We thus implemented Algorithm HermiteReduce (original version) in [21, Section 2.2] (carefully avoiding redundant extended gcd calculations), and noted that it performs better. We then implemented a variant of Algorithm HermiteTelescoping in Figure 3.2, using HermiteReduce in place of HermiteEvalInterp, and including the optimisation at the end of § 3.4.1, refined by additional modular calculations. The improved Almkvist-Zeilberger’s algorithm RationalAZ is also implemented. Those functions form the module RationalCT.

> eval(RationalCT);

module()

exportSquareFreeParFrac, HermiteReduce, HermiteTelescoping, RationalAZ ; optionpackage; description“Computing the minimal telescoper for a rational function”; end module

The function SquareFreeParFrac is used to compute the partial fraction decomposition of a rational function with respect to the squarefree factorization of its denominator.

> f :=1/(x^3+5*x^2+8*x+4);

f := \frac{1}{x^{3} + 5 x^{2} + 8 x + 4}

> SquareFreeParFrac(f, x, 'pfd');

1, [x, 0, [[x + 2, [1, −1], [1, −1]], [x + 1, [1, 1]]]]

> pfd;

−\frac{1}{x + 2} − \frac{1}{(x + 2)^{2}} + \frac{1}{x + 1}

The function HermiteReduce returns the additive decomposition of a rational function with respect to the specific variable.

> f := 1/(-y+y^2+x)^2;

f := \frac{1}{(−y + y^{2} + x)^{2}}

> HermiteReduce(f, y);

[\frac{−1 + 2 y}{(4 x − 1) (−y + y^{2} + x)}, \frac{2}{(4 x − 1) (−y + y^{2} + x)}]

In Maple, the function for computing minimal telescopers is DEtools[Zeilberger], which works for any hyperexponential inputs.

> f := 1/(-y+y^2+x);

f := \frac{1}{−y + y^{2} + x}

> DEtools[Zeilberger](f, x, y, Dx);

[2 + Dx (4 x − 1), −\frac{−1 + 2 y}{−y + y^{2} + x}]

Our implementation of the algorithms HermiteTelescoping and RationalAZ is as follows.

   > HermiteTelescoping(f, x, y, Dx);

[2 + Dx (4 x − 1), −\frac{−1 + 2 y}{−y + y^{2} + x}]

   > HermiteTelescoping(f, x, y, Dx, 'No_Certificate');

2 + Dx (4 x − 1)

   > RationalAZ(f, x, y, Dx);

[2 + Dx (4 x − 1), −\frac{−1 + 2 y}{−y + y^{2} + x}]

3.6.2 Experimental results

Now, we show some timings by our implementation and others. Random rational functions We use the function randpoly in Maple to generate the numerator P and the denominator Q of a testing rational function with P and Q having the same degree in x and y. We test the following set of rational functions:

f = \frac{P}{Q}, d = deg_{x}(P) = deg_{y}(P) ∈ \{1, 2, …, 7\}.

For brevity, we denote • AZ: the function DEtools[Zeilberger] in Maple 13;

• Hermite: Hermite-reduction based method;

• RatAZ: improved Almkvist–Zeilberger algorithm.

d 1 2 3 4 5 6 7
AZ  0.054 0.158 2.731 64.75 619.0 > hr > hr
RatAZ 0.019 0.059 0.402 4.461 34.13 220.5 792.1
Hermite 0.016 0.057 0.398 2.664 18.80 106.2 422.5

Table 3.2: Timings on random examples (in seconds)

Differential equations for algebraic functions

Let α(x) be a univariate algebraic function over k(x) with minimal polynomial P in k[x, y]. The following lemma shows that one can compute differential equations for algebraic functions via rational-function telescoping. For the proof, see [20, Proposition 2].

Lemma 3.6.1. If L is a telescoper for yD_{y}(P)/P with respect to y, then L(α(x)) = 0.

In the next table, we show the timings on a set of polynomials in ℤ[x, y], which are also generated by the function randpoly. Let d denote the total degree of a polynomial in x and y. The column gfun (Maple 13) shows the timing by the function algeqtodiffeq in Maple package gfun and gfun (Algo) by the one in the updated Algolib 13.0 of INRIA Algorithms project.

d  4 5 6 7 8  9  10
RatAZ  0.30 1.05 4.90 21.6 69.5 237. 846.
Hermite  0.21 0.94 4.53 20.5 84.7 231. 864.
gfun (Maple 13) 0.14 0.75 6.92 79.6 1661 > hr > hr
gfun (Algo)  0.10 0.46 2.44 12.2 52.7 157. 464.

Table 3.3: Timings on computing differential equations for algebraic functions (in seconds)

Differential equations for diagonals

For a bivariate rational power series

f = ∑_{i,j≥0} f_{i,j}x^{i}y^{j} ∈ k(x, y) ∩ k[[x, y]],

define its diagonal by

diag(f) := ∑_{i=0}^{∞} f_{i,i}x^{i}.

The following lemma shows that one can compute differential equations for diagonals of bivariate rational functions via rational-function telescoping. For the proof, see [67].

Lemma 3.6.2. If L is a telescoper for f(y, x/y)/y with respect to y, then L(diag(f)) = 0.

First, we compare the various algorithms on an example borrowed from [74]:

f = \frac{1}{1 − x − y − xy(1 − x^{d})}, where d ∈ ℕ.

d  8 9 10 11 12 13 14 15
AZ  3.53 6.33 13.6 38.5 68.1 145. 263. 368.
RatAZ 5.27 4.63 8.72 16.9 36.1 55.4 99.4 352.
Hermite 2.33 4.52 8.71 18.6 36.1 65.3 121. 169.

Table 3.4: Timing on Pemantle and Wilson’s example (in seconds)

The next test is on an example from Stanley’s book [83] related to the plane walks. Let S_{d} = {(i, j) ∈ ℕ^{2}|i + j = d}, consider

f(x, y, d) = \frac{1}{1 − ∑_{(i,j)∈S_{d}} x^{i}y^{j}} , where 11 ≤ d ≤ 20.

From the table 3.5, we observed that the order of minimal telescopers is d/2 when d is even.

d  11 12 13 14 15 16 17 18 19 20
AZ  48.7 5.72 144. 12.4 400. 23.9 1016. 46.7 > hr. 81.2
RatAZ 43.8 5.61 129. 11.8 269. 27.9 663.4 45.8 2976. 88.4
Hermite 11.7 2.55 31.9 5.71 91.3 12.8 227.8 21.1 617.9 40.3
Order 11 6 13 7 15 8 17 9 19 10

Table 3.5: Timings on the plane walk examples (in seconds)

We conjecture that this pattern is true for all even d and expect a combinatorial explanation.

50

Structure of Multivariate

Hyperexponential-Hypergeometric

Functions

4.1 Introduction

Multivariate hyperexponential-hypergeometric functions are a generalization of usual exponential functions and hypergeometric terms. In their paper [92], Wilf and Zeilberger observed that a large class of identities on special functions are related to the integrals or sums of proper hyperexponential-hypergeometric functions. Moreover, they have shown that elegant and computer-constructible proofs for those identities can be obtained efficiently via Zeilberger’s method of creative telescoping. Hyperexponential-hypergeometric functions also play an important role in factoring modules over Laurent-Ore algebras [93].

Motivation. It is useful to decompose a hyperexponential-hypergeometric function as a product of “simpler” functions of the same kind. For example, Christopher [29] decomposes a hyperexponential function and uses its multiplicative form to compute Liouvillian first integrals [29]. The Ore-Sato theorem [70, 81], which describes the multiplicative structure of hypergeometric terms, is rediscovered and used by Payne [73], Hou [54], Abramov, and Petkovšek [14] to prove that Wilf and Zeilberger’s conjecture holds for hypergeometric terms. Feng, Singer, and Wu present a multiplicative form for bivariate hyperexponential-hypergeometric function to study Liouvillian solutions of linear difference-differential equations [43, 42]. We will decompose a multivariate hyperexponential-hypergeometric function in the hope of studying the general Wilf and Zeilberger’s conjecture.

Our contribution. We decompose a multivariate hyperexponential-hypergeometric function into a product of a rational function, a hyperexponential function, a hypergeometric term and some “simpler” functions. Our result is a generalization of a result by Feng, Singer and Wu [42, Proposition 5]. Combining this result with multivariate extension of Christopher’s theorem [29, 97] and the Ore-Sato theorem, we present a structure theorem for multivariate hyperexponentialhypergeometric functions.

The rest of this chapter is organized as follows. We describe in Section 4.2 an algebraic setting for hyperexponential-hypergeometric functions. Two rational normal forms are reviewed and a new one is introduced in Section 4.3. We describe a structure of compatible rational functions in Section 4.4. Based on this structure, we obtain a multiplicative form of hyperexponentialhypergeometric functions in Section 4.5.

An earlier version of this chapter can be found in [27], which is a joint work (in progress) with Ziming Li.

4.2 Algebraic setting

We describe an algebraic setting for hyperexponential-hypergeometric functions, which is introduced and used in [62, 22, 93, 66]. We will regard hyperexponential-hypergeometric functions as elements in some differential-difference extensions over the field of multivariate rational functions.

4.2.1 Hyperexponential-hypergeometric functions

Let A be a commutative ring. Recall that a derivation δ on A is an additive map from A to itself satisfying the Leibniz rule

δ(ab) = aδ(b) + δ(a)b for all a, b ∈ A.

Let Δ be a finite set of derivations and automorphisms from A to itself. The pair (A, Δ) is called a differential-difference ring, or Δ-ring for short. If the maps in Δ commute pairwise, then (A, Δ) is said to be orthogonal [62].

Let k be a field of characteristic zero. For brevity, we let x stand for the continuous variables (x_{1}, …, x_{m}) and y for the discrete variables (y_{1}, …, y_{n}). Let k(x, y) be the field of rational functions in x_{1}, …, x_{m} and y_{1}, …, y_{n} over k. On the field k(x, y), the derivations δ_{i} (1 ≤ i ≤ m) and shift operators σ_{j} (1 ≤ j ≤ n) are defined for all f ∈ k(x, y) by

δ_{i}(f) = \frac{∂f}{∂x_{i}} and σ_{j}(f) = f(x, y_{1}, …, y_{j−1}, y_{j} + 1, y_{j+1}, …, y_{n}).

Put Δ = {δ_{1}, …, δ_{m}, σ_{1}, …, σ_{n}} and note that the elements of Δ commute pairwise over k(x, y). So the pair (k(x, y), Δ) is a Δ-field.

A ring R is called a Δ-extension of the field k(x, y) if R contains k(x, y), all derivations can be extended to R, all the shift operators can be extended to monomorphisms to R, and the extended operators also commute pairwise over R. The set of all extended derivations and shift operators is still denoted by Δ. The ring R is said to be simple if it contains no ideal closed under all maps in Δ except the zero ideal and the whole ring. An element c of R is called a constant with respect to a derivation δ ∈ Δ if δ(c) = 0 and a constant with respect to a shift operator σ ∈ Δ if σ(c) = c. An element c of R is called a constant if it is a constant with respect to all operators in Δ. The set of all constants of R, denoted by C_{R}, is a subring of R. The ring C_{R} is a subfield if R is a field.

Definition 4.2.1. Let R be a Δ-extension of k(x, y). A nonzero element h of R is said to be hyperexponential-hypergeometric over k(x, y) if there exist a_{1}, …, a_{m}, b_{1}, …, b_{n} ∈ k(x, y) such that b_{1} ⋯ b_{n} ≠ 0 and

δ_{1}(h) = a_{1}h, …, δ_{m}(h) = a_{m}h, and σ_{1}(h) = b_{1}h, …, σ_{n}(h) = b_{n}h.  (4.1)

The rational functions a_{i} and b_{j} above are called the certificates of h with respect to x_{i} and y_{j}, respectively.

According to the definition above, it is easy to see that derivatives, shifts, and products of hyperexponential-hypergeometric functions are also hyperexponential-hypergeometric. If R is simple, then any hyperexponential-hypergeometric element in R is invertible [66, Lemma 2.1].

Since all elements in Δ commute pairwise, the certificates a_{i}’s and b_{j}’s in Definition 4.2.1 satisfy three sets of integrability conditions:

δ_{i}(a_{j}) = δ_{j}(a_{i}) for 1 ≤ i < j ≤ m,  (4.2)
σ_{i}(b_{j})b_{i} = σ_{j}(b_{i})b_{j} for 1 ≤ i < j ≤ n,  (4.3)
\frac{δ_{i}(b_{j})}{b_{j}} = σ_{j}(a_{i}) − a_{i} for 1 ≤ i ≤ m and 1 ≤ j ≤ n.  (4.4)

In the discrete case, rational functions b_{1}, …, b_{n} ∈ k(y) in [14] are said to be compatible if they satisfy the integrability conditions (4.3). We follow this and call a_{1}, …, a_{m}, b_{1}, …, b_{n} in k(x, y) compatible if they satisfy the integrability conditions (4.2), (4.3), and (4.4). As opposed to the treatment in [14], we regard the discrete variables y_{1}, …, y_{n} as indeterminates, and, thus, the function h = |y_{1} − y_{2}| is excluded in the thesis.

4.2.2 First-order fully integrable systems

A Δ-extension of k(x, y) may not be an integral domain in general. For this reason, we need to construct a simple Δ-ring, in which any hyperexponential-hypergeometric element is invertible. We specialize the construction in [22, Section 3] for finitely many first-order systems.

Definition 4.2.2. Let a_{1}, …, a_{m}, b_{1}, …, b_{n} be rational functions in k(x, y). A first-order system

δ_{1}(z) = a_{1}z, …, δ_{m}(z) = a_{m}z, σ_{1}(z) = b_{1}z, …, σ_{n}(z) = b_{n}z,  (4.5)

is said to be fully integrable over k(x, y) if b_{1} ⋯ b_{n} ≠ 0 and a_{1}, …, a_{m}, b_{1}, …, b_{n} are compatible.

According to Theorem 2 in [22], given a finite number of first-order fully integrable systems, there exists a simple differential-difference extension R of k(x, n) such that R contains a nonzero solution of each system. Moreover, the subring of constants in R is equal to k if k is algebraically closed. Let H(a, b) denote the solution space of the first-order fully integrable system (4.5) in such an extension R. By Theorem 2 in [22], H(a, b) is one-dimensional over k if k is algebraically closed. In the rest of this chapter, hyperexponential-hypergeometric functions will be regarded as elements in such an extension R. Therefore, it is legitimate to add, multiply and invert hyperexponential-hypergeometric functions. Following Definition 2 in [14], two hyperexponential-hypergeometric functions are said to be conjugate if they have the same

Function  Expression  Conjugates given by H(a, b)-Notation
exponential  exp(f)  H(δ_{1}(f), …, δ_{m}(f), 1, …, 1)
constant power β^{λ}  H(λ\frac{δ_{1}β}{β} , …, λ\frac{δ_{m}β}{β} , 1, …, 1)
symbolic power β^{y_{j}}  H (y_{j} \frac{δ_{1}β}{β} , …, y_{j} \frac{δ_{m}β}{β} , 1, …, 1, β, 1, …, 1)
factorial  (λ)_{e·y}  H (0, …, 0, ∏_{ℓ=0}^{e_{1}−1}(e · y + λ + ℓ), …, ∏_{ℓ=0}^{e_{n}−1}(e · y + λ + ℓ))
(In the table above, f, β ∈ k(x), λ ∈ k, and e ∈ ℤ^{n}.)

Table 4.1: Familiar functions and their H(a, b)-representations

certificates. Two conjugate functions can only differ by a nonzero multiplicative element in k if k is algebraically closed.

The following lemma shows some basic properties of the space H(a, b). They are derived directly from the definition.

Lemma 4.2.3. (i) For any hyperexponential-hypergeometric function g over k(x, y), we have

gH(a_{1}, …, a_{m}, b_{1}, …, b_{n}) = H (a_{1} + \frac{δ_{1}(g)}{g}, …, a_{m} + \frac{δ_{m}(g)}{g}, b_{1} \frac{σ_{1}(g)}{g}, …, b_{n} \frac{σ_{n}(g)}{g}) .

(ii) H(a, b)H(ã, b̃) = H(a + ã, bb̃), where a + ã and bb̃ are defined termwise.

(iii) δ_{i}(H(a, b)) = a_{i}H(a, b) for 1 ≤ i ≤ m and σ_{j}(H(a, b)) = b_{j}H(a, b) for 1 ≤ j ≤ n.

Some examples for hyperexponential-hypergeometric functions are listed in Table 4.1.

4.3 Rational normal forms

In this section, we review two normal forms for univariate rational functions in [13, 10, 47]. Those normal forms have played an important role in the minimal decompositions of hyperexponential and hypergeometric functions [47, 12, 13, 10]. In order to study the compatible multivariate rational functions in the next sections, we introduce a new kind of rational normal forms for the continuous-discrete rational functions in k(x, y).

4.3.1 Differential and shift rational normal forms

Let F be a field of characteristic zero. The field F (z) of univariate rational functions is equipped with both differential and difference structures by

δ(f(z)) = \frac{d(f(z))}{dz} and σ(f(z)) = f(z + 1) for all f ∈ F (z) .

A polynomial P ∈ F [z] is said to be squarefree with respect to z over F if gcd(P, δ(P)) = 1. It is said to be shift-free if gcd(P, σ^{i}(P)) = 1 for all i ∈ ℤ \ {0}. In other words, any two roots of a shift-free polynomial have a non-integer distance. In the following, we collect basic facts concerning rational functions whose denominators are squarefree or shift-free, which has appeared implicitly in the literature [2, 71, 52, 53]. Lemma 4.3.1. Let f = P/Q be in F (z) with gcd(P, Q) = 1 and deg(P) < deg(Q). Then

(i) If Q is squarefree and f = δ(g) for some g ∈ F (z), then f = 0.

(ii) If Q is shift-free and f = σ(g) − g for some g ∈ F (z), then f = 0. Definition 4.3.2 (Differential-reduced, shift-reduced). A rational function f = P/Q ∈ F (z) with gcd(P, Q) = 1 is said to be differential-reduced with respect to z over F if

gcd(Q, P − iδ(Q)) = 1, for all i ∈ ℤ.

It is said to be shift-reduced with respect to z over F if

gcd(P, σ^{i}(Q)) = 1, for all i ∈ ℤ.

The following lemmas show basic properties of differential-reduced and shift-reduced rational functions, respectively. Those results are seemingly classical in [51, 15, 47, 13], but we still present their proofs for completeness. Lemma 4.3.3. Let f = P/Q ∈ F (z) be differential-reduced with respect to z, where P, Q ∈ F [z] and gcd(P, Q) = 1. If f = δ(g)/g for some g ∈ F (z), then g ∈ F and f = 0. Proof. Suppose that g ∈ F (z) \ F. Then

f = \frac{δ(g)}{g} = ∑_{i=1}^{s} \frac{m_{i}}{z − α_{i}} , where m_{i} ∈ ℤ and α_{i} ∈ \overline{F} for all i with 1 ≤ i ≤ s.

By Theorem 2.3.8, the m_{i}’s are roots of RT_{z}(f). So gcd(Q, P − m_{i}δ(Q)) ≠ 1, which contradicts the assumption that f is differential-reduced. So g ∈ F and therefore f = 0.

Lemma 4.3.4. Let f ∈ F (z) be shift-reduced with respect to z over F. If f = σ(g)/g for some g ∈ F (z), then g ∈ F and f = 1.

Proof. Let g = a/b with a, b ∈ F [z] and gcd(a, b) = 1. Suppose that g is not in F. Then either a or b has a positive degree. Assume that the degree of a is positive. Then there exists an element α ∈ \overline{F} such that α is a root of a, but α − 1 is not. Moreover, there exists a nonnegative integer ℓ such that α + ℓ is a root of a but α + ℓ + 1 is not.

Set β = α + ℓ. Since f = σ(g)/g,

f = \frac{σ(a)b}{aσ(b)}.

We have that α − 1 is a root of σ(a), but not a root of aσ(b) by the definition of α and since gcd(σ(a), σ(b)) = 1. Similarly, β is a root of a, but not a root of σ(a)b by the definition of β and since gcd(a, b) = 1. It follows that α−1 is a root of the numerator of f, while β is a root of the denominator of f. Since β − (α − 1) is an integer, f is not shift-reduced, a contradiction. Hence, a belongs to F. In the same vein, b belongs to F.

Definition 4.3.5 (Differential rational normal form). For f ∈ F (z), call a pair (K, S) ∈ F (z)^{2} a differential rational normal form (abbreviated as DRNF) of f if f = K + δ(S)/S and K is differential-reduced with respect to z over F. If, in addition, the denominators of K and S are coprime, then the pair (K, S) is said to be a strict DRNF of f.

Definition 4.3.6 (Shift rational normal form). For f ∈ F (z), call a pair (K, S) ∈ F (z)×F (z) a shift rational normal form (abbreviated as SRNF) of f if f = K · σ(S)/S and K is shift-reduced with respect to z over F. Let k_{1} = num(K), k_{2} = den(K), s_{1} = num(S) and s_{2} = den(S). If, furthermore,

gcd(k_{1}, σ(s_{2})s_{1}) = gcd(k_{2}, σ(s_{1})s_{2}) = 1,

then the pair (K, S) is said to be a strict SRNF of f.

Both strict DRNF’s and strict SRNF’s exist and can be computed efficiently [47, 13]. Let h(z) be a hyperexponential or hypergeometric function over F (z) with certificate f in F (z). Any differential or shift rational normal form (K, S) of f leads to a multiplicative factorization of h in the form

h(z) = Sh′, where the certificate of h′ is K.

From this factorization, one can perform two kinds of minimal decompositions on h, which will be reviewed in the next chapter.

4.3.2 𝕐-rational normal forms

We introduce a new normal form for rational functions in k(x, y), which will help us derive the structure of certificates in the continuous-discrete setting.

For brevity, we set F_{i} = k(x_{1}, …, x_{i−1}, x_{i+1}, … x_{m}, y) for all i ∈ {1, …, m}, and abbreviate “irreducible partial fraction decomposition” as IPFD.

For a nonzero element r ∈ k(x, y), its IPFD with respect to x_{i} can be written as

r = p_{0} + ∑_{i=1}^{ℓ} ∑_{j=1}^{n_{i}} \frac{p_{ij}}{q_{i}^{j}}

where p_{0}, p_{ij}, q_{i} are in F_{i}[x_{i}] such that deg_{x_{i}}(p_{ij}) < deg_{x_{i}}(q_{i}) and q_{i} is monic and irreducible.

Let 𝕐 be the additive group generated by 1, y_{1}, …, y_{n}. Let r and its IPFD with respect to x_{i} be given above. We define r_{S} to be the sum of all fractions of the form p_{i,1}/q_{i} with the following two properties

(i) p_{i,1} = yδ_{i}(q_{i}) for some nonzero element y ∈ 𝕐, and

(ii) n_{i} = 1, that is, q_{i}^{2} does not divide the denominator of r.

Furthermore, set r_{K} = r − r_{S}. We call r_{K} the 𝕐-kernel, and r_{S} the 𝕐-shell of r with respect to x_{i}. For completeness, the 𝕐-kernel and 𝕐-shell of zero are both defined to be zero.

Note that a 𝕐-shell with respect to x_{i} can be written as

r_{S} = \frac{δ_{i}(f)}{f} + ∑_{j=1}^{n} y_{j} \frac{δ_{i}(u_{j})}{u_{j}}

for some f, u_{j} ∈ k(x, y). Moreover, the denominators of the 𝕐-kernel and 𝕐-shell of a nonzero rational function are coprime.

Lemma 4.3.7. Let r be a nonzero rational function in k(x, y) with 𝕐-kernel r_{K} and 𝕐-shell r_{S} with respect to x_{i}. Let g = yδ_{i}(q)/q, where y is a nonzero element of 𝕐, and q is a monic and irreducible polynomial in F_{i}[x_{i}] with positive degree.

(i) If q is a divisor of the denominator of r_{K}, then the 𝕐-kernel and 𝕐-shell of r +g are r_{K} +g and r_{S}, respectively.

(ii) If q is not a divisor of the denominator of r_{K}, then the 𝕐-kernel and 𝕐-shell of r +g are r_{K} and r_{S} + g, respectively.

Proof. Since the denominator q of g is irreducible, it suffices to look at q-expansions of r instead of its IPFD. The q-adic expansion of r with respect to x_{i} can be written as:

r = \frac{p_{ℓ}}{q^{ℓ}} + ⋯ + \frac{p_{2}}{q^{2}} + \frac{p_{1}}{q} + higher terms,  (4.6)

where p_{ℓ}, …, p_{2}, p_{1} ∈ F_{i}[x_{i}] with degrees less than the degree of q in x_{i}. It follows that the q-adic expansion of r + g with respect to x_{i} is

r + g = \frac{p_{ℓ}}{q^{ℓ}} + ⋯ + \frac{p_{2}}{q^{2}} + \frac{p_{1} + yδ(q)}{q} + higher terms.  (4.7)

Assume that q is a divisor of the denominator of r_{K}. Then one of p_{1}, p_{2}, …, p_{ℓ} is nonzero. If one of p_{2}, …, p_{ℓ} is nonzero, then there is no fraction with denominator q appearing in the IPFD of the 𝕐-shell of r + g. The first assertion holds. If p_{2}, …, p_{ℓ} are all equal to zero, then p_{1} is nonzero. By the definition of 𝕐-shells, p_{1} is not a product of δ(q) and a nonzero element in 𝕐, and nor is p_{1} + yδ_{i}(q). Hence, the fraction (p_{1} + yδ_{i}(q))/q appears in the IPFD of the 𝕐-kernel of r + g. The first assertion holds again.

Assume that q is not a divisor of the denominator of r_{K}. Then p_{1} in (4.6) is of the form y′δ_{i}(q) for some y′ ∈ 𝕐 (noting that y′ may possibly be zero). Moreover, all p_{2}, …, p_{ℓ} are equal to zero. Hence, (4.7) becomes

r + g = (y′ + y)\frac{δ_{i}(q)}{q} + higher terms.

It follows that the 𝕐-shell of r + g is equal to r_{S} + g.

We now define five additive subgroups of k(x, y) in order to suppress complicated expressions that would appear. Let

L_{i} = \{\frac{δ_{i}(f)}{f} | f ∈ k(x, y), f ≠ 0\} ,

M_{i} = \{∑_{ℓ=1}^{n} y_{ℓ} \frac{δ_{i}(u_{ℓ})}{u_{ℓ}} | u_{ℓ} ∈ k(x), u_{ℓ} ≠ 0\} ,

and

M_{i,j} = \{∑_{ℓ=1}^{j−1} y_{ℓ} \frac{δ_{i}(v_{ℓ})}{v_{ℓ}} + ∑_{ℓ=j+1}^{n} y_{ℓ} \frac{δ_{i}(v_{ℓ})}{v_{ℓ}} | v_{ℓ} ∈ k(x, y_{j}), v_{ℓ} ≠ 0\} .

Moreover, let N_{i} = L_{i} + M_{i} + k(x) and N_{i,j} = L_{i} + M_{i,j} + k(x, y_{j}).

Remark 4.3.8. Let G be one of the subgroups L_{i}, M_{i}, N_{i}, M_{i,j}, N_{i,j}. If r belongs to G, so does every fraction appearing in the IPFD of r with respect to x_{i}.

The next lemma is rather technical. It tells us how the 𝕐-shell and 𝕐-kernel of an element in N_{i,j} look like.

Lemma 4.3.9. Let i ∈ {1, …, m} and j ∈ {1, …, n}. Let a be a nonzero element of N_{i,j} with 𝕐-kernel a_{K} and 𝕐-shell a_{S} with respect to x_{i}. Then there exists u ∈ k(x, y_{j}) such that

a_{K} ≡ 0 mod M_{i,j} + k(x, y_{j}) and a_{S} ≡ y_{j} \frac{δ_{i}(u)}{u} mod L_{i} + M_{i,j}.

Proof. Without loss of generality, assume that i = 1 and j = n. All the IPFD’s, 𝕐-kernels and 𝕐-shells in the proof are with respect to x_{1}. Our goal is to show that

a_{K} ≡ 0 mod M_{1,n} + k(x, y_{n}) and a_{S} ≡ y_{n} \frac{δ_{1}(u)}{u} mod L_{1} + M_{1,n}  (4.8)

for some u ∈ k(x, y_{n}).

Since a ∈ N_{1,n}, there exists r ∈ k(x, y_{n}) such that

a ≡ r mod L_{1} + M_{1,n}.

Assume that r_{K} and r_{S} are the 𝕐-kernel and 𝕐-shell of r, respectively. Then

r_{K} ∈ k(x, y_{n}) and r_{S} = \frac{δ_{1}(f)}{f} + y_{n} \frac{δ_{1}(u)}{u}.

where f, u ∈ k(x, y_{n}). Using congruences, we have

r_{K} ≡ 0 mod M_{1,n} + k(x, y_{n}) and r_{S} ≡ y_{n} \frac{δ_{1}(u)}{u} mod L_{1} + M_{1,n}  (4.9)

To proceed, we denote by p the denominator of r_{K}, which is viewed as a polynomial in x_{1} over k(x_{2}, …, x_{m}, y_{n}). Let b = a − r, which is in L_{1} + M_{1,n}. If b = 0, then there is nothing to prove. Assume that b ≠ 0, and that g is a fraction appearing in the IPFD of b. By the definitions of L_{1} and M_{1,n}, we have that g = yδ_{1}(q)/q, where y is a nonzero element in 𝕐 free of y_{n}, and q is a monic and irreducible polynomial in x_{1} over k(x_{2}, …, x_{m}, y). We now make a case distinction. Case 1. If q is a divisor of p, then we set r′_{K} = r_{K} + g and r′_{S} = r_{S}. Case 2. If q is not a divisor of p, then we set r′_{K} = r_{K} and r′_{S} = r_{S} + g.

Let r′ = r′_{K} +r′_{S}. By Lemma 4.3.7, r′_{K} and r′_{S} are the 𝕐-kernel and 𝕐-shell of r′, respectively. In Case 1, q is a factor of p. So q is in k(x, y_{n}). Accordingly, g is in M_{1,n} + k(x, y_{n}). It follows from (4.9) that

r′_{K} ≡ 0 mod M_{1,n} + k(x, y_{n}) and r′_{S} ≡ y_{n} \frac{δ_{1}(u)}{u} mod L_{1} + M_{1,n}.

By Remark 4.3.8, g is in L_{1} + M_{1,n}, which, together with (4.9), implies that the above two congruences also hold in Case 2.

Since a − r′ = a − r − g and g ∈ L_{1} + M_{1,n}, we have

a ≡ r′ mod L_{1} + M_{1,n}.

Let b′ = a − r′. Then the number of fractions in the IPFD of b′ is less than that in the IPFD of b. Repeating the same argument on b′, we see that (4.8) follows from an easy induction.

The lemma below is useful to generalize Proposition 5 in [42].

Lemma 4.3.10. Let i ∈ {1, …, m} and p, q ∈ {1, …, n} with p ≠ q. Then

N_{i,p} ∩ N_{i,q} = N_{i}.

Proof. Without loss of generality, we assume that i = 1, p = 1, q = n and n > 1. Since

M_{1} ⊂ M_{1,1} + k(x, y_{1}) and M_{1} ⊂ M_{1,n} + k(x, y_{n}),

we have N_{1} is included in the intersection N_{1,1} and N_{1,n} by definition. It suffices to show the opposite inclusion. Assume that a is a nonzero element in the intersection of N_{1,1} and N_{1,n}. Let a_{S} and a_{K} be, respectively, the 𝕐-shell and 𝕐-kernel of a with respect to x_{1}. By Lemma 4.3.9, the 𝕐-shell of a with respect to x_{1} is of the form

a_{S} = \frac{δ_{1}(f_{1})}{f_{1}} + ∑_{j=1}^{n} y_{j} \frac{δ_{1}(v_{j})}{v_{j}} = \frac{δ_{1}(f′_{1})}{f′_{1}} + ∑_{j=1}^{n} y_{j} \frac{δ_{1}(v′_{j})}{v′_{j}},

where f_{1}, f′_{1} ∈ k(x, y), v_{j} ∈ k(x, y_{1}) and v′_{j} ∈ k(x, y_{n}) for all j with 1 ≤ j ≤ n. Since the elements 1, y_{1}, …, y_{n} are ℤ-linearly independent constants with respect to δ_{1}, Corollary 2.3.5 implies that

\frac{δ_{1}(f_{1})}{f_{1}} = \frac{δ_{1}(f′_{1})}{f′_{1}} and \frac{δ_{1}(v_{j})}{v_{j}} = \frac{δ_{1}(v′_{j})}{v′_{j}}

for all j with 1 ≤ j ≤ n. Since the differential equation

δ_{1}(z) = \frac{δ_{1}(v′_{j})}{v′_{j}}z

has a solution v_{j} ∈ k(x, y_{1}) and its coefficients belong to k(x, y_{n}), it must have a solution u_{j} in k(x). Hence,

a_{S} = \frac{δ_{1}(f_{1})}{f_{1}} + ∑_{j=1}^{n} y_{j} \frac{δ_{1}(u_{j})}{u_{j}} ∈ N_{1}.

It remains to show that the 𝕐-kernel of a with respect to x_{1} belongs to N_{1}. By Lemma 4.3.9, the 𝕐-kernel of a is of the form

a_{K} = y_{n} \frac{δ_{1}(w_{n})}{w_{n}} + ∑_{ℓ=2}^{n−1} y_{ℓ} \frac{δ_{1}(w_{ℓ})}{w_{ℓ}} + s = y_{1} \frac{δ_{1}(w′_{1})}{w′_{1}} + ∑_{ℓ=2}^{n−1} y_{ℓ} \frac{δ_{1}(w′_{ℓ})}{w′_{ℓ}} + s′,

where w_{2}, …, w_{n}, s ∈ k(x, y_{1}) and w′_{1}, …, w′_{n−1}, s′ ∈ k(x, y_{n}). Note that a_{K} is a polynomial in y_{1} over k(x, y_{2}, …, y_{n}) whose degree (in y_{1}) is less than or equal to 1 by the second equality given above. So w_{2}, …, w_{n} in the first equality can be taken as elements in k(x). For, otherwise, the denominator of a_{K} would involve y_{1}, a contradiction. It follows that

s = cy_{1} + d,

where c, d ∈ k(x). In the same vein, we have

a_{K} = y_{1} \frac{δ_{1}(w′_{1})}{w′_{1}} + ∑_{ℓ=2}^{n−1} y_{ℓ} \frac{δ_{1}(w′_{ℓ})}{w′_{ℓ}} + c′y_{n} + d′,

where w′_{1}, …, w′_{n−1}, and c′, d′ are chosen to be in k(x). It follows that c = δ_{1}(w′_{1})/w′_{1}. Consequently,

a_{K} = y_{1} \frac{δ_{1}(w′_{1})}{w′_{1}} + ∑_{ℓ=2}^{n} y_{ℓ} \frac{δ_{1}(w_{ℓ})}{w_{ℓ}} + d,

which belongs to N_{1}. We have proved that both a_{S} and a_{K} are in N_{1}, and so is a.

The lemma below can be viewed as a discrete analogue of Lemma 4.3.10.

Lemma 4.3.11. Let p, q ∈ {1, …, m} with p ≠ q, and j ∈ {1, …, n}. If a nonzero rational function b ∈ k(x, y) can be written as

b = \frac{σ_{j}(f_{p})}{f_{p}}β_{p}α_{p} = \frac{σ_{j}(f_{q})}{f_{q}}β_{q}α_{q}  (4.10)

for some f_{p}, f_{q} ∈ k(x, y), β_{p}, β_{q} ∈ k(x), α_{p} ∈ k(x_{p}, y) and α_{q} ∈ k(x_{q}, y), then there exist f ∈ k(x, y), α ∈ k(y), and β ∈ k(x) such that

b = \frac{σ_{j}(f)}{f}βα.

Proof. Without loss of generality, we may assume that both α_{p} and α_{q} are shift-reduced with respect to y_{j}. Suppose that

\frac{α_{p}}{α_{q}} = \frac{σ_{j}(g)}{g} \frac{1}{u},  (4.11)

in which (g, 1/u) is a shift rational normal form of α_{p}/α_{q} with respect to y_{j}. Then g belongs to k(y), because x_{p} and x_{q} are two distinct indeterminates, and α_{p}, α_{q} are shift-reduced. It follows from (4.10) and (4.11) that

\frac{σ_{j}(w)}{w} = \frac{β_{q}}{β_{p}} u,

where w = gf_{p}/f_{q}. Since u is shift-reduced with respect to y_{j}, and β_{p}, β_{q} ∈ k(x), the above equation and Lemma 4.3.4 imply that σ_{j}(w) = w, and

\frac{β_{q}}{β_{p}} u = 1,

and, hence, u ∈ k(x). By (4.11), α_{p} ∈ k(x_{p}, y) is a product of an element in k(y), an element in k(x), and α_{q} in k(x_{q}, y). By the uniqueness of the factorization of rational functions, α_{p} can be written as a product cd for some c ∈ k(x_{p}) and d ∈ k(y). Setting f = f_{p}, β = β_{p}c, and α_{p} = d yields the lemma.

4.4 Structure of compatible rational functions

The goal of this section is to show the following theorem, which describes the structure of compatible rational functions. To this end, let us recall a useful notation from Payne’s thesis [73].

For any integers a, b ∈ ℤ and a sequence of expressions A_{i}, define

\begin{matrix}∏ & {b \atop a} \\ i \end{matrix} A_{i} = \{ \begin{matrix}∏_{i=a}^{b−1} A_{i}, & if b > a; \\ 1, & if b = a; \\ 1/ ∏_{i=b}^{a−1} A_{i}, & if b < a. \end{matrix}

Theorem 4.4.1. Assume that a_{1}, …, a_{m}, b_{1}, …, b_{n} ∈ k(x, y) are m + n rational functions such that b_{1} ⋯ b_{n} ≠ 0, and that all the equalities in (4.2), (4.3) and (4.4) hold. Then there exist f ∈ k(x, y), univariate rational functions r_{v} ∈ k(z) for each v in a finite set V ⊂ ℤ^{n}, c_{1}, …, c_{L} ∈ k̄, g_{0}, β_{1}, …, β_{n} ∈ k(x), and g_{1}, …, g_{L} ∈ k̄(x) such that

a_{i} = δ_{i}(g_{0}) + \frac{δ_{i}(f)}{f} + ∑_{ℓ=1}^{L} c_{ℓ} \frac{δ_{i}(g_{ℓ})}{g_{ℓ}} + ∑_{j=1}^{n} y_{j} \frac{δ_{i}(β_{j})}{β_{j}} for all i with 1 ≤ i ≤ m,  (4.12)

and

b_{j} = \frac{σ_{j}(f)}{f}β_{j} ∏_{v∈V} \begin{matrix}∏ & {v_{j} \atop 0} \\ p \end{matrix} r_{v}(y · v + p) for all j with 1 ≤ j ≤ n.  (4.13)

where y · v denotes the inner product y_{1}v_{1} + … + y_{n}v_{n}.

Before proving this theorem, we review two special cases in Sections 4.4.1 and 4.4.2.

4.4.1 The Ore-Sato theorem

The structure of rational solutions of the recurrence equation

R_{1}(m, n + 1)R_{2}(m, n) = R_{1}(m, n)R_{2}(m + 1, n)

has been described by Ore [70]. The multivariate extension of Ore’s theorem was obtained by Sato [81] in the 1960s when he developed the theory of prehomogeneous vector spaces. In the process of proving the discrete case of Wilf and Zeilberger’s conjecture on holonomic hypergeometric terms, the Ore-Sato theorem was discovered and proved again by Payne in his Ph.D. thesis [73] and independently by Abramov and Petkovšek [14]. In particular, the case of two variables has also been shown by Hou [54, 55] and by Abramov and Petkovšek [12]. Their results reveal a multiplicative structure of the certificates b_{1}, …, b_{n} of a hyperexponential-hypergeometric function from its integrability conditions (4.3). The following presentation of the Ore-Sato theorem is taken from Payne’s thesis [73, Theorem 2.8.4].

Theorem 4.4.2 (Ore-Sato theorem). Let b_{1}, …, b_{n} ∈ k(y) be nonzero rational functions such that

b_{i}σ_{i}(b_{j}) = b_{j}σ_{j}(b_{i}), for 1 ≤ i < j ≤ n.

Then there exist a rational function f ∈ k(y), a finite set V ⊂ ℤ^{n}, and univariate rational functions r_{v} in k(z) for each v ∈ V such that for all j with 1 ≤ j ≤ n,

b_{j} = \frac{σ_{j}(f)}{f} ∏_{v∈V} \begin{matrix}∏ & {v_{j} \atop 0} \\ p \end{matrix} r_{v}(y · v + p),

where y · v denotes the inner product y_{1}v_{1} + … + y_{n}v_{n}.

4.4.2 Alternative proof of multivariate Christopher’s theorem

A multivariate hyperexponential function h(x) is a nonzero solution of the first-order fully integrable system

δ_{1}(z) = a_{1}z, …, δ_{m}(z) = a_{m}z,

where a_{1}, …, a_{m} are compatible rational functions in k(x), i.e., they satisfy the integrability conditions (4.3). In his refinement of Singer’s theorem [82] on Liouvillian first integrals, Christopher [29] has described a possible form for those compatible a_{i}’s in the bivariate case.

Theorem 4.4.3 (Christopher, 1999). Any bivariate hyperexponential function h(x_{1}, x_{2}) over the field ℂ(x_{1}, x_{2}) can be written as

exp(f) ∏_{ℓ=1}^{L} g_{ℓ}^{c_{ℓ}}, where f, g_{ℓ} ∈ ℂ(x_{1}, x_{2}) and c_{ℓ} ∈ ℂ.

Consequently, the two certificates a_{1} and a_{2} of h(x_{1}, x_{2}) can be written as

a_{i} = δ_{i}(f) + ∑_{ℓ=1}^{L} c_{ℓ} \frac{δ_{i}(g_{ℓ})}{g_{ℓ}}, for i = 1, 2.

Zoladek [97] has shown the multivariate extension of Christopher’s theorem by using a result in Cerveau and Mattei’s book [25]. We offer an alternative proof involving only rational functions.

Theorem 4.4.4 (Multivariate Christopher’s theorem). Let a_{1}, …, a_{m} ∈ k(x) be rational functions such that

δ_{i}(a_{j}) = δ_{j}(a_{i}), for 1 ≤ i < j ≤ m.

Then there exist f ∈ k(x), nonzero elements c_{ℓ} ∈ \overline{k}, and g_{ℓ} ∈ k(c_{ℓ})(x) for 1 ≤ ℓ ≤ n such that

a_{i} = δ_{i}(f) + ∑_{ℓ=1}^{n} c_{ℓ} \frac{δ_{i}(g_{ℓ})}{g_{ℓ}}, for i = 1, …, m.

The proof proceeds by induction on m. To this end, we first present a lemma for the induction step. In the proof of the lemma below, we use the formula:

δ_{i} (\frac{δ_{j}(f)}{f}) = δ_{j} (\frac{δ_{i}(f)}{f}) for all nonzero f ∈ k(x) and 1 ≤ i < j ≤ m.  (4.14)

Lemma 4.4.5. Let K denote k(x_{2}, x_{3}, …, x_{m}) and a_{1}, …, a_{m} ∈ k(x) be rational functions such that

δ_{i}(a_{j}) = δ_{j}(a_{i}) for all i, j with 1 ≤ i < j ≤ m.

Then there exist f ∈ k(x), A_{ℓ} ∈ K with ℓ = 2, …, m, nonzero elements c_{j} ∈ \overline{k} and p_{j} ∈ K(c_{j})[x_{1}] \ K(c_{j}) with 1 ≤ j ≤ n for some finite n ∈ ℕ such that δ_{i}(A_{j}) = δ_{j}(A_{i}) for all i, j with 2 ≤ i < j ≤ m, and

a_{1} = δ_{1}(f) + ∑_{j=1}^{n} c_{j} \frac{δ_{1}(p_{j})}{p_{j}},
a_{ℓ} = δ_{ℓ}(f) + ∑_{j=1}^{n} c_{j} \frac{δ_{ℓ}(p_{j})}{p_{j}} + A_{ℓ} for ℓ = 2, …, m.

Moreover, the p_{j}’s are pairwise coprime polynomials over K(c_{1}, …, c_{n}).

Proof. Lemma 2.3.9 asserts that the theorem holds for m = 1. Applying Lemma 2.3.9 to a_{1} ∈ K(x_{1}) yields that there exist f ∈ k(x), nonzero elements c_{j} ∈ \overline{K} and p_{j} ∈ K(c_{j})[x_{1}] \ K(c_{j}) with 1 ≤ j ≤ n for some n ∈ ℕ such that

a_{1} = δ_{1}(f) + ∑_{j=1}^{n} c_{j} \frac{δ_{1}(p_{j})}{p_{j}}.  (4.15)

Moreover, the c_{j}’s are the distinct roots of the Rothstein-Trager resultant of the logarithmic part A/D of a_{1} with respect to x_{1} and p_{j} = gcd(D, A − c_{j}δ_{1}(D)). So the p_{j}’s are pairwise coprime polynomials over K(c_{1}, …, c_{n}) by Lemma 2.3.7.

We show that all the c_{j}’s are in \overline{k}. For all ℓ with 2 ≤ ℓ ≤ m, the commutative formula (4.14) implies

δ_{ℓ}(a_{1}) = δ_{1}(δ_{ℓ}(f)) + δ_{1} (∑_{j=1}^{n} c_{j} \frac{δ_{ℓ}(p_{j})}{p_{j}} ) + ∑_{j=1}^{n} δ_{ℓ}(c_{j})\frac{δ_{1}(p_{j})}{p_{j}}.

Now, it follows from the integrability condition δ_{ℓ}(a_{1}) = δ_{1}(a_{ℓ}) that

δ_{1} (a_{ℓ} − δ_{ℓ}(f) − ∑_{j=1}^{n} c_{j} \frac{δ_{ℓ}(p_{j})}{p_{j}} ) = ∑_{j=1}^{n} δ_{ℓ}(c_{j})\frac{δ_{1}(p_{j})}{p_{j}}.

Since p_{j} ∈ K(c_{j})[x_{1}] \ K(c_{j}), we have δ_{1}(p_{j}) ≠ 0. By Lemma 2.3.4, for all ℓ with 2 ≤ ℓ ≤ m, we have δ_{ℓ}(c_{j}) = 0 and

a_{ℓ} = δ_{ℓ}(f) + ∑_{j=1}^{n} c_{j} \frac{δ_{ℓ}(p_{j})}{p_{j}} + A_{ℓ}, for some A_{ℓ} ∈ K(c_{1}, …, c_{n}).

So the c_{j}’s are in \overline{k} by Lemma 3.3.2 in [21]. Moreover, for all ℓ with 2 ≤ ℓ ≤ m, the sum

∑_{j=1}^{n} c_{j} \frac{δ_{ℓ}(p_{j})}{p_{j}}

is in k(x), because the c_{j}’s are distinct roots of some Rothstein-Trager resultant with coefficients in k(x_{2}, …, x_{m}) and the sum above is invariant under every permutation of the c_{j}’s. This implies the A_{ℓ}’s are in K. The integrability conditions δ_{i}(a_{j}) = δ_{j}(a_{i}) imply that δ_{i}(A_{j}) = δ_{j}(A_{i}) for all i, j with 2 ≤ i < j ≤ m. This completes the proof.

Now, we present the proof of Theorem 4.4.4.

Proof. We proceed by induction on m. Lemma 2.3.9 asserts the theorem holds for the base case m = 1. We assume that m ≥ 2 and that the theorem holds for m − 1. Let K denote k(x_{2}, x_{3}, …, x_{m}). By Lemma 4.4.5, there exist f ∈ k(x), A_{ℓ} ∈ K for ℓ = 2, …, m, nonzero elements c_{j} ∈ \overline{k} and p_{j} ∈ K(c_{j})[x_{1}] \ K(c_{j}) with 1 ≤ j ≤ n for some finite n ∈ ℕ such that δ_{i}(A_{j}) = δ_{j}(A_{i}), for all i, j with 2 ≤ i < j ≤ m, and

a_{1} = δ_{1}(f) + ∑_{j=1}^{n} c_{j} \frac{δ_{1}(p_{j})}{p_{j}},
a_{ℓ} = δ_{ℓ}(f) + ∑_{j=1}^{n} c_{j} \frac{δ_{ℓ}(p_{j})}{p_{j}} + A_{ℓ}, for ℓ = 2, …, m.

By the induction hypothesis, for the m − 1 rational functions A_{ℓ}, there exist f̄ ∈ K, nonzero elements c̄_{j} ∈ \overline{k} and p̄_{j} ∈ K(c_{j}) for j = 1, …, n̄ such that

A_{ℓ} = δ_{ℓ}(f̄) + ∑_{j=1}^{n̄} c̄_{j} \frac{δ_{ℓ}(p̄_{j})}{p̄_{j}}, for all ℓ with 2 ≤ ℓ ≤ m.

Since the f̄ and p̄_{j}’s are free of x_{1}, we have

a_{i} = δ_{i}(f + f̄) + ∑_{j=1}^{n} c_{j} \frac{δ_{i}(p_{j})}{p_{j}} + ∑_{j=1}^{n̄} c̄_{j} \frac{δ_{i}(p̄_{j})}{p̄_{j}}, for all i with 1 ≤ i ≤ m.

This completes the proof.

4.4.3 Multivariate extension of Feng-Singer-Wu’s lemma

In the bivariate continuous-discrete case, Feng, Singer, and Wu [42, Proposition 5] have shown that there exist f ∈ k(x_{1}, y_{1}), α ∈ k(y_{1}), and β, γ ∈ k(x_{1}) such that

a_{1} = \frac{δ_{1}(f)}{f} + y_{1} \frac{δ_{1}(β)}{β} + γ and b_{1} = \frac{σ_{1}(f)}{f}βα.

These two equalities are used to develop algorithms for computing Liouvillian solutions of primeorder linear difference-differential equations in [42]. We generalize Feng-Singer-Wu’s result to the multivariate case. Theorem 4.4.6. Assume that a_{1}, …, a_{m}, b_{1}, …, b_{n} are in k(x, y) with b_{1} ⋯ b_{n} ≠ 0. If all the equalities in (4.2), (4.3) and (4.4) hold, then there exist f ∈ k(x, y), β_{1}, …, β_{n} ∈ k(x), γ_{1}, …, γ_{m} ∈ k(x) and α_{1}, …, α_{n} ∈ k(y) such that

a_{i} = \frac{δ_{i}(f)}{f} + ∑_{j=1}^{n} y_{j} \frac{δ_{i}(β_{j})}{β_{j}} + γ_{i} for all i with 1 ≤ i ≤ m,  (4.16)

and

b_{j} = \frac{σ_{j}(f)}{f}β_{j}α_{j} for all j with 1 ≤ j ≤ n.  (4.17)

Moreover, all the equalities in (4.2), (4.3) and (4.4) remain valid if a_{i} is replaced by γ_{i} and b_{j} is replaced by α_{j} for all i with 1 ≤ i ≤ m and j with 1 ≤ j ≤ n.

Before presenting a proof of the theorem, we look at a few special cases. If n = 0, then we may set f = 1 and γ_{i} = a_{i} for all i with 1 ≤ i ≤ m. Similarly, if m = 0, then we set f = 1, β_{j} = 1 and b_{j} = α_{j} for all j with 1 ≤ j ≤ n. Hence, the theorem holds if either n = 0 or m = 0. If m = n = 1, then the theorem holds by Proposition 5 in [42]. Note that k is assumed to be algebraically closed in [42]. But the proposition holds without this assumption when one reads their proof carefully. Hence, it suffices to show the theorem when m > 1 and n ≥ 1 or when m ≥ 1 and n > 1.

The proof proceeds by induction on m and n. To this end, we first show that Theorem 4.4.6 holds for m = 1 and n arbitrary. Lemma 4.4.7. Let a, b_{1}, …, b_{n} be rational functions in k(x, y) such that b_{1}, …, b_{n} satisfy the integrability condition (4.3), b_{1} ⋯ b_{n} ≠ 0 and

\frac{δ(b_{j})}{b_{j}} = σ_{j}(a) − a for all j with 1 ≤ j ≤ n,  (4.18)

where δ = d/dx. Then there exist f ∈ k(x, y), α_{1}, …, α_{n} ∈ k(y), β_{1}, …, β_{n} ∈ k(x), and γ ∈ k(x) such that

a = \frac{δ(f)}{f} + ∑_{j=1}^{n} y_{j} \frac{δ(β_{j})}{β_{j}} + γ and b_{j} = \frac{σ_{j}(f)}{f}β_{j}α_{j} for 1 ≤ j ≤ n.  (4.19)

Moreover, the α_{i}’s also satisfy (4.3). Proof. We proceed by induction on n. If n = 1, then the lemma follows from Proposition 5 in [42]. Assume that n is greater than 1, and the claim holds for the values lower than n. Regard a, b_{1}, …, b_{n} as rational functions in x, y_{1}, …, y_{n−1} over k(y_{n}). By the induction hypothesis, there exist f̃ in k(x, y), β̃_{1}, …, β̃_{n−1} in k(x, y_{n}) and γ̃ in k(x, y_{n}) such that

a = \frac{δ(f̃)}{f̃} + y_{1} \frac{δ(β̃_{1})}{β̃_{1}} + ⋯ + y_{n−1} \frac{δ(β̃_{n−1})}{β̃_{n−1}} + γ̃.

Similarly, regarding a, b_{1}, …, b_{n} as rational functions in x, y_{2}, …, y_{n} over k(y_{1}), we have

a = \frac{δ(f′)}{f′} + y_{2} \frac{δ(β′_{2})}{β′_{2}} + ⋯ + y_{n} \frac{δ(β′_{n})}{β′_{n}} + γ′,

where f′ is in k(x, y), β′_{2}, …, β′_{n} are in k(x, y_{1}), and γ′ is in k(x, y_{1}). It follows from Lemma 4.3.10 that there exist f ∈ k(x, y) and β_{1}, …, β_{n}, γ ∈ k(x) such that

a = \frac{δ(f)}{f} + ∑_{j=1}^{n} y_{j} \frac{δ(β_{j})}{β_{j}} + γ.  (4.20)

Assume that, for j with 1 ≤ j ≤ n,

b_{j} = \frac{σ_{j}(f)}{f}β_{j}α_{j}  (4.21)

for some α_{j} ∈ k(x, y). By (4.20),

σ_{j}(a) − a = \frac{σ_{j}(δ(f))}{σ_{j}(f)} − \frac{δ(f)}{f} + \frac{δ_{j}(β)}{β},

and, by (4.21),

\frac{δ(b_{j})}{b_{j}} = \frac{δ(σ_{j}(f))}{σ_{j}(f)} − \frac{δ(f)}{f} + \frac{δ(β_{j})}{β_{j}} + \frac{δ(α_{j})}{α_{j}}.

These two equalities, the integrability condition (4.18), and the commutativity of δ and σ_{j} imply

\frac{δ(α_{j})}{α_{j}} = 0.

Consequently, α_{j} belongs to k(y). By (4.3), the α_{j}’s in the above proof satisfy

\frac{σ_{q}(α_{j})}{α_{j}} = \frac{σ_{j}(α_{q})}{α_{q}}

for all j, q with 1 ≤ j < q ≤ n. This completes the proof.

Now, we present a proof of Theorem 4.4.6.

Proof. We proceed by induction on m with Lemma 4.4.7 as our induction base. Let n be a fixed positive integer, and assume that the theorem holds for the values lower than m. Consider the continuous variables x_{1}, …, x_{m−1} and the discrete ones y_{1}, …, y_{n} over k(x_{m}). We have, for all j with 1 ≤ j ≤ n, there exist f̃_{j} ∈ k(x, y), β̃_{j} ∈ k(x) and α̃_{j} ∈ k(x_{m}, y) such that

b_{j} = \frac{σ_{j}(f̃_{j})}{f̃_{j}} β̃_{j}α̃_{j}.

Consider the continuous variables x_{2}, …, x_{m} and the discrete ones y_{1}, …, y_{n} over k(x_{1}). We have that for all j with 1 ≤ j ≤ n, there exist f̂_{j} ∈ k(x, y), β̂_{j} ∈ k(x) and α̂_{j} ∈ k(x_{1}, y) such that

b_{j} = \frac{σ_{j}(f̂_{j})}{f̂_{j}} β̂_{j}α̂_{j}.

By the above two equalities and Lemma 4.3.11, there exist f′_{j} ∈ k(x, y), β′_{j} ∈ k(x), α′_{j} ∈ k(y) such that

b_{j} = \frac{σ_{j}(f′_{j})}{f′_{j}}β′_{j}α′_{j} for all j with 1 ≤ j ≤ n.  (4.22)

Now, we claim that there exist f ∈ k(x, y), β_{1}, …, β_{n} ∈ k(x) and α_{1}, …, α_{n} ∈ k(y) such that

b_{j} = \frac{σ_{j}(f)}{f}β_{j}α_{j} for all j with 1 ≤ j ≤ n.

We prove the claim by induction on n. If n = 1, then the claim is true by (4.22). Assume that the claim holds for ℓ < n, and that

b_{ℓ+1} = \frac{σ_{ℓ+1}(f)}{f}u  (4.23)

for some u ∈ k(x, y). By (4.3),

σ_{1}(u) = \frac{σ_{ℓ+1}(α_{1})}{α_{1}}u, …, σ_{ℓ}(u) = \frac{σ_{ℓ+1}(α_{ℓ})}{α_{ℓ}}u.

Since α_{1}, …, α_{ℓ} ∈ k(y),

u = vw,

where v ∈ k(x, y_{ℓ+1}, …, y_{n}) and w ∈ k(y). By (4.22) and (4.23),

\frac{σ_{ℓ+1}(f)}{f}vw = \frac{σ_{ℓ+1}(f′_{ℓ+1})}{f′_{ℓ+1}}β′_{ℓ+1}α′_{ℓ+1}.  (4.24)

From the equality (4.24), it follows that

v = \frac{σ_{ℓ+1}(p)}{p}β′_{ℓ+1}q,  (4.25)

where p = f′/f and q = α′_{ℓ+1}/w. Since k is of characteristic zero, there exist c_{1}, …, c_{ℓ} ∈ k such that the denominators and numerators of p and q evaluated at y_{1} = c_{1}, …, y_{ℓ} = c_{ℓ} are all nonzero. Substituting c_{1}, …, c_{ℓ} for y_{1}, …, y_{ℓ}, respectively, into (4.25), we see that

v = \frac{σ_{ℓ+1}(g)}{g}β′_{ℓ+1}r,

where g ∈ k(x, y_{ℓ+1}, …, y_{n}) and r ∈ k(y_{ℓ+1}, …, y_{n}). Setting β_{ℓ+1} = β′_{ℓ+1} and α_{ℓ+1} = rw, we see that

b_{ℓ+1} = \frac{σ_{ℓ+1}(fg)}{fg}β_{ℓ+1}α_{ℓ+1}.

Moreover, since g is a constant with respect to σ_{1}, …, σ_{ℓ},

b_{j} = \frac{σ_{j}(fg)}{fg}β_{j}α_{j} for all j with 1 ≤ j ≤ ℓ.

The claim holds for ℓ + 1 when we replace f by fg.

It remains to obtain the form for the a_{i}’s. Assume that, for all i with 1 ≤ i ≤ m,

a_{i} = \frac{δ_{i}(f)}{f} + ∑_{j=1}^{m} y_{i} \frac{δ_{i}(β_{j})}{β_{j}} + γ_{i},

where γ_{i} ∈ k(x, y) is to be determined. Using the integrability conditions in (4.4), and performing a calculation similar to that below (4.21) in the proof of Lemma 4.4.7, we find that

σ_{j}(γ_{i}) = γ_{i} for all j with 1 ≤ j ≤ n,

which implies that γ_{i} belongs to k(x). Moreover, for all 1 ≤ i < p ≤ m,

δ_{i}(γ_{p}) = δ_{p}(γ_{i}).

By the integrability conditions (4.3)

\frac{σ_{j}(α_{q})}{α_{q}} = \frac{σ_{q}(α_{j})}{α_{j}}

for all j, q with 1 ≤ j < q ≤ n. This completes the proof of Theorem 4.4.6.

Remark 4.4.8. Applying Theorem 4.4.2 and Theorem 4.4.4 to the γ_{i}’s and α_{j}’s in Theorem 4.4.6, respectively, we readily see that the structure Theorem 4.4.1 holds.

4.5 Multiplicative structure

In this section, we assume that k is algebraically closed. We derive a multiplicative form of a hyperexponential-hypergeometric function from the structure of its certificates.

Let h(x, y) be a hyperexponential-hypergeometric function over k(x, y) with m+n certificates a_{1}, …, a_{m}, b_{1}, …, b_{n} such that the product b_{1} ⋯ b_{n} is nonzero. Then (4.2), (4.3) and (4.4) hold. By Theorem 4.4.6, there exist f ∈ k(x, y), β_{1}, …, β_{n} ∈ k(x) \ {0}, γ_{1}, …, γ_{m} ∈ k(x) and α_{1}, …, α_{n} ∈ k(y) such that

a_{i} = \frac{δ_{i}(f)}{f} + ∑_{j=1}^{n} y_{j} \frac{δ_{i}(β_{j})}{β_{j}} + γ_{i}

for all i with 1 ≤ i ≤ m, and

b_{j} = \frac{σ_{j}(f)}{f}β_{j}α_{j}

for all j with 1 ≤ j ≤ n. Moreover,

γ_{1}, …, γ_{m}, α_{1}, …, α_{n}

satisfy the integrability conditions (4.2), (4.3) and (4.4). By Theorem 2 in [22], there exists a simple differential-difference extension R of k(x, y) containing functions h, β^{y_{i}}, and any h′ ∈ H(γ_{1}, …, γ_{m}, α_{1}, …, α_{n}). Every hyperexponential-hypergeometric element in R is invertible, since R is simple. Moreover, two hyperexponential-hypergeometric elements having the same certificates differ by a multiplicative constant in k. It is straightforward to verify that h and fβ_{1}^{y_{1}} ⋯ β_{n}^{y_{n}}h′ have the same certificates. Therefore,

h = cfβ_{1}^{y_{1}} ⋯ β_{n}^{y_{n}}h′,

where c ∈ k and the certificates of h′ are γ_{1}, …, γ_{m}, α_{1}, …, α_{n}. Note that the γ_{i}’s are functions in x, and the α_{j}’s are functions in y. Then we have the factorization

H(γ_{1}, …, γ_{m}, α_{1}, …, α_{n}) = H(γ_{1}, …, γ_{m}, 1, …, 1)H(0, …, 0, α_{1}, …, α_{n}).

By Theorem 4.4.4, any element h_{e} ∈ H(γ_{1}, …, γ_{m}, 1, …, 1) is conjugate to

exp(g_{0}) ∏_{ℓ=1}^{L} g_{ℓ}^{c_{ℓ}},

where c_{1}, …, c_{L} ∈ k and g_{0}, g_{1}, …, g_{L} ∈ k(x). By Corollary 3.7.3 in [73] or Corollary 4 in [14], any element h_{g} ∈ H(0, …, 0, α_{1}, …, α_{n}) is conjugate to

f̃(y)T (y),

where f̃ ∈ k(y) and T (y) is a factorial term in y (see [73, Definition 3.5.1] or [14, Definition 5]). We summarize the above discussion by the following theorem.

Theorem 4.5.1. If k is algebraically closed, then any hyperexponential-hypergeometric function h over k(x, y) is conjugate to a function

r(x, y) exp(g_{0}(x)) ∏_{ℓ=1}^{L} g_{ℓ}(x)^{c_{ℓ}} ∏_{j=1}^{n} β_{j}(x)^{y_{j}}T (y)  (4.26)

where r ∈ k(x, y), g_{0}, g_{1}, …, g_{L}, β_{1}, …, β_{n} ∈ k(x), c_{l}, …, c_{L} ∈ k, and T (y) is a factorial term.

Example 4.5.2 (Jacobi polynomials).

P_{n}^{(m_{1},m_{2})}(x) = ∑_{ℓ} h(x, ℓ, n, m_{1}, m_{2}),

where h is a hyperexponential-hypergeometric function over k(x, ℓ, n, m_{1}, m_{2}) with certificates:

a = \frac{2ℓ + nx − n}{x^{2} − 1}, b_{1} = \frac{(x − 1)(n − ℓ)(n + m_{2} − ℓ)}{(m_{1} + ℓ + 1)(ℓ + 1)(x + 1)},
b_{2} = \frac{(x + 1)(n + m_{2} + 1)(n + m_{1} + 1)}{2(n − ℓ + 1)(n + m_{2} − ℓ + 1)}, b_{3} = \frac{n + m_{1} + 1}{m_{1} + ℓ + 1}, b_{4} = \frac{n + m_{2} + 1}{n + m_{2} − ℓ + 1}

According to Theorem 4.4.6, we write (a, b_{1}, b_{2}, b_{3}) as

a = ℓ\frac{δ(\frac{x−1}{x+1})}{\frac{x−1}{x+1}} + n\frac{δ(x + 1)}{x + 1}, b_{1} = \frac{x − 1}{x + 1} · \frac{(n − ℓ)(n + m_{2} − ℓ)}{(m_{1} + ℓ + 1)(ℓ + 1)},
b_{2} = (x + 1) · \frac{1}{2} · \frac{(n + m_{2} + 1)(n + m_{1} + 1)}{(n − ℓ + 1)(n + m_{2} − ℓ + 1)}, b_{3} = \frac{n + m_{1} + 1}{m_{1} + ℓ + 1}, b_{4} = \frac{n + m_{2} + 1}{n + m_{2} − ℓ + 1}

Then h is a hyperexponential-hypergeometric function with a multiplicative form

h = (\frac{x − 1}{x + 1})^{ℓ} (x + 1)^{n}h′, where h′ = \frac{1}{2^{n}} \binom{n + m_{1}}{n − ℓ}\binom{n + m_{2}}{ℓ}.

Definition 4.5.3 (Properness). A hyperexponential-hypergeometric function over k(x, y) is said to be proper if it is conjugate to a function

p(x, y) exp(g_{0}(x)) ∏_{ℓ=1}^{L} g_{ℓ}(x)^{c_{ℓ}} ∏_{j=1}^{n} β_{j}(x)^{y_{j}}T (y)  (4.27)

where p is a polynomial in k[x, y], g_{0}, g_{1}, …, g_{L}, β_{1}, …, β_{n} ∈ k(x), c_{1}, …, c_{L} ∈ k, and T (y) is a factorial term over k(y).

By Definition 4.5.3 and Theorem 4.4.4, we obtain the following result. Corollary 4.5.4. Any multivariate hyperexponential function H(x) over k(x) is proper.

In their paper [92], Wilf and Zeilberger have listed several familiar special functions which are proper hyperexponential-hypergeometric functions.

Existence of Telescopers for

Hyperexponential-Hypergeometric

Functions

5.1 Introduction

The termination problem of Zeilberger’s algorithms for creative telescoping has been extensively studied in the last two decades. The algorithms terminate if the existence of telescopers is guaranteed. Zeilberger [95] has shown that holonomicity guarantees the success of his algorithms. In particular, Wilf and Zeilberger have provided an elementary proof, based on the ideas of Fasenmyer [40] and Verbaeten [89], that telescopers always exist for proper hyperexponentialhypergeometric functions [95, 92]. However, holonomicity is only a sufficient condition, i.e., there are cases in which the input functions are not holonomic but Zeilberger’s algorithm still terminates [34]. Therefore, one challenge is to find theoretical criteria which enable us to algorithmically detect the existence of telescopers.

In view of the theoretical difficulty, special attention has been focused on the subclass of hyperexponential-hypergeometric functions. In the continuous case, the work by Bernstein [17] and Lipshitz [67] shows that every hyperexponential function has a telescoper. This implies that Zeilberger’s algorithm always succeeds on those inputs. However, the situation in other cases turns out to be more involved. In the discrete case and its q-analogue, the first complete solution to the termination problem is Abramov and Le’s criterion [64, 9], which decides whether telescopers exist for a given bivariate rational function in the discrete variables m and n. According to their criterion, the rational function

f = \frac{1}{m^{2} + n^{2}}

has no telescoper. Soon, the criterion was extended to the general case of bivariate hypergeometric terms by Abramov [5, 6]. Basically, Abramov proves that a hypergeometric term can be written as a sum of a hypergeometric-summable term and a proper one if it has a telescoper [6, Theorem 10]. Similar results have been obtained in the general q-shift case by Chen, Hou and Mu [28]. These results are fundamental for detecting the termination of Zeilberger’s algorithm.

The continuous-discrete analogue of Zeilberger’s algorithm is presented by Almkvist and Zeilberger [15]. This analogue has been shown to be useful in the study of orthogonal polynomials [60, Chapters 10–13].

Motivation. In the bivariate mixed case, not all hyperexponential-hypergeometric functions have telescopers. For example, we will show in this chapter that the simple rational function

f = \frac{1}{x + n}

has no telescoper with respect to either the continuous variable x or the discrete variable n. Therefore, an Abramov-like criterion is also needed in this mixed case.

Our contribution. The main contribution in this chapter is two criteria for the existence of telescopers for bivariate hyperexponential-hypergeometric functions. Our criteria are of the same type as Abramov’s, which says that a hyperexponential-hypergeometric function can be written as a sum of a hypergeometric-summable (integrable) function and a proper one if it has a telescoper with respect to the discrete (continuous) variable. The key ingredients of our proof are: standard representations of bivariate hyperexponential-hypergeometric functions in [42], and two additive decompositions: one is the Abramov-Petkovšek algorithm [11, 13] with respect to the discrete variable; and the other is the Hermite-like reduction [47] with respect to the continuous one.

The rest of this chapter is organized as follows. We describe an algebraic setting for bivariate hyperexponential-hypergeometric functions in Section 5.2. We discuss on the existence of telescopers in Section 5.3. A standard representation of hyperexponential-hypergeometric functions is introduced in Section 5.4. The Abramov-Petkovšek decomposition and the Hermite-like reduction are adapted to the continuous-discrete setting in Section 5.5. Our criteria for the existence of telescopers are presented in Section 5.6.

An earlier version of this chapter can be found in [26], which is a joint work (in progress) with Frédéric Chyzak, Ruyong Feng, and Ziming Li.

5.2 Preliminaries

In his papers [77, 78], Risch based his algorithm on Ritt’s theory of differential fields and their extensions. From then on, differential algebra [79, 56, 61] has been the algebraic foundation of symbolic integration [21]. The analog investigation [57, 58] has been carried out for summation problems via difference algebra [36]. In this section, we describe an algebraic setting for hyperexponential-hypergeometric functions, which permits us to deal with the telescoping problem in an algebraic way.

In the rest of this chapter, let k be an algebraically closed field of characteristic zero, and k(x, n) be the field of rational functions in x and n over k. We view x as the continuous variable and n as the discrete variable. On the field k(x, n), the derivation δ and the shift operator σ are defined as

δ(f(x, n)) = \frac{∂f}{∂x} and σ(f(x, n)) = f(x, n + 1)

for all f ∈ k(x, n), respectively.

5.2.1 Ring of sequences

We recall the ring of sequences defined in [41, 43] and regard a hyperexponential-hypergeometric function as an element in that ring.

Let K be the universal Picard-Vessiot field [88, Chapter 10] of the field k(x). Then (K, δ) is a differential field extension of (k(x), δ). Furthermore, the constants in K with respect to δ are exactly the elements of k, since k is algebraically closed. Let K^{ℕ} be the ring of infinite sequences over K, where the addition and multiplication of sequences are defined termwise. For a sequence

s = (s_{0}, s_{1}, s_{2}, …)

in K^{ℕ}, define

δ′(s) = (δ(s_{0}), δ(s_{1}), δ(s_{2}), …).

Then δ′ is a derivation operator on K^{ℕ}. Let I ⊂ K^{ℕ} be the ideal consisting of sequences with finitely many nonzero terms. Since I is closed under δ′, the quotient ring S = K^{ℕ}/I is a differential ring whose derivation δ̄ is induced by δ′. In other words:

δ̄(s + I) = δ′(s) + I.

We define a map σ̄ : S → S by sending

s + I ↦ (s_{1}, s_{2}, s_{3}, …) + I.

It is easy to verify that σ̄ is a well-defined automorphism, and δ̄ and σ̄ commute with each other.

We now embed the field k(x, n) of rational functions into S. For an element f ∈ k(x, n), there exists m ∈ ℕ such that f(x, i) is well-defined when i ≥ m. Define a map ϕ : k(x, n) → S by sending

f(x, n) ↦ (\underbrace{0, …, 0}_{m}, f(x, m), f(x, m + 1), f(x, m + 2), …) + I.

By the definition of I, ϕ is a well-defined monomorphism such that ϕ∘δ = δ̄∘ϕ and ϕ∘σ = σ̄ ∘ϕ. Let us identify k(x, n) with ϕ(k(x, n)). Then δ̄ and σ̄ can be viewed as extensions of δ and σ, respectively. Hence, we may also identify δ̄ with δ, and σ̄ with σ. Consequently, S is a differentialdifference ring extension of k(x, n), and k is the set of constants in S.

Lemma 5.2.1. Let a and b be in k(x, n) with b ≠ 0 and

\frac{δ(b)}{b} = σ(a) − a.  (5.1)

Then there exists an invertible hyperexponential-hypergeometric function in S such that its xcertificate is equal to a and its n-certificate is equal to b.

Proof. View a and b as rational functions in n. Let N be a sufficiently large integer such that for any j ≥ N, a(x, j) and b(x, j) are well-defined and b(x, j) ≠ 0. There exists a nonzero element v ∈ K such that

δ(v) = a(x, N)v.

Let h_{i} = 0 for all i with 0 ≤ i ≤ N − 1, h_{N} = v, and h_{i+1} = b(x, i)h_{i} for all i > N. We claim that

h = (h_{0}, h_{1}, …, h_{N}, h_{N+1}, …) + I

is an invertible element. Furthermore, its x-certificate and n-certificate are a and b, respectively.

Invertibility follows from the fact h_{i} ≠ 0 for all i ≥ N. We now verify that σ(h) = bh by the following calculation:

σ(h) = (\underbrace{0, …, 0}_{N−1}, h_{N}, h_{N+1}, h_{N+2}, ⋯) + I

= (\underbrace{0, …, 0}_{N−1}, h_{N}, b(x, N)h_{N}, b(x, N + 1)h_{N+1}, ⋯) + I

= (\underbrace{0, …, 1}_{N}, b(x, N), b(x, N + 1), ⋯)(\underbrace{0, …, 0}_{N−1}, h_{N}, h_{N}, h_{N+1}, ⋯) + I

= bh.

Next, we verify that δ(h_{i}) = a(x, i)h_{i} for i ≥ N, which implies δ(h) = ah. When i = N, the equality holds by the definition of h_{N}. Assume that the equality holds for all i with N ≤ i ≤ ℓ. By (5.1), we have

δ(b(x, ℓ)) = a(x, ℓ + 1)b(x, ℓ) − a(x, ℓ)b(x, ℓ).

It follows that

δ(h_{ℓ+1}) = δ(b(x, ℓ)h_{ℓ}) = δ(b(x, ℓ))h_{ℓ} + b(x, ℓ)δ(h_{ℓ})
= (a(x, ℓ + 1)b(x, ℓ) − a(x, ℓ)b(x, ℓ))h_{ℓ} + b(x, ℓ)a(x, ℓ)h_{ℓ}
= a(x, ℓ + 1)b(x, ℓ)h_{ℓ} = a(x, ℓ + 1)h_{ℓ+1}.

The claim is proved for ℓ + 1, and so is the lemma by induction.

Example 5.2.2. Assume that β(x) is a nonzero element of k(x). Then the pair (n\frac{δ(β)}{β}, β) satisfies (5.1). It follows that

(1, β, β^{2}, …) + I

is a hyperexponential-hypergeometric function in S with certificates nδ(β)/β and β. For later convenience, this element is denoted as β^{n}.

Let a and b be the same as in Lemma 5.2.1, and assume that (5.1) holds. Recall that

H(a, b) ≜ {h ∈ S | δ(h) = ah and σ(h) = bh}.

Proposition 5.2.3. The set H(a, b) is a one-dimensional vector space over k, in which every nonzero element is invertible in S.

Proof. By Lemma 5.2.1, H(a, b) contains an invertible element, say h. Let g be another element of H(a, b). Then δ(gh^{−1}) = 0 and σ(gh^{−1}) = gh^{−1} by a straightforward calculation. Hence, gh^{−1} is a constant. The lemma follows from the fact that k is the set of constants in S.

This proposition implies that every first-order fully integrable system over k(x, n) has a onedimensional solution space over k. For later use, we specialize the three formulas in Lemma 4.2.3 for H(a, b).

Corollary 5.2.4. (i) For any hyperexponential-hypergeometric function g ∈ S, we have

gH(a, b) = H (a + \frac{δ(g)}{g}, b\frac{σ(g)}{g});

(ii) For all a, b, a′, b′ ∈ k(x, n) with bb′ ≠ 0 such that both (a, b) and (a′, b′) satisfy (5.1), we

have

H(a, b)H(a′, b′) = H(a + a′, bb′);

(iii) δ(H(a, b)) = aH(a, b) and σ(H(a, b)) = bH(a, b).

From now on, we assume that all hyperexponential-hypergeometric functions belong to S, unless otherwise stated.

5.2.2 Ring of differential-recurrence operators

Let k(x, n)⟨D_{x}, S_{n}⟩ be the ring of differential-recurrence operators whose commutation rules are

S_{n}D_{x} = D_{x}S_{n}, D_{x}f = fD_{x} + δ(f), and S_{n}f = σ(f)S_{n},

for all f ∈ k(x, n). In this ring, we denote Δ_{n} the forward difference operator S_{n} − 1. For every s ∈ S, define two actions: D_{x}(s) = δ(s) and S_{n}(s) = σ(s). Then S becomes a left module over k(x, n)⟨D_{x}, S_{n}⟩. Assume that h is a hyperexponential-hypergeometric function with x-certificate a and n-certificate b. Then the integrability condition (5.1) can be rewritten as S_{n}D_{x}(h) = D_{x}S_{n}(h).

Two hyperexponential-hypergeometric functions h_{1} and h_{2} are said to be similar if the ratio h_{1}/h_{2} is in k(x, n). It is easy to verify that similarity is an equivalence relation among hyperexponential-hypergeometric functions.

Lemma 5.2.5. Let h_{1} and h_{2} be two hyperexponential-hypergeometric functions over k(x, n). If h_{1} and h_{2} are similar, then

(i) h_{1} + h_{2} is either equal to zero or similar to h_{1};

(ii) L(h_{1}) is either equal to zero or similar to h_{1} for any L ∈ k(x, n)⟨D_{x}, S_{n}⟩.

Proof. Let r ∈ k(x, n) be the nonzero ratio h_{1}/h_{2}. Then the first assertion follows from the equality h_{1} + h_{2} = (1 + 1/r)h_{1}. Since h_{1} is hyperexponential-hypergeometric over k(x, n), successive derivatives and shifts of h_{1} are all similar to h_{1}. By the first assertion, the second one follows.

5.2.3 Split polynomials

For a nonzero element f ∈ k(x, n), the denominator and numerator of f are denoted by den(f) and num(f), respectively. den(f) and num(f) are two coprime polynomials in k[x, n]. We introduce the definition below to describe a standard form of rational functions in k(x, n).

Definition 5.2.6. A polynomial p ∈ k[x, n] is said to be split if it is of the form p_{1}(x)p_{2}(n) with p_{1} ∈ k[x] and p_{2} ∈ k[n].

A rational function f can be decomposed as f_{1}(x)f_{2}(n)f_{3}(x, n), where f_{1} ∈ k(x), f_{2} ∈ k(n) and both num(f_{3}) and den(f_{3}) have no split factors. The product f_{1}f_{2} is called the split part of f.

Algorithms for computing rational solutions of linear differential and difference equations with polynomial coefficients have been presented by Abramov [3, 8, 4]. The lemma below concerns the shape of denominators of rational solutions of linear differential and difference equations with the leading coefficients split.

Lemma 5.2.7. Let L be an operator either in k[x, n]⟨D_{x}⟩ or in k[x, n]⟨S_{n}⟩ and let p be a polynomial in k[x, n]. If the leading coefficient of L is split, so is the denominator of rational solutions of L(y(x, n)) = p.

Proof. Let f ∈ k(x, n) be any rational solution of L(y(x, n)) = p and q ∈ k[x, n] be the leading coefficient of L. If L is in k[x, n]⟨D_{x}⟩, the splitness of den(f) follows from the well-known fact that any root of den(f) in \overline{k(n)} is also a root of q, see [3]. If L is in k[x, n]⟨S_{n}⟩, the denominator of f is a divisor of the product Q(x, n) := q(x, n)q(x, n − 1) ⋯ q(x, n − d) for some finite d ∈ ℕ according to Abramov’s algorithm for computing a denominator bound in [4]. If q is split, so is the product Q. Therefore, den(f) is split in both cases.

5.3 Existence problems

The fact that differential and difference cases can be treated in similar-looking ways has been observed for a long time in mathematics. Once a result is obtained in one side, it can often be analogized in the other one. In the sequel, our presentation will reflect the symmetry between the two counterparts of differential and difference variables.

The method of differencing under the integral sign was first formulated by Almkvist and Zeilberger in [15], which can be used to find a linear recurrence equation for the integral (if exists)

H(n) := ∫_{0}^{+∞} h(x, n) dx,

where h(x, n) is a hyperexponential-hypergeometric function over k(x, n). Here, we suppose integrals are well-defined over k, say ℂ. The key step of Almkvist and Zeilberger’s algorithm tries to find a nonzero linear recurrence operator L(n, S_{n}) in k(n)⟨S_{n}⟩ such that

L(n, S_{n})(h) = D_{x}(g),  (5.2)

for some hyperexponential-hypergeometric function g over k(x, n).

Definition 5.3.1. Let h be a hyperexponential-hypergeometric function over k(x, n). A nonzero linear recurrence operator L(n, S_{n}) ∈ k(n)⟨S_{n}⟩ is called a telescoper for h with respect to x if there exists another hyperexponential-hypergeometric function g such that (5.2) holds.

After that, applying the integral sign to both sides of (5.2) yields

L(n, S_{n})(H(n)) = g(+∞, n) − g(0, n).

This further implies that L(n, S_{n}) is indeed the recurrence relation satisfied by H(n) under certain nice boundary condition, say g(+∞, n) = g(0, n). For example, consider the integral

A(n) = ∫_{0}^{+∞} x^{n−1} exp(−x) dx.

Almkvist and Zeilberger’s algorithm returns a pair (L, g) with

L = S_{n} − n and g = −x^{n} exp(−x).

Note that g(+∞, n) = g(0, n) = 0, which implies that L(A(n)) = A(n + 1) − nA(n) = 0. So we recognize that A(n) = (n − 1)! since the value A(1) is equal to 1.

The differential analogue of Almkvist and Zeilberger’s algorithm tries to construct a linear differential equation L(x, D_{x}) ∈ k(x)⟨D_{x}⟩ such that

L(x, D_{x})(h) = Δ_{n}(g),  (5.3)

for some hyperexponential-hypergeometric function g over k(x, n).

Definition 5.3.2. Let h be a hyperexponential-hypergeometric function over k(x, n). A nonzero linear differential operator L(x, D_{x}) ∈ k(x)⟨D_{x}⟩ is called a telescoper for h with respect to n if there exists another hyperexponential-hypergeometric function g such that (5.3) holds.

Both the g’s in (5.2) and (5.3) are either zero or similar to h by Lemma 5.2.5. This analogue can be used to construct a linear differential equation for the sum

H(x) := ∑_{n=0}^{+∞} h(x, n).

For more interesting examples, see the appendix of [15] or Koepf’s book [60, Chapters 10–13].

As in the discrete case, not all hyperexponential-hypergeometric functions have telescopers with respect to x or n. That is, Almkvist and Zeilberger’s algorithm does not terminate on all hyperexponential-hypergeometric inputs.

Example 5.3.3. The non-existence phenomenon is well illustrated by the rational function

f = \frac{1}{x + n}.

Let us show that f has no telescoper with respect to either x or n.

Differential case: Suppose that f has a telescoper with respect to x, i.e., there exists a nonzero operator L in k(n)⟨S_{n}⟩ such that L(f) = D_{x}(g) for some g ∈ k(x, n). Write L = ∑_{i=0}^{ρ} ℓ_{i}(n)S_{n}^{i} with ρ ∈ ℕ and ℓ_{i} ∈ k(n). Then we have

L(f) = ∑_{i=0}^{ρ} ℓ_{i}(n)S_{n}^{i} (\frac{1}{x + n}) = ∑_{i=0}^{ρ} (\frac{ℓ_{i}(n)}{x + n + i}) = \frac{A}{D} ,

where D divides the product (x + n)(x + n + 1) ⋯ (x + n + ρ), which is squarefree with respect to x, and A ∈ k(n)[x] with deg_{x}(A) < deg_{x}(D) and gcd(A, D) = 1 as polynomials in x. Since A/D = D_{x}(g), A = 0 by Lemma 4.3.1 (i). So f is a rational solution of the linear difference equation L(z) = 0. By Lemma 5.2.7, the denominator of f should be split in k[x, n] since the leading coefficient of L is free of x, but x + n is not split, a contradiction.

Difference case: Suppose that f has a telescoper with respect to n, i.e., there exists a nonzero operator L in k(x)⟨D_{x}⟩ such that L(f) = Δ_{n}(g) for some g ∈ k(x, n). Write L = ∑_{i=0}^{ρ} ℓ_{i}(x)D_{x}^{i} with ρ ∈ ℕ and ℓ_{i} ∈ k(x). Then we have

L(f) = ∑_{i=0}^{ρ} ℓ_{i}(x)D_{x}^{i} (\frac{1}{x + n}) = ∑_{i=0}^{ρ} (\frac{(−1)^{i}ℓ_{i}(x)}{(x + n)^{i+1}}) = \frac{A}{D} ,

where D is a factor of (x + n)^{ρ+1}, which is shift-free with respect to n, and A ∈ k(x)[n] with deg_{n}(A) < deg_{n}(D) and gcd(A, D) = 1 as polynomials in n. Since A/D = Δ_{n}(g), A = 0 by Lemma 4.3.1 (ii). So f is a rational solution of the linear differential equation L(z) = 0. By Lemma 5.2.7, the denominator of f should be split in k[x, n] since the leading coefficient of L is free of n, but x + n is not split, a contradiction.

In this chapter, we solve the following problem, which is equivalent to the termination problem of Almkvist and Zeilberger’s algorithm.

Problem 5.3.4. Given a hyperexponential-hypergeometric function h(x, n) over k(x, n), decide whether h has a telescoper with respect to x or n.

5.4 Standard representations

We derive a standard representation of hyperexponential-hypergeometric functions, based on [42, Proposition 5], i.e., the bivariate case of Theorem 4.4.6.

Lemma 5.4.1. Let a and b be two nonzero rational functions in k(x, n) such that the integrability condition (5.1) holds. Then there exist f ∈ k(x, n) whose split part is trivial, α ∈ k(n), and β, γ ∈ k(x) such that

a = \frac{δ(f)}{f} + \frac{δ(β(x))}{β(x)}n + γ(x) and b = \frac{σ(f)}{f}α(n)β(x).  (5.4)

Proof. By Proposition 5 in [42], there exist f ∈ k(x, n), α ∈ k(n), and β, γ ∈ k(x) such that (5.4) holds. Write f as f_{1}f_{2}f_{3} where f_{1} ∈ k(x), f_{2} ∈ k(n), and f_{3} ∈ k(x, n) whose split part is trivial. Replacing f by f_{3}, γ by γ + δ(f_{1})/f_{1}, and α by ασ(f_{2})/f_{2} in (5.4) yields the lemma.

By Lemma 5.4.1, any hyperexponential-hypergeometric function can be written in a unique multiplicative form.

Proposition 5.4.2. Let h be a hyperexponential-hypergeometric function. Then there exist f ∈ k(x, n) with trivial split part, α ∈ k(n) with monic numerator and denominator, and β, γ ∈ k(x) such that

h ∈ f(x, n)β(x)^{n}H(γ(x), α(n)).  (5.5)

Moreover, if

h ∈ f′(x, n)β′(x)^{n}H(γ′(x), α′(n))

where f′ ∈ k(x, n) with trivial split part, α′ ∈ k(n) with monic numerator and denominator, and β′, γ′ ∈ k(x), then f/f′ ∈ k, α = α′, β = β′, and γ = γ′.

Proof. Let a and b be the x-certificate and n-certificate of h, respectively. By Lemma 5.4.1, there exist f ∈ k(x, n) with trivial split part, α ∈ k(n), and β, γ ∈ k(x) such that (5.4). By a straightforward calculation, h and a nonzero element of fβ^{n}H(γ, α) have the same certificates. Hence, (5.5) holds by Proposition 5.2.3. Assume that h belongs to f′(x, n)β′(x)^{n}H(γ′(x), α′(n)), as given in the statement of the proposition. Then the x-certificate and n-certificate of every nonzero element in f′(x, n)β′(x)^{n}H(γ′(x), α′(n)) are a and b, respectively. It follows that

\frac{δ(f)}{f} + \frac{δ(β(x))}{β(x)}n + γ(x) = \frac{δ(f′)}{f′} + \frac{δ(β′(x))}{β′(x)}n + γ′(x),  (5.6)

and

\frac{σ(f)}{f}α(n)β(x) = \frac{σ(f′)}{f′}α′(n)β′(x).  (5.7)

Consequently,

\frac{δ(f)}{f} = \frac{δ(f′)}{f′}, \frac{δ(β(x))}{β(x)} = \frac{δ(β′(x))}{β′(x)}, and γ(x) = γ′(x)

by (5.6), the uniqueness of the partial fraction decomposition of a, and the assumption that both f and f′ have trivial split parts. Hence, f = cf′ and β = c′β′ for some c, c′ ∈ k(n). Again, by the assumption on f and f′, c belongs to k. Since β and β′ are free of n, so is c′. It follows from (5.7) that α(n) = c′α(n). Thus, c′ is equal to 1 since the numerators and denominators of α(n) and α′(n) are assumed to be monic.

Definition 5.4.3. Let h(x, n) ∈ H(a, b) be a hyperexponential-hypergeometric function. We say that a quadruple (f(x, n), α(n), β(x), γ(x)) ∈ k(x, n)^{4} is the standard representation of h if f ∈ k(x, n) with trivial split part, α ∈ k(n) with monic numerator and denominator, β, γ ∈ k(x), and (5.5) holds.

By the definition in [92] or Definition 4.5.3 in Chapter 4, a bivariate hyperexponentialhypergeometric function is proper if its standard representation is

(p(x, n), α(n), β(x), γ(x)),  (5.8)

where p ∈ k[x, n], α ∈ k(n), and β, γ ∈ k(x).

Proposition 5.4.4. Let h ∈ f(x, n)β(x)^{n}H(γ(x), α(n)) with f ∈ k(x, n), α ∈ k(n), and β, γ ∈ k(x). If the denominator of f is split, then h is proper.

Proof. Write f = p/q with p, q ∈ k[x, n] and gcd(p, q) = 1. Since q is split, q = A(x)B(n) for some A ∈ k[x] and B ∈ k[n]. By Corollary 5.2.4 (i), we have

h ∈ p(x, n)β(x)^{n}H (γ(x) − \frac{δ(A(x))}{A(x)}, α(n)\frac{B(n)}{σ(B(n))}) .

Therefore, h is proper by definition.

5.5 Two additive decompositions

For a bivariate hyperexponential-hypergeometric function h(x, n) in the continuous variable x and the discrete one n, we can perform two kinds of complete additive decompositions: one is with respect to the continuous variable, obtaining h = D_{x}(h_{1}) + h_{2}, and the other is with respect to the discrete variable, obtaining h = Δ_{n}(h_{1}) + h_{2}. The h_{1}’s and h_{2}’s above are hyperexponentialhypergeometric functions with h_{2} minimal in some sense. Algorithms for complete additive decompositions generalize the capability of Gosper’s algorithm [51] and its differential analogue [15]. That is, it returns h_{2} = 0 if h = D_{x}(g) or h = Δ_{n}(g) for some g hyperexponential-hypergeometric. In this section, we adapt the two additive decompositions to hyperexponential-hypergeometric functions described by their standard representations. Notation. In what follows, let F be either k(n) if the decomposition is with respect to x or k(x) if it is with respect to n.

5.5.1 Additive decomposition with respect to x

In Chapter 3, Hermite reduction or Ostrogradsky and Horowitz’s method decomposes a rational function f ∈ F (x) into f = D_{x}(g) + r, where r = a/b is minimal in the sense that deg_{x}(a) < deg_{x}(b) and b is squarefree, that is, the factors of b have minimal multiplicity. Hermite reduction was later extended by Davenport [39] and by Geddes, Le and Li [47] with two different motivations: one is to solve the Risch differential equation

D_{x}(y) + fy = g, where f and g are in F (x),

and the other is to solve the decomposition problem of hyperexponential functions. We will use the latter to compute additive decompositions of hyperexponential-hypergeometric functions with respect to x. To this end, we first recall some terminology from [47]. Definition 5.5.1 (GLL triple). Let u, v ∈ F [x] and w ∈ F (x). We call (u, v, w) a Geddes–Le–Li triple, in short GLL triple, if the following conditions are satisfied:

(i) gcd(u, v) = 1; (ii) u is squarefree with respect to x; (iii) w is differential-reduced with respect to x; (iv) gcd(u, den(w)) = 1. Definition 5.5.2 (Hyperexponential-integrable). A hyperexponential function h over F (x) is said to be hyperexponential-integrable if there exists a hyperexponential function g such that h = D_{x}(g).

Note that the function g in the above definition is similar to h if it exists. Definition 5.5.3 (Additive decomposition with respect to x). Let h be in H(a, b) with a, b ∈ F (x). An additive decomposition of h with respect to x is of the form

h = D_{x}(h_{1}) + h_{2},

where h_{1} is hyperexponential-hypergeometric over F (x), and either h_{2} is zero or

h_{2} ∈ \frac{v}{u} · H(a_{2}, b_{2}), where b_{2} ∈ F (x) and (u, v, a_{2}) is a GLL triple.

Such a decomposition is said to be complete if h_{2} = 0 whenever h is hyperexponential-integrable. Remark 5.5.4. In the sequel we never use complete additive decompositions. In fact, the h_{1} and h_{2} above are similar to h by the same argument as in [75, Proposition 5.6.2].

Geddes–Le–Li’s algorithm

We recall the algorithm from [47] to perform additive decompositions with respect to x. View a hyperexponential-hypergeometric function h ∈ H(a, b) as a hyperexponential function with xcertificate a ∈ F (x). Let (K, S) ∈ F (x) × F (x) be a strict DRNF of a, that is, K is differentialreduced with respect to x and

a = K + \frac{δ(S)}{S}, where gcd(den(K), den(S)) = 1.

Immediately, the additive form above for a leads to a multiplicative form for h, that is,

h = S · h̃, where h̃ ∈ H (K, b\frac{S}{σ(S)}) .

Now comes the Hermite-like reduction: Theorem 2 in [47] shows that there exists S_{1} ∈ F (x) such that

S = δ(S_{1}) + KS_{1} + \frac{v}{u · den(K)^{i}} ,  (5.9)

where i ∈ {0, 1} and (u, v, K) is a GLL triple. Since D_{x}(h̃) = Kh̃, the decomposition for S in (5.9) further implies that

h = D_{x}(S_{1}h̃) + \frac{v}{u · den(K)^{i}} h̃,

In order to get a complete additive decomposition for h, we still need to decide whether

h_{2} = \frac{v}{u · den(K)^{i}} h̃  (5.10)

is hyperexponential-integrable or not. A necessary condition is as follows.

Lemma 5.5.5 (Geddes–Le–Li, 2004). Let h be a hyperexponential function over F (x). If h is hyperexponential-integrable, and

\frac{δ(h)}{h} = w + \frac{δ(v/u)}{v/u},

where v, u ∈ F [x], w ∈ F (x), and (u, v, w) is a GLL triple, then u is in F.

Proof. For the proof, see [47, Theorem 4].

Additive decompositions with respect to x from standard representations

Let (f(x, n), α(n), β(x), γ(x)) be the standard representation of h ∈ H(a, b), that is,

h ∈ f(x, n)β(x)^{n}H(γ(x), α(n)),

where f ∈ k(x, n) with trivial split part, α ∈ k(n) with monic numerator and denominator, and β, γ ∈ k(x). We show how to obtain an additive decomposition of h with respect to x from its standard representation.

First, we compute a strict DRNF for a from the standard representation of h. The following lemma helps us decide whether certain rational functions are differential-reduced.

Lemma 5.5.6. Let

g = \frac{p}{q} + n\frac{s}{t} + \frac{u}{v},

where p, q, s, t, u, v ∈ k[x] with gcd(p, q) = gcd(s, t) = gcd(u, v) = 1. Assume that v | t. Then

(i) the denominator of g is the least common multiple of q and t;

(ii) if p/q is differential-reduced, so is g;

(iii) if p/q + ns/t is differential-reduced, so is g.

Proof. Let M = lcm(q, t), w_{1} = M/q and w_{2} = M/t. We first show the following claim.

Claim. If w is a nontrivial irreducible factor of

gcd(M, sw_{2}n + r), where r is an element in k[x],

then w divides w_{2}, r and q, but does not divide w_{1}.

Proof of the Claim. Since n is an indeterminate over k, w is a factor of sw_{2}. Suppose that w is not a factor of w_{2}. Then it is a factor of s. Since w is a factor of M, it is a factor of w_{2}t, and, thus, a factor of t, a contradiction to the assumption that gcd(s, t) = 1. We have concluded that w is a factor of w_{2}. Consequently, it is a factor of r because n is an indeterminate, but it is not a factor of w_{1} because gcd(w_{1}, w_{2}) = 1. Since w divides qw_{1}, it divides q. The claim is proved.

To prove the first assertion, we let w_{3} = M/v, which is a multiple of w_{2} since t is a multiple of v. Put W = pw_{1} + sw_{2}n + uw_{3}. We have g = W/M. It remains to show that gcd(M, W) = 1. Set r = pw_{1} + uw_{3} and assume that w is an irreducible factor of gcd(M, W). By the claim, w divides r and w_{2}, hence, it divides pw_{1}. The claim then implies that w divides p. But w also divides q by the same claim. Thus, w ∈ k since gcd(p, q) = 1.

To prove the second assertion, we only need to show that gcd(q, p − ℓδ(q)) = 1 implies that gcd(M, W − ℓδ(M)) = 1. Setting r = −pw_{1} − uw_{2} − ℓδ(M), we have

gcd(M, W − ℓδ(M)) = gcd(M, sw_{2}n + r).

Assume that w is an irreducible factor of gcd(M, W −ℓδ(M)). By the claim, w divides r and w_{2}, therefore, it divides pw_{1} − ℓδ(M), which is equal to pw_{1} − ℓw_{1}δ(q) − ℓδ(w_{1})q. Since w divides q, it divides (p−ℓδ(q))w_{1}. Again, by the claim, w divides p−ℓδ(q). On the other hand, w divides q because w | M and w ∤ w_{1}. Hence, w is a factor of gcd(q, p − ℓδ(q)), which is equal to 1, because p/q is differential-reduced. Consequently, w is in k.

To prove the third assertion, it suffices to show that gcd(M, pw_{1} + sw_{2}n − ℓδ(M)) = 1 implies gcd(M, W − ℓδ(M)) = 1. Let r and w be defined as above. By the claim, w divides w_{3} since w_{2} divides w_{3}. It follows from the definition of r that w divides pw_{1} + sw_{2}n − ℓδ(M), and, hence, w is in k.

For later use, we recall the definition of pumps, given in [13, page 533].

Definition 5.5.7 (Pump). Let p, q ∈ F [x] be such that p | q. A factor p̃ of q is said to be the pump of p in q if gcd(q/p̃, p̃) = 1, p | p̃, and any irreducible factor of p̃ divides p.

Example 5.5.8. Let p_{1}, p_{2}, p_{3} be in F [x] with positive degrees. Assume that they are pairwise coprime and p = p_{1}p_{2}^{2}p_{3}^{3}. Then the pump of p_{1}p_{3}^{2} in p is equal to p_{1}p_{3}^{3}.

In other words, a pump of a factor p of q has the same squarefree part as p’s, but each of its irreducible factors has the maximal multiplicity in q. For a rational function f ∈ k(x, n), denote by f^{∗} the squarefree part of num(f) · den(f). Now, we show how to compute a strict DRNF of a from that of γ.

Lemma 5.5.9. Let (f(x, n), α(n), β(x), γ(x)) be a standard representation of a nonzero element in H(a, b). Let (K̃, S̃) be a strict DRNF of γ. Let u = gcd(den(S̃), β^{∗}) and T the pump of u in den(S̃). If

K = K̃ + n\frac{δ(β)}{β} − \frac{δ(T)}{T} and S = fS̃T,  (5.11)

then (K, S) is a strict DRNF of a with respect to x.

Proof. Let (K̃, S̃) be a strict DRNF of γ(x). Then by definition, we have

a = \frac{δ(f)}{f} + n\frac{δ(β)}{β} + γ = \frac{δ(f)}{f} + n\frac{δ(β)}{β} + \frac{δ(S̃)}{S̃} + K̃.

By Lemma 2.3.3, the denominators of δ(β)/β and δ(T)/T are β^{∗} and T^{∗}, respectively. By the definition of pump and u | β^{∗}, u is the squarefree part of T, which implies T^{∗} | β^{∗}. Hence, K is differential-reduced by Lemma 5.5.6 (ii).

It remains to show gcd(den(K), den(S)) = 1. By the definition of pumps, den(S̃) = T w, where w ∈ k[x] and gcd(T, w) = 1. On the other hand,

den(S) = den(f)den(S̃T) = den(f)w,

since num(f)den(f) has only trivial factors in k[x]. It follows that

gcd(den(K), den(S)) = gcd(den(K), den(f)w),

which is equal to gcd(den(K), w) since den(K) ∈ k[x]. So it suffices to show that den(K) and w are coprime. By Lemma 5.5.6 (i), den(K) = lcm(den(K̃), β^{∗}). So it suffices to show that

gcd(den(K̃), w) = 1 and gcd(β^{∗}, w) = 1.

The first equality is immediate from the assumptions that (K̃, S̃) is a strict DRNF and that w is a factor of den(S̃). To prove the second, we assume that p is an irreducible factor of gcd(β^{∗}, w). Then p is a factor of u by the definitions of u and w. Hence, it is a factor of T. Consequently, p is a common factor of T and w. We conclude that p is in k.

We apply Geddes–Le–Li’s algorithm to a hyperexponential-hypergeometric function described by its standard representation and obtain a special form of its additive decomposition with respect to x.

Theorem 5.5.10. Let (f(x, n), α(n), β(x), γ(x)) be the standard representation of a nonzero element h in H(a, b). Then one can find a hyperexponential-hypergeometric function h_{1} and γ̃(x) ∈ k(x) such that

h − D_{x}(h_{1}) ∈ \frac{v}{u} · H (γ̃ + n\frac{δ(β)}{β}, βα)  (5.12)

where u, v are in k(n)[x] and (u, v, γ̃ + nδ(β)/β) is a GLL triple.

Proof. By Lemma 5.5.9, the x-certificate a of h has a strict DRNF (K, S) of the form

(K̃ + n · \frac{δ(β)}{β} − \frac{δ(T)}{T}, f · S̃ · T ) ,

where (K̃, S̃) is a strict DRNF of γ(x) and T ∈ k[x]. A direct calculation verifies that

h = Sh̃ for some h̃ ∈ H (K, βα).

By (5.9), Hermite-like reduction rewrites S into

S = δ(S_{1}) + S_{1}K + \frac{v}{uden(K)^{i}} ,

where i ∈ {0, 1}, S_{1} ∈ k(x, n) and (u, v, K) is a GLL triple. It follows that

h = D_{x}(S_{1}h̃) + \frac{v}{uden(K)^{i}} h̃.

Since den(K) is in k[x] and by Corollary 5.2.4

\frac{v}{uden(K)^{i}} H(K, βα) = \frac{v}{u}H (K − \frac{δ(den(K)^{i})}{den(K)^{i}}, βα\frac{den(K)^{i}}{σ(den(K)^{i})})

= \frac{v}{u}H (K − \frac{δ(den(K)^{i})}{den(K)^{i}}, βα) .

Setting h_{1} = S_{1}h̃ and γ̃ = K̃ − \frac{δ(T)}{T} − \frac{δ(den(K)^{i}}{den(K)^{i}} yields (5.12). It remains to show that the triple (u, v, γ̃ + n\frac{δ(β)}{β}) is a GLL triple. Note that

g := γ̃ + n\frac{δ(β)}{β} = K − \frac{δ(den(K)^{i})}{den(K)^{i}}.

Since K is differential-reduced, so is g by definition. Moreover, u and den(g) are coprime with respect to x, as u and den(K) are coprime. It follows that (u, v, g) is a GLL triple.

Remark 5.5.11. An additive decomposition with respect to x can be computed by the algorithm

ReducedHyperexp in [47].

5.5.2 Additive decomposition with respect to n

In the discrete case, the algorithms in [2, 4, 76, 72] decompose a rational function f ∈ F (n) into f = Δ_{n}(g) + r, where r = a/b is minimal in the sense that deg_{n}(a) < deg_{n}(b) and b is shift-free, that is b has minimal dispersion. Those decomposition algorithms are extended to the case of hypergeometric terms by Abramov and Petkovšek [11, 13]. This extension also generalizes the capability of the well-known Gosper algorithm [51], which decides whether the indefinite sum of a hypergeometric term is hypergeometric. More precisely, for a given hypergeometric term H(n) over F (n), Abramov and Petkovšek’s algorithm computes two hypergeometric terms H_{1}(n) and H_{2}(n) such that H = Δ_{n}(H_{1}) + H_{2}, where H_{2} is minimal in some sense. In particular, H_{2} is identically zero if H(n) has a hypergeometric indefinite sum. Here, we will use Abramov and Petkovšek’s algorithm to compute additive decompositions of hyperexponentialhypergeometric functions with respect to the discrete variable n. To this end, we recall some terminology from [13].

Definition 5.5.12 (AP triple). Let u, v ∈ F [n] and w ∈ F (n). We call (u, v, w) an Abramov– Petkovšek triple, in short AP triple, if the following conditions are satisfied:

(i) gcd(u, v) = 1;

(ii) u is shift-free with respect to n;

(iii) w is shift-reduced with respect to n;

(iv) gcd(u, σ^{−ℓ}(num(w))) = gcd(u, σ^{ℓ}(den(w))) = 1 for all ℓ ∈ ℕ.

Definition 5.5.13 (Hypergeometric-summable). A hypergeometric term h over F (n) is said to be hypergeometric-summable if there exists a hypergeometric function g such that h = Δ_{n}(g).

Note that the g in the above definition is similar to h if it exists.

Definition 5.5.14 (Additive decomposition with respect to n). Let h be in H(a, b) with a, b ∈ F (n). An additive decomposition of h with respect to n is of the form

h = Δ_{n}(h_{1}) + h_{2},

where h_{1} is hyperexponential-hypergeometric over F (n), and either h_{2} is zero or

h_{2} ∈ \frac{v}{u} · H(a_{2}, b_{2}), where a_{2} ∈ F (n) and (u, v, b_{2}) is an AP triple.

Such a decomposition is said to be complete if h_{2} = 0 whenever h is hypergeometric-summable.

Remark 5.5.15. In the sequel we never use complete additive decompositions. In fact, the h_{1} and h_{2} above are similar to h by the same argument in [75, Proposition 5.6.2].

Abramov and Petkovšek’s algorithm

We recall the algorithm from [11, 13] to perform additive decomposition with respect to n. View a hyperexponential-hypergeometric function h ∈ H(a, b) as a hypergeometric term with ncertificate b ∈ F (n). Let (K, S) ∈ F (n) × F (n) be a strict SRNF of b, that is, K is shift-reduced with respect to n and

b = K · \frac{σ(S)}{S},

where gcd(num(K), σ(den(S)) · num(S)) = gcd(den(K), σ(num(S)) · den(S)) = 1. Immediately, the multiplicative form above for b leads to a multiplicative form for h, that is,

h = S · h̃, where h̃ ∈ H (a − \frac{δ(S)}{S}, K) .

Now comes Abramov and Petkovšek’s reduction: Lemma 9 in [13] implies that there exists S_{1} ∈ F (n) such that

S = σ(S_{1})K − S_{1} + \frac{v}{u · (σ^{−1}(k_{1}))^{i} · k_{2}^{j}} ,  (5.13)

where i, j ∈ {0, 1}, k_{1} and k_{2} are the numerator and denominator of K, respectively, and (u, v, K) is an AP triple. Since S_{n}(h̃) = Kh̃, the decomposition for S in (5.13) further implies that

h = Δ_{n}(S_{1}h̃) + \frac{v}{u · (σ^{−1}(k_{1}))^{i} · k_{2}^{j}} h̃.

In order to get a complete additive decomposition of h, we still need to decide whether

h_{2} = \frac{v}{u · (σ^{−1}(k_{1}))^{i} · k_{2}^{j}} h̃  (5.14)

is hypergeometric-summable or not. A necessary condition is as follows.

Lemma 5.5.16 (Abramov–Petkovšek, 2002). Let h be a hypergeometric function over F (n). If h is hypergeometric-summable, and

\frac{σ(h)}{h} = \frac{σ(v/u)}{v/u}w,

where v, u ∈ F [n], w ∈ F (n), and (u, v, w) is an AP triple, then u is in F.

Proof. For the proof, see [13, Theorem 11].

Additive decompositions with respect to n from standard representations

Let (f(x, n), α(n), β(x), γ(x)) be the standard representation of h ∈ H(a, b). We show how to obtain an additive decomposition of h with respect to n from its standard representation.

The next lemma relates a strict SRNF of the n-certificate b of h to its standard representation.

Lemma 5.5.17. Let (f(x, n), α(n), β(x), γ(x)) be the standard representation of a nonzero element in h ∈ H(a, b). If (K̃, S̃) is a strict SRNF of α(n) with respect to n, then (βK̃, fS̃) is a strict SRNF of the n-certificate a of h.

Proof. Let K = βK̃ and S = fS̃. Since (K̃, S̃) is a strict SRNF of α(n), α = K̃σ(S̃)/S̃. It follows that

b = \frac{σ(h)}{h} = βα\frac{σ(f)}{f} = β (K̃ \frac{σ(S̃)}{S̃}) \frac{σ(f)}{f} = βK̃ · \frac{σ(fS̃)}{fS̃}.

The product βK̃ is shift-reduced with respect to n, because K̃ is shift-reduced and β is in k(x). Write K = k_{1}/k_{2} and S = s_{1}/s_{2} with gcd(k_{1}, k_{2}) = gcd(s_{1}, s_{2}) = 1 in k[x, n]. It remains to verify the GCD condition

gcd(k_{1}, σ(s_{2})s_{1}) = gcd(k_{2}, σ(s_{1})s_{2}) = 1  (5.15)

for (K, S). Since β ∈ k(x) and K̃ ∈ k(n), we have

k_{1} = num(β) · num(K̃) and k_{2} = den(β) · den(K̃).

Similarly, since S̃ ∈ k(n) and f has a trivial split part, we have

s_{1} = num(f) · num(S̃), and s_{2} = den(f) · den(S̃).

Hence, the condition (5.15) holds, because (K̃, S̃) is a strict SRNF, and every nonzero polynomial in k[n] is coprime with both num(f) and den(f).

We apply Abramov–Petkovšek’s algorithm to a hyperexponential-hypergeometric function described by standard representation and obtain a special form of its additive decomposition with respect to n.

Theorem 5.5.18. Let (f(x, n), α(n), β(x), γ(x)) be the standard representation of a nonzero element in h ∈ H(a, b). Then one can find γ̃(x) ∈ k(x), α̃ ∈ k(n), and a hyperexponentialhypergeometric function h_{1} such that

h − Δ_{n}(h_{1}) ∈ \frac{v}{u} · H (γ̃ + n\frac{δ(β)}{β}, βα̃)

where u, v ∈ k(x)[n] and (u, v, βα̃) is an AP triple.

Proof. By Lemma 5.5.17, the n-certificate b of h has a strict SRNF (K, S) of the form

(βK̃, fS̃) ,

where (K̃, S̃) is a strict SRNF of α(n). A direct calculation verifies that

h = Sh̃ for some h̃ ∈ H (γ + n\frac{δ(β)}{β}, K).

Write K = k_{1}/k_{2} and S = s_{1}/s_{2} with gcd(k_{1}, k_{2}) = gcd(s_{1}, s_{2}) = 1 in k[x, n]. Now, Abramov– Petkovšek’s reduction [13] decomposes S into

S = δ(S_{1})K − S_{1} + \frac{v}{u(σ^{−1}(k_{1}))^{i}k_{2}^{j}} ,

where i, j ∈ {0, 1}, S_{1} ∈ k(x, n) and (u, v, K) is an AP triple. It follows that

h = Δ_{n}(S_{1}h̃) + \frac{v}{u(σ^{−1}(k_{1}))^{i}k_{2}^{j}} h̃.

Set T = (σ^{−1}(k_{1}))^{i}k_{2}^{j} ∈ k[x, n]. Since K is split in k[x, n], T is split in k[x, n]. By Corollary 5.2.4, we have

\frac{v}{u · T} h̃ ∈ \frac{v}{u · T} H (γ + n\frac{δ(β)}{β}, K) = \frac{v}{u}H (γ + n\frac{δ(β)}{β} − \frac{δ(T)}{T}, K \frac{T}{σ(T)}) .

Set h_{1} = S_{1}h̃, γ̃ = γ − \frac{δ(T)}{T}, and α̃ = K̃ \frac{T}{σ(T)}. Since T is split, γ ∈ k(x), and K̃ ∈ k(n), we have that γ̃ ∈ k(x) and α̃ ∈ k(n) hold. It remains to show that (u, v, βα̃) is an AP triple.

We only need to verify the last two conditions in Definition 5.5.12. Since i, j ∈ {0, 1}, βα̃ is in the set

\{\frac{k_{1}}{k_{2}} , \frac{σ^{−1}(k_{1})}{k_{2}}, \frac{k_{1}}{σ(k_{2})}, \frac{σ^{−1}(k_{1})}{σ(k_{2})}\} .

Then βα̃ is shift-reduced, because K = k_{1}/k_{2} is. Note that the numerator and denominator of βα̃ are just some shifts of that of K, respectively. Then the last gcd condition follows from the assumptions on u and k_{1}, k_{2}.

Remark 5.5.19. An additive decomposition with respect to n can be computed by the algorithm dterm in [13].

5.5.3 Applying operators to additive decompositions

In this section, we apply an operator in k(x)⟨D_{x}⟩, resp. in k(n)⟨S_{n}⟩, to an additive decomposition with respect to n, resp. x, and prove that the application results in another additive decomposition. To this end, we need a simple lemma.

Lemma 5.5.20. Let u and v be in k(x)[n] and let w be any factor of u^{m} with m ∈ ℕ. Then

(i) if u is shift-free with respect to n, so is w;

(ii) for all i ∈ ℤ, if gcd(u, σ^{i}(v)) = 1, then gcd(w, σ^{i}(v)) = 1.

Proof. A polynomial in k(x)[n] is shift-free with respect to n if and only if the largest integer distance of its roots is equal to zero [2]. This fact implies the first assertion. For any i ∈ ℤ, if gcd(u, σ^{i}(v)) = 1, then gcd(u^{m}, σ^{i}(v)) = 1 for all m ∈ ℕ, which implies gcd(w, σ^{i}(v)) = 1. The second assertion holds.

Proposition 5.5.21. With the assumptions and notation introduced in Theorem 5.5.10, assume further that h_{2} = h − D_{x}h_{1}, and that L belongs to k(n)⟨S_{n}⟩ such that L(h_{2}) is nonzero. Set

A_{2} = γ̃ + n\frac{δ(β)}{β} − ρ\frac{δ(den(β))}{den(β)} and B_{2} = βα,

where ρ denotes the order of L. Then

L(h_{2}) ∈ \frac{V}{U} H(A_{2}, B_{2})

for some U, V ∈ k(n)[x], and (U, V, A_{2}) is a GLL triple.

Proof. Setting a_{2} = γ̃ + nδ(β)/β. we get A_{2} = a_{2} − ρδ(den(β))/den(β).

By Theorem 5.5.10, h_{2} ∈ \frac{v}{u} · H (a_{2}, B_{2}) . An easy calculation implies

L(h_{2}) = \frac{v′}{u′den(β)^{ρ}} · H (a_{2}, B_{2}) ,

where v′ ∈ k(n)[x] and u′ is the least common multiple of u, S_{n}(u), …, S_{n}^{ρ}(u) in k(n)[x]. Since u is squarefree with respect to x, so are its shifts. Consequently, u′ is squarefree with respect to x. Set

U = den(v′/u′) and V = num(v′/u′).

Then L(h_{2}) = (V/U) · H(A_{2}, B_{2}).

It remains to show that (U, V, A_{2}) is a GLL triple. First, U and V are coprime by their definition. Second, U is squarefree with respect to x, as u′ is squarefree with respect to x. To verify the third condition in Definition 5.5.1, we observe that the squarefree part of den(β)num(β) is the denominator of δ(β)/β, and that den(β) is the denominator of ρδ(den(β))/den(β), which divides the former denominator. This observation enables us to apply Lemma 5.5.6. We conclude that A_{2} is differential-reduced by Lemma 5.5.6 (iii) and since γ̃ + n\frac{δ(β)}{β} is differential-reduced. At last, we verify that U and den(A_{2}) are coprime over k(n). By Lemma 5.5.6 (i), den(a_{2}) = den(A_{2}). Since (u, v, a_{2}) is a GLL triple, gcd(u, den(a_{2})) = 1. It follows from the fact den(a_{2}) ∈ k[x] that

gcd (S_{n}^{i}(u), den(a_{2})) = 1

for all i ∈ ℕ, and, consequently, gcd(u′, den(a_{2})) = 1. This implies that gcd(U, den(a_{2})) = 1, that is, gcd(U, den(A_{2})) = 1.

Proposition 5.5.22. With the assumptions and notation introduced in Theorem 5.5.18, assume further that h_{2} = h − Δ_{n}(h_{1}), and that L belongs to k(x)⟨D_{x}⟩ such that L(h_{2}) is nonzero. Set

A_{2} = γ̃ + n\frac{δ(β)}{β} and B_{2} = βα̃.

Then

L(x, D_{x})(h_{2}) ∈ \frac{V}{U} · H (A_{2}, B_{2})

for some U, V ∈ k(x)[n], and (U, V, B_{2}) is an AP triple.

Proof. By Theorem 5.5.18,

h_{2} ∈ \frac{v}{u}H(A_{2}, B_{2})

for some u, v ∈ k(x)[n] and (u, v, B_{2}) is an AP triple. A straightforward calculation leads to

L(x, D_{x})(h_{2}) ∈ (V/U) · H(A_{2}, B_{2})

for some U, V ∈ k(x)[n] with gcd(U, V) = 1. By Theorem 5.5.18, den(A_{2}) is in k[x]. It follows that U is a factor of a power of u in k(x)[n]. Thus, (U, V, B_{2}) is an AP triple by Lemma 5.5.20 and since (u, v, B_{2}) is an AP triple, as stated in Theorem 5.5.18.

5.6 Two existence criteria

The Fundamental Theorem in [92] shows that telescopers exist for proper hyperexponentialhypergeometric functions. However, properness is just a sufficient condition. For instance, the rational function f = 1/(x + n)^{2} is not proper, but it does have a telescoper with respect to x since f = D_{x}(−1/(x + n)). On the other hand, Example 5.3.3 indicates that not all hyperexponential-hypergeometric functions have telescopers. In this section, we present sufficient and necessary conditions in the following theorems for the existence of telescopers.

Theorem 5.6.1. Let h be a hyperexponential-hypergeometric function over k(x, n) and h = D_{x}(h_{1}) + h_{2} be an additive decomposition of h with respect to x. Then h has a telescoper with respect to x if and only if h_{2} is either zero or proper.

Theorem 5.6.2. Let h be a hyperexponential-hypergeometric function over k(x, n) and h = Δ_{n}(h_{1}) + h_{2} be an additive decomposition of h with respect to n. Then h has a telescoper with respect to n if and only if h_{2} is either zero or proper.

Our proofs of these two theorems will be divided into several steps.

5.6.1 Existence of telescopers implies properness

The lemma below shows that the existence problem for h is equivalent to that for h_{2}. The proofs take the same argument as in [9, page 3].

Lemma 5.6.3. Let h, h_{1} and h_{2} be three hyperexponential-hypergeometric functions in the ring S of sequences.

(i) Assume that h = D_{x}(h_{1}) + h_{2}. Then L is a telescoper for h with respect to x if and only if it is a telescoper for h_{2} with respect to x.

(ii) Assume that h = Δ_{n}(h_{1}) + h_{2}. Then L is a telescoper for h with respect to n if and only if it is a telescoper for h_{2} with respect to n.

Proof. By the commutativity between L(n, S_{n}) ∈ k(n)⟨S_{n}⟩ and D_{x}, we have

L(n, S_{n})(h) = D_{x}(g) ⇔ L(n, S_{n})(h_{2}) = D_{x}(g − L(h_{1})).

This implies the first assertion. Similarly, the second one follows from the commutativity between L(x, D_{x}) ∈ k(x)⟨D_{x}⟩ and Δ_{n}.

Proof of the necessity in Theorem 5.6.1

Proof. Assume that h has a telescoper with respect to x, and that h_{2} is nonzero. Our goal is to show that h_{2} is proper. By Theorem 5.5.10,

h_{2} = \frac{v}{u}H for some H ∈ H (a_{2}, b_{2}) ,

where b_{2} = βα with β ∈ k(x) and α ∈ k(n), and (u, v, a_{2}) is a GLL triple. Without loss of generality, we further assume that u is in k[x, n]. To prove that h_{2} is proper, it suffices to prove that u is split by Proposition 5.4.4. Let L ∈ k(n)⟨S_{n}⟩ be a telescoper for h with respect to x. By Lemma 5.6.3, L is also a telescoper for h_{2} with respect to x. Then

L(h_{2}) = L (\frac{v}{u}H) = D_{x}(g),  (5.16)

for some hyperexponential-hypergeometric function g. Since b_{2} is split, a direct calculation yields

L(h_{2}) = M (\frac{v}{u}) H  (5.17)

for some nonzero element M ∈ k(n, x)⟨S_{n}⟩ whose coefficients are split. On the other hand, Proposition 5.5.21 implies

L(h_{2}) ∈ \frac{V}{U} H(A_{2}, B_{2}).  (5.18)

where U, V ∈ k(n)[x] and (U, V, A_{2}) is a GLL triple. Since L(h_{2}) is hyperexponential-integrable by (5.16), it follows from Lemma 5.5.5 that U is free of x.

By Proposition 5.5.21, A_{2} = a_{2} − ρδ(den(β))/den(β) and B_{2} = b_{2}, where ρ denotes the order of L. It follows from (5.18) and den(β) ∈ k[x] that

L(h_{2}) ∈ \frac{V}{Uden(β)^{ρ}} H(a_{2}, b_{2}),

which, together with (5.17) and H(a_{2}, b_{2}) = {cH | c ∈ k}, implies that

M (\frac{v}{u}) = \frac{cV}{Uden(β)^{ρ}} for some c ∈ k.

Hence, there exist a nonzero operator M′ ∈ k[x, n]⟨S_{n}⟩ with leading coefficient in k[n] and a rational function g ∈ k(x) such that M′(gv/u) belongs to k[x, n]. Consequently, u is split by Lemma 5.2.7, and, hence, h_{2} is proper.

Proof of the necessity in Theorem 5.6.2

Proof. Assume that h has a telescoper with respect to n, and that h_{2} is nonzero. Our goal is to show that h_{2} is proper. By Theorem 5.5.18,

h_{2} = \frac{v}{u}H for some H ∈ H (a_{2}, b_{2}) ,

where a_{2} = γ̃ +n\frac{δ(β)}{β} for some γ̃ ∈ k(x), and (u, v, b_{2}) is an AP triple. Without loss of generality, we further assume that u is in k[x, n]. To prove that h_{2} is proper, it suffices to prove that u is split by Proposition 5.4.4.

Let L ∈ k(x)⟨D_{x}⟩ be a telescoper for h with respect to n. By Lemma 5.6.3, L is also a telescoper for h_{2} with respect to n. Then

L(h_{2}) = L (\frac{v}{u}H) = Δ_{n}(g),  (5.19)

for some hyperexponential-hypergeometric function g. For all i ∈ ℕ, let f_{i} be the rational function in k(x, n) such that

D_{x}^{i}(H) = f_{i}H.

Then den(f_{i}) is in k[x], because den(a_{2}) is in k[x]. It follows from Leibniz’s formula that there exists M_{i} ∈ k(x)[n]⟨D_{x}⟩ such that

D_{x}^{i}(h_{2}) = M_{i} (\frac{v}{u}) H.

Moreover, the leading coefficient of M_{i} is in k. Consequently, there exists M ∈ k(x)[n]⟨D_{x}⟩ with leading coefficient in k(x) such that

L(h_{2}) = M (\frac{v}{u}) H.  (5.20)

On the other hand, Proposition 5.5.22 and H(a_{2}, b_{2}) = {cH | c ∈ k} imply,

L(h_{2}) = \frac{V}{U} H.  (5.21)

where U and V are in k(x)[n], and (U, V, b_{2}) is an AP triple. Since L(h_{2}) is hypergeometricsummable by (5.19), it follows from Lemma 5.5.16 that U is free of n.

Comparing (5.20) with (5.21), we have that M(v/u) = V/U, which is in k(x)[n]. It follows that there exists a nonzero operator M′ ∈ k[x, n]⟨D_{x}⟩ with leading coefficient in k[x] such that M′(v/u) belongs to k[x, n]. By Lemma 5.2.7, u is split, and then h_{2} is proper.

5.6.2 Properness implies the existence of telescopers

Wilf and Zeilberger present an elementary proof of the existence of telescopers for proper multivariate hypergeometric terms [92, Theorem 3.1], and indicate that their argument should be also valid in the continuous-discrete setting. For the sake of completeness, we elaborate on the proof of the following theorem.

Theorem 5.6.4. If h is a proper hyperexponential-hypergeometric function over k(x, n), then h has a telescoper with respect to n and a telescoper with respect to x.

Preparation lemmas

Before the proof of the theorem above, we shall show some lemmas.

Lemma 5.6.5. Let h ∈ H(a, b) be a hyperexponential-hypergeometric function over k(x, n).

(i) If h has a telescoper with respect to n, so does the product Hh for any hyperexponentialhypergeometric function H with Δ_{n}(H) = 0.

(ii) If h has a telescoper with respect to x, so does the product Hh for any hyperexponentialhypergeometric function H with D_{x}(H) = 0.

Proof. Assume that H is a nonzero element of H(A, B). If Δ_{n}(H) = 0, then B = 1 and a ∈ k(x) by the integrability condition σ(A)−A = δ(B)/B. It is straightforward to verify that HD_{x}(h) = (D_{x} − A)Hh, which, together with an easy induction, implies that

HD_{x}^{i}(h) = (D_{x} − A)^{i}(Hh) for all i ∈ ℕ.

Let L ∈ k(x)⟨D_{x}⟩ be a telescoper of h. By the equality above,

HL(h) = M(Hh)

where M is obtained by replacing D_{x}^{i} in L by (D_{x} − A)^{i}. Since there exists a hypergeometric function g such that L(h) = Δ_{n}(g) and Δ_{n}(H) = 0, the equality above implies that

M(Hh) = HΔ_{n}(g) = Δ_{n}(Hg).

The first assertion is proved.

To prove the second assertion, assume that H is a nonzero element in H(A, B). If D_{x}(H) is equal to 0, then A = 0 and B ∈ k(n) by the integrability condition σ(A) − A = δ(B)/B. Since HS_{n}(h) = S_{n}(Hh)/B, we have that, for all i ∈ ℕ,

HS_{n}^{i}(h) = B_{i}S_{n}^{i}(Hh) for some B_{i} ∈ k(n).

Let L ∈ k(n)⟨S_{n}⟩ be a telescoper of h. By the equality above,

HL(h) = M(Hh)

where M is obtained by replacing S_{n}^{i} in L by B_{i}S_{n}^{i}. Since there exists a hypergeometric function g such that L(h) = D_{x}(g) and D_{x}(H) = 0, the equality above implies that

M(Hh) = HD_{x}(g) = D_{x}(Hg).

The second assertion is proved.

Lemma 5.6.6. Let h_{1} and h_{2} be two similar hyperexponential-hypergeometric functions over the field k(x, n). Let 𝕍_{n} be the vector space spanned by {Δ_{n}^{i}(h_{j}) | i ∈ ℕ, j ∈ {1, 2}} over k(x), and 𝕍_{x} be the vector space spanned by {D_{x}^{i}(h_{j}) | i ∈ ℕ, j ∈ {1, 2}} over k(n).

1. If both h_{1} and h_{2} have telescopers with respect to n, so does any nonzero element in 𝕍_{n}.

2. If both h_{1} and h_{2} have telescopers with respect to x, so does any nonzero element in 𝕍_{x}. Proof. Assume that there exists a nonzero L_{j} ∈ k(x)⟨D_{x}⟩ such that L_{j}(h_{j}) = Δ_{n}(g_{j}) for some hyperexponential-hypergeometric function g_{j}, and for j = 1, 2. Then

L_{j}(Δ_{n}^{i}(h_{j})) = Δ_{n}(Δ_{n}^{i}(g_{j})).

Hence, every element of {Δ_{n}^{i}(h_{j}) | i ∈ ℕ, j ∈ {1, 2}} has a telescoper. If two hypergeometric functions have telescopers with respect to n, say T_{1} and T_{2}, then their sum, assumed to be nonzero, also has a telescoper with respect to n, which can be taken as a common left multiple of T_{1} and T_{2}. The first assertion then follows from Lemma 5.6.5 (ii). The second assertion can be proved in the same manner.

In the proof of the Fundamental Theorem [92], Wilf and Zeilberger neglected the case that the remainder of a nonzero m-free operator P (n, S_{n}, S_{m}) ∈ k(n)⟨S_{n}, S_{m}⟩ by S_{m} − 1 may be zero. Wegschaider used a “noncommutative trick” to deal with this case. Here, we will apply Wegschaider’s trick and its differential analogue to transform an x-free or n-free operator to a telescoper with respect to n or x. Lemma 5.6.7 (Wegschaider’s trick). Let h be a hypergeometric function over k(x, n).

(i) If there exists a nonzero operator A ∈ k(x)⟨S_{n}, D_{x}⟩ such that A(h) = 0, then h has a telescoper with respect to n.

(ii) If there exists a nonzero operator A ∈ k(n)⟨S_{n}, D_{x}⟩ such that A(h) = 0, then h has a

telescoper with respect to x. Proof. Assume that there exists a nonzero operator A in k(x)⟨S_{n}, D_{x}⟩ such that A(h) = 0. Since Δ_{n} = S_{n} − 1, we can write A = Δ_{n}^{m}(L(x, D_{x}) + Δ_{n}M), where m is in ℕ, L is a nonzero operator in k(x)⟨D_{x}⟩ and M is in k(x)⟨S_{n}, D_{x}⟩. By Wegschaider’s trick in [91, Theorem 3.2], there exist w ∈ k(n) and r ∈ k with r ≠ 0 such that

wΔ_{n}^{m} = Δ_{n}Q + r  (5.22)

for some Q ∈ k(n)⟨S_{n}⟩. In particular, r = (−1)^{m}m! if we take w = n^{m}. Using the fact rΔ_{n} = Δ_{n}r and (5.22), we find

\frac{w}{r} A = L + Δ_{n}N for some N ∈ k(x, n)⟨D_{x}, S_{n}⟩.

Hence, L is a telescoper for h with respect to n.

The second assertion can be proved in a similar way. Instead of (5.22), we need to find w ∈ k(x) and r ∈ k \ {0} such that

wD_{x}^{m} = D_{x}Q + r

for some Q ∈ k(x)⟨D_{x}⟩. In particular, r = (−1)^{m}m! if we take w = x^{m}.

Proof of Theorem 5.6.4

Now, we present the proof of Theorem 5.6.4.

Proof. Assume that h is a proper hyperexponential-hypergeometric function. Then the standard representation of h is of the form

(p(x, n), β(x), γ(x), α(n)),

where p ∈ k[x, n], β, γ ∈ k(x), and α ∈ k(n). Writing p as ∑_{i=0}^{m} p_{i}(x)n^{i} with p_{i} ∈ k[x] yields

pβ^{n}H(γ, α) = ∑_{i=0}^{m} p_{i}n^{i}β^{n}H(γ, α) = ∑_{i=0}^{m} β^{n}H (γ + \frac{δ(p_{i})}{p_{i}}, α(n)\frac{σ(n^{i})}{n^{i}}) .

It follows that

h = ∑_{i=0}^{m} β(x)^{n}G_{i}H_{i}  (5.23)

for some G_{i} ∈ H (γ(x) + \frac{δ(p_{i})}{p_{i}}, 1) and H_{i} ∈ H (0, α(n)\frac{σ(n^{i})}{n^{i}}) .

First, we show that h has a telescoper with respect to n. Note that each G_{i} in (5.23) has n-certificate equal to 1. By Lemmas 5.6.5 and 5.6.6, it suffices to show that a hypergeometric function ĥ ∈ β(x)^{n}·H (0, g(n)) with g ∈ k(n) has a telescoper with respect to n. Let s = num(β), t = den(β), a = num(g) and b = den(g). Moreover, let v = num(δ(β)/β) and w = den(δ(β)/β). A straightforward calculation yields that the x-certificate and n-certificate of ĥ are, respectively,

\frac{nv}{w} and \frac{sa}{tb}.

Note that s, t, v, w are in k[x], and a, b in k[n].

We claim that there exists a nonzero operator A in k(x)⟨S_{n}, D_{x}⟩ such that A(ĥ) = 0. The claim will be proved by a well-known argument used in [67]. The first assertion will then follow from Lemma 5.6.7.

Let F_{N} be the linear space spanned by {S_{n}^{i}D_{x}^{j} | i+j ≤ N} over k(x). Let μ be the maximum of degrees of a and b in n, and let

W_{N} = span_{k(x)} \{\frac{n^{i}ĥ}{b(n + N − 1) ⋯ b(n)} | i ≤ (μ + 1)N\} .

An easy induction on i and j yields

S_{n}^{i}D_{x}^{j}(ĥ) = \frac{q(n, x)ĥ}{b(n + i − 1) ⋯ b(n)}, where deg_{n}(q) ≤ iμ + j.

Hence, S_{n}^{i}D_{x}^{j}(ĥ) belongs to W_{N} if i+j ≤ N. Accordingly, there is a k(x)-linear map ϕ_{N} from F_{N} to W_{N} that sends L to L(ĥ) for all L ∈ F_{N}. Since the dimension of F_{N} over k(x) is \binom{N+2}{2}, while that of W_{N} is (μ + 1)N + 1, the kernel of ϕ_{N} is nontrivial when N is sufficiently large. Any nonzero element in the kernel annihilates ĥ. The claim is established.

Second, we show that h has a telescoper with respect to x. Note that each nonzero element in H (0, α(n)\frac{σ(n^{i})}{n^{i}}) in (5.23) has x-certificate equal to zero. By Lemmas 5.6.5 and 5.6.6, it suffices to show that any hypergeometric function ĥ(x, n) in the set

β(x)^{n} · H (g(x), 1) with g ∈ k(x)

has a telescoper with respect to x. Let s = num(β), t = den(β),

v = num (n\frac{δ(β)}{β} + g) and w = den (n\frac{δ(β)}{β} + g) .

A straightforward calculation yields that the x-certificate and n-certificate of ĥ are, respectively,

\frac{v(x, n)}{w(x)} and \frac{s(x)}{t(x)}.

Note that s, t, and w are in k[x] and v is in k[x, n] with deg_{n}(v) = 1.

We claim that there exists a nonzero operator A in k(n)⟨S_{n}, D_{x}⟩ such that A(ĥ) = 0. The second assertion then follows from Lemma 5.6.7.

Following an argument similar to the one above, we now consider a linear space F_{N} spanned by {S_{n}^{i}D_{x}^{j} | i + j ≤ N} over k(n). Let μ be the maximum of the degrees in x of s, t, v and w, and let

W_{N} = span_{k(n)} \{\frac{x^{i}ĥ}{(tw)^{N}} | i ≤ 2μN\} .

Algorithm ExistenceContinuous Input: A hyperexponential-hypergeometric function h ∈ H(a, b) with a, b ∈ k(x, n). Output: True, if h has a telescoper with respect to x; false, otherwise.

1. Compute a strict DRNF (K, S) of a with respect to x;

2. Apply ReduceCert in [47] to (K, S) to get the decomposition

S = δ(S_{1}) + S_{1}K + \frac{v}{u den(K)^{i}} .

3. If u is split, then return true, otherwise, return false.

Figure 5.1: Algorithm for deciding the existence of telescopers with respect to x.

An easy induction on i and j yields

S_{n}^{i}D_{x}^{j}(ĥ) = \frac{q(n, x)}{t^{i}w^{j}} ĥ, where deg_{x}(q) ≤ (i + j)μ.

Hence, S_{n}^{i}D_{x}^{j}(ĥ) belongs to W_{N} if i+j ≤ N. Accordingly, there is a k(n)-linear map ϕ_{N} from F_{N} to W_{N} that sends L to L(ĥ) for all L ∈ F_{N}. Since the dimension of F_{N} over k(n) is \binom{N+2}{2}, while that of W_{N} is 2μN + 1, the kernel of ϕ_{N} is nontrivial when N is sufficiently large. Any nonzero element in the kernel annihilates ĥ. The claim is established.

Remark 5.6.8. Theorems 5.6.1 and 5.6.2 hold by the proofs of necessity in Theorems 5.6.1 and 5.6.2, and by Theorem 5.6.4

5.7 Algorithms and examples

5.7.1 Algorithms

According to our criteria, we only use the continuous, resp. discrete, certificate for deciding the existence of telescopers with respect to the continuous, resp. discrete, variable. Algorithms in Figures 5.1 and 5.2 are based on strict DRNF and strict SRNF computation, two additive decompositions, and our criteria.

Algorithm ExistenceDiscrete Input: A hyperexponential-hypergeometric function h ∈ H(a, b) with a, b ∈ k(x, n). Output: True, if h has a telescoper with respect to n; false, otherwise.

1. Compute a strict SRNF (K, S) of b with respect to n;

2. Apply dcert in [13] to (K, S) to get the decomposition

S = σ(S_{1})K − S_{1} + \frac{v}{u (σ^{−1}(num(K)))^{i}den(K)^{j}} .

3. If u is split, then return true, otherwise, return false. Figure 5.2: Algorithm for deciding the existence of telescopers with respect to n.

5.7.2 Examples

Example 5.7.1. It is possible that a hyperexponential-hypergeometric function has a telescoper with respect to one variable but no telescoper with respect to the other variable. Consider the rational function

h = \frac{1}{(x + n)^{2}} .

Applying Hermite reduction to h with respect to x yields

h = D_{x} (\frac{−1}{x + n}) ,

which implies that 1 is a telescoper for h with respect to x. However, h has no telescoper with respect to n, because (x + n)^{2} is shift-free with respect to n and it is not split. Similarly, consider the rational function

h = \frac{1}{(x + n)(x + n + 1)}.

Since h = Δ_{n}(−1/(x+n)), 1 is a telescoper for h with respect to n. However, h has no telescoper with respect to x, because (x + n)(x + n + 1) is squarefree with respect to x and it is not split.

Example 5.7.2. As we know, the properness is only a sufficient condition for the existence of telescopers. We show a concrete example to illustrate this fact. Consider the hyperexponentialhypergeometric function

h(x, n) = \frac{−x + 2nx + 2n^{2}}{(x + n)^{2}x} · x^{n} · e^{−x}.

Though h is not proper, it still has a telescoper with respect to x since h can be decomposed into

h = D_{x} (\frac{1}{x + n} · x^{n} · e^{−x}) + x^{n−1} · e^{−x},

where x^{n−1} · e^{−x} is proper. In fact, our criteria indicate that all nonproper examples are of this form.

Conclusion and Perspectives

In this thesis, we have focused on the practical efficiency, theoretical complexity, and the termination of creative-telescoping algorithms. For bivariate rational functions, we have managed to blend the general method of creative telescoping with the Hermite reduction. According to our complexity analyses and experiments, the Hermite reduction based method has been proved to outperform the existing methods. In order to study the general case of Wilf and Zeilberger’s conjecture, we have presented a structure theorem for multivariate hyperexponentialhypergeometric functions. With the help of standard representations and two additive decompositions, we have derived two existence criteria for deciding whether a bivariate hyperexponentialhypergeometric function has a telescoper with respect to the continuous variable or the discrete one. Our criteria solve the termination problem of Zeilberger-style algorithms for bivariate hyperexponential-hypergeometric functions. Still, there are many further studies to be carried out. In what follows, we would like to mention a few directions for future work.

6.1 Extensions of the Hermite-reduction based method

For the moment, the Hermite-reduction based method only works for bivariate rational-function inputs. As we have already seen in Chapter 5, Hermite reduction can be extended to hyperexponential or hypergeometric cases. So it is natural to wonder whether the approach of the Hermite-reduction based method can be used to compute the minimal telescopers for bivariate hyperexponential functions or hypergeometric terms. We have considered two possible extensions: one is for hyperexponential functions, especially in simple radical extensions; and the other is for multivariate rational functions. In the multivariate rational case, the corresponding telescoping problem is as follows.

Problem 6.1.1. Given a rational function f ∈ k(x, y_{1}, …, y_{n}), construct a nonzero operator L ∈ k(x)⟨D_{x}⟩ such that

L(x, D_{x})(f) = ∑_{i=1}^{n} D_{y_{i}}(g_{i}) for some g_{i} ∈ k(x, y_{1}, …, y_{n}).

In order to solve the above problem, we have tried to extend the Hermite reduction to the multivariate case. That is, we would like to decompose a multivariate rational function f in k(y_{1}, …, y_{n}) into

f = D_{y_{1}}(g_{1}) + ⋯ + D_{y_{n}}(g_{n}) + r,

where r is minimal in some sense. For a special class of rational functions, we show this can be obtained.

Let f = P/Q^{m} be a rational function in k(y_{1}, …, y_{n}) with m > 1, P, Q ∈ k[y_{1}, …, y_{n}], and gcd(P, Q) = 1. Let I be the ideal generated by Q, D_{y_{1}}(Q), …, D_{y_{n}}(Q) in k[y_{1}, …, y_{n}]. We observe that the ideal I is equal to the ring k[y_{1}, …, y_{n}] in certain applications, see examples gathered in [74]. In the sequel, we proceed the discussion under the above constraints on Q. Then we can always express P in the form

P = AQ + ∑_{i=1}^{n} B_{i}D_{y_{i}}(Q),

where A, B_{1}, …, B_{n} are polynomials in k(y_{1}, …, y_{n}). Now, integration by parts yields

\frac{P}{Q^{m}} = \frac{AQ + ∑_{i=1}^{n} B_{i}D_{y_{i}}(Q)}{Q^{m}} = ∑_{i=1}^{n} D_{y_{i}} (\frac{−B_{i}}{mQ^{m−1}} ) + \frac{A + m^{−1} ∑_{i=1}^{n} D_{y_{i}}(B_{i})}{Q^{m−1}}.

Repeating the above process at most m − 1 times, we get

\frac{P}{Q^{m}} = ∑_{i=1}^{n} D_{y_{i}}(g_{i}) + \frac{P̃}{Q},  (6.1)

where g_{i} ∈ k(y_{1}, …, y_{n}) and P̃ ∈ k[y_{1}, …, y_{n}]. In the univariate case when Q is squarefree, this is exactly Hermite reduction. When n > 1, the essential difference from the univariate case is that one cannot make the remaining fraction P̃/Q proper with respect to all the variables. The uniqueness of the decomposition in (6.1) is also not well studied. In the future, we would like to understand this decomposition better and expect it to be useful to solve the telescoping problem for multivariate rational functions.

6.2 Existence criteria for telescopers in the q-setting

There are nine termination problems of Zeilberger’s method for bivariate hyperexponentialhypergeometric functions if telescopers could involve the q-shift operator, denoted by Q_{n}. In Table 6.1, we mark the solved instances by ✓ and unsolved by ?. In order to solve the remaining termination problems, we would like to extend Theorem 4.4.1 to the q-shift setting. Moreover, we expect a unified way to deal with all of the nine cases.

\begin{matrix}(L, g) & D_{y} & S_{n} − 1 & Q_{n} − 1 \\ L(x, D_{x}) & ✓ & ✓ & ? \\ L(n, S_{n}) & ✓ & ✓ & ? \\ L(n, Q_{n}) & ? & ? & ✓ \end{matrix}

Table 6.1: Solved and unsolved termination problems

6.3 Wilf and Zeilberger’s conjecture

We present some remarks on the general case of Wilf and Zeilberger’s conjecture. First, we recall basic notation about holonomic functions and modules following the presentation in Coutinho’s book [38]. Let A_{m} denote the m-th Weyl algebra k[x_{1}, …, x_{m}]⟨D_{1}, …, D_{m}⟩, in which we have the following multiplication rule: for all 1 ≤ i, j ≤ m,

x_{i}x_{j} = x_{j}x_{i}, D_{i}D_{j} = D_{j}D_{i}, and D_{i}x_{j} − x_{j}D_{i} = δ_{ij},

where δ_{ij} is the Kronecker delta symbol: it equals one if i = j and zero otherwise. A finitely generated left A_{m}-module is holonomic if it is zero, or if it has Hilbert dimension m [38, Chapter 5]. By Bernstein’s inequality, any finitely generated nonzero left A_{m}-module has dimension at least m. So holonomicity indicates the minimality of dimension for nonzero A_{m}-modules.

Let H(x) be a function in a left A_{m}-module of functions. Then one can define the annihilating ideal of H(x) in A_{m} by

ann_{A_{m}}(H) = {p ∈ A_{m} | p · H = 0}.

The function H(x) is said to be holonomic if its annihilating ideal is holonomic as a left A_{m}- module. When the function H(x) can be viewed as an element in a left module over the ring k(x_{1}, …, x_{m})⟨D_{1}, …, D_{m}⟩ of differential operators, H(x) is said to be D-finite if the vector space generated by all derivatives D_{1}^{i_{1}} ⋯ D_{m}^{i_{m}}(H) over k(x) is finite dimensional over k(x) [67, 68]. In the continuous case, an elementary proof of Kashiwara’s equivalence [59] between holonomicity and D-finiteness was presented by Takayama in [85]. For this reason, we choose the following definition, used in [24, page 59], for holonomicity of continuous-discrete functions.

Definition 6.3.1. A function h(x, y) is said to be holonomic if its associated generating function

H(x, z) = ∑_{y_{1},...,y_{n}∈ℕ}h(x, y)z_{1}^{y_{1}} ⋯ z_{n}^{y_{n}}

is a holonomic function in a left module over the ring k(x, z)⟨D_{x_{1}}, …, D_{x_{m}}, D_{z_{1}}, …, D_{z_{n}}⟩.

Now, we can describe the conjecture by Wilf and Zeilberger [92, page 585] on holonomic hyperexponential-hypergeometric functions.

Conjecture 6.3.2 (Wilf and Zeilberger). Let h(x, y) be a hyperexponential-hypergeometric function over k(x, y). Then h is holonomic if and only if it is proper.

In the case of several discrete variables, a slightly modified version of the conjecture was proved independently by Payne in his Ph.D. thesis [73] and by Abramov and Petkovšek [14]. In particular, the case of two variables has also been shown by Hou [54, 55] and by Abramov and Petkovšek [12]. In the continuous case, any multivariate hyperexponential function is D-finite, and then holonomic by Kashiwara’s equivalence. By Corollary 4.5.4, it is also proper. Thus, Wilf and Zeilberger’s conjecture holds naturally in this case.

In the discrete case, the structure theorem on multivariate hypergeometric terms helps us reduce Wilf and Zeilberger’s conjecture to rational case [73, 14], where one only need to show that any holonomic rational functions are proper. For the moment, we have not obtained any complete proof of the general continuous-discrete case of the conjecture. In the future, we would like to prove the following conjecture on holonomic rational functions in the mixed setting.

Conjecture 6.3.3. Let f(x, y) be a rational function in k(x, y). Then f is holonomic if and only if the denominator of f splits into the form

A(x) ∏_{ℓ=1}^{L}(v_{ℓ} · y + λ_{ℓ})

where A ∈ k[x], λ_{1}, …, λ_{L} ∈ k̄, v_{ℓ} ∈ ℤ^{n} for all ℓ with 1 ≤ ℓ ≤ L, and the v_{ℓ} · y means the inner product of vectors v_{ℓ} and y.

116

Bibliography

[1] N. H. Abel. Œuvres Complètes de Niels Henrik Abel. Tome II. Imprimerie de Grøndahl & Son, Christiania, 1981. Contenant les mémoirs posthumes d’Abel. [Containing the posthumous memoirs of Abel], Edited and with notes by L. Sylow and S. Lie.

[2] S. A. Abramov. The rational component of the solution of a first order linear recurrence relation with rational right hand side. Ž. Vyčisl. Mat. i Mat. Fiz., 15(4):1035–1039, 1090, 1975.

[3] S. A. Abramov. Rational solutions of linear differential and difference equations with polynomial coefficients. USSR Comput. Math. Math. Phys., 29(6):7–12, 1989.

[4] S. A. Abramov. Rational solutions of linear difference and q-difference equations with polynomial coefficients. In ISSAC ’95: Proceedings of the 1995 International Symposium on Symbolic and Algebraic Computation, pages 285–289, New York, NY, USA, 1995. ACM.

[5] S. A. Abramov. Applicability of Zeilberger’s algorithm to hypergeometric terms. In IS-SAC’02: Proceedings of the 2002 International Symposium on Symbolic and Algebraic Computation, pages 1–7, New York, 2002. ACM.

[6] S. A. Abramov. When does Zeilberger’s algorithm succeed? Adv. in Appl. Math., 30(3):424– 441, 2003.

[7] S. A. Abramov, J. J. Carette, K. O. Geddes, and H. Q. Le. Telescoping in the context of symbolic summation in Maple. J. Symbolic Comput., 38(4):1303–1326, 2004.

[8] S. A. Abramov and K. Y. Kvasenko. Fast algorithms to search for the rational solutions of linear differential equations with polynomial coefficients. In ISSAC ’91: Proceedings of the 1991 International Symposium on Symbolic and Algebraic Computation, pages 267–270, New York, NY, USA, 1991. ACM.

[9] S. A. Abramov and H. Q. Le. A criterion for the applicability of Zeilberger’s algorithm to rational functions. Discrete Math., 259(1-3):1–17, 2002.

[10] S. A. Abramov, H. Q. Le, and M. Petkovšek. Rational canonical forms and efficient representations of hypergeometric terms. In ISSAC’03: Proceedings of the 2003 International Symposium on Symbolic and Algebraic Computation, pages 7–14, New York, 2003. ACM.

[11] S. A. Abramov and M. Petkovšek. Minimal decomposition of indefinite hypergeometric sums. In ISSAC’01: Proceedings of the 2001 International Symposium on Symbolic and Algebraic Computation, pages 7–14 (electronic), New York, 2001. ACM.

[12] S. A. Abramov and M. Petkovšek. Proof of a conjecture of Wilf and Zeilberger, 2001. Preprints Series of the Institute of Mathematics, Physics and Mechanics 39(748)(2001), Ljubljana, March 9, 2001.

[13] S. A. Abramov and M. Petkovšek. Rational normal forms and minimal decompositions of hypergeometric terms. J. Symbolic Comput., 33(5):521–543, 2002.

[14] S. A. Abramov and M. Petkovšek. On the structure of multivariate hypergeometric terms. Adv. Appl. Math., 29(3):386–411, 2002.

[15] G. Almkvist and D. Zeilberger. The method of differentiating under the integral sign. J. Symbolic Comput., 10:571–591, 1990.

[16] M. Apagodu and D. Zeilberger. Multi-variable Zeilberger and Almkvist-Zeilberger algorithms and the sharpening of Wilf- Zeilberger theory. Adv. in Appl. Math., 37(2):139–152, 2006.

[17] J. Bernšteĭn. Modules over a ring of differential operators. An investigation of the fundamental solutions of equations with constant coefficients. Funkcional. Anal. i Priložen., 5(2):1–16, 1971.

[18] A. Bostan, S. Chen, F. Chyzak, and Z. Li. Complexity of creative telescoping for bivariate rational functions. In ISSAC ’10: Proceedings of the 2010 International Symposium on Symbolic and Algebraic Computation, pages 203–210, New York, NY, USA, 2010. ACM.

[19] A. Bostan, F. Chyzak, T. Cluzeau, and B. Salvy. Low complexity algorithms for linear recurrences. In J.-G. Dumas, editor, ISSAC’06: Proceedings of the 2006 International Symposium on Symbolic and Algebraic Computation, pages 31–38. ACM Press, 2006.

[20] A. Bostan, F. Chyzak, B. Salvy, G. Lecerf, and É. Schost. Differential equations for algebraic functions. In ISSAC’07: Proceedings of the 2007 International Symposium on Symbolic and Algebraic Computation, pages 25–32. ACM, New York, 2007.

[21] M. Bronstein. Symbolic Integration I: Transcendental Functions, volume 1 of Algorithms and Computation in Mathematics. Springer-Verlag, Berlin, second edition, 2005.

[22] M. Bronstein, Z. Li, and M. Wu. Picard–Vessiot extensions for linear functional systems. In ISSAC ’05: Proceedings of the 2005 International Symposium on Symbolic and Algebraic Computation, pages 68–75, New York, USA, 2005. ACM.

[23] M. Bronstein and M. Petkovšek. An introduction to pseudo-linear algebra. Theoretical Computer Science, 157:3–33, 1996.

[24] P. Cartier. Démonstration “automatique” d’identités et fonctions hypergéométriques (d’après D. Zeilberger). Astérisque, (206):Exp. No. 746, 3, 41–91, 1992. Séminaire Bourbaki, Vol. 1991/92.

[25] D. Cerveau and J.-F. Mattei. Formes Intégrables Holomorphes Singulières, volume 97 of Astérisque. Société Mathématique de France, Paris, 1982. With an English summary.

[26] S. Chen, F. Chyzak, R. Feng, and Z. Li. The existence of telescopers for hyperexponentialhypergeometric functions, 2010. MM-Res. Preprints (2010) No. 29, 239-267.

[27] S. Chen and Z. Li. A multiplicative form of multivariate hyperexponential-hypergeometric functions, 2010. MM-Res. Preprints (2010) No. 29, 25-35.

[28] W. Y. C. Chen, Q.-H. Hou, and Y.-P. Mu. Applicability of the q-analogue of Zeilberger’s algorithm. J. Symbolic Comput., 39(2):155–170, 2005.

[29] C. Christopher. Liouvillian first integrals of second order polynomial differential equations. Electron. J. Differential Equations, pages No. 49, 7 pp. (electronic), 1999.

[30] D. V. Chudnovsky and G. V. Chudnovsky. On expansion of algebraic functions in power and Puiseux series. I. J. Complexity, 2(4):271–294, 1986.

[31] D. V. Chudnovsky and G. V. Chudnovsky. On expansion of algebraic functions in power and Puiseux series. II. J. Complexity, 3(1):1–25, 1987.

[32] F. Chyzak. Fonctions Holonomes en Calcul Formel. PhD thesis, École Polytechnique, Paris, France, 1998. INRIA, TU 0531. 227 pages.

[33] F. Chyzak. An extension of Zeilberger’s fast algorithm to general holonomic functions. Discrete Math., 217(1-3):115–134, 2000.

[34] F. Chyzak, M. Kauers, and B. Salvy. A non-holonomic systems approach to special function identities. In ISSAC ’09: Proceedings of the 2009 International Symposium on Symbolic and Algebraic Computation, pages 111–118, New York, NY, USA, 2009. ACM.

[35] F. Chyzak and B. Salvy. Non-commutative elimination in Ore algebras proves multivariate identities. J. Symbolic Comput., 26(2):187–227, 1998.

[36] R. M. Cohn. Difference Algebra. Interscience Publishers John Wiley & Sons, New York-London-Sydeny, 1965.

[37] O. Cormier, M. F. Singer, B. M. Trager, and F. Ulmer. Linear differential operators for polynomial equations. J. Symbolic Comput., 34(5):355–398, 2002.

[38] S. C. Coutinho. A Primer of Algebraic D-modules, volume 33 of London Mathematical Society Student Texts. Cambridge University Press, Cambridge, 1995.

[39] J. H. Davenport. The Risch differential equation problem. SIAM J. Comput., 15(4):903–918, 1986.

[40] M. C. Fasenmyer. Some generalized hypergeometric polynomials. Bull. Amer. Math. Soc., 53:806–812, 1947.

[41] R. Feng. A note on the Galois group of linear differential-difference equations, 2010. MM-Res. Preprints (2010) No. 29, 1–8.

[42] R. Feng, M. F. Singer, and M. Wu. An algorithm to compute Liouvillian solutions of prime order linear difference-differential equations. J. Symbolic Comput., 45(3):306–323, 2010.

[43] R. Feng, M. F. Singer, and M. Wu. Liouvillian solutions of linear difference-differential equations. J. Symbolic Comput., 45(3):287–305, 2010.

[44] R. V. Gamkrelidze, editor. Analysis. I, volume 13 of Encyclopaedia of Mathematical Sciences. Springer-Verlag, Berlin, 1989. Integral representations and asymptotic methods, A translation of Sovremennye problemy matematiki. Fundamentalnye napravleniya, Tom 13, Akad. Nauk SSSR, Vsesoyuz. Inst. Nauchn. i Tekhn. Inform., Moscow, 1986, Translation by D. Newton, Translation edited by R. V. Gamkrelidze.

[45] K. O. Geddes, S. R. Czapor, and G. Labahn. Algorithms for Computer Algebra. Kluwer Academic Publishers, Boston, MA, 1992.

[46] K. O. Geddes and H. Q. Le. An algorithm to compute the minimal telescopers for rational functions (differential-integral case). In Mathematical Software, pages 453–463. WSP, 2002.

[47] K. O. Geddes, H. Q. Le, and Z. Li. Differential rational normal forms and a reduction algorithm for hyperexponential functions. In ISSAC’04: Proceedings of the 2004 International Symposium on Symbolic and Algebraic Computation, pages 183–190, New York, USA, 2004. ACM.

[48] J. Gerhard. Fast modular algorithms for squarefree factorization and Hermite integration. Appl. Algebra Engrg. Comm. Comput., 11(3):203–226, 2001.

[49] J. Gerhard. Modular Algorithms in Symbolic Summation and Symbolic Integration (LNCS). SpringerVerlag, 2004.

[50] J. Gerhard, M. Giesbrecht, A. Storjohann, and E. V. Zima. Shiftless decomposition and polynomial-time rational summation. In Proceedings of the 2003 International Symposium on Symbolic and Algebraic Computation, ISSAC ’03, pages 119–126, New York, NY, USA, 2003. ACM.

[51] R. W. Gosper, Jr. Decision procedure for indefinite hypergeometric summation. Proc. Nat. Acad. Sci. U.S.A., 75(1):40–42, 1978.

[52] C. Hermite. Sur l’intégration des fractions rationnelles. Ann. Sci. École Norm. Sup. (2), 1:215–218, 1872.

[53] E. Horowitz. Algorithms for partial fraction decomposition and rational function integration. In SYMSAC’71, pages 441–457, New York, USA, 1971. ACM.

[54] Q.-H. Hou. Algebraic Method in Combinatorics. PhD thesis, Nankai University, Tianjin, China, 2001.

[55] Q.-H. Hou. k-free recurrences of double hypergeometric terms. Adv. in Appl. Math., 32(3):468–484, 2004.

[56] I. Kaplansky. An Introduction to Differential Algebra. Actualités Sci. Ind., No. 1251, Publ. Inst. Math. Univ. Nancago, No. 5. Hermann, Paris, 1957.

[57] M. Karr. Summation in finite terms. J. Assoc. Comput. Mach., 28(2):305–350, 1981.

[58] M. Karr. Theory of summation in finite terms. J. Symbolic Comput., 1(3):303–315, 1985.

[59] M. Kashiwara. On the holonomic systems of linear differential equations. II. Invent. Math., 49(2):121–135, 1978.

[60] W. Koepf. Hypergeometric Summation. Advanced Lectures in Mathematics. Friedr. Vieweg & Sohn, Braunschweig, 1998. An algorithmic approach to summation and special function identities.

[61] E. R. Kolchin. Differential Algebra and Algebraic Groups. Academic Press, New York, 1973. Pure and Applied Mathematics, Vol. 54.

[62] G. Labahn and Z. Li. Hyperexponential solutions of finite-rank ideals in orthogonal Ore rings. In ISSAC’04: Proceedings of the 2004 International Symposium on Symbolic and Algebraic Computation, pages 213–220. ACM, New York, 2004.

[63] H. Q. Le. On the differential-integral analogue of Zeilberger’s algorithm to rational functions. In Computer mathematics (Chiang Mai, 2000), volume 8 of Lecture Notes Ser. Comput., pages 204–213. World Sci. Publ., River Edge, NJ, 2000.

[64] H. Q. Le. On the q-analogue of Zeilberger’s algorithm to rational functions. Program. Comput. Softw., 27(1):35–42, 2001.

[65] H. Q. Le. Algorithms for the construction of the minimal telescopers. PhD thesis, University of Waterloo, Waterloo, Ont., Canada, 2003. AAINQ83003.

[66] Z. Li and D. Zheng. Determining whether a multivariate hyperexponential function is algebraic. J. Syst. Sci. Complex., 19(3):352–364, 2006.

[67] L. M. Lipshitz. The diagonal of a D-finite power series is D-finite. J. Algebra, 113(2):373–378, 1988.

[68] L. M. Lipshitz. D-finite power series. J. Algebra, 122(2):353–373, 1989.

[69] D. Mack. On rational integration. Research Report UCP-38, 1975.

[70] O. Ore. Sur la forme des fonctions hypergéométriques de plusieurs variables. J. Math. Pures Appl. (9), 9(4):311–326, 1930.

[71] M. V. Ostrogradskiĭ. De l’intégration des fractions rationnelles. Bull. de la classe physicomathématique de l’Acad. Impériale des Sciences de Saint-Pétersbourg, 4:145–167, 286–300, 1845.

[72] P. Paule. Greatest factorial factorization and symbolic summation. J. Symbolic Comput., 20(3):235–268, 1995.

[73] G. H. Payne. Multivariate Hypergeometric Terms. PhD thesis, Pennsylvania State University, Pennsylvania, USA, 1997.

[74] R. Pemantle and M. C. Wilson. Twenty combinatorial examples of asymptotics derived from multivariate generating functions. SIAM Rev., 50(2):199–272, 2008.

[75] M. Petkovšek, H. S. Wilf, and D. Zeilberger. A = B. A K Peters Ltd., Wellesley, MA, 1996. With a foreword by Donald E. Knuth.

[76] R. Pirastu and V. Strehl. Rational summation and Gosper-Petkovšek representation. J. Symbolic Comput., 20(5-6):617–635, 1995. Symbolic computation in combinatorics (Ithaca, NY, 1993).

[77] R. H. Risch. The problem of integration in finite terms. Trans. Amer. Math. Soc., 139:167– 189, 1969.

[78] R. H. Risch. The solution of the problem of integration in finite terms. Bull. Amer. Math. Soc., 76:605–608, 1970.

[79] J. F. Ritt. Differential Algebra. American Mathematical Society Colloquium Publications, Vol. XXXIII. American Mathematical Society, New York, N. Y., 1950.

[80] M. Rothstein. A new algorithm for integration of exponential and logarithmic functions. In Proceedings of the 1977 MACSYMA Users Conference (Berkeley, CA), pages 263–274, Washington DC, 1977. NASA.

[81] M. Sato. Theory of prehomogeneous vector spaces (algebraic part)—the English translation of Sato’s lecture from Shintani’s note. Nagoya Math. J., 120:1–34, 1990. Notes by Takuro Shintani, Translated from the Japanese by Masakazu Muro.

[82] M. F. Singer. Liouvillian first integrals of differential equations. Trans. Amer. Math. Soc., 333(2):673–688, 1992.

[83] R. P. Stanley. Enumerative Combinatorics. Vol. 2, volume 62 of Cambridge Studies in Advanced Mathematics. Cambridge University Press, 1999.

[84] A. Storjohann and G. Villard. Computing the rank and a small nullspace basis of a polynomial matrix. In ISSAC’05: Proceedings of the 2005 International Symposium on Symbolic and Algebraic Computation, pages 309–316. ACM, New York, 2005.

[85] N. Takayama. An approach to the zero recognition problem by Buchberger algorithm. J. Symbolic Comput., 14(2-3):265–282, 1992.

[86] N. Takayama. An algorithm for finding recurrence relations of binomial sums and its complexity. J. Symbolic Comput., 20(5-6):637–651, 1995. Symbolic computation in combinatorics Δ_{1} (Ithaca, NY, 1993).

[87] B. M. Trager. Algebraic factoring and rational function integration. In SYMSAC’76: Proceedings of the Third ACM Symposium on Symbolic and Algebraic Computation, pages 219– 226, New York, NY, USA, 1976. ACM.

[88] M. van der Put and M. F. Singer. Galois Theory of Linear Differential Equations, volume 328 of Grundlehren der Mathematischen Wissenschaften. Springer-Verlag, Berlin, 2003.

[89] P. Verbaeten. The automatic construction of pure recurrence relations. SIGSAM Bull., 8(3):96–98, 1974.

[90] J. von zur Gathen and J. Gerhard. Modern Computer Algebra. Cambridge University Press, Cambridge, second edition, 2003.

[91] K. Wegschaider. Computer Generated Proofs of Binomial Multi-Sum Identities. Master’s thesis, RISC, J. Kepler University, May 1997.

[92] H. S. Wilf and D. Zeilberger. An algorithmic proof theory for hypergeometric (ordinary and “q”) multisum/integral identities. Invent. Math., 108(3):575–633, 1992.

[93] M. Wu. On Solutions of Linear Functional Systems and Factorization of Modules over Laurent-Ore Algebras. PhD thesis, Chinese Academy of Sciences and Université de Nice-Sophia Antipolis, Beijing, China and Nice, France, 2005.

[94] D. Y. Y. Yun. Fast algorithm for rational function integration. In Information processing 77 (Proc. IFIP Congr., Toronto, Ont., 1977), pages 493–498. IFIP Congr. Ser., Vol. 7. North-Holland, Amsterdam, 1977.

[95] D. Zeilberger. A holonomic systems approach to special functions identities. J. Comput. Appl. Math., 32:321–368, 1990.

[96] D. Zeilberger. The method of creative telescoping. J. Symbolic Comput., 11(3):195–204, 1991.

[97] H. Zoladek. The extended monodromy group and Liouvillian first integrals. J. Dynam. Control Systems, 4(1):1–28, 1998.

126

Index

Abramov–Petkovšek triple, 94 additive decomposition, 18

bidegree, 30

certificate, 10 compatible, 54 complete, 88 conjugate, 54 constant, 15 creative telescoping, 9

D-finite, 114 derivation, 15 difference field, 16 difference ring, 16 differential field, 15 differential Gosper form, 40 differential rational normal form, 57 differential ring, 15 differential-difference ring, 16 differential-reduced, 56

existence criteria, 100 fully integrable, 54

Geddes–Le–Li triple, 88

Hermite reduction, 17, 28

holonomic, 113

hyperexponential-hypergeometric, 51 hyperexponential-integrable, 88 hypergeometric-summable, 94

integrability conditions, 54

logarithmic part, 17

lucky, 32

minimal telescoper, 34

orthogonal, 16

Ostrogradsky–Horowitz system, 22

properness, 73 pump, 91

rational part, 17

residue, 22

ring of Ore polynomials, 16 Rothstein–Trager resultant, 23 shift rational normal form, 57 shift-reduced, 56 similar, 81 split, 81 squarefree, 56 squarefree factorization, 29 squarefree part, 29 standard representation, 86 strict DRNF, 57 strict SRNF, 57 telescoper, 10 Zeilberger’s method, 9

Curriculum Vitae

Personal Data • Name: Shaoshi Chen

• Gender: Male

• Nationality: Chinese

• Date and place of birth: July 17, 1983, Zhejiang Province, China

• Marital state: Married, no children Contact

• Email: schen@amss.ac.cn Education

2001-2005, Studies in mathematics at Jiangsu University, China

2005-2007, Studies in symbolic computation at Key Laboratory of Mathematics Mechanization, Chinese Academy of Sciences. Supervisor: Ziming Li 2007-2010, Doctorate studies in symbolic computation at Key Laboratory of Mathematics Mechanization (Chinese Academy of Sciences, China) and Algorithms Project (INRIA-Rocquencourt, France). Co-supervisors: Frédéric Chyzak and Ziming Li

Published Papers • A. Bostan, S. Chen, F. Chyzak, and Z. Li. Complexity of creative telescoping for bivariate

rational functions. In Stephen M. Watt, editor, ISSAC 2010 Proc. 2010 Internat. Symp. Symbolic Algebraic Comput., New York, N. Y., 2010. ACM Press pp. 203-210.

Posters

• S. Chen and Z. Li. On the structure of multivariate hyperexponential-hypergeometric functions. Poster at the conference DART IV (Beijing, China), 2010.

• A. Bostan, S. Chen, F. Chyzak, and Z. Li. Rational-function telescoping: blending creative telescoping with Hermite reduction. Poster at the conference ISSAC’09 (Seoul, South Korea), 2009.

• S. Chen and Z. Li. A speed-up of the Hermite reduction for rational functions. Poster at the conference ICMM’09 (Beijing, China), 2009.

Preprints

• S. Chen, F. Chyzak, R. Feng and Z. Li. The existence of telescopers for hyperexponentialhypergeometric functions. MM Research Preprints, KLMM, AMSS, Academia Sinica, Vol. 29, p. 239-267, July 2010

• S. Chen and Z. Li. A note on Ostrogradsky and Horowitz’s method. MM Research Preprints, KLMM, AMSS, Academia Sinica, Vol. 29, p. 36–47, July 2010.

• S. Chen and Z. Li. A multiplicative form of multivariate hyperexponential-hypergeometric functions. MM Research Preprints KLMM, AMSS, Academia Sinica, Vol. 29, p. 25–35, July 2010.

• S. Chen, R. Feng, Z. Li and H. Wang. An exercise on real elementary functions in the book Symbolic Integration I (second edition). MM Research Preprints KLMM, AMSS, Academia Sinica, Vol. 26, p. 115–125, January 2008.