Linear-Time and Constant-Space Algorithms to compute Multi-Sequences that arise in Enumerative Combinatorics (and Elsewhere)

Shalosh B. EKHAD and Doron ZEILBERGER

Abstract: How many ways, exactly, can a Chess King, always moving forward, i.e. with the set of steps {[1, 0], [0, 1], [1, 1]}, walk from the origin, (0, 0), to the lattice point (100000, 200000)? Thanks to the amazing Apagodu-Zeilberger extension of the Almkvist-Zeilberger algorithm, adapted in this article for combinatorial applications, this 104492-digit number, viewable from

http://www.math.rutgers.edu/~zeilberg/tokhniot/oPureRecRat2.txt, can be computed in

less than 33 seconds. But not just this particular number! Many other numbers that come up in enumerative combinatorics, and elsewhere, can be computed just as efficiently.

Maple Packages

This article is accompanied by two Maple packages, PureRec.txt, and PureRecRat.txt, available from

https://sites.math.rutgers.edu/~zeilberg/tokhniot/PureRec.txt and
https://sites.math.rutgers.edu/~zeilberg/tokhniot/PureRecRat.txt .

The web-page of this article,

https://sites.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/pure.html ,

contains input and output files, referred to in this paper.

The purpose of this short article is to be a pointer to the much longer Maple packages (written by DZ) and numerous output files (generated by SBE) given in the above mentioned web-page. Recall that the Almkvist-Zeilberger algorithm [AlZ] (see [D] for an engaging exposition) produces a pure linear recurrence equation with polynomial coefficients, satisfied by the Taylor coefficients of many kinds of generating functions. All we need is that R′(x)/R(x) is a rational function of x. This includes all rational functions, their powers, R(x)^{α} (either for symbolic or numeric exponent α) and all functions of the form e^{P (x)} where P (x) is a polynomial, and many other functions besides. This follows from the fact, that thanks to Cauchy, it can be written as a contour integral \frac{1}{2πi} ∫_{|x|=1} \frac{R(x)}{x^{n+1}} dx. In fact, it can do much more, it can handle any expressions of the form ∫_{|x|=1} R(x, n) dx, where R(x, n) satisfies the condition that both \frac{\frac{d}{dx} R(n,x)}{R(n,x)} and \frac{R(n+1,x)}{R(n,x)} are rational functions of the discrete variable n and the continuous variable x.

Since the recurrence is finite order (often fairly low), this enables a very fast linear-time and constant-memory computation of such sequences. If the order of the recurrence is L, one only has to ‘remember’ the last L values. Of course, since the numbers (usually) grow exponentially, the bit-size memory is linear, but still it is very good!

After dividing by the leading term, one can write the recurrence for the sequence, let’s call is a(n) as

a(n) = ∑_{i=1}^{L} f_{i}(n)a(n − i) ,

where f_{i}(n) are rational functions of n. There is a potential issue of ‘dividing by zero’, but if n_{0} is the largest positive integer root of the denominators, one can directly compute the first n_{0} values and start the recurrence there.

The amazing Multi-variable Almkvist-Zeilberger algorithm, designed by Moa Apagodu and Doron Zeilberger [ApZ] extends this to several variables. If a(n_{1}, …, n_{d}) is such a multi-sequence, gotten, for example. from the Taylor coefficients of a multi-variable rational function in d variables (and many other kinds of multi-sequences) , then the algorithm finds d pure recurrences, one for each discrete axis:

a(n_{1}, …, n_{d}) = ∑_{i=1}^{L_{j}} f_{i}^{(j)}(n)a(n_{1}, …, n_{j−1}, n_{j} − i, n_{j+1}, …, n_{d}) ,

where f_{i}^{(j)}(n_{1}, …, n_{d}) are rational functions of n_{1}, …, n_{d}. Then to compute a specific value, a(n_{1}, …, n_{d}), one constructs a discrete-path from the origin [0, …, 0] to the lattice point (a_{1}, …, a_{d}) with d segments, one for each direction. One such path does not suffice, since in order to say, climbup in the last axis n_{d}, we need L_{d} initial values, so we need recursively to take care of them L_{d} ‘initial conditions’, each requiring its own path in (d − 1)-discrete space, but one only needs finitely many paths, making it linear time and constant memory.

