A Case Study in Meta-AUTOMATION:

AUTOMATIC Generation of Congruence AUTOMATA For Combinatorial Sequences

Eric ROWLAND and Doron ZEILBERGER

Abstract: In this paper, that may be considered a sequel to a recent article by Eric Rowland and Reem Yassawi, we present yet another approach for the automatic generation of automata (and an extension that we call Congruence Linear Schemes) for the fast (log-time) determination of congruence properties, modulo small (and not so small!) prime powers, for a wide class of combinatorial sequences. Even more interesting than the new results that could be obtained, is the illustrated methodology, that of designing ‘meta-algorithms’ that enable the computer to develop algorithms, that it (or another computer) can then proceed to use to actually prove (potentially!) infinitely many new results. This paper is accompanied by a Maple package, AutoSquared, and numerous sample input and output files, that readers can use as templates for generating their own, thereby proving many new ‘theorems’ about congruence properties of many famous (and, of course, obscure) combinatorial sequences.

Very Important: This article is accompanied by the general Maple package

http://www.math.rutgers.edu/~zeilberg/tokhniot/AutoSquared ,

and several other specific ones, and numerous input and output files that are obtainable, by one click, from the webpage (“front”) of this article

http://www.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/meta.html  .

They could (and should!) be used as templates for generating as many input files that the human would care to type, and that the computer would agree to run.

Prologue: What are the Last Three (Decimal) Digits of the Googol-th Catalan, Motzkin, and Delannoy Numbers?

We will never know all the (decimal) digits of the Googol-th terms of the famous Catalan, Motzkin,

and (Central) Delannoy sequences [http://oeis.org/A000108, http://oeis.org/A001006, and

http://oeis.org/A001850, respectively], if nothing else because our computers are not big enough to store them!

But thanks to the Maple packages accompanying this article, we know for sure that the last three digits are 000, 187, and 281 respectively. These packages can compute in logarithmic time (i.e. linear in the number of digits of the input) the values of the n-th term modulo many different m (but alas, not too big!). These fast algorithms were generated by a meta-algorithm implemented in the main Maple package, AutoSquared.

Fast Exponentiation

E-commerce is possible (via RSA) thanks to the fact that it is very easy (for computers!) to compute

a^{n} mod m ,

for a and m several-hundred-digits long, and large n. Reminding you that a^{n} is shorthand for the sequence, let’s call it x_{n} defined by the linear recurrence equation with constant coefficients, of order one:

x_{n} − ax_{n−1} = 0 , x_{0} = 1 .

In order to compute x_{10^{100}} mod m, you don’t compute all the 10^{100} previous terms, but use the implied recurrences

x_{2n} = x_{n}^{2} mod m , x_{2n+1} = ax_{n}^{2} mod m .

This takes only log_{2} 10^{100} operations!

What about sequences defined by higher-order recurrences, but still with constant coefficients? For example, what are the last three decimal digits of the googol-th Fibonacci number, F_{10^{100}}? You would get the answer, 875, in 0.008 seconds!

All you need is type

Fnm(10**100, 1000);  ,

once you typed (or copied-and-pasted) the following short code into a Maple session:

Fnm:=proc(n, m) option remember;
if n = 1 or n = 2 then 1
elif n mod 2 = 0 then Fnm(1/2*n, m)*(Fnm(1/2*n + 1, m) + Fnm(1/2*n - 1, m)) mod m
else Fnm(1/2*n - 1/2, m)**2 + Fnm(1/2*n + 1/2, m)**2 mod m
fi:
end:

It implements the (nonlinear) recurrence scheme

F_{2n} = F_{n}(F_{n−1} + F_{n+1}) , F_{2n+1} = F_{n}^{2} + F_{n+1}^{2} , F_{1} = 1 , F_{2} = 1 ,

and of course takes it modulo m at every step.

Another way is to take the (1, 2) entry of the matrix

( \begin{matrix}1 & 1 \\ 1 & 0 \end{matrix} )^{10^{100}} mod 1000 ,

and use the ‘iterated-squaring’ trick applied to matrix (rather than scalar) exponentiation.

Both these simple methods are applicable for the fast (linear-in-bit-size) computation of the terms, modulo any m, of any integer sequence defined in terms of a linear recurrence equation with constant coefficients (aka C-finite integer sequences).

But what about sequences that are defined via linear recurrence equations with polynomial coefficients, aka P-recursive sequences, aka holonomic sequences?

In a beautiful and deep paper ([KKM]), dedicated to one of us (DZ) on the occasion of his 60^{th} birthday, Manuel Kauers, Christian Krattenthaler, and Thomas Müller developed a deep and ingenious method for the determination of holonomic sequences modulo powers of 2 (with applications to group theory!). This has been extended to powers of 3 in [KM1], and a different method for obtaining congruences was developed in [KM2].

An important subclass of the class of holonomic integer sequences is the class of integer sequences whose (ordinary) generating function, let’s call it f(x), satisfies an algebraic equation of the form P (f(x), x) = 0 where P is a polynomial of two variables with integer coefficients. For this class, and an even wider class, the sequences arising from the diagonals of rational functions of several variables, Rowland and Yassawi ([RY]) developed a very general method for computing finite automata for the fast computation (once the automaton is found, of course) of the congruence behavior modulo prime powers. Of course, as the primes and/or their powers get larger, the automata get larger too, but if the automaton is precomputed once and for all (and saved!), it is logarithmic time (i.e. linear in the bit-size). Of course, the implied constant in the O(log n) computation times gets larger with the moduli.

History

Many papers, in the past, proved isolated results about congruence properties for specific sequences and for specific moduli. We refer the reader to [RY] for many references, that we will not repeat here.

The Present Method: Using Constant Terms

Most (perhaps all) of the combinatorial sequences treated in [RY] can be written in the form

a_{n} := ConstantT ermOf [ P (x)^{n}Q(x) ] ,

where both P (x) and Q(x) are Laurent polynomials with integer coefficients, where x is either a single variable or a multi-variable x = (x_{1}, …, x_{m}), and ConstantTermOf means “coefficient of x^{0} ”, or “the coefficient of x_{1}^{0} ⋯ x_{m}^{0} ”.

For example, the arguably second-most famous combinatorial sequence (after the Fibonacci sequence), is the sequence of the Catalan Numbers (http://oeis.org/A000108), that may be defined by

C_{n} := ConstantT ermOf [ (\frac{1}{x} + 2 + x)^{n}(1 − x) ] .

Not as famous, but also popular, are the Motzkin numbers (http://oeis.org/A001006), that may be defined by

M_{n} := ConstantT ermOf [ (\frac{1}{x} + 1 + x)^{n}(1 − x^{2}) ] ,

and also fairly famous are the Central Delannoy Numbers (http://oeis.org/A001850), that may be defined by

D_{n} := ConstantT ermOf [ (\frac{1}{x} + 3 + 2x)^{n} ] .

So far, we got away with a single variable.

Another celebrated sequence is the sequence of Apéry Numbers, that were famously used by 64year-old Roger Apéry (in 1978) to prove the irrationality of ζ(3). These are defined in terms of a binomial coefficient sum

A(n) := ∑_{k=0}^{n} \binom{n}{k}^{2}\binom{n + k}{k}^{2} .

These may be equivalently defined (see below) as

A(n) := ConstantT ermOf [(\frac{(1 + x_{1}) (1 + x_{2}) (1 + x_{3}) (1 + x_{2} + x_{3} + x_{2}x_{3} + x_{1}x_{2}x_{3})}{x_{1}x_{2}x_{3}})^{n}] .

How to convert ANY Binomial Coefficient Sum into a Constant Term Expression?

Before describing our new method, let us indicate how any binomial coefficient sum of the form

A(n) = ∑_{k=0}^{n} \binom{n}{k}g^{k} ∏_{i=1}^{m} \binom{a_{i}n + b_{i}k + c_{i}}{d_{i}n + e_{i}k + f_{i}} ,

where all the a_{i}, b_{i}, c_{i}, d_{i}, e_{i}, f_{i} and g are integers, can be made into a constant term expression. (This is essentially Georgy Petrovich EGORYCHEV’s celebrated method of coefficients). We introduce m variables x_{1}, …, x_{m} and use the fact, that, by definition

\binom{a_{i}n + b_{i}k + c_{i}}{d_{i}n + e_{i}k + f_{i}} = ConstantT ermOf_{x_{i}} [ \frac{(1 + x_{i})^{a_{i}n+b_{i}k+c_{i}}}{x_{i}^{d_{i}n+e_{i}k+f_{i}}} ]

Hence

A(n) = ∑_{k=0}^{n} \binom{n}{k}g^{k} ∏_{i=1}^{m} \binom{a_{i}n + b_{i}k + c_{i}}{d_{i}n + e_{i}k + f_{i}}

= ∑_{k=0}^{n} \binom{n}{k}g^{k} ∏_{i=1}^{m} ConstantT ermOf_{x_{i}} [ \frac{(1 + x_{i})^{a_{i}n+b_{i}k+c_{i}}}{x_{i}^{d_{i}n+e_{i}k+f_{i}}} ]

= ConstantT ermOf_{x_{1},...,x_{m}} [ ∑_{k=0}^{n} \binom{n}{k}g^{k} ∏_{i=1}^{m} \frac{(1 + x_{i})^{a_{i}n+b_{i}k+c_{i}}}{x_{i}^{d_{i}n+e_{i}k+f_{i}}} ]

= ConstantT ermOf_{x_{1},...,x_{m}} [ (∏_{i=1}^{m} \frac{(1 + x_{i})^{a_{i}n+c_{i}}}{x_{i}^{d_{i}n+f_{i}}} ) ∑_{k=0}^{n} \binom{n}{k}g^{k} ∏_{i=1}^{m} (\frac{(1 + x_{i})^{b_{i}k}}{x_{i}^{e_{i}k}} ) ]

= ConstantT ermOf_{x_{1},...,x_{m}} [ ( ∏_{i=1}^{m} \frac{(1 + x_{i})^{a_{i}n+c_{i}}}{x_{i}^{d_{i}n+f_{i}}} ) ∑_{k=0}^{n} \binom{n}{k}g^{k} ∏_{i=1}^{m} (\frac{(1 + x_{i})^{b_{i}}}{x_{i}^{e_{i}}} )^{k} ]

= ConstantT ermOf_{x_{1},...,x_{m}} [ ∏_{i=1}^{m} (\frac{(1 + x_{i})^{a_{i}n+c_{i}}}{x_{i}^{d_{i}n+f_{i}}} ) (1 + g ∏_{i=1}^{m} \frac{(1 + x_{i})^{b_{i}}}{x_{i}^{e_{i}}} )^{n} ] .

= ConstantT ermOf_{x_{1},...,x_{m}} [ ∏_{i=1}^{m} (\frac{(1 + x_{i})^{c_{i}}}{x_{i}^{f_{i}}} ) (\frac{(1 + x_{i})^{a_{i}}}{x_{i}^{d_{i}}} )^{n} (1 + g ∏_{i=1}^{m} \frac{(1 + x_{i})^{b_{i}}}{x_{i}^{e_{i}}} )^{n} ] .

= ConstantT ermOf_{x_{1},...,x_{m}} [ ∏_{i=1}^{m} \frac{(1 + x_{i})^{c_{i}}}{x_{i}^{f_{i}}} ((∏_{i=1}^{m} \frac{(1 + x_{i})^{a_{i}}}{x_{i}^{d_{i}}} ) (1 + g ∏_{i=1}^{m} \frac{(1 + x_{i})^{b_{i}}}{x_{i}^{e_{i}}} ))^{n} ] .

This is implemented in procedure BinToCT(L,x,a) in our Maple package AutoSquared. For example, we got the above constant-term rendition of the Apéry numbers by typing:

BinToCT([ [[1,0,0],[0,1,0]], [[1,1,0],[0,1,0]]$2 ],x,1); .

Illustrating the Constant Term Approach In Terms of the Simplest-Not-Entirely-Trivial Example

Recall from above that the Catalan numbers may be defined by the constant-term formula

C_{n} := ConstantT ermOf [ (\frac{1}{x} + 2 + x)^{n}(1 − x) ] .

We are interested in the mod 2 behavior of C_{n}, in other words we want to have a quick way of computing C_{n} modulo 2. So let’s define

A_{1}(n) := C_{n} mod 2 .

Using the above formula for C_{n}, and taking it modulo 2, we have:

A_{1}(n) := ConstantT ermOf [ (1 + x)(\frac{1}{x} + x)^{n} ] .

We will try to find a constant-term expression for A_{1}(2n) .

A_{1}(2n) = ConstantT ermOf [ (1+x)(\frac{1}{x}+x)^{2n} ] mod 2 = ConstantT ermOf [ (1+x)((\frac{1}{x}+x)^{2})^{n} ] mod 2

= ConstantT ermOf [ (1+x)(\frac{1}{x^{2}} +2+x^{2})^{n} ] mod 2 = ConstantT ermOf [ (1+x)(\frac{1}{x^{2}} +x^{2})^{n} ] mod 2

But

ConstantT ermOf [ (1 + x)(\frac{1}{x^{2}} + x^{2})^{n} ] = ConstantT ermOf [ 1 · (\frac{1}{x^{2}} + x^{2})^{n} ] ,

since, obviously,

ConstantT ermOf [ x · (\frac{1}{x^{2}} + x^{2})^{n} ] = 0 .

Since the constant-termand of

ConstantT ermOf [ 1 · (\frac{1}{x^{2}} + x^{2})^{n} ] ,

only depends on x^{2}, we can replace x^{2} by x, implying that

A_{1}(2n) = ConstantT ermOf [ (\frac{1}{x} + x)^{n} ] mod 2 .

This forces us to put up with a new kid on the block, let’s call it A_{2}(n):

A_{2}(n) := ConstantT ermOf [ (\frac{1}{x} + x)^{n} ] mod 2 ,

and we got the recurrence

A_{1}(2n) = A_{2}(n) .

We will handle A_{2}(n) in due course, but first let’s consider A_{1}(2n + 1). We have

A_{1}(2n + 1) = ConstantT ermOf [ (1 + x)(\frac{1}{x} + x)^{2n+1} ] mod 2

= ConstantT ermOf [ (1 + x)(\frac{1}{x} + x)((\frac{1}{x} + x)^{2})^{n} ] mod 2

= ConstantT ermOf [ (\frac{1}{x} + x + 1 + x^{2})(\frac{1}{x^{2}} + 2 + x^{2})^{n} ] mod 2

= ConstantT ermOf [ (\frac{1}{x} + x + 1 + x^{2})(\frac{1}{x^{2}} + x^{2})^{n} ] mod 2 .

But

ConstantT ermOf [ (\frac{1}{x} + x + 1 + x^{2})(\frac{1}{x^{2}} + x^{2})^{n} ] = ConstantT ermOf [ (1 + x^{2}) · (\frac{1}{x^{2}} + x^{2})^{n} ] ,

since, obviously

ConstantT ermOf [ (\frac{1}{x} + x) · (\frac{1}{x^{2}} + x^{2})^{n} ] = 0 .

Since the constant-termand of

ConstantT ermOf [ (1 + x^{2}) · (\frac{1}{x^{2}} + x^{2})^{n} ] ,

only depends on x^{2}, we can replace x^{2} by x, implying that

A_{1}(2n + 1) = ConstantT ermOf [ (1 + x)(\frac{1}{x} + x)^{n} ] mod 2 .

But this looks familiar! It is good-old A_{1}(n), so we have established, so far, that

A_{1}(2n) = A_{2}(n) , A_{1}(2n + 1) = A_{1}(n) .

But in order to establish a recurrence scheme, we need to handle A_{2}(n). A priori, this may force us to introduce yet more discrete functions, and that would be OK, as long as we would finally stop, after finitely many steps, getting a scheme with finitely many discrete functions, that would enable the fast (logarithmic time) computation of our initial function A_{1}(n). We will see that this would always be the case, no matter how complicated P (x) and Q(x) are (and even with many variables). Alas, as P (x) gets more complicated, the ‘finite’ gets bigger and bigger, so eventually the ‘logarithmic time’ in n would be impractical, since the implied constant would be eeeeeeeeeeeeenormous.

But in this toy example, don’t worry! The ‘finitely many discrete functions’, is only two! As we will shortly see, all we need is A_{2}(n), in addition to A_{1}(n).

Recall that

A_{2}(n) := ConstantT ermOf [ (\frac{1}{x} + x)^{n} ] mod 2 .

Let’s try to find a constant-term expression for A_{2}(2n) .

A_{2}(2n) = ConstantT ermOf [ (\frac{1}{x} + x)^{2n} ] mod 2 = ConstantT ermOf [ ((\frac{1}{x} + x)^{2})^{n} ] mod 2

= ConstantT ermOf [ (\frac{1}{x^{2}} + 2 + x^{2})^{n} ] mod 2 = ConstantT ermOf [ (\frac{1}{x^{2}} + x^{2})^{n} ] mod 2

Since the constant-termand only depends on x^{2}, we can replace x^{2} by x, implying that

A_{2}(2n) = ConstantT ermOf [ (\frac{1}{x} + x)^{n} ] mod 2 .

But that’s exactly A_{2}(n), so we have found out that

A_{2}(2n) = A_{2}(n) .

What about A_{2}(2n + 1)? Here goes:

A_{2}(2n+1) = ConstantT ermOf [ (\frac{1}{x}+x)^{2n+1} ] mod 2 = ConstantT ermOf [ (\frac{1}{x}+x)((\frac{1}{x}+x)^{2})^{n} ] mod 2

= ConstantT ermOf [ (\frac{1}{x}+x)(\frac{1}{x^{2}} +2+x^{2})^{n} ] mod 2 = ConstantT ermOf [ (\frac{1}{x}+x)(\frac{1}{x^{2}} +x^{2})^{n} ] mod 2 .

But the constant-termand now only has odd powers, so the coefficient of x^{0}, alias the constant term, is 0. We have just established, the fast recurrence scheme:

A_{1}(2n) = A_{2}(n) , A_{1}(2n + 1) = A_{1}(n) ,

A_{2}(2n) = A_{2}(n) , A_{2}(2n + 1) = 0 ,

subject to the initial conditions

A_{1}(0) = 1 , A_{2}(0) = 1 .

[The above human-generated scheme can be also done (much faster) by the Maple package

AutoSquared. Having downloaded http://www.math.rutgers.edu/~zeilberg/tokhniot/AutoSquared

into your computer, that has Maple installed, you stay in the same directory, and you type:

read AutoSquared:  CA([1/x+2+x,1-x],x,2,1,2)[1];  ,

and you would get (in 0 seconds!), the output

[[[2, 1], [2, 0]], [1, 1]]  ,

which is our package’s way of encoding the above ‘scheme’.

Another way of describing the scheme is via the binary representation of n (for some k ≥ 1)

n = ∑_{i=1}^{k} α_{i}2^{k−i} ,

where α_{i} ∈ {0, 1}, α_{1} = 1, and it is abbreviated, in the positional notation, as a word, of length k, in the alphabet {0, 1}

α_{1} ⋯ α_{k} .

Phrased in terms of such ‘words’, the above scheme can be written, (where w is any word in the alphabet {0, 1})

A_{1}(w0) = A_{2}(w) , A_{1}(w1) = A_{1}(w) .
A_{2}(w0) = A_{2}(w) , A_{2}(w1) = 0 ,

subject to the initial conditions (here ϕ is the empty word):

A_{1}(ϕ) = 1 , A_{2}(ϕ) = 1 .

Let’s revert to post-fix notation for representing functions, and omit the parentheses, i.e. write wA_{1} instead of A_{1}(w) and wA_{2} instead of A_{2}(w). This will not cause any ambiguity, since the alphabet of function names, {A_{1}, A_{2}} is disjoint from the alphabet of letters, {0, 1}. The above scheme becomes

w0A_{1} = wA_{2} , w1A_{1} = wA_{1}
w0A_{2} = wA_{2} , w1A_{2} = 0 ,

subject to the initial conditions

ϕA_{1} = 1 , ϕA_{2} = 1 .

Let’s try to find A_{1}(30), alias, A_{1}(11110_{2}), alias, with our new convention, 11110A_{1}. We get in two steps

11110A_{1} = 1111A_{2} = 0 .

This only took two steps due to a premature exit to an output gate. The default number of steps is the length of the word, that keeps traveling until it becomes the empty word, and then it is forced to move to an output gate.

It is readily seen that if the input word has a zero in it, the output would be 0. Hence the only words that output 1 are those given by the regular expression

1^{∗} .

Equivalently, the only integers n for which the Catalan number C_{n} is odd are those of the form n = 2^{k} − 1 for k = 0, 1, 2, . . ..

The words in the alphabet {0, 1} that output 0 (i.e. those words that have at least one 0 in their binary representation) are the complement ‘language’, whose regular expression rendition is

1 {0, 1}^{∗} 0 {0, 1}^{∗} .

What we have here is a finite automaton with output. The set of states is {A_{1}, A_{2}} while the alphabet is the set {0, 1}. There are 2 directed edges coming out of each state, one for each letter of the alphabet, leading to another (possibly the same) state, or possibly to an output gate (in our case always 0, via ‘exit edges’ that prematurely end the journey. You have a starting state (in this example, A_{1}) and an input word, and you travel along the automaton, according to the current state and the current rightmost letter, until you run out of letters, i.e. have the empty word, or wind-up in the output 0 prematurely, since some states have edges that lead directly to 0. (In our example when you are at state A_{2} and the rightmost letter is 1 you immediately output 0.)

Yet another way of describing it is via a type-three grammar (aka regular grammar) in the famous Chomsky hierarchy (see e.g. [R]). For each possible output (in this example, 0 and 1, NOT TO BE CONFUSED WITH THE LETTERS OF THE ALPHABET), there is a regular grammar describing the language (set of words) that yield that output.

In this example, the set of non-terminal symbols is {A_{1}, A_{2}} and the set of terminal symbols is {0, 1}. For a grammar for the language yielding 1 (i.e. the binary representations of the integers n for which C_{n} is odd) the non-terminal symbol A_{2} is not needed (is superfluous), and the grammar is extremely simple

A_{1} → ϕ , A_{1} → 1A_{1} .

We leave it to the interested reader to write down the only slightly more complicated grammar for the language of binary representations of integers n for which C_{n} is even.

It is well-known that the notions of finite automata, regular expressions, and regular grammars are equivalent (as far as the generated languages), and there are easy algorithms for going between them.

These are all very nice, but for the present formulation, it is more convenient not to write the input integers n in base 2 (or more generally, base p, if the desired modulus is a power of a prime p), but stick to integers (as inputs). Let’s make the following formal definition.

Definition: Let ℕ be the set of non-negative integers, let p be a positive integer, and let E be any set. An automatic p-scheme for a function f : ℕ → E is a set of finitely many (say r) auxiliary functions A_{1}(n), …, A_{r}(n), where f(n) = A_{1}(n) and there is a function

σ : {0, …, p − 1} × {1, …, r} → {1, …, r} ,

such that, for each 1 ≤ i ≤ r and 0 ≤ α ≤ p − 1, we have the recurrence

A_{i}(pn + α) = A_{σ(α,i)}(n) .

We also have initial conditions

A_{i}(0) = a_{i} ,

for some a_{i} ∈ E , 1 ≤ i ≤ r.

Note: In the application to schemes for congruence properties of combinatorial sequences modulo prime powers p^{a}, treated in the present article, p will always be a prime, and the output set, E, would be

{0, 1, …, p^{a} − 1} .

Teaching the Computer How to Create Automatic p-schemes

All the tricks described above, in excruciating detail, for finding the scheme for determining the mod 2 behavior of the Catalan numbers

C_{n} := ConstantT ermOf [ (\frac{1}{x} + 2 + x)^{n}(1 − x) ] ,

can be taught to the computer (in our case using the symbolic programming language Maple), to find without human touch, an automatic p-scheme for determining the mod p^{a} behavior, for any prime p, and any power a, for any combinatorial sequence defined by

A(n) := ConstantT ermOf [ P (x_{1}, …, x_{m})^{n}Q(x_{1}, …, x_{n}) ] mod p^{a} ,

for any polynomials with integer coefficients, P (x_{1}, …, x_{m}) and Q(x_{1}, …, x_{m}), for any number of variables.

We will associate A(n) with the pair [P, Q].

We first rename A(n), A_{1}(n), and [P, Q], [P_{1}, Q_{1}]. We then try to find constant-term expressions for A_{1}(np), A_{1}(np + 1), …, A_{1}(np + p − 1). After using the multinomial theorem and doing it mod p^{a}, we would get, e.g.,

A_{1}(pn) = ConstantT ermOf [ P_{1}(x_{1}, …, x_{m})^{np}Q_{1}(x_{1}, …, x_{m}) ] mod p^{a}
= ConstantT ermOf [ (P_{1}(x_{1}, …, x_{m})^{p})^{n}Q_{1}(x_{1}, …, x_{m}) ] mod p^{a} ,

that after simplification (expanding, taking modulo p^{a}, and, if applicable, replacing x^{p} by x) will force us to put up with a brand-new discrete function, let’s call it A_{2}(n), given by

A_{2}(n) = ConstantT ermOf [ P_{2}(x_{1}, …, x_{m})^{n}Q_{2}(x_{1}, …, x_{m}) ] mod p^{a} ,

So A_{2} corresponds to a brand-new pair [P_{2}, Q_{2}]. We do likewise for A_{1}(pn + 1), all the way to A_{1}(pn + p − 1), getting (at the beginning) new pairs. Then we do the same for A_{2}(pn) through A_{2}(pn + p − 1). After awhile, by the pigeonhole principle, we will get old friends, and eventually there won’t be any ‘new guys’, and we get a finite (alas, often very large!) automatic p-scheme. The proof is as follows. If P (x) is a Laurent polynomial in x_{1}^{p}, …, x_{m}^{p}, let Λ(P (x)) denote the Laurent polynomial obtained by replacing each x_{j}^{p} by x_{j}. Since Λ commutes with raising to the pth power, the first component of each pair [P_{i}, Q_{i}] after a iterations is Λ^{k}(P (x)^{p^{a}}) for some k ≥ 0. The only terms of P (x)^{p^{a}} whose coefficients are non-zero modulo p^{a} are those in which the exponent of each x_{j} is a multiple of p; therefore k ≥ 1. It is not too difficult to see (for example, using Proposition 1.9 in [RY]) that

Λ (P (x)^{p^{a}}) ≡ P (x)^{p^{a−1}} (mod p^{a}).

From this it follows that

Λ^{k} (P (x)^{p^{a}}) ≡ Λ^{k−1} (P (x)^{p^{a−1}}) (mod p^{a}).

On the next iteration, we raise this polynomial to the pth power and apply Λ; this gives

Λ ((Λ^{k−1}(P (x)^{p^{a−1}}))^{p}) mod p^{a} = Λ^{k} (P (x)^{p^{a}}) mod p^{a} = Λ^{k−1} (P (x)^{p^{a−1}}) mod p^{a},

so the first component of [P_{i}, Q_{i}] stays the same after a iterations. There are only finitely many possibilities for the second component as well, since after the first component stabilizes then we can apply Λ to both P and (after deleting some terms) Q at each iteration, and this puts bounds on the degree and low-degree of Q.

All of this is implemented in AutoSquared by procedure CA for single-variable polynomials P and Q and by procedure CAmul for multivariate P and Q (of course, CAmul can handle also a single variable, but we kept CA both for old-time-sake and because it may be a bit faster for this special case).

The syntax is

CA(Z,x,p,a,K):  ,

where Z is a pair of single-variable functions [P, Q], in the variable x, p is a prime, a is a positive integer, and K is a (usually large) positive integer, stating the maximum number of ‘states’ (auxiliary functions) that you are willing to put up with. (It returns FAIL if the number of states exceeds K.) For example, to get an automatic 2-scheme for the Motzkin numbers, modulo 2, (if you are willing to tolerate a scheme with up to 30 members),

you type: gu:=CA([1/x+1+x,1-x**2],x,2,1,30):  .

The output (that we call gu) has two parts. The second part, gu[2], that is not needed for the application for the fast determination of the sequence modulo 2 (and in general modulo p^{a}) consists in the ‘definition’, in terms of constant term expressions A_{i}(n) := ConstantT erm[P_{i}(x)^{n}Q_{i}(x)], of the various auxiliary functions. So, in this example, gu[2] is

[[1/x+1+x, 1+x**2], [1/x+1+x, 1+x], [1/x+1+x, 1], [1/x+1+x, x]] ,

meaning that A_{1}(n) = ConstantT erm[(1/x+1+x)^{n}(1+x^{2})] , A_{2}(n) = ConstantT erm[(1/x+1+x)^{n}(1+x)]

A_{3}(n) = ConstantT erm[(1/x + 1 + x)^{n}] , A_{4}(n) = ConstantT erm[(1/x + 1 + x)^{n} · x)] . The more interesting part, the one needed for the actual fast computation, is gu[1]. Typing : lprint(gu[1]) in the same Maple session, gives

[[[2, 2], [3, 4], [3, 3], [0, 2]], [1, 1, 1, 0]] ,

that in humanese means the 2-scheme A_{1}(2n) = A_{2}(n) , A_{1}(2n + 1) = A_{2}(n) ,

A_{2}(2n) = A_{3}(n) , A_{2}(2n + 1) = A_{4}(n) ,
A_{3}(2n) = A_{3}(n) , A_{3}(2n + 1) = A_{3}(n) ,
A_{4}(2n) = 0 , A_{4}(2n + 1) = A_{2}(n) .

The initial conditions are A_{1}(0) = 1 , A_{2}(0) = 1 , A_{3}(0) = 1 , A_{4}(0) = 0 .

Moving right along, to get an automatic 2-scheme for the Motzkin numbers mod 4 (let’s tolerate from now on systems up to 10000 states):

gu:=CA([1/x+1+x,1-x**2],x,2,2,10000):  ,

getting (by typing nops(gu[2]) (or nops(gu[1][1]))) a scheme with 24 states. To get an automatic 2-scheme for the Motzkin numbers mod 8 (still with ≤ 10000 states, if possible), you type

gu:=CA([1/x+1+x,1-x**2],x,2,3,10000):  ,

getting a certain scheme with 128 states. For mod 16, we type

gu:=CA([1/x+1+x,1-x**2],x,2,4,10000):  ,

getting a certain scheme with 801 states. For mod 32, we type

gu:=CA([1/x+1+x,1-x**2],x,2,5,10000):  ,

getting a certain scheme with 5093 states. For mod 64, we type

gu:=CA([1/x+1+x,1-x**2],x,2,6,10000); ,

getting the output FAIL, meaning that the number of needed states exceeds our ‘cap’, 10000. Fast Evaluation mod p^{a} Once an automatic p-scheme, S, is found for a combinatorial sequence modulo p^{a}, AutoSquared can find very fast the N^{th} term of the sequence modulo p^{a}, for very large N, using the procedure EvalLS(Z,N,i,p), with i = 1. For example, after first finding an automatic 5-scheme for the Motzkin numbers modulo 25, by typing

gu:=CA([1/x+1+x,1-x**2],x,5,2,1000)[1]: ,

to get the remainder upon dividing M_{10^{100}} by 25, you should type:

EvalCA(gu,10**100,1,5);

getting 12. To get the first N terms of the sequence (modulo p^{a}), once a scheme, S, has been computed, type:

SeqCA(S,N,p);

For example, with the above scheme (that we called gu) (for the Motzkin numbers modulo 25)

SeqCA(gu,100000,5);

takes 2.36 seconds to give you the first 100000 terms, and getting the first million terms, by typing “SeqCA(gu,10**6,5);”, only takes 30 seconds.

Congruence Linear Schemes

The notion of automatic p-scheme defined above is conceptually attractive, since it can be modeled by a finite automaton with output. But, as can be seen by the above example, the number of ‘states’ (auxiliary functions) grows very fast. But note that the space of polynomials modulo p^{a} is a nice module over the ring ℤ/(p^{a}ℤ), and it is a shame to not take advantage of it. So rather than waiting until no new pairs [P (x), Q(x)] show up among the “children”, it may be a good idea, whenever a new pair comes along, to see whether it can be expressed as a linear combination of previously encountered pairs with the same P (x) (which we already know stays the same after a iterations, and only the Q(x)’s change).

One can get away with many fewer auxiliary functions (‘states’) with the following notion.

Definition: Let ℕ be the set of non-negative integers, and let p be a prime, a a positive integer, and let M be a module over the ring of integers modulo p^{a}, ℤ/(p^{a}ℤ). A linear p-scheme for a function f : ℕ → M is a set of finitely many (say r) auxiliary functions A_{1}(n), …, A_{r}(n), where f(n) = A_{1}(n), and such that for each i (1 ≤ i ≤ r), and each α (0 ≤ α < p), there exists a linear combination

A_{i}(pn + α) = ∑_{j=1}^{r} C_{i,j}^{(α)}A_{j}(n) ,

for some C_{i,j}^{(α)} ∈ {0, 1, …, p^{a} − 1}, and there are initial conditions:

A_{i}(0) = a_{i} .

Note that the previous notion of automatic p-scheme is the very special case, where for each α and i, there is exactly one j (that equals σ(α, i)) such that C_{i,j}^{(α)} is non-zero, and it has to be a 1.

Finding Linear p-Schemes in AutoSquared

This is implemented, in AutoSquared, by procedure LS for single-variable P and Q and by procedure LSmul for multivariate P and Q (of course, LSmul can handle also a single variable, and we kept LS both for old-time-sake and because it may be a bit faster for this special case).

The syntax for LS is

LS(Z,x,p,a,A,K):

where Z is a pair of single-variable functions [P, Q], x is the (single) variable name x that serves as the argument of P and Q, p is a prime, a is a positive integer, A is a symbol for expressing the linear expressions (where A[i] means our humanese A_{i}), and K is (usually fairly large) positive integer, stating the maximum number of ‘states’ (auxiliary functions) that you are willing to put up with. (It returns FAIL if the number of states exceeds K.)

For example, to get a Linear 2-scheme for the Motzkin numbers, modulo 2, (if you are willing to tolerate a scheme with up to 30 members), you type

gu:=LS([1/x+1+x,1-x**2],x,2,1,A,30):  .

getting

[[[A[2], A[2]], [A[3], A[4]], [A[3], A[3]], [0, A[2]]], [1, 1, 1, 0]]  ,

which is the same as the automatic 2-scheme, spelled-out above, except it is phrased more verbosely. If you type:

LS([1/x+1+x,1-x**2],x,2,2,A,30)[1];

you would get the following linear 2-scheme with 8 states:

[[[A[2], A[8]], [A[3], A[7]], [A[4], A[5]], [A[4], A[6]],
[A[4], 2*A[3]+2*A[4]+3*A[5]], [3*A[4], 2*A[3]+2*A[4]+A[5]], [A[3]+A[4], A[2]+A[3]+A[4]],
[A[3], 3*A[2]+A[3]+A[4]+3*A[5]]],
[1, 1, 1, 1, 1, 3, 2, 1]]  ,

that means that

A_{1}(2n) = A_{2}(n) , A_{1}(2n + 1) = A_{8}(n) , …,
A_{8}(2n) = A_{3}(n) , A_{8}(2n + 1) = 3A_{2}(n) + A_{3}(n) + A_{4}(n) + 3A_{5}(n) mod 4 .

The corresponding automatic 2-scheme has 24 states. For modulo 8 we get 18 states, compared to 128 for the automatic 2-scheme. For modulo 16 we get 43 states, compared to 801 states, and for modulo 32 we get 96 states, compared to 5093 states. Having gotten a scheme, S, phrased in terms of A, to get the first N terms of the sequence (modulo p^{a}), type

SeqLS(S,N,p,a,A) ;

Other Highlights of AutoSquared Procedures BinCA and BinLS find automatic p-schemes and linear p-schemes respectively for any binomial coefficient sum. See the on-line help. As mentioned at the beginning, there are quite a few sample input and output files linked to from the front of this article

http://www.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/meta.html  .

What about congruences modulo integers that are NOT primes or prime powers?

The Chinese Remainder Theorem comes to the rescue!

One first constructs as many automatic p-schemes, or linear p-schemes, for as many prime powers as one could afford, or care about, and then one can very fast find the congruence class modulo any integer involving these primes up to the given power.

The Maple packages CatalanLS, MotzkinLS, DelannoyLS

Using the main package AutoSquared, our computer precomputed schemes for quite a few prime powers, that enables us to find the remainder upon dividing by m, for many m, in particular, m = 1000, getting the last three digits of the Catalan, Motzkin, and Delannoy numbers given at the prologue.

See the on-line help in these packages.

Disclaimer

Both the automatic p-schemes and the linear p-schemes that our Maple package outputs are not guaranteed to be minimal. Of course the size does not change the fact that they run in logarithmic time in the input, but the ‘implied constants’ in the O(log n) algorithms are most probably not best-possible.

Conclusion

The present project is yet another case study in teaching computers to do research all by themselves, once they were taught (programmed) the human tricks. Once the computer mastered them, it can reproduce, in a few seconds, all the previous results accomplished by humans, and go on to output much deeper results, that no human, by himself, or herself, would be able to do, hence getting, much deeper results. So the fact that the last three decimal digits of M_{googol} are 187, may not be as interesting as Fermat’s Last Theorem, but is, in some sense, much deeper!

Acknowledgment: The second-named author (DZ) was supported in part by a grant from the National Science Foundation of the United States of America. The authors thank Shalosh B. Ekhad for its many diligent and tedious computations and proofs!

References

[KKM] Manuel Kauers, Christian Krattenthaler, and Thomas W. Müller, A method for determining the mod-2^{k} behaviour of recursive sequences, with applications to subgroup counting, Electron. J.

Combin. 18(2) (2012), Article P37. http://arxiv.org/abs/1107.2015

[KM1] Christian Krattenthaler and Thomas W. Müller,A method for determining the mod-3^{k} be-

haviour of recursive sequences, preprint. http://www.mat.univie.ac.at/~kratt/artikel/3psl2z.html

[KM2] Christian Krattenthaler and Thomas W. Müller, A Riccati differential equation and free subgroup numbers for lifts of P SL_{2}(ℤ) modulo powers of primes, J. Combin. Theory Ser. A 120

(2013), 2039–2063. http://arxiv.org/abs/1211.2947

[R] György E. Révész, “Introduction to Formal Languages”, Dover, 1991. [Originally published by McGraw-Hill, 1983].

[RY] Eric Rowland and Reem Yassawi, Automatic congruences for diagonals of rational functions

http://arxiv.org/abs/1310.8635 .

Eric Rowland, Université du Québec à Montréal, Montréal, Canada.

rowland at lacim dot ca

Doron Zeilberger, Mathematics Department, Rutgers University (New Brunswick), Piscataway, NJ 08854, USA.

zeilberg at math dot rutgers dot edu

Nov. 18, 2013