A Case Study in Meta-AUTOMATION: AUTOMATIC Generation of Congruence AUTOMATA For Combinatorial Sequences
aa r X i v : . [ m a t h . C O ] N ov 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 Rowlandand 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) determinationof congruence properties, modulo small (and not so small!) prime powers, for a wide class ofcombinatorial sequences. Even more interesting than the new results that could be obtained, is theillustrated methodology , that of designing ‘meta-algorithms’ that enable the computer to developalgorithms, 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 , andnumerous sample input and output files, that readers can use as templates for generating theirown, thereby proving many new ‘theorems’ about congruence properties of many famous (and, ofcourse, obscure) combinatorial sequences.
Very Important : This article is accompanied by the general Maple package ,and several other specific ones, and numerous input and output files that are obtainable, by one click, from the webpage (“front”) of this article .They could (and should!) be used as templates for generating as many input files that the humanwould 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 enoughto store them!But thanks to the Maple packages accompanying this article, we know for sure that the last threedigits 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 inthe main Maple package, AutoSquared . 1 ast 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 − = 0 , x = 1 . In order to compute x mod m , you don’t compute all the 10 previous terms, but use the implied recurrences x n = x n mod m , x n +1 = ax n mod m . This takes only log operations!What about sequences defined by higher-order recurrences, but still with constant coefficients? Forexample, what are the last three decimal digits of the googol-th Fibonacci number, F ? Youwould get the answer, 875, in 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 1elif n mod 2 = 0 then Fnm(1/2*n, m)*(Fnm(1/2*n + 1, m) + Fnm(1/2*n - 1, m)) mod melse Fnm(1/2*n - 1/2, m)**2 + Fnm(1/2*n + 1/2, m)**2 mod mfi:end:
It implements the (nonlinear) recurrence scheme F n = F n ( F n − + F n +1 ) , F n +1 = F n + F n +1 , F = 1 , F = 1 , and of course takes it modulo m at every step.Another way is to take the (1 ,
2) entry of the matrix (cid:18) (cid:19) mod 1000 , and use the ‘iterated-squaring’ trick applied to matrix (rather than scalar) exponentiation.2oth 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 con-stant coefficients (aka C -finite integer sequences).But what about sequences that are defined via linear recurrence equations with polynomial coef-ficients, 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¨uller developed a deep and in-genious group-theoretical method for the determination of holonomic sequences modulo powers of2. This has been extended to powers of 3 in [KM1], and further developed in [KM2].An important subclass of the class of holonomic integer sequences is the class of integer sequenceswhose (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 severalvariables, Rowland and Yassawi ([RY]) developed a very general method for computing finiteautomata for the fast computation (once the automaton is found, of course) of the congruencebehavior modulo prime powers. Of course, as the primes and/or their powers get larger, theautomata get larger too, but if the automaton is precomputed once and for all (and saved!), itis 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 sequencesand for specific moduli. We refer the reader to [RY] for many references, that we will not repeathere.
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 asingle variable or a multi-variable x = ( x , . . . , x m ), and ConstantTermOf means “coefficient of x ”, or “the coefficient of x · · · x m ”.For example, the arguably second-most famous combinatorial sequence (after the Fibonacci se-quence), is the sequence of the Catalan Numbers ( http://oeis.org/A000108 ), that may be definedby C n := ConstantT ermOf [ ( 1 x + 2 + x ) n (1 − x ) ] . Not as famous, but also popular, are the
Motzkin numbers ( http://oeis.org/A001006 ), that3ay be defined by M n := ConstantT ermOf [ ( 1 x + 1 + x ) n (1 − x ) ] , and also fairly famous are the Central Delannoy Numbers ( http://oeis.org/A001850 ), thatmay be defined by D n := ConstantT ermOf [ ( 1 x + 3 + 2 x ) n ] . So far, we got away with a single variable.Another celebrated sequence is the sequence of
Ap´ery Numbers , that were famously used by 64-year-old Roger Ap´ery (in 1978) to prove the irrationality of ζ (3). These are defined in terms of a binomial coefficient sum A ( n ) := n X k =0 (cid:18) nk (cid:19) (cid:18) n + kk (cid:19) . These may be equivalently defined (see below) as A ( n ) := ConstantT ermOf (cid:20)(cid:18) (1 + x ) (1 + x ) (1 + x ) (1 + x + x + x x + x x x ) x x x (cid:19) n (cid:21) . 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 ) = n X k =0 (cid:18) nk (cid:19) g k m Y i =1 (cid:18) a i n + b i k + c i d i n + e i k + f i (cid:19) , 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 in-troduce m variables x , . . . , x m and use the fact, that, by definition (cid:18) a i n + b i k + c i d i n + e i k + f i (cid:19) = ConstantT ermOf x i " (1 + x i ) a i n + b i k + c i x d i n + e i k + f i i Hence A ( n ) = n X k =0 (cid:18) nk (cid:19) g k m Y i =1 (cid:18) a i n + b i k + c i d i n + e i k + f i (cid:19) = n X k =0 (cid:18) nk (cid:19) g k m Y i =1 ConstantT ermOf x i [ (1 + x i ) a i n + b i k + c i x d i n + e i k + f i i ]= ConstantT ermOf x ,...,x m " n X k =0 (cid:18) nk (cid:19) g k m Y i =1 (1 + x i ) a i n + b i k + c i x d i n + e i k + f i i = ConstantT ermOf x ,...,x m " m Y i =1 (1 + x i ) a i n + c i x d i n + f i i ! n X k =0 (cid:18) nk (cid:19) g k m Y i =1 (1 + x i ) b i k x e i ki ! ConstantT ermOf x ,...,x m " m Y i =1 (1 + x i ) a i n + c i x d i n + f i i ! n X k =0 (cid:18) nk (cid:19) g k m Y i =1 (cid:18) (1 + x i ) b i x e i i (cid:19) k = ConstantT ermOf x ,...,x m " m Y i =1 (1 + x i ) a i n + c i x d i n + f i i ! g m Y i =1 (1 + x i ) b i x e i i ! n . = ConstantT ermOf x ,...,x m " m Y i =1 (1 + x i ) c i x f i i ! (1 + x i ) a i x d i i ! n g m Y i =1 (1 + x i ) b i x e i i ! n . = ConstantT ermOf x ,...,x m " m Y i =1 (1 + x i ) c i x f i i m Y i =1 (1 + x i ) a i x d i i ! g m Y i =1 (1 + x i ) b i x e i i !! n . This is implemented in procedure
BinToCT(L,x,a) in our Maple package
AutoSquared . For exam-ple, we got the above constant-term rendition of the Ap´ery 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 [ ( 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 ofcomputing C n modulo 2. So let’s define A ( n ) := C n mod 2 . Using the above formula for C n , and taking it modulo 2, we have: A ( n ) := ConstantT ermOf [ (1 + x )( 1 x + x ) n ] . We will try to find a constant-term expression for A (2 n ) . A (2 n ) = ConstantT ermOf [ (1+ x )( 1 x + x ) n ] mod 2 = ConstantT ermOf [ (1+ x )(( 1 x + x ) ) n ] mod 2= ConstantT ermOf [ (1+ x )( 1 x +2+ x ) n ] mod 2 = ConstantT ermOf [ (1+ x )( 1 x + x ) n ] mod 2But ConstantT ermOf [ (1 + x )( 1 x + x ) n ] = ConstantT ermOf [ 1 · ( 1 x + x ) n ] , since, obviously, ConstantT ermOf [ x · ( 1 x + x ) n ] = 0 . ConstantT ermOf [ 1 · ( 1 x + x ) n ] , only depends on x , we can replace x by x , implying that A (2 n ) = ConstantT ermOf [ ( 1 x + x ) n ] mod 2 . This forces us to put up with a new kid on the block , let’s call it A ( n ): A ( n ) := ConstantT ermOf [ ( 1 x + x ) n ] mod 2 , and we got the recurrence A (2 n ) = A ( n ) . We will handle A ( n ) in due course, but first let’s consider A (2 n + 1).We have A (2 n + 1) = ConstantT ermOf [ (1 + x )( 1 x + x ) n +1 ] mod 2= ConstantT ermOf [ (1 + x )( 1 x + x )(( 1 x + x ) ) n ] mod 2= ConstantT ermOf [ ( 1 x + x + 1 + x )( 1 x + 2 + x ) n ] mod 2= ConstantT ermOf [ ( 1 x + x + 1 + x )( 1 x + x ) n ] mod 2 . But
ConstantT ermOf [ ( 1 x + x + 1 + x )( 1 x + x ) n ] = ConstantT ermOf [ (1 + x ) · ( 1 x + x ) n ] , since, obviously ConstantT ermOf [ ( 1 x + x ) · ( 1 x + x ) n ] = 0 . Since the constant-termand of
ConstantT ermOf [ (1 + x ) · ( 1 x + x ) n ] , only depends on x , we can replace x by x , implying that A (2 n + 1) = ConstantT ermOf [ (1 + x )( 1 x + x ) n ] mod 2 . But this looks familiar! It is good-old A ( n ), so we have established, so far, that A (2 n ) = A ( n ) , A (2 n + 1) = A ( n ) . recurrence scheme , we need to handle A ( n ). A priori, this mayforce 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, thatwould enable the fast (logarithmic time) computation of our initial function A ( n ). We will seethat this would always be the case, no matter how complicated P ( x ) and Q ( x ) are (and evenwith many variables). Alas, as P ( x ) gets more complicated, the ‘finite’ gets bigger and bigger, soeventually the ‘logarithmic time’ in n would be impractical, since the implied constant would beeeeeeeeeeeeeenormous.But in this toy example , don’t worry! The ‘finitely many discrete functions’, is only two! As wewill shortly see, all we need is A ( n ), in addition to A ( n ).Recall that A ( n ) := ConstantT ermOf [ ( 1 x + x ) n ] mod 2 . Let’s try to find a constant-term expression for A (2 n ) . A (2 n ) = ConstantT ermOf [ ( 1 x + x ) n ] mod 2 = ConstantT ermOf [ (( 1 x + x ) ) n ] mod 2= ConstantT ermOf [ ( 1 x + 2 + x ) n ] mod 2 = ConstantT ermOf [ ( 1 x + x ) n ] mod 2Since the constant-termand only depends on x , we can replace x by x , implying that A (2 n ) = ConstantT ermOf [ ( 1 x + x ) n ] mod 2 . But that’s exactly A ( n ), so we have found out that A (2 n ) = A ( n ) . What about A (2 n + 1)? Here goes: A (2 n +1) = ConstantT ermOf [ ( 1 x + x ) n +1 ] mod 2 = ConstantT ermOf [ ( 1 x + x )(( 1 x + x ) ) n ] mod 2= ConstantT ermOf [ ( 1 x + x )( 1 x +2+ x ) n ] mod 2 = ConstantT ermOf [ ( 1 x + x )( 1 x + x ) n ] mod 2 . But the constant-termand now only has odd powers, so the coefficient of x , alias the constantterm , is 0. We have just established, the fast recurrence scheme : A (2 n ) = A ( n ) , A (2 n + 1) = A ( n ) ,A (2 n ) = A ( n ) , A (2 n + 1) = 0 , subject to the initial conditions A (0) = 1 , A (0) = 1 . much faster ) by the Maple package AutoSquared . Having downloaded 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 ≥ n = k X i =1 α i k − i , where α i ∈ { , } , α = 1, and it is abbreviated, in the positional notation, as a word , of length k , in the alphabet { , } α · · · α k . Phrased in terms of such ‘words’, the above scheme can be written, (where w is any word in thealphabet { , } ) A ( w
0) = A ( w ) , A ( w
1) = A ( w ) .A ( w
0) = A ( w ) , A ( w
1) = 0 , subject to the initial conditions (here φ is the empty word): A ( φ ) = 1 , A ( φ ) = 1 . Let’s revert to post-fix notation for representing functions, and omit the parentheses, i.e. write wA instead of A ( w ) and wA instead of A ( w ). This will not cause any ambiguity, since the alphabetof function names , { A , A } is disjoint from the alphabet of letters , { , } . The above schemebecomes w A = wA , w A = wA w A = wA , w A = 0 , subject to the initial conditions φA = 1 , φA = 1 . Let’s try to find A (30), alias, A (11110 ), alias, with our new convention, 11110 A . We get in two steps 11110 A = 1111 A = 0 . 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 onlywords that output 1 are those given by the regular expression ∗ . Equivalently, the only integers n for which the Catalan number C n is odd are those of the form n = 2 k − k = 0 , , , . . . .The words in the alphabet { , } that output 0 (i.e. those words that have at least one 0 in theirbinary representation) are the complement ‘language’, whose regular expression rendition is1 { , } ∗ { , } ∗ . What we have here is a finite automaton with output . The set of states is { A , A } while the alphabet is the set { , } . There are 2 directed edges coming out of each state, one for each letterof the alphabet, leading to another (possibly the same) state, or possibly to an output gate (in ourcase always 0, via ‘exit edges’ that prematurely end the journey. You have a starting state (in thisexample, A ) and an input word, and you travel along the automaton, according to the currentstate and the current rightmost letter, until you run out of letters, i.e. have the empty word, orwind-up in the output 0 prematurely, since some states have edges that lead directly to 0. (In ourexample when you are at state A 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 thefamous 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 aregular grammar describing the language (set of words) that yield that output.In this example, the set of non-terminal symbols is { A , A } and the set of terminal symbols is { , } . 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 is not needed (is superfluous), and the grammaris extremely simple A → φ , A → A . We leave it to the interested reader to write down the only slightly more complicated grammar forthe 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 areequivalent (as far as the generated languages), and there are easy algorithms for going betweenthem. 9hese are all very nice, but for the present formulation, it is more convenient not to write theinput 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 N be the set of non-negative integers, let p be a positive integer, and let E be anyset. An automatic p -scheme for a function f : N → E is a set of finitely many (say r ) auxiliaryfunctions A ( n ) , . . . , A r ( n ), where f ( n ) = A ( n ) and there is a function σ : { , . . . , p − } × { , . . . , r } → { , . . . , 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 , ≤ i ≤ r . Note : In the application to schemes for congruence properties of combinatorial sequences moduloprime powers p a , treated in the present article, p will always be a prime, and the output set, A ,would be { , , . . . , p a − } . Teaching the Computer How to Create Automatic p -schemes All the tricks described above, in excruciating detail, for finding the scheme for determining themod 2 behavior of the Catalan numbers C n := ConstantT ermOf [ ( 1 x + 2 + x ) n (1 − x ) ] , can be taught to the computer (in our case using the symbolic programming language Maple), tofind 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 , . . . , x m ) n Q ( x , . . . , x n ) ] mod p a , for any polynomials with integer coefficients, P ( x , . . . , x m ) and Q ( x , . . . , x m ), for any numberof variables.We will associate A ( n ) with the pair [ P, Q ].We first rename A ( n ), A ( n ), and [ P, Q ], [ P , Q ]. We then try to find constant-term expressionsfor A ( np ) , A ( np + 1) , . . . , A ( np + p − p a , we would get, e.g., A ( pn ) = ConstantT ermOf [ P ( x , . . . , x m ) np Q ( x , . . . , x m ) ] mod p a ConstantT ermOf [ ( P ( x , . . . , x m ) p ) n Q ( x , . . . , x m ) ] mod p a , that after simplification (expanding, taking modulo p a , and, if applicable, replacing x p by x ) willforce us to put up with a brand-new discrete function, let’s call it A ( n ), given by A ( n ) = ConstantT ermOf [ P ( x , . . . , x m ) n Q ( x , . . . , x m ) ] mod p a , So A corresponds to a brand-new pair [ P , Q ]. We do likewise for A ( pn + 1), all the way to A ( pn + p − A ( pn ) through A ( pn + p − pigeonhole principle , we will get old friends, and eventuallythere 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 p , . . . , x pm , let Λ( P ( x )) denote theLaurent polynomial obtained by replacing each x pj by x j . Since Λ commutes with raising to the p thpower, the first component of each pair [ P i , Q i ] after a iterations is Λ k ( P ( x ) p a ) for some k ≥
0. Theonly terms of P ( x ) p a whose coefficients are non-zero modulo p a are those in which the exponentof each x j is a multiple of p ; therefore k ≥
1. It is not too difficult to see (for example, usingProposition 1.9 in [RY]) that Λ (cid:16) P ( x ) p a (cid:17) ≡ P ( x ) p a − (mod p a ) . From this it follows that Λ k (cid:16) P ( x ) p a (cid:17) ≡ Λ k − (cid:16) P ( x ) p a − (cid:17) (mod p a ) . On the next iteration, we raise this polynomial to the p th power and apply Λ; this givesΛ (cid:16)(cid:16) Λ k − ( P ( x ) p a − ) (cid:17) p (cid:17) mod p a = Λ k (cid:16) P ( x ) p a (cid:17) mod p a = Λ k − (cid:16) P ( x ) p a − (cid:17) mod p a , so the first component of [ P i , Q i ] stays the same after a iterations. There are only finitely manypossibilities for the second component as well, since after the first component stabilizes then wecan apply Λ to both P and (after deleting some terms) Q at each iteration, and this puts boundson 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 singlevariable, but we kept CA both for old-time-sake and because it may be a bit faster for this specialcase).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 positiveinteger, and K is a (usually large) positive integer, stating the maximum number of ‘states’ (auxiliaryfunctions) that you are willing to put up with. (It returns FAIL if the number of states exceeds K .)11or example, to get an automatic 2-scheme for the Motzkin numbers, modulo 2, (if you are willingto 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 theapplication for the fast determination of the sequence modulo 2 (and in general modulo p a ) consistsin the ‘definition’, in terms of constant term expressions A i ( n ) := ConstantT erm [ P i ( x ) n Q i ( x )], ofthe 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 ( n ) = ConstantT erm [(1 /x + 1 + x ) n (1 + x )] , A ( n ) = ConstantT erm [(1 /x + 1 + x ) n (1 + x )] A ( n ) = ConstantT erm [(1 /x + 1 + x ) n ] , A ( 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 (2 n ) = A ( n ) , A (2 n + 1) = A ( n ) ,A (2 n ) = A ( n ) , A (2 n + 1) = A ( n ) ,A (2 n ) = A ( n ) , A (2 n + 1) = A ( n ) ,A (2 n ) = 0 , A (2 n + 1) = A ( n ) . The initial conditions are A (0) = 1 , A (0) = 1 , A (0) = 1 , A (0) = 0 . Moving right along, to get an automatic 2-scheme for the Motzkin numbers mod 4 (let’s toleratefrom 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 ≤ u:=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 theMotzkin 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 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 beencomputed, 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 seconds to give you the first 100000 terms, and getting the first million terms, by typing“
SeqCA(gu,10**6,5) ;”, only takes seconds. 13 ongruence Linear Schemes The notion of automatic p -scheme defined above is conceptually attractive, since it can be modeledby 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 Z / ( p a Z ), and it is a shame to not take advantage of it. So ratherthan 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 ofpreviously 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 N 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 , Z / ( p a Z ). A linear p -scheme for afunction f : N → M is a set of finitely many (say r ) auxiliary functions A ( n ) , . . . , A r ( n ), where f ( n ) = A ( n ), and such that for each i (1 ≤ i ≤ r ), and each α (0 ≤ α < p ), there exists a linearcombination A i ( pn + α ) = r X j =1 C ( α ) i,j A j ( n ) , for some C ( α ) i,j ∈ { , , . . . , p a − } , 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 servesas the argument of P and Q , p is a prime, a is a positive integer, A is a symbol for expressingthe linear expressions (where A[i] means our humanese A i ), and K is (usually fairly large) positiveinteger, stating the maximum number of ‘states’ (auxiliary functions) that you are willing to putup 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 totolerate a scheme with up to 30 members), 14ou 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 (2 n ) = A ( n ) , A (2 n + 1) = A ( n ) , . . . ,A (2 n ) = A ( n ) , A (2 n + 1) = 3 A ( n ) + A ( n ) + A ( n ) + 3 A ( 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 get43 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 fromthe front of this article .15 hat 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 powersas one could afford, or care about, and then one can very fast find the congruence class moduloany 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 primepowers, 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 atthe prologue.See the on-line help in these packages. Disclaimer
Both the automatic p -schemes and the linear p -schemes that our Maple package output are notguaranteed to be minimal. Of course the size does not change the fact that they run in logarithmictime 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 them-selves, 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 onto output much deeper results, that no human, by himself, or herself, would be able to do, hencegetting, much deeper results. So the fact that the last three decimal digits of M googol are 187, maynot 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 theNational Science Foundation of the United States of America. The authors thank Shalosh B. Ekhadfor its many diligent and tedious computations and proofs!
References [KKM] Manuel Kauers, Christian Krattenthaler, and Thomas W. M¨uller,
A method for determiningthe mod- k behaviour of recursive sequences, with applications to subgroup counting , Electron. J.Combin. (2) (2012), Article P37. http://arxiv.org/abs/1107.2015 [KM1] Christian Krattenthaler and Thomas W. M¨uller, A method for determining the mod- k be-haviour of recursive sequences , preprint. A Riccati differential equation and freesubgroup numbers for lifts of
P SL ( Z ) modulo powers of primes , J. Combin. Theory Ser. A (2013), 2039–2063. http://arxiv.org/abs/1211.2947 [R] Gy¨orgy E. R´ev´esz, “Introduction to Formal Languages , Dover, 1991. [Originally published byMcGraw-Hill, 1983].[RY] Eric Rowland and Reem Yassawi, Automatic congruences for diagonals of rational functions http://arxiv.org/abs/1310.8635 .Eric Rowland, Universit´e du Qu´ebec `a Montr´eal, Montr´eal, Canada. rowland at lacim dot ca
Doron Zeilberger, Mathematics Department, Rutgers University (New Brunswick), Piscataway, NJ08854, USA. zeilberg at math dot rutgers dot edu