Alas, now the division by zero is not so easy to overcome. We chose to only use the most obvious generic paths (see the Maple code for procedure EvalScheme in the Maple package

https://sites.math.rutgers.edu/~zeilberg/tokhniot/PureRec.txt ), and our program re-

turns FAIL if one encounters a division by 0. However, we believe that this can be overcome by finding other paths that stay-away from dividing by 0.

Sample Output

For five random examples of sequences that are Taylor coefficients of one-variable functions of the form 1/P (x)^{\frac{1}{3}} , where P (x) is a randomly generated polynomial, see:

https://sites.math.rutgers.edu/~zeilberg/tokhniot/oPureRec1.txt .

For five random examples of sequences that are Taylor coefficients of one-variable functions that are exponentials of polynomials, see:

https://sites.math.rutgers.edu/~zeilberg/tokhniot/oPureRec1a.txt .

For five random examples of bi-sequences that are Taylor coefficients of two-variable functions of the form 1/P (x_{1}, x_{2})^{\frac{1}{3}} , where P (x_{1}, x_{2}) is a randomly generated polynomial of x_{1}, x_{2}, see:

https://sites.math.rutgers.edu/~zeilberg/tokhniot/oPureRec2.txt .

For five random examples of bi-sequences that are Taylor coefficients of functions that are exponentials of polynomials of two variables, see

https://sites.math.rutgers.edu/~zeilberg/tokhniot/oPureRec2a.txt .

For three random examples of triple-sequences that are Taylor coefficients of a three-variable function of the form 1/P (x_{1}, x_{2}, x_{3})^{\frac{1}{3}} , where P (x_{1}, x_{2}, x_{3}) is a randomly generated polynomial of x_{1}, x_{2}, x_{3}, see:

https://sites.math.rutgers.edu/~zeilberg/tokhniot/oPureRec3.txt .

Now it takes much longer to generate the scheme, but once found, it is still very fast to compute specific values, even those very far from the origin, For three random examples of tri-sequences that are Taylor coefficients of functions that are exponentials of polynomials of three variables, see

https://sites.math.rutgers.edu/~zeilberg/tokhniot/oPureRec3a.txt .

Here we encountered some singularities, so our (incomplete) implementation sometimes (but fairly rarely) returns FAIL. Applications to the Enumeration of Lattice Path in the Plane with Many choices of Atomic Steps Let St be an arbitrary set of positive steps, e.g. for the forward-moving Chess King it is {[1, 0], [0, 1], [1, 1]}. How can we compute fast the number of ways of walking from the origin to [a, b] using the members of St? The lattice-path oriented companion Maple package

https://sites.math.rutgers.edu/~zeilberg/tokhniot/PureRecRat.txt ,

has a procedure

Walk2Dpaper(St,K) ,

that inputs an arbitrary set of positive steps in 2-dimensions, and a large positive integer K (for illustrative purposes) and outputs the scheme for the bi-variate function F (a, b), defined as the number of walks from [0, 0] to [a, b] using the steps in St, and if no singularity is encountered, outputs F (K, 2K). It also outputs the one-variable linear recurrence for the diagonal sequence {F (a, a)} and uses it to compute F (2K, 2K) (now there are never singularities). For the sake of Neil Sloane, it also spits out the first 30 terms of the diagonal sequence.

You are welcome to try out procedure Walk2Dpaper(St,K) with St = {[0, 1], [1, 0], [1, 1]} and K = 100000. getting the 104922-digit number mentioned in the abstract.

For 29 such theorems see:

https://sites.math.rutgers.edu/~zeilberg/tokhniot/oPureRecRat1.txt .

Christoph Koutschan’s even better program (alas, in Mathematica)

A Mathematica implementation of the multi-variate Almkvist-Zeilberger algorithm is contained, inter alia, in Christoph Koutschan’s very efficient and versatile Mathematica program [K]. It would be nice if he would use the simple ideas in our paper to interface them and create such a multisequence calculator.

References

[AlZ] Gert Almkvist and Doron Zeilberger, The method of differentiating Under The integral sign, J. Symbolic Computation 10 (1990), 571-591.

https://sites.math.rutgers.edu/~zeilberg/mamarim/mamarimPDF/duis.pdf .

[ApZ] Moa Apagodu and Doron Zeilberger, Multi-Variable Zeilberger and Almkvist-Zeilberger Algorithms and the Sharpening of Wilf-Zeilberger Theory, Adv. Appl. Math. 37 (2006), 139-152. [Special issue in honor of Amitai Regev]

https://sites.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/multiZ.html .

[D] Robert Dougherty-Bliss, Integral recurrences from A to Z, American Mathematical Monthly, to appear.

https://arxiv.org/abs/2102.10170 .

[K] Christoph Koutschan, Advanced applications of the holonomic systems approach, PhD thesis, Research Institute for Symbolic Computation (RISC), Johannes Kepler University, Linz, Austria, 2009.

http://www.koutschan.de/publ/Koutschan09/thesisKoutschan.pdf,

Shalosh B. Ekhad and Doron Zeilberger, Department of Mathematics, Rutgers University (New Brunswick), Hill Center-Busch Campus, 110 Frelinghuysen Rd., Piscataway, NJ 08854-8019, USA.

Email: [ShaloshBEkhad, DoronZeil] at gmail dot com .

Exclusively published in the Personal Journal of Shalosh B. Ekhad and Doron Zeilberger and arxiv.org

March 9, 2022