Séminaire Lotharingien de Combinatoire 66 (2011), Article B66a

EXPLICIT FORMULA FOR THE GENERATING SERIES OF DIAGONAL 3D ROOK PATHS

ALIN BOSTAN^{†}, FRÉDÉRIC CHYZAK^{†}, MARK VAN HOEIJ^{‡}, AND LUCIEN PECH^{†}

Abstract. Let a_{n} denote the number of ways in which a chess rook can move from a corner cell to the opposite corner cell of an n × n × n three-dimensional chessboard, assuming that the piece moves closer to the goal cell at each step. We describe the computer-driven discovery and proof of the fact that the generating series G(x) = ∑ _{n≥0} a_{n}x^{n} admits the following explicit expression in terms of a Gaussian hypergeometric function:

G(x) = 1 + 6 · ∫_{0}^{x} \frac{_{2}F_{1}( \begin{matrix}1/3 & & 2/3 \\ & 2 \end{matrix} | \frac{27w(2−3w)}{(1−4w)^{3}} )}{(1 − 4w)(1 − 64w)} dw.

1. Introduction

In this article, we solve a problem of enumerative combinatorics addressed and left open in [16]. The initial question is formulated in terms of paths on an infinite three-dimensional chessboard. The 3D chessboard being identified with ℕ^{3}, a 3D rook is a piece which is allowed to move parallelly to one of the three axes. The general objective is to count paths (i.e., finite sequences of moves) of a 3D rook on the 3D chessboard. Following [16], we further restrict to 3D rook paths starting from the cell (0, 0, 0), whose steps are positive integer multiples of (1, 0, 0), (0, 1, 0), or (0, 0, 1). In other words, we assume that the piece moves closer to the goal cell at each step. For example, one such rook path is

(0, 0, 0) → (5, 0, 0) → (5, 0, 6) → (10, 0, 6) → (10, 2, 6).

Let r_{i,j,k} denote the number of rook paths from (0, 0, 0) to (i, j, k) ∈ ℕ^{3}, and let a_{n} = r_{n,n,n} be the number of diagonal rook paths from (0, 0, 0) to (n, n, n). The sequence (a_{n}) is called the 3D diagonal rook sequence; its first terms are

Key words and phrases. Enumerative combinatorics, generating functions, lattice paths, algebraic functions, hypergeometric functions, computer algebra, automated guessing, creative telescoping, fast algorithms.

† Supported in part by the Microsoft Research – INRIA Joint Centre.

‡ Supported by NSF 1017880.

1, 6, 222, 9918, 486924, 25267236, 1359631776, 75059524392, 4223303759148, …

This is Sequence A144045 in Sloane’s on-line encyclopedia of integer sequences at http://oeis.org/A144045. It was conjectured in [16] that the sequence (a_{n}) satisfies the fourth-order recurrence relation

(1) 2n^{2}(n − 1)a_{n} − (n − 1)(121n^{2} − 91n − 6)a_{n−1}

−(n−2)(475n^{2}−2512n+2829)a_{n−2} +18(n−3)(97n^{2}−519n+702)a_{n−3}

− 1152(n − 3)(n − 4)^{2}a_{n−4} = 0, for n ≥ 4.

The aim of the present article is not only to give a computer-driven proof of this conjecture, but also to show that the generating function G(x) = ∑ _{n≥0} a_{n}x^{n} ∈ ℚ[[x]] admits the following explicit expression:

(2)  G(x) = 1 + 6 · ∫_{0}^{x} \frac{_{2}F_{1}( \begin{matrix}1/3 & & 2/3 \\ & 2 \end{matrix} | \frac{27w(2−3w)}{(1−4w)^{3}} )}{(1 − 4w)(1 − 64w)} dw,

where the Gaussian hypergeometric series _{2}F_{1} with parameters a, b, c ∈ ℂ, −c ∉ ℕ, is defined by

2F_{1}( \begin{matrix}a b \\ c \end{matrix}| z) = ∑_{n=0}^{∞} \frac{(a)_{n}(b)_{n}}{(c)_{n}} \frac{z^{n}}{n!},

using the notation (a)_{n} for the Pochhammer symbol (a)_{n} = a(a + 1) ⋯ (a + n − 1).

The core of this article is Section 2, in which we complete the enumerative study of 3D rooks, proving in particular recurrence (1) and the explicit form (2). Section 3 then presents further comments, variant proofs and calculations, and discussions of a possible extension.

In more detail, Section 2 proceeds as follows. The question of the enumeration of a constrained combinatorial walk is re-expressed in Section 2.1 in terms of a diagonal of a rational series. Then, this diagonal is encoded as a residue in Section 2.2, leading to the problem of finding a special ODE it satisfies. In the literature, an effective approach to finding this special ODE was developed by Lipshitz. In Section 2.3, we show that, while the approach applies in theory, it leads to unfeasible calculations in practice. A faster approach, Zeilberger’s creative telescoping and its evolution by Chyzak, is summarised and applied in Section 2.4. The obtained ODE allows us to prove the conjectured recurrence (1) and even a lower-order one, (15), in Section 2.5. We then derive the _{2}F_{1} explicit form in Section 2.6.

When possible, we have attempted to describe the mathematical phenomena, ideas, and techniques on a general level in Sections 2.1– 2.6. When the complete proofs of the results we state require a more involved calculation, and also when they are based on more involved notions or on the fine tuning of algorithms, we have postponed some of the technicalities to Appendix A.

2. Main steps of the solution

2.1. From combinatorics to algebra. The combinatorial problem translates into a question on generating series, best expressed in terms of diagonals: given ϕ(s, t, u) = ∑ _{i,j,k≥0} c_{i,j,k} s^{i}t^{j}u^{k} in ℚ[[s, t, u]], its diagonal is defined as the univariate power series Diag(ϕ) = ∑ _{n≥0} c_{n,n,n} x^{n} in ℚ[[x]].

Question. Define the trivariate rational function f(s, t, u) in ℚ(s, t, u) as

( 1 − ∑_{n≥1} s^{n} − ∑_{n≥1} t^{n} − ∑_{n≥1} u^{n} ) ^{−1} =

\frac{(1 − s)(1 − t)(1 − u)}{1 − 2(s + t + u) + 3(st + tu + us) − 4stu} ,

which is seen to be at the same time in ℚ[[s, t, u]]. What can be said about its diagonal Diag(f) = ∑ _{n≥0} a_{n}x^{n}?

This translation is based on the following straightforward extension of Proposition 6.3.7 in [31], which we shall use with the choice 𝔖 = { (n, 0, 0) : n ≥ 1 } ∪ { (0, n, 0) : n ≥ 1 } ∪ { (0, 0, n) : n ≥ 1 }.

Lemma 1. Given a finite set of allowed steps 𝔖 ⊆ ℕ^{3} \ {(0, 0, 0)}, let N_{𝔖}(m, n, p) denote the number of paths from (0, 0, 0) to (m, n, p), with path steps taken from 𝔖. Then, the generating series

G_{𝔖}(s, t, u) = ∑_{m,n,p≥0} N_{𝔖}(m, n, p)s^{m}t^{n}u^{p}

is rational, and has the form

G_{𝔖}(s, t, u) = \frac{1}{1 − ∑ _{(i,j,k)∈𝔖} s^{i}t^{j}u^{k}} .

The same result holds for infinite step sets 𝔖 provided the power series ∑ _{(i,j,k)∈𝔖} s^{i}t^{j}u^{k} is rational.

2.2. The key equation. A first, theoretical answer to our question comes from the main result in Lipshitz’s article [24]: the diagonal Diag(f) is a D-finite series, i.e., it satisfies a linear differential equation in ∂_{x} = d/dx with polynomial coefficients in x. Therefore, the sequence (a_{n}) of its coefficients satisfies a linear recurrence with polynomial coefficients in n.

Moreover, let F be the trivariate rational function

F = \frac{1}{st} · f(s, t/s, x/t) =

\frac{(s − t)(s − 1)(t − x)}{st(−st + 2s^{2}t + 2t^{2} + 2xs − 3st^{2} − 3xt − 3xs^{2} + 4xst)} .

Then, Remark 3 in [24] shows that

Lemma 2. If the equation

(3)  P (F) = \frac{∂S}{∂s} + \frac{∂T}{∂t}

admits a solution (P, S, T) where P (x, ∂_{x}) is a non-zero linear differential operator with polynomial coefficients in ℚ[x], and S and T are two rational functions in ℚ(x, s, t), then P (x, ∂_{x}) annihilates Diag(f).

Lemma 2 is the heart of the result in [24], and the basis of all methods for effectively computing differential equations satisfied by diagonals. For completeness, we give here a self-contained proof, which closely follows [24].

Proof. Although F is not a formal (Laurent) power series in x, s, t, it does come from a formal power series after a change of variables, and thus it can be seen as an element of the differential ℚ[x, s, t]-module of all formal sums of the form

B = _{∑_{α∈ℕ,β∈ℤ,γ∈ℤ} ; α+β≥−μ,β+γ≥−ν} c_{α,β,γ} s^{α}t^{β}x^{γ}

for some (μ, ν) ∈ ℕ^{2} depending on B. As a consequence, it makes sense to speak of the series expansion of F and of P (F), but in later calculations, these series may be multiplied solely by elements of ℚ[x, s, t].

The diagonal Diag(f) is, by construction, equal to the coefficient of s^{−1} t^{−1} in F. On the other hand, equation (3) implies that the coefficient of s^{−1} t^{−1} in P (F) is zero. The conclusion then follows from the sequence of equalities

P (Diag(f)) = P ([s^{−1} t^{−1} ]F) = [s^{−1} t^{−1} ]P (F) = 0,

the second equality being a consequence of the fact that P is free of ∂_{s} and ∂_{t}. □

The question remaining is how to determine a non-trivial triple (P, S, T) solution of (3).

2.3. Lipshitz’s approach. Lipshitz [24, Lemma 3] uses a counting argument, similar to Fasenmyer’s [17, 39], to show that there exists a non-zero linear differential operator L(x, ∂_{x}, ∂_{s}, ∂_{t}), free of s and t, that annihilates F. Then, writing this operator under the form

L(x, ∂_{x}, ∂_{s}, ∂_{t}) = P (x, ∂_{x})+ (higher-order terms in ∂_{s} and ∂_{t}, free of s and t),

it follows from Lemma 2 that P annihilates Diag(f).

Let us sketch Lipshitz’s argument on our concrete example; as we shall see, while this proof is algorithmic and ultimately reduces the existence of P to a linear-algebra argument, it produces unreasonably large systems. The key observation is that, by Leibniz’s rule, for any nonnegative integer N, the \binom{N+4}{4} rational functions

x^{i}∂_{x}^{j}∂_{s}^{k}∂_{t}^{ℓ}(F), 0 ≤ i + j + k + ℓ ≤ N,

belong to the ℚ-vector space of dimension at most 18(N + 1)^{3} spanned by the elements

(4) \frac{x^{i}s^{j}t^{k}}{q^{N+1}} , 0 ≤ i ≤ 2N + 1, 0 ≤ j ≤ 3N + 2, 0 ≤ k ≤ 3N + 2,

where q is the denominator st(−st+2s^{2}t+2t^{2}+2xs−3st^{2}−3xt−3xs^{2}+ 4xst) of F. Thus, if N is an integer such that \binom{N+4}{4} > 18(N + 1)^{3}, then there exists a non-zero linear differential operator L(x, ∂_{x}, ∂_{s}, ∂_{t}) of total degree at most N in x, ∂_{x}, ∂_{s}, and ∂_{t}, such that L(F) = 0.

Now, it appears that N = 425 is the smallest integer satisfying the inequality \binom{N+4}{4} > 18(N + 1)^{3}; therefore, explicitly finding the operator L would require solving a homogeneous linear system with some 1,391,641,251 unknowns and some 1,391,557,968 equations!

It is possible to refine the previous argument in various ways, and thus to reduce the search of L to a linear algebra problem of smaller size. Such methods will be discussed and compared thoroughly in a later work. In our case, the best reduction in size that we could achieve uses a variation of the main idea in [9, §3.2].

First, as before, Leibniz’s rule shows that the u_{N} := \binom{N+3}{3} rational functions

∂_{x}^{j}∂_{s}^{k}∂_{t}^{ℓ}(F), 0 ≤ j + k + ℓ ≤ N,

belong to the ℚ(x)-vector space spanned by the elements

(5)  \frac{s^{j}t^{k}}{q^{N+1}} , 0 ≤ j ≤ 3N + 2, 0 ≤ k ≤ 3N + 2,

whose dimension d_{N} is at most 9(N + 1)^{2}.

A further idea is to exploit the sparse feature of the monomial supports of the numerator and denominator of F. By studying the behaviour of these supports under differentiations with respect to x, s, and t, one shows that the terms (5) will not be involved for the \binom{2N+2}{2} + \binom{N+2}{2} pairs (j, k) ∈ ℕ^{2} in the set

{ (j, k) : j + k ≤ 2N } ∪ { (j, k) : 2N + 2 ≤ j, k ≤ 3N + 2, j + k ≥ 5N + 4 }.

This implies that d_{N} is bounded above by e_{N} := \frac{1}{2}(13N +14)(N +1), and shows that the search of L boils down to a homogeneous linear system over ℚ(x), with u_{N} unknowns and e_{N} equations. The smallest value of N such that u_{N} > e_{N} is N = 36, leading to the task of computing the (right) kernel of a matrix of size 8917 × 9139, with entries in ℚ[x] of degree at most 37.

Despite the spectacular reduction in dimension compared to Lipshitz’s original argument, these sizes are still too high to be dealt with in practice. Of course, the value e_{N} is only an upper bound on d_{N}, leaving room for a magical rank drop occurring before the predicted value N = 36. That being said, our implementation of the method constructs the linear system for increasing values of N, but does not find any non-trivial kernel for N ≤ 18, and is not able to finish the computation for N = 19 before consuming all the available memory.^{1}

The conclusion is that Lipshitz’s approach and its improvements are not sufficient to obtain effectively an equation for the diagonal of f(x, s, t).

2.4. The creative-telescoping approach. An alternative approach to solving (3) was popularised among the computer-algebra community by Zeilberger under the name of creative telescoping, and it is the one that we shall use to effectively construct an equation for our diagonal. This method, first introduced for hypergeometric summation [40, 42] and hyperexponential integration [3], was later generalised by Chyzak [11] to more general integrands involving a larger class of special functions called holonomic. Zeilberger’s definition of this class [41, Section 2.2.4] has its roots in D-module theory [5, 6] and is not needed to understand the present paper. See also [12] for recent extensions of creative telescoping to more general special functions.

1Our implementation of an optimisation of Lipshitz’s algorithm, available at http://algo.inria.fr/chyzak/Rooks/OptimisedLipshitz.mpl, is nevertheless able to predict that, with high probability, a non-trivial kernel does exist for N = 19.

Prior to summarising the method, we show its behaviour and its result on our concrete case, postponing the explicit details of intermediate steps to Appendix A.

2.4.1. Obtaining the result. In our case, Chyzak’s algorithm, as implemented by the package Mgfun^{2} for Maple, produces a solution (P, S, T) of (3), with P = P_{2}∂_{x}, where

(6) P_{2} = x(x − 1)(64x − 1)(3x − 2)(6x + 1)∂_{x}^{2}

+ (4608x^{4} − 6372x^{3} + 813x^{2} + 514x − 4)∂_{x}
+ 4(576x^{3} − 801x^{2} − 108x + 74).

Let q = st · q_{1} with q_{1} irreducible in ℚ[x, s, t] of degree (1, 2, 2) in (x, s, t). Then, S and T are of the very special forms:

S = \frac{(s − t) · U}{2st · q_{1}^{2} · disc_{t}(q_{1})} , T = \frac{(s − t) · V}{2s^{2} · q_{1}^{3} · disc_{t}(q_{1})^{2}},

where U and V are (irreducible) polynomials in ℚ[x, s, t] of respective degrees (5, 8, 3) and (8, 14, 5) in (x, s, t). Here, disc_{t}(q_{1}) denotes the discriminant of q_{1} with respect to t. Observe the very moderate degrees in comparison to those pessimistically predicted by Lipshitz’s method.

Maple commands to produce P, S, and T are provided online at

http://algo.inria.fr/chyzak/Rooks/CreativeTelescoping.mpl ;

they take less than 10 seconds on an average laptop. The rational functions S and T are not displayed here, but can be found online at

http://algo.inria.fr/chyzak/Rooks/Certificates.mpl .

2.4.2. Sketch of the algorithm. We proceed to sketch the algorithm of [11] for the creative telescoping of multiple “integrals.” By this, we actually mean in this article answering the algebraic question of how to solve the key equation, whether in the form (7) below for “simple integrals,” or in the form (3) for “double integrals.”

Presenting Chyzak’s algorithm directly in its generality would prove too arid, so we go instead progressively from the simpler case of rational simple integrals to the integration of more general functions, then to double integrals. In effect, this section is mostly bibliographic, with the exception that, in the case of multiple integrals, we provide an extension of the algorithm in [11] beyond the so-called “natural boundaries.”

2version 4.0, available as part of Algolib 13.0 from http://algo.inria.fr/libraries/. Initially developed by F. Chyzak. Specific version used here reimplemented by L. Pech.

For a single integral of a bivariate rational function F (u, v), the creative-telescoping method finds a linear differential operator P (u, ∂_{u}) and a rational function S in ℚ(u, v) such that (cf. equation (3))

(7)  P (F) = \frac{∂S}{∂v}.

To this end, the algorithm in [3] writes P and S in undetermined form

(8)  P = η_{r}(u)∂_{u}^{r} + ⋯+ η_{0}(u) and S = ϕ(u, v)F

for some tentative order r and unknown rational functions η_{i} from ℚ(u) and ϕ from ℚ(u, v). To get ϕ, the original algorithm of [3] exploits the specific nature of F (hyperexponential). We prefer describing here the less optimal variant of [11], as it will extend to the more general class of inputs we need: (7) is rewritten into an auxiliary ordinary linear differential equation on ϕ, before one uses known algorithms for solving a linear differential equation with polynomial coefficients for its rational solutions. More precisely, such algorithms, like the one in [1], have to be refined to solve for the η_{i} as well; this is easily taken into account just by having more unknowns in the crucial step of linear algebra over ℚ(u) in this rational solving.

For reasons to be given shortly, the double integration of trivariate rational functions is based on integration of more general functions F, leading to a more general form for S in (7). Specifically, the case dealt with in [11] is that of functions F for which there exists a finite maximal set of derivatives ∂_{u}^{a}∂_{v}^{b}(F) that are linearly independent over ℚ(u, v). Such functions are known as differentiably finite in the literature, or in short D-finite [30, 25]. However, note that the notions of holonomic and differentiably finite functions, which look technically very different, cover in fact exactly the same sets of functions. This fact is explained in elementary terms in [32, Appendix]. Still, the more elaborate holonomic point of view is needed in proofs that require a subtler counting argument.

Once such a set of pairs (a, b) is fixed for a given D-finite function F, equation (7) is changed by setting S to the undetermined form

(9)  S = ∑_{(a,b)}ϕ_{a,b}(u, v)∂_{u}^{a}∂_{v}^{b}(F)

for unknown rational functions ϕ_{a,b} from ℚ(u, v). The only change then is that the auxiliary equation in ϕ becomes an auxiliary ordinary linear differential system in the ϕ_{a,b}’s. Solving is performed either by uncoupling and using the algorithm mentioned previously several times, or by a direct approach as in [4].

The case of double integrals leads to generalising (7) into (3), but the solving does not generalise so smoothly: attempting to write S = ϕ_{1}(s, t, x)F and T = ϕ_{2}(s, t, x)F for undetermined rational functions ϕ_{1} and ϕ_{2} from ℚ(s, t, x) yields a linear partial differential equation relating these functions with ∂ϕ_{1}/∂s and ∂ϕ_{2}/∂t. To the best of our knowledge, although this overdetermined equation has a very specific form, no algorithm is available to solve it for its rational solutions.

Therefore, instead of a direct approach, Chyzak developed in [11] a cascading approach which we now summarise. Noting that the dependency of P on a single derivation ∂_{u} in (7) is inessential, the same approach is possible for the creative telescoping with respect to the (single) variable v of a trivariate rational function F from ℚ(u_{1}, u_{2}, v). Indeed, setting P to the undetermined form

(10)  P = ∑_{0≤i+j≤r} η_{i,j} (u_{1}, u_{2})∂_{u_{1}}^{i}∂_{u_{2}}^{j}

for some tentative total order r and unknown rational functions η_{i,j} from ℚ(u_{1}, u_{2}) and performing the same solving as previously, now relying on linear algebra over ℚ(u_{1}, u_{2}), leads to a basis of P^{(α)}’s of total order at most r for which there exists a rational function ϕ^{(α)}(u_{1}, u_{2}, v) satisfying

P^{(α)}(F) = \frac{∂}{∂v} (ϕ^{(α)}F ) .

The theory (as developed, e.g., by Zeilberger in [41]) guarantees that the function F is D-finite and that the set of P^{(α)}’s obtained for sufficiently large r can be used to determine the finite set needed as an input to the algorithm of [11] and in (9).

Finally, a double integration algorithm is obtained by continuing the approach used for natural boundaries in [11] (Stages A and B below) by a suitable recombination of the outputs (Stage C below). The resulting treatment of multiple integrals over non-natural boundaries is an extension over [11], and the corresponding algorithm is as follows:

• Stage A: First iteration of creative telescoping. Using the univariate algorithm for trivariate rational functions and variables (s, x, t) in place of (u_{1}, u_{2}, v) delivers identities

(11)  P^{(α)}(s, x, ∂_{s}, ∂_{x})(F) = \frac{∂}{∂t} (ϕ^{(α)}(s, t, x)F ) .

• Stage B: Second iteration of creative telescoping. Considering a function F̂ of (s, x) that is annihilated by all P^{(α)} and using the univariate algorithm for general functions and variables (s, x) in place of (u, v) delivers an identity

(12)  P (x, ∂_{x})(F̂) = \frac{∂}{∂s} (Q(s, x, ∂_{s}, ∂_{x})(F̂)).

• Stage C: Recombination. By the theory of linear-differential-operators ideals,^{3} the calculations of the algorithm can be interpreted as a proof of existence of operators L^{(α)}(s, x, ∂_{s}, ∂_{x}) satisfying

(13) P (x, ∂_{x})−∂_{s}Q(s, x, ∂_{s}, ∂_{x}) = ∑_{α} L^{(α)}(s, x, ∂_{s}, ∂_{x})P^{(α)}(s, x, ∂_{s}, ∂_{x}).

These L^{(α)} can effectively be obtained either by following the calculations step by step or (less efficiently) by a postprocessing (non-commutative multivariate division). Hence, defining

(14)

S = Q(s, x, ∂_{s}, ∂_{x})(F) and T = ∑_{α} L^{(α)}(s, x, ∂_{s}, ∂_{x})(ϕ^{(α)}(s, t, x)F )

leads to a solution (P, S, T) of (3). Note that this two-stage process inherently introduces a dissymmetry

in the treatment of the variables s and t: the output from the first iteration tends to be larger than its input; in turn, the output from the second is larger than the output from the first. As a consequence, the order we deal with the variables may have an impact on the running time.

2.5. Proving the recurrence. Translating the differential equation P (Diag(f)) = 0 into a recurrence satisfied by the coefficients a_{n} of Diag(f) yields the fourth-order recurrence relation (1).

It is possible to prove that the sequence (a_{n}) satisfies the even shorter recurrence

(15) 2(n−1)(35n−52)n^{2}a_{n}−(n−1)(4655n^{3}−11781n^{2}+8494n−1776)a_{n−1}

+ (n − 2)(11305n^{3} − 41856n^{2} + 46487n − 13128)a_{n−2}

− 192(n − 3)^{2}(35n − 17)(n − 2)a_{n−3} = 0, for n ≥ 3.

Recurrence (15) can be guessed (e.g., by linear algebra) using, say, the first 25 terms of the sequence (a_{n}). Once it is conjectured, proving its correctness is easy: if (b_{n}) denotes the sequence defined by the lefthand side of (15), then recurrence (1) implies that (b_{n}) satisfies the recurrence

3This observation was made during a discussion between F. Chyzak, M. Kauers, and Ch. Koutschan on July 17, 2008, in the “UFO” room at RISC (Linz, Austria).

b_{n} + 6b_{n−1} = (35n − 52)( 2(n − 1)n^{2}a_{n} − (n − 1)(121n^{2} − 91n − 6)a_{n−1}

−(n−2)(475n^{2}−2512n+2829)a_{n−2} +18(n−3)(97n^{2}−519n+702)a_{n−3}

− 1152(n − 3)(n − 4)^{2}a_{n−4} ) = 0, for n ≥ 3,

which, combined with the initial condition b_{0} = −54864 · 6 + 2 · 43362 · 222 − 2 · 2 · 53 · 3^{2} · 9918 = 0, yields b_{n} = 0 for all n ≥ 0.

Maple code to obtain the first terms of the sequence and derive the recurrences can be obtained at

http://algo.inria.fr/chyzak/Rooks/Recurrence.mpl.

2.6. Solution in explicit form. Operator P_{2} from equation (6) has a power series solution with integer coefficients, namely ∂_{x}(Diag(f)). By a speculation of Dwork’s [15], it appears plausible that P_{2} has a 2F_{1}-type solution. We shall briefly sketch how, for a second-order operator L, one might find such solutions. A similar problem (for Bessel instead of _{2}F_{1}) was solved in [14, 36], and many of the details discussed in those papers are needed for the _{2}F_{1}-case as well (in particular the first part of Section 3 in [14]). The result of our computation is:

(16) ∂_{x}(Diag(f)) = \frac{6}{(1 − 4x)(1 − 64x)} ·_{2}F_{1}( \begin{matrix}1/3 & & 2/3 \\ & 2 \end{matrix} | \frac{27x(2 − 3x)}{(1 − 4x)^{3}} ) .

2.6.1. Formal asymptotic-series solutions of linear differential operators. We recall the necessary facts from the classical theory of asymptotic-series solutions of linear differential operators [22, 38], and specialise them here to our need: the order 2.

With respect to a given linear differential operator L of order r with polynomial coefficients, each point p from ℂ ∪ {∞} falls in one of two categories:

• it is called an ordinary point of L if the coefficient of the highest derivative in the equation does not vanish at p, in which case a solution f is fully determined by the initial conditions f(p), …, f^{(r−1)} (p) (Cauchy theorem);

• it is called a singular point of L otherwise. Among singular points, the more tractable situation is that of a reg-

ular singular point (of L). At such a point, the growth of solutions is bounded (in any small sector) by an algebraic function. In the neighbourhood of a regular singular point p, any solution can be written as a linear combination of products of a (possibly multivalued) term of the form (x−p)^{e} with e ∈ ℂ, logarithms ln(x−p), (respectively (1/x)^{e} and ln(1/x) if p = ∞) and analytic (univalued) functions. The constants e are called the exponents of L at x = p. They are the zeroes of the so-called indicial equation of L at x = p, an algebraic equation that is defined for any L and p, and can be obtained algorithmically.

Another simple situation occurs at a singular point p of L when there exist a function u (possibly singular at p) and functions y_{1}, …, y_{r} that are analytic in a neighbourhood of p, such that uy_{1}, …, uy_{r} constitute a basis of the solution space of L. Then, we say that p is a removable singularity of L. Otherwise, we say that it is a non-removable singularity.

Assume from now on that L is irreducible and of order 2. Then, the indicial equation at x = p has degree at most 2; in other words, there are at most 2 exponents at x = p (counting with multiplicity). We shall also assume that L is regular singular at x = p, which implies that the degree of the indicial equation is equal to 2 and that there are two exponents, e_{1} and e_{2} (counting with multiplicity). The exponent difference d at x = p is defined up to sign as ±(e_{1} − e_{2}).

The numbers e_{1}, e_{2}, and d describe a lot of the asymptotic behaviour of the solutions of L when x tends to p. Let us be more specific about this. For notational convenience, consider first the case p = 0. Then e is an exponent at x = 0 if and only if L has a formal solution u for which u/x^{e} ∈ R \ xR, where R := ℂ[[x]][ln(x)]. If ln(x) actually occurs in u/x^{e} then the point p is called logarithmic. If ln(x) does not occur, then our assumption that L is regular singular can be used to prove that u/x^{e} is analytic at x = 0, so that x^{e} describes the local behaviour of u at x = 0. Similar statements hold for p ≠ 0, after replacing x by x − p if p is finite, or x by 1/x if p = ∞. The following lemma summarises the link between the exponent difference of a second-order operator and the asymptotic behaviour of its solutions.

Lemma 3. Let L be a linear difference operator of order 2 with exponent difference d at some point p ∈ ℂ ∪ {∞}. Assume p to be either an ordinary or a regular singular point of L. Then:

(1) If p is ordinary, then d = 1.

(2) If p is a removable singularity, then d ∈ ℤ \ {0}.

(3) If p is logarithmic, then d ∈ ℤ.

(4) If d = 0, then p is logarithmic.

(5) If d ∈ ℤ \ {0}, then p is either logarithmic or a removable singularity.

For proofs, see [14, Theorem 1], or see [13] for more details.

2.6.2. Transporting the singularities of the Gauss hypergeometric equation. Let L_{c}^{a,b} denote the differential operator of the Gauss hypergeometric equation,

x(1 − x)\frac{d^{2}w}{dx^{2}} + (c − (a + b + 1)x) \frac{dw}{dx} − abw = 0.

Of its two linearly independent classical solutions, we shall only need one: _{2}F_{1}(a, b; c; x). It has 3 singular points: 0, 1, and ∞. The respective exponent differences are (e_{0}, e_{1}, e_{∞} ) = (±(1−c), ±(c−a−b), ±(a−b)).

Let f be a rational function, and denote by L_{c;f}^{a,b} the operator obtained from L_{c}^{a,b} by the change of variables x = f, so that it has _{2}F_{1}(a, b; c; f) as a solution.

The following lemma explains how the exponent differences of the Gauss differential equation get transported by the composition under f. The proof is the same as in [14, Theorem 2(i)] (more details are given in [13]).

Lemma 4. Let d be the exponent difference of L_{c;f}^{a,b} at x = q. Then: (1) If q is a root of f(x) with multiplicity m, then d = m · e_{0}.

(2) If q is a root of f(x) − 1 with multiplicity m, then d = m · e_{1}.

(3) If q is a pole of f(x) of order m, then d = m · e_{∞} .

Recall that L is irreducible and of order 2. Let y = _{2}F_{1}(a, b; c; f) where f is some rational function. By solving L in terms of y, we mean finding three rational functions r, r_{0}, and r_{1} such that

(17)  exp( ∫ r) · (r_{0}y + r_{1}y′)

is a solution of L (see [14, Definition 2, parts (ii),(iii)]). A necessary condition for such a solution to exist is that L has the same set of non-removable singularities as L_{c;f}^{a,b}, and moreover, that the exponent differences match modulo ℤ (see [14, Lemma 4]). The expression exp(∫ r) is usually a radical function, i.e., a function whose rth power is rational for some integer r.

Consider the operator P_{2} in (6). Its non-removable singularities are 0, 1, 1/64, 2/3, and ∞, and all these singular points are logarithmic; on the other hand, the root x = −1/6 of the leading coefficient of P_{2} is a removable singularity.

Suppose that we choose a, b, c so that (e_{0}, e_{1}, e_{∞} ) = (0, 1/3, 0) (this choice will be explained in Section 2.6.3), and want to solve P_{2} in terms of y := _{2}F_{1}(a, b; c; f). By the discussion above, if p is a logarithmic point for P_{2} then it must be a logarithmic point for L_{c;f}^{a,b} and vice versa. Since x = 0 and x = ∞ are the logarithmic points of L_{c}^{a,b} one sees that the union of the set of roots and the set of poles of f must be exactly {0, 1, 1/64, 2/3, ∞}. For each possible degree of f, this leaves finitely many candidate integer values for n_{0}, n_{1}, n_{1/64}, and n_{2/3}, such that f is of the form c(x − 0)^{n_{0}}(x − 1)^{n_{1}}(x − 1/64)^{n_{1/64}} (x − 2/3)^{n_{2/3}} for an unknown constant c. For an exhaustive search, we therefore try successive potential degrees 2, 3, 4, … for f.

The exponent difference of L_{c;f}^{a,b} at x = q is m·e_{1} if q is a root of f −1 with multiplicity m. So m · e_{1} must match (modulo ℤ) an exponent difference of P_{2}. Since P_{2} has only integer exponent differences, and e_{1} = 1/3, it follows that 3 divides m and so the numerator of f −1 must be a cube. Thus, we try various n_{0}, n_{1}, … ∈ ℤ, and for each, we check if there exists c ∈ ℂ for which the numerator of c(x − 0)^{n_{0}}(x − 1)^{n_{1}}(x − 1/64)^{n_{1/64}} (x−2/3)^{n_{2/3}} −1 is a cube. A computer implementation quickly finds

n_{0} = 1, n_{1} = −2, n_{2/3} = 1, n_{1/64} = −1, c = −81/64,

and we get

f = \frac{−81x(x − 2/3)}{64(x − 1)^{2}(x − 1/64)} .

Next, we compute a, b, c for which the exponent differences e_{0}, e_{1}, e_{∞} are 0, 1/3, 0 and find L_{c}^{a,b}. Applying a change of variables produces L_{c;f}^{a,b}. Then a call to the program equiv [34] provides a map of the form y ↦ exp(∫ r)(r_{0}y + r_{1}y′) from the solutions of L_{c;f}^{a,b} to the solutions of P_{2}. This way we find a solution of P_{2} in the form (17).

2.6.3. Simplifying the obtained solution. The question now is how to make this solution more compact. First of all, the numerator of f −1 is a cube. Interchanging e_{1} and e_{∞} makes the denominator of f a cube. That makes the expression for f slightly smaller. Now, composition with x/(x − 1) interchanges 1, ∞ while keeping 0 fixed, so we replace f by f/(f − 1) and obtain a new f,

f = \frac{27x(2 − 3x)}{(1 − 4x)^{3}} ,

which corresponds to the new (e_{0}, e_{1}, e_{∞} ) = (0, 0, 1/3). Next, we try to change e_{0}, e_{1}, e_{∞} by some integers in such a way that the exponent differences of L_{c;f}^{a,b} not only match those of P_{2} modulo ℤ, but are also as close as possible, because this tends to reduce the size of r_{0}, r_{1}. Surprisingly, one can even match the exponent differences exactly! That allows r_{1} to become 0. We get an exact match for (e_{0}, e_{1}, e_{∞} ) = (1, 1, 1/3), corresponding to a = 1/3, b = 2/3, c = 2, and end up with the following solution

\frac{6}{(1 − 4x)(1 − 64x)} ·_{2}F_{1}( \begin{matrix}1/3 & & 2/3 \\ & 2 \end{matrix} | \frac{27x(2 − 3x)}{(1 − 4x)^{3}} ) .

This expression is not unique because we can interchange (e_{0}, e_{1}, e_{∞} ), updating f accordingly. We can also change e_{0}, e_{1}, e_{∞} by integers (as long as those integers add up to an even number) but doing that would make r_{1} non-zero again. Moreover, there are other exponent differences that would have solved P_{2} as well.

The set of (e_{0}, e_{1}, e_{∞} )’s (modulo ℤ) we should try is given in [33, Table (I), Section 4]. To make the expression size as small as possible, the (e_{0}, e_{1}, e_{∞} ) that leads to the lowest degree for f should be tried first. The table shows how the degrees are related to each other; to minimise the degree, we should check (0, 0, 0) first, and (0, 1/2, 1/3) last. The first case (0, 0, 0) did not lead to a solution. The second case we tried was (0, 0, 1/3), and since that led to a solution, we did not need to try more cases. The table in [33] shows some of the other cases must also lead to a solution; one sees in the table that (0, 0, 1/3) (and hence P_{2}) is solvable in terms of (0, 1/2, 1/6) and (0, 1/2, 1/3). We can also read from the table that if we had used those exponent differences for our _{2}F_{1}, then the degree of f would have been 2, respectively 4, times higher.

A Maple session to perform all operations of this section is given

online at http://algo.inria.fr/chyzak/Rooks/Compute2F1.mpl.

3. Further comments

3.1. Asymptotics. It was pointed out in [16] that general results from [27] on the asymptotics of diagonal coefficients of multivariate generating functions imply

(18)  a_{n} ∼ \frac{9\sqrt{3}}{40π} · \frac{64^{n}}{n}.

Here we show that the same result can be derived in a simpler way (using univariate, instead of multivariate, singularity analysis), based on the explicit expression (16).

Indeed, basic singularity analysis [19, 20] shows that, since the dominant singularity of ∂_{x}(Diag(f)) is x = 1/64, with residue

r = \frac{6}{(1 − \frac{4}{64})} ·_{2}F_{1}( \begin{matrix}1/3 & & 2/3 \\ & 2 \end{matrix} | \frac{\frac{27}{64}(2 − \frac{3}{64})}{(1 − \frac{4}{64})^{3}} ) = \frac{32}{5} ·_{2}F_{1}( \begin{matrix}1/3 & & 2/3 \\ & 2 \end{matrix} | 1) ,

the asymptotic behaviour of the coefficient h_{n} of ∂_{x}(Diag(f)) is h_{n} ∼

r · 64^{n}, and thus a_{n} ∼ \frac{r}{64} · \frac{64^{n}}{n}.

The value of r is found using Gauss’s formula

2F_{1}( \begin{matrix}a b \\ c \end{matrix}| 1) = \frac{Γ(c)Γ(c − a − b)}{Γ(c − a)Γ(c − b)}

and the multiplication formula

Γ(z) · Γ ( z + \frac{1}{m}) ⋯ Γ ( z + \frac{m − 1}{m} ) = (2π)^{(m−1)/2} m^{\frac{1}{2} −mz} Γ(mz)

with m = 3 and z = 1. This gives r = 72\sqrt{3}/(5π), and proves (18).

3.2. Non-algebraicity and simpler formulas. One may wonder whether the function Diag(f) is algebraic. This is not the case, and can be seen in several ways: (i) as a consequence of the asymptotic behaviour of the coefficients a_{n} (generating series of sequences asymptotically equivalent to n^{α}ρ^{n} with α ∈ ℤ_{−} are necessarily transcendental [18, Thm D]); (ii) by doing a local analysis of P near its singularities (the presence of logarithmic terms in a local expansion excludes algebraicity); (iii) by applying Schwarz’s classification [28] of algebraic _{2}F_{1}’s.

It is also legitimate to seek for an even simpler closed form, not involving any integration sign. Such a reduction is not possible: if it were, the algorithm given in [2] should have found it. The sequence (a_{n}) satisfies the third-order recurrence (15), but no second-order recurrence with polynomial coefficients, for otherwise (15) or a closely related adjoint recurrence would need to possess a hypergeometric solution, which is easily sorted out by Petkovšek’s algorithm [26].

3.3. An alternative hypergeometric expression. Using a completely different method, Frits Beukers [8] independently found that the diagonal 3D rook path series G(x) satisfies

(19)  G′(x) = \frac{1 − x}{2(1 + 6x)} ( (1 − 4x)H′(x) − 4H(x)) ,

where

H(x) = \frac{1}{g_{2}^{1/4}} ·_{2}F_{1}( \begin{matrix}1/12 & & 5/12 \\ & 1 \end{matrix}| \frac{1}{J(x)}) ,

with

g_{2} = (1 − 4x)(1 − 60x + 120x^{2} − 64x^{3}),

J(x) = \frac{g_{2}^{3}}{1728(1 − x)^{2}x^{3}(2 − 3x)^{3}(1 − 64x)} .

The function J(x) is obtained, by a method similar to the one in [7], as the J-invariant of a family of elliptic curves naturally associated to the integral interpretation of the diagonal problem.

Of course, once obtained, it is an easy task to check that the new hypergeometric expression coincides with ours: the difference of the two hypergeometric expressions satisfies a simple linear differential equation derived from the defining linear differential equation for the hypergeometric _{2}F_{1}, and it is easily expanded as a Taylor series, showing that enough first terms are zero, so that the difference is zero. This calculation is available as Maple code at

http://algo.inria.fr/chyzak/Rooks/CompareFormulas.mpl .

Interestingly, one can algebraically reduce one expression to the other by combining only two classical hypergeometric identities:

(a) Gauss’ contiguity relation

(20)  \frac{9(1 − x)}{2} ·_{2}F_{1}( \begin{matrix}1/3 & & 2/3 \\ & 1 \end{matrix} | x) ′= _{2}F_{1}( \begin{matrix}1/3 & & 2/3 \\ & 2 \end{matrix} | x) ,

obtained by evaluating at a = \frac{1}{3}, b = \frac{2}{3}, and c = 1 the more general contiguous transformation

c(1 − x) · \frac{d}{dx}_{2}F_{1}( \begin{matrix}a b \\ c \end{matrix}| x) =

(c − a)(c − b) ·_{2}F_{1}( \begin{matrix}a b \\ c + 1 \end{matrix} | x) + c(a + b − c) ·_{2}F_{1}( \begin{matrix}a b \\ c \end{matrix}| x) ,

(b) Goursat’s quartic relation (21)

2F_{1}( \begin{matrix}1/3 & & 2/3 \\ & 1 \end{matrix} | x) = _{2}F_{1}( \begin{matrix}1/12 & & 5/12 \\ & 1 \end{matrix}| \frac{64x^{3}(1 − x)}{(9 − 8x)^{3}} ) · ( 1 − \frac{8x}{9}) ^{−1/4} ,

obtained by taking α = \frac{1}{12} in the identity [21, Eq. (126)] (see also [37, Eq. (25)]):

2F_{1}( \begin{matrix}4α & 4α + \frac{1}{3} \\ 6α + \frac{1}{2} \end{matrix} | x) = _{2}F_{1}( \begin{matrix}α α + \frac{1}{3} \\ 2α + \frac{5}{6} \end{matrix} | \frac{64x^{3}(1 − x)}{(9 − 8x)^{3}} ) · ( 1 − \frac{8x}{9}) ^{3α}.

The detailed reduction is as follows: Equations (20) and (21) show that

2F_{1}( \begin{matrix}1/3 & & 2/3 \\ & 2 \end{matrix} | x) = \frac{9(1 − x)}{2} · ( _{2}F_{1}( \begin{matrix}1/12 & & 5/12 \\ & 1 \end{matrix}| ψ(x)) · χ(x)) ′,

where ψ(x) = 64x^{3}(1 − x) /(9 − 8x)^{3} and χ(x) = (1 − \frac{8x}{9} )^{−1/4} . By substitution, the latter equality implies

2F_{1}( \begin{matrix}1/3 & & 2/3 \\ & 2 \end{matrix} | φ(x)) =

\frac{9}{2} · \frac{1 − φ(x)}{φ′(x)} · ( _{2}F_{1}( \begin{matrix}1/12 & & 5/12 \\ & 1 \end{matrix}| ψ(φ(x))) · χ(φ(x))) ′,

with φ(x) = 27x(2 − 3x) /(1 − 4x)^{3} . The right-hand side term rewrites

\frac{(1− 64x)(1− x)(1− 4x)}{12(1 + 6x)} · ( _{2}F_{1}( \begin{matrix}1/12 & & 5/12 \\ & 1 \end{matrix}| \frac{1}{J(x)}) · ( \frac{g_{2}}{(1 − 4x)^{4}} ) ^{−1/4} ) ′,

and this, in combination with equation (16), finally implies equation (19).

3.4. The case of 3D queens. A closely related problem is to find the generating function for the diagonal sequence for 3D queen paths. This is Sequence A143734 in Sloane’s encyclopedia

http://oeis.org/A143734, whose first terms read

1, 13, 638, 41476, 3015296, 232878412, 18691183682, 1540840801552, …

Computationally, it turns out that this problem is more difficult than the one on 3D rooks. By Lemma 1, it boils down to finding the complete diagonal Q(x) of the rational function

( 1− \frac{s}{1 − s} − \frac{t}{1 − t} − \frac{x}{1 − x} − \frac{st}{1 − st} − \frac{tx}{1 − tx} − \frac{xs}{1 − xs} − \frac{stx}{1 − stx} ) ^{−1} .

We empirically found (using guessing routines based on Hermite– Padé approximants) that Q(x) satisfies a differential equation L_{Q}(x, ∂_{x}) of order 6 and polynomial coefficients (in x) of degree at most 71. The corresponding sequence satisfies a recurrence of order 14 with polynomial coefficients (in n) of degree at most 52. Moreover, it is very likely that these are the minimal-order equations. They are available online at

     http://algo.inria.fr/chyzak/Rooks/QueensRecEq.mpl

for the recurrence and

    http://algo.inria.fr/chyzak/Rooks/QueensDiffEq.mpl

for the differential equation.

The conjectured equations already have a large size, so that proving their correctness by the creative-telescoping approach becomes a challenge; we have not yet succeeded in proving this in a fully algorithmic fashion.

However, we are convinced that our guessed equations are correct: they pass all the various check tests described in [10]. In addition, the degree-71 leading coefficient of the differential equation factors as

x^{2} · (5x − 4) · (x + 1) · (x − 1)^{6} · (36x^{2} − 40x + 9) ·

(512x^{4} − 661x^{3} + 84x^{2} + 95x − 1) ·

(14063x^{6} − 15940x^{5} − 5918x^{4} + 12063x^{3} − 4118x^{2} + 575x − 40) · P (x)

where P (x) ∈ ℚ[x] is irreducible of degree 49.

Note the known fact [27, 16] that the nth coefficient of Q(x) is asymptotically equal to κc^{−3n} /n, for some κ ∈ ℝ and c ≈ 0.2185. More specifically, c is the smallest positive real root of the equation

1 − \binom{3}{1} \frac{x}{x − 1} − \binom{3}{2} \frac{x^{2}}{1 − x^{2}} − \binom{3}{3} \frac{x^{3}}{1 − x^{3}} = 0,

which is consistent with the presence of the factor 512x^{4} − 661x^{3} + 84x^{2} + 95x − 1 having c^{3} ≈ 0.0104 as a root (in fact, its root that is closest to zero).

The sixth-order differential operator L_{Q}(x, ∂_{x}) can be factored as L_{Q} = L_{4} · R_{1}^{(a)} · R_{1}^{(b)}, where L_{4} has order 4, and both R_{1}^{(a)} and R_{1}^{(b)} have order 1. By analogy with our result on rook paths, a natural question is whether L_{4} is solvable in terms of indefinite integrals, algebraic functions, and _{2}F_{1}’s expressions. If it were, then it would be solvable in terms of solutions of second-order equations, and by [29] at least one of the following should then be true (see [35, Sections 1.1 and 1.2] for definitions and details):

(1) L_{4} is reducible,

(2) L_{4} is gauge equivalent (under Maple’s DEtools[Homomorphisms]) to the third symmetric power of a second-order operator,

(3) L_{4} is gauge equivalent to the symmetric product of two secondorder operators.

Regarding option (1), the operator does not factor (DFactor in Maple). Regarding options (2) and (3), [35, Section 3] shows how to find order reductions of such type (if they exist). We found that no such reduction exists. Option (2) can be ruled out quickly by local data (if a third symmetric power has a logarithmic solution at x = p, then the cube of a logarithm should appear as well). For option (3), we computed the eigenring of the second exterior power of L_{4}. To conclude, the fourthorder factor L_{4} is not solvable (in the sense of Singer’s definition [29]) in terms of solutions of second-order equations; in particular, this means that, in this definition, there is no expression in terms of _{2}F_{1} functions for the generating function Q(x) of diagonal 3D queen paths.

It appears that L_{4} has no _{4}F_{3} type solution either, but we have not obtained a proof of this fact yet.

3.5. Faster but heuristic approach to getting the certificates. As already mentioned in the introduction, we have chosen in Section 2 to confine ourselves to calculations that follow an algorithmic approach, that is, that would provably (theoretically) succeed on any other rational input. Beside this, heuristics are available to cleverly guess a solution to the key equation, including certificates, and often lead to faster computations when they succeed, although they can fail in theory and sometimes do in practice. One of the most advanced strategies in this spirit is due to Koutschan [23]. Note that by rational-function normalisation, verifying the key equation a posteriori proves the obtained results.

We have implemented such a strategy in the rational case, and we could get a solution (P, S′, T′), where P is the same operator as in Section 2.4.1, but for other certificates S′and T′. In effect, this heuristically solves a partial differential equation for (some of) its rational solutions. The calculation is 4 times faster than in the algorithmic way. Interestingly enough, the denominators of S′ and T′ only involve q_{1}, and not its discriminant disc_{t}(q_{1}).

Acknowledgements. We wish to thank the referee for his quick and serious job and for his constructive comments on the manuscript.

References

[1] S. A. Abramov and K. Yu. Kvashenko. Fast algorithms for the search of the rational solutions of linear differential equations with polynomial coefficients. In Stephen M. Watt, editor, Symbolic and Algebraic Computation, pages 267– 270, New York, 1991. ACM Press. Proc. ISSAC’91, Bonn, Germany.

[2] S.A. Abramov and M. van Hoeij. Integration of solutions of linear functional equations. Integral Transforms and Special Functions, 8:3–12, 1999.

[3] G. Almkvist and D. Zeilberger. The method of differentiating under the integral sign. J. Symbolic Comput., 10(6):571–591, 1990.

[4] M. A. Barkatou. On rational solutions of systems of linear differential equations. J. Symbolic Comput., 28(4-5):547–567, 1999. Differential algebra and differential equations.

[5] I. N. Bernšteĭn. Modules over a ring of differential operators. Study of the fundamental solutions of equations with constant coefficients. Funct. Anal. Appl., 5(2):1–16 (Russian); 89–101 (English transl.), 1971.

[6] I. N. Bernšteĭn. The analytic continuation of generalized functions with respect to a parameter. Funct. Anal. Appl., 6(4):26–40 (Russian); 273–285 (English transl.), 1972.

[7] F. Beukers. Irrationality of π^{2}, periods of an elliptic curve and Γ_{1}(5). In Diophantine approximations and transcendental numbers (Luminy, 1982), volume 31 of Progr. Math., pages 47–66. Birkhäuser Boston, Mass., 1983.

[8] F. Beukers, 2010. Private communication.

[9] A. Bostan, F. Chyzak, B. Salvy, G. Lecerf, and É. Schost. Differential equations for algebraic functions. In ISSAC 2007, pages 25–32. ACM, New York, 2007.

[10] A. Bostan and M. Kauers. Automatic classification of restricted lattice walks. In DMTCS Proceedings of the 21st International Conference on Formal Power Series and Algebraic Combinatorics (FPSAC’09), Hagenberg, Austria, pages 203–217, 2009.

[11] F. Chyzak. An extension of Zeilberger’s fast algorithm to general holonomic functions. Discrete Math., 217(1-3):115–134, 2000. Formal power series and algebraic combinatorics (Vienna, 1997).

[12] F. Chyzak, M. Kauers, and B. Salvy. A non-holonomic systems approach to special function identities. In John May, editor, ISSAC ’09: Proceedings of the twenty-second international symposium on Symbolic and algebraic computation, pages 111–118, 2009.

[13] R. Debeerst. Solving differential equations in terms of Bessel functions. Master’s thesis, Universität Kassel, 2007.

[14] R. Debeerst, M. van Hoeij, and W. Koepf. Solving differential equations in terms of Bessel functions. ISSAC’08 Proceedings, pages 39–46, 2008.

[15] B. Dwork. Differential operators with nilpotent p-curvature. Amer. J. Math., 112(5):749–786, 1990.

[16] M. Erickson, S. Fernando, and K. Tran. Enumerating Rook and Queen Paths. Bull. Inst. Combin. Appl., 60:37–48, 2010.

[17] M. C. Fasenmyer. Some generalized hypergeometric polynomials. Bull. Amer. Math. Soc., 53:806–812, 1947.

[18] P. Flajolet. Analytic models and ambiguity of context-free languages. Theoret. Comput. Sci., 49(2-3):283–309, 1987. Twelfth international colloquium on automata, languages and programming (Nafplion, 1985).

[19] P. Flajolet and A. Odlyzko. Singularity analysis of generating functions. SIAM J. Discrete Math., 3(2):216–240, 1990.

[20] P. Flajolet and R. Sedgewick. Analytic combinatorics. Cambridge University Press, Cambridge, 2009.

[21] É. Goursat. Sur l’équation différentielle linéaire qui admet pour intégrale la série hypergéométrique. Ann. Sci. École Norm. Sup. (2), 10:3–142, 1881.

[22] E. L. Ince. Ordinary Differential Equations. Dover, 1956. Reprint of the 1926 edition.

[23] C. Koutschan. A fast approach to creative telescoping. Math. Comput. Sci., 4(2–3):259–266, 2010.

[24] L. Lipshitz. The diagonal of a D-finite power series is D-finite. J. Algebra, 113(2):373–378, 1988.

[25] L. Lipshitz. D-finite power series. J. Algebra, 122(2):353–373, 1989.

[26] M. Petkovšek. Hypergeometric solutions of linear recurrences with polynomial coefficients. J. Symbolic Comput., 14(2-3):243–264, 1992.

[27] A. Raichev and M. C. Wilson. A new method for computing asymptotics of diagonal coefficients of multivariate generating functions. In 2007 Conference on Analysis of Algorithms, AofA 07, Discrete Math. Theor. Comput. Sci. Proc., AH, pages 439–449. Assoc. Discrete Math. Theor. Comput. Sci., Nancy, 2007.

[28] H. A. Schwarz. Ueber diejenigen Fälle, in welchen die Gaussische hypergeometrische Reihe eine algebraische Function ihres vierten Elementes darstellt. J. reine angew. Math., 75:292–335, 1873.

[29] M. F. Singer. Solving homogeneous linear differential equations in terms of second order linear differential equations. Amer. J. Math., 107(3):663–696, 1985.

[30] R. P. Stanley. Differentiably finite power series. European J. Combin., 1:175– 188, 1980.

[31] R. P. Stanley. Enumerative combinatorics. Vol. 2, volume 62 of Cambridge Studies in Advanced Mathematics. Cambridge University Press, Cambridge, 1999. With a foreword by Gian-Carlo Rota and appendix 1 by Sergey Fomin.

[32] N. Takayama. An approach to the zero recognition problem by Buchberger algorithm. J. Symbolic Comput., 14(2-3):265–282, 1992.

[33] K. Takeuchi. Commensurability classes of arithmetic triangle groups. J. Fac. Sci. Univ. Tokyo, 24:201–212, 1977.

[34] M. van Hoeij, 2001. Implementation:

http://www.math.fsu.edu/~hoeij/files/equiv.

[35] M. van Hoeij. Solving third order linear differential equations in terms of second order equations. ISSAC’07 Proceedings, pages 355–360, 2007.

[36] M. van Hoeij and Q. Yuan. Finding all Bessel type solutions for linear differential equations with rational function coefficients. ISSAC’10 Proceedings, pages 37–44, 2010.

[37] R. Vidūnas. Algebraic transformations of Gauss hypergeometric functions. Funkcial. Ekvac., 52(2):139–180, 2009.

[38] W. Wasow. Asymptotic expansions for ordinary differential equations. Dover, New York, 1987. Reprint of the 1976 edition.

[39] D. Zeilberger. Sister Celine’s technique and its generalizations. J. Math. Anal. Appl., 85(1):114–145, 1982.

[40] D. Zeilberger. A fast algorithm for proving terminating hypergeometric identities. Discrete Math., 80(2):207–211, 1990.

[41] D. Zeilberger. A holonomic systems approach to special functions identities. J. Comput. Appl. Math., 32(3):321–368, 1990.

[42] D. Zeilberger. The method of creative telescoping. J. Symbolic Comput., 11(3):195–204, 1991.

Appendix A. Details of calculations of the iterated

creative telescoping for the 3D rook generating series

In this appendix, the calculation sketched in §2.4.2 is explained in more detail, and we provide explicit values of intermediate computations. We start by introducing the notions needed for our description.

Differential operators. Whereas the sketch of the algorithm was given with references to the function F to be integrated and its derivatives, here we prefer presenting the calculations as they are done in actual implementations: a primal role is played by differential operators that represent linear combinations of derivatives of a given function.

As certain operators map the function under consideration to zero, operators have to be viewed up to addition of such annihilating operators. This explains the use of congruences in the calculations below, denoted L_{1} ≡ L_{2} mod I, when L_{1}(F) = L_{2}(F) and I is the (left) ideal of all operators that annihilate F.

Representing a D-finite function. Recall from [30, 25] the definition that a D-finite function of variables u and v is a function F for which there exists a finite maximal set of derivatives ∂_{u}^{a}∂_{v}^{b}(F) that are linearly independent over ℚ(u, v). The natural question is how a function is algorithmically recognised and represented as D-finite.

The best algorithmic answer would be in terms of non-commutative Gröbner bases (and ideals of Hilbert dimension 0) [12], but here, to keep the sophistication to a minimum, we shall use the alternate definition that a function is D-finite if there are non-zero operators A(u, v, ∂_{u}) and B(u, v, ∂_{v}) that annihilate it, one for each derivative. Following [12], the pair (A, B) of annihilating operators will be called a rectangular system.

Rational solutions of parametrised systems of linear differential equations. Chyzak’s algorithm relies on an algorithm that, given matrices A ∈ M_{n}(𝕂(t)) and B ∈ M_{n×d} (𝕂(t)) for some field 𝕂 that depends on the application, solves

(22)  \frac{∂y}{∂t} + Ay = Be

for y ∈ 𝕂(t)^{n} and e ∈ 𝕂^{d}. In our uses, the field 𝕂 is a rational function field in the other variables of the problem. We do not describe this algorithm here, and refer the reader for instance to [1] (and the description of its use in [11]).

For convenience, let us recall the notation of §2.4.1. Let q be the denominator of F, then q = stq_{1} with q_{1} = −st + 2s^{2}t + 2t^{2} + 2xs − 3st^{2} − 3xt − 3xs^{2} + 4xst an irreducible polynomial of ℚ[x, s, t]. Let disc_{t}(q_{1}) = (x−s)(16xs^{2}−4s^{3}−24xs+4s^{2}+9x−s) be the discriminant of q_{1} with respect to t.

A.1. Stage A: First iteration of Chyzak’s algorithm. Stage A in §2.4.2 looks for operators P^{(α)} ∈ ℚ(x, s)⟨∂_{x}, ∂_{s}⟩ (according to the ansatz (10)) and rational functions ϕ^{(α)} ∈ ℚ(x, s, t) such that P^{(α)}(F) = \frac{∂}{∂t}(ϕ^{(α)}F) (cf. (11)). A property of the left ideal ⟨P^{(α)}⟩ thus obtained in ℚ(x, s)⟨∂_{x}, ∂_{s}⟩ is that it contains a rectangular system.

Ansatz with r = 1. We first look for operators P^{(α)} of total degree r = 1. Let η_{0,0}, η_{1,0}, η_{0,1} ∈ ℚ(x, s) be unknowns such that P = η_{0,0} + η_{1,0}∂_{x}+η_{0,1}∂_{s}, and let ϕ ∈ ℚ(x, s, t) be another unknown. The equation P (F) = \frac{∂}{∂t}(ϕF) rewrites

η_{0,0}F + η_{0,1} \frac{∂F}{∂s} + η_{1,0} \frac{∂F}{∂x} = \frac{∂ϕ}{∂t}F + ϕ\frac{∂F}{∂t} ,

which is an equation in ϕ and η_{0,0}, η_{0,1}, η_{1,0}, of the form (22) with 𝕂 = ℚ(x, s). All solutions are proportional to the one given by (23)

\{ \begin{matrix}η_{0,0} & = 2(s − 1)(3s^{2} − 6s + 2), \\ η_{1,0} & = 6xs^{3} − 2s^{3} − 10xs^{2} + s^{2} − 4x^{2}s + 10xs + 3x^{2} − 4x, \\ η_{0,1} & = 2s(3s − 2)(s − 1)^{2}, \\ ϕ & = \frac{−t(s − 1)(−6xs^{2} + 6s^{2}t − s^{2} + 11xs − 9st − 4x − xt + 4t)}{t − x}  . \end{matrix}

Set P_{1} to the operator P obtained for these specific values. The ideal generated by it obviously does not contain any rectangular system in ℚ(x, s)⟨∂_{x}, ∂_{s}⟩.

Ansatz with r = 2. Since all the monomials ∂_{s}^{i} can be reduced by P_{1}, we set r = 2 and look for P = η_{0}+η_{1}∂_{x}+η_{2}∂_{x}^{2}, with η_{0}, η_{1}, η_{2} ∈ ℚ(x, s), and ϕ ∈ ℚ(x, s, t) such that P (F) = \frac{∂}{∂t}(ϕF). (This ansatz is not strictly of the form (10) that would involve ∂_{s}; this will require additional manipulations when obtaining a rectangular system.) The equation on F becomes

η_{0}F + η_{1} \frac{∂F}{∂x} + η_{2} \frac{∂^{2}F}{∂x^{2}} = \frac{∂ϕ}{∂t}F + ϕ\frac{∂F}{∂t} ,

which is a parametrised equation in ϕ and η_{0}, η_{1}, η_{2}, with 𝕂 = ℚ(x, s). All solutions are proportional to the one given by:

(24) \{ \begin{matrix}η_{0} & = 0, \\ η_{1} & = −2(−19s^{2} − 9x + 13s^{3} + 7s − 16xs^{2} + 24xs), \\ η_{2} & = disc_{t}(q_{1}), \\ ϕ & = \frac{−t(3s − 2)(s − 1)(2s^{2} − 4st − s + 3t)(s − t)^{2}}{(t − x)q_{1}}  , \end{matrix}

and we call P_{2} the operator P obtained for these specific values. The ideal ⟨P_{1}, P_{2}⟩ now contains a rectangular system: eliminating ∂_{x} from P_{1} and P_{2} results in an operator P′_{2} = η′_{0}+η′_{1}∂_{s} +η′_{2}∂_{s}^{2}, with coefficients of respective degrees (4, 9), (4, 10), and (4, 11) in (x, s). This ends the calculation.

The Maple code to compute P′_{2} can be obtained at

http://algo.inria.fr/chyzak/Rooks/RectangularSystem.mpl .

A.2. Stage B: Second iteration of Chyzak’s algorithm. To perform Stage B in §2.4.2, we now work modulo the ideal I = ⟨P_{1}, P_{2}⟩ of ℚ(x, s)⟨∂_{x}, ∂_{s}⟩. The vector space V = ℚ(x, s)⟨∂_{x}, ∂_{s}⟩/I is 2-dimension-al over ℚ(x, s), with basis (1 + I, ∂_{x} + I). In other words, we consider a generic solution F̂ of P_{1} and P_{2}, and represent it by 1 + I.

To simplify the notation in what follows, we normalise the operators into the form P_{1}/(2s(3s−2)(s−1)^{2}) = ∂_{s}−p_{1}∂_{x}−p_{0} and P_{2}/disc_{t}(q_{1}) = ∂_{x}^{2} − q∂_{x}, for coefficients p_{0}, p_{1}, q ∈ ℚ(x, s).

The algorithm now starts from ansatz (8) for P, but with the form (9) for S, in order to find a relation of the form (12). For increasing values of d, it looks for an operator P = ∑ _{i=0}^{d} η_{i}∂_{x}^{i}, such that P ≡ ∂_{s}Q mod I, with η_{i} ∈ ℚ(x), ϕ_{0}, ϕ_{1} ∈ ℚ(x, s), and Q = ϕ_{0} + ϕ_{1}∂_{x}.

There is no solution for d < 3, and we give details for d = 3 only. First, the algorithm simplifies (7). Modulo I, we have ∂_{x}^{2} ≡ q∂_{x}, and therefore

∂_{x}^{3} ≡ \frac{∂q}{∂x}∂_{x} + q∂_{x}^{2} ≡ ( \frac{∂q}{∂x} + q^{2}) ∂_{x}.

Additionally, ∂_{s} ≡ p_{0} + p_{1}∂_{x}, so

∂_{x}∂_{s} ≡ \frac{∂p_{0}}{∂x} + ( p_{0} + \frac{∂p_{1}}{∂x}) ∂_{x} + p_{1}∂_{x}^{2} ≡ \frac{∂p_{0}}{∂x} + ( p_{0} + \frac{∂p_{1}}{∂x} + p_{1}q) ∂_{x}.

As to the right-hand side:

∂_{s}Q = \frac{∂ϕ_{0}}{∂s} + ϕ_{0}∂_{s} + \frac{∂ϕ_{1}}{∂s} ∂_{x} + ϕ_{1}∂_{s}∂_{x}

≡ \frac{∂ϕ_{0}}{∂s} +(p_{0} + p_{1}∂_{x}) ϕ_{0}+ \frac{∂ϕ_{1}}{∂s} ∂_{x} +ϕ_{1} ( \frac{∂p_{0}}{∂x} + ( p_{0} + \frac{∂p_{1}}{∂x} + p_{1}q) ∂_{x}) .

Since P ≡ η_{0} + (η_{1} + qη_{2} + ( \frac{∂q}{∂x} + q^{2}) η_{3}) ∂_{x}, the congruence P ≡ ∂_{s}Q mod I expressed in the basis (1 + I, ∂_{x} + I) is equivalent to the differential system

(25)

\{ \begin{matrix}\frac{∂ϕ_{0}}{∂s} & = −p_{0}ϕ_{0} − \frac{∂p_{0}}{∂x}ϕ_{1} + η_{0}, \\ \frac{∂ϕ_{1}}{∂s} & = −p_{1}ϕ_{0} − ( p_{0} + \frac{∂p_{1}}{∂x} + p_{1}q) ϕ_{1} + η_{1} + qη_{2} + ( \frac{∂q}{∂x} + q^{2}) η_{3}. \end{matrix}

Solving this parametrised system of linear differential equations in ℚ(x) results in a 1-dimensional vector space, with basis given by the solution

(26)  \{ \begin{matrix}η_{0} & = 0, \\ η_{1} & = 4(576x^{3} − 801x^{2} − 108x + 74), \\ η_{2} & = 4608x^{4} + 813x^{2} − 6372x^{3} + 514x − 4, \\ η_{3} & = x(x − 1)(64x − 1)(3x − 2)(6x + 1), \\ ϕ_{0} & = \frac{(53 + 108x)(3s − 2)s}{s − 1} , \\ ϕ_{1} & = \frac{γ}{2(s − 1)^{2}disc_{t}(q_{1})} , \end{matrix}

where γ is an irreducible polynomial in ℚ[x, s] of degree (5, 7) in (x, s).

The Maple code to compute γ can be obtained at

http://algo.inria.fr/chyzak/Rooks/SecondIteration.mpl .

A.3. Stage C: Reconstruction. We finally derive an explicit identity of the form (13) in Stage C in §2.4.2 and obtain the corresponding values for S and T given by (14), thus leading to the solution (P, S, T) of (3).

The previous calculation of P_{1} and P_{2} can be reinterpreted as follows: Introducing the left ideal J that annihilates F, that is, the ideal J generated by

F ∂_{x} − \frac{dF}{dx}, F ∂_{s} − \frac{dF}{ds}, F ∂_{t} − \frac{dF}{dt},

we have obtained the congruences

P_{1} ≡ ∂_{t}ψ_{1} mod J and P_{2} ≡ ∂_{t}ψ_{2} mod J,

for rational functions ψ_{1} and ψ_{2} that are the functions ϕ in (23) and (24), respectively. Additionally, P ≡ ∂_{s}Q mod I for P = η_{0} + η_{1}∂_{x} + η_{2}∂_{x}^{2} and Q = ϕ_{0} + ϕ_{1}∂_{x} described by (26), and I = ⟨P_{1}, P_{2}⟩.

To express P − ∂_{s}Q ∈ ⟨P_{1}, P_{2}⟩ more explicitly, we look for A_{1}, A_{2} ∈ ℚ(x, s)⟨∂_{x}, ∂_{s}⟩ such that P − ∂_{s}Q = A_{1}P_{1} + A_{2}P_{2}. This division is made algorithmic by observing that (P_{1}, P_{2}) is a Gröbner basis, and by performing a reduction by this Gröbner basis. In simple terms, the calculation is as follows: We first use P_{1} to eliminate ∂_{s} from P − ∂_{s}Q, then use P_{2} to divide the remainder. Collecting the quotients, we then find:

A_{1} = −\frac{108x + 53}{2(s − 1)^{3}} + \frac{γ_{1}}{4s(3s − 2)(s − 1)^{4}disc_{t}(q_{1})} ∂_{x},
A_{2} = \frac{γ_{2}}{4s(3s − 2)(s − 1)^{4}disc_{t}(q_{1})^{2}} +

\frac{x(x − 1)(64x − 1)(3x − 2)(6x + 1)}{disc_{t}(q_{1})} ∂_{x},

with γ_{1}, γ_{2} irreducible polynomials in ℚ[x, s] of respective degrees (5, 7) and (7, 10) in (x, s).

Since P_{1}(F) = ∂_{t}(ψ_{1}F) and P_{2}(F) = ∂_{t}(ψ_{2}F), we get the equality (P −∂_{s}Q)(F) = (A_{1}∂_{t}ψ_{1}+A_{2}∂_{t}ψ_{2})(F). Since A_{1} and A_{2} do not depend on t, we get P (F) = (∂_{s}Q + ∂_{t}(A_{1}ψ_{1} + A_{2}ψ_{2}))(F).

Let S = Q(F) and T = (A_{1}ψ_{1} + A_{2}ψ_{2})(F). Then S, T ∈ ℚ(x, s, t) and P (F) = \frac{∂S}{∂s} + \frac{∂T}{∂t}. It turns out that S and T are of the very special form

S = \frac{(s − t)U}{2stq_{1}^{2}disc_{t}(q_{1})}, T = \frac{(s − t)V}{2s^{2}q_{1}^{3}disc_{t}(q_{1})^{2}},

with U, V irreducible polynomials in ℚ[x, s, t] of respective degrees (5, 8, 3) and (8, 14, 5) in (x, s, t).

The Maple code for this reconstruction of S and T can be obtained

at http://algo.inria.fr/chyzak/Rooks/Reconstruction.mpl .

INRIA (France)

  E-mail address: alin.bostan@inria.fr

INRIA (France)

  E-mail address: frederic.chyzak@inria.fr

Florida State University, Department of Mathematics (Tallahassee, FL, USA)

  E-mail address: hoeij@math.fsu.edu

École Normale Supérieure (France)

  E-mail address: lucien.pech@gmail.com