A Meta-Algorithm for Creating Fast Algorithms for Counting ON Cells in Odd-Rule Cellular Automata

By Shalosh B. EKHAD, N. J. A. SLOANE, and Doron ZEILBERGER

Abstract: By using the methods of Rowland and Zeilberger (2014), we develop a meta-algorithm that, given a polynomial (in one or more variables), and a prime p, produces a fast (logarithmic time) algorithm that takes a positive integer n and outputs the number of times each residue class modulo p appears as a coefficient when the polynomial is raised to the power n and the coefficients are read modulo p. When p = 2, this has applications to counting the ON cells in certain “Odd-Rule” cellular automata. (This article is accompanied by a Maple package, CAcount, as well as numerous examples of input and output files, all of which can be obtained from the web page for this article:

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

Preface

The number of ON cells in the nth generation of an “Odd-Rule” cellular automaton is found by raising the defining polynomial (in which the number of variables is equal to the dimension of the ambient space) to the nth power, reading the coefficients modulo 2, and counting the remaining monomials—or equivalently, setting all the variables equal to 1 (see [Sl] for a detailed discussion).

The purpose of this article is to describe a meta-algorithm, inspired by a recent paper of Eric Rowland and Doron Zeilberger [RZ], that takes such a polynomial as input, and outputs a recurrence scheme that enables the fast (logarithmic time) computation of terms of the sequence giving the number of ON cells at time n. This provides an alternative, computer proof of Theorems 4 and 5 of [Sl].

A toy example

Following the Gelfand Principle, let’s illustrate the method with a simple example that can be done by hand. We will later describe how this method can be ‘taught’ to a computer, which will then be able to do far more complicated cases, impossible for humans.

Consider the sequence

a_{1}(n) := (1 + x + x^{2})^{n} mod 2 |_{x=1} ,

(sequence A071053 in [OEIS]), and suppose we wish to compute a_{1}(10^{100}), or a_{1}(n) for any very large n.

Of course, direct computation is hopeless, even if we reduce modulo 2 at each step and use the repeated squaring trick that makes RSA possible (P^{n} = (P^{n/2})^{2} if n is even, P^{n} = P P^{n−1} if n is odd), since the polynomials, before we set x = 1, are far too big for our modest universe. What we will do is adapt this trick so that we can also make the substitution x = 1 at intermediate steps.

First let’s try to relate a_{1}(2n) to a_{1}(n), using the Freshman’s Dream identity P (x)^{p} ≡ P (x^{p}) mod p:

a_{1}(2n) = (1 + x + x^{2})^{2n} mod 2 |_{x=1}= ((1 + x + x^{2})^{2})^{n} mod 2 |_{x=1}
= (1 + x^{2} + x^{4})^{n} mod 2 |_{x=1}= (1 + x + x^{2})^{n} mod 2 |_{x=1}  (EvenCase1)

(replacing x^{2} by x). Hence

a_{1}(2n) = a_{1}(n) .  (Recurrence1even)

Now we do the same thing for a_{1}(2n + 1):

a_{1}(2n + 1) = (1 + x + x^{2})^{2n+1} mod 2 |_{x=1}= (1 + x + x^{2}) ((1 + x + x^{2})^{2})^{n} mod 2 |_{x=1}
= (1 + x + x^{2}) (1 + x^{2} + x^{4})^{n} mod 2 |_{x=1}
= (1 + x^{2}) (1 + x^{2} + x^{4})^{n} mod 2 |_{x=1} + x (1 + x^{2} + x^{4})^{n} mod 2 |_{x=1} .  (OddCase1)

In the first term, once again, we can replace x^{2} by x, getting an uninvited guest, a_{2}(n), say:

a_{2}(n) := (1 + x) (1 + x + x^{2})^{n} mod 2 |_{x=1} .

As for the second term of Eq. (OddCase1), multiplying by x does not change anything, so this is equal to (1 + x^{2} + x^{4})^{n} mod 2 |_{x=1}, which, again replacing x^{2} by x, is our old friend a_{1}(n). Hence

a_{1}(2n + 1) = a_{2}(n) + a_{1}(n) .  (Recurrence1odd)

But this pair of recurrences is useless unless we can handle a_{2}(n). So let’s try the same technique on it. A priori, this may force us to introduce terms a_{3}(n), a_{4}(n), etc., and lead us into an infinite regression, also known as a Ponzi scheme, but let’s hope for the best.

Again we start with a_{2}(2n). Using the Freshman’s Dream, and the fact that multiplying a polynomial by x (or any other monomial) does not affect the result if we are going to read it modulo 2 and set x = 1, we have

a_{2}(2n) = (1 + x) (1 + x + x^{2})^{2n} mod 2 |_{x=1} = (1 + x) · ((1 + x + x^{2})^{2})^{n} mod 2 |_{x=1}

= (1 + x) · (1 + x^{2} + x^{4})^{n} mod 2 |_{x=1} = 1 · (1 + x^{2} + x^{4})^{n} mod 2 |_{x=1} + x (1 + x^{2} + x^{4})^{n} mod 2 |_{x=1}

= 2 (1 + x^{2} + x^{4})^{n} mod 2 |_{x=1} = 2 (1 + x + x^{2})^{n} mod 2 |_{x=1} = 2a_{1}(n) .

Hence

a_{2}(2n) = 2a_{1}(n) .  (Recurrence2even)

Now for a_{2}(2n + 1). We have

a_{2}(2n + 1) = (1 + x) · (1 + x + x^{2})^{2n+1} mod 2 |_{x=1}

= ((1 + x) · (1 + x + x^{2})) · ((1 + x + x^{2})^{2})^{n} mod 2 |_{x=1}
= (1 + 2x + 2x^{2} + x^{3}) · (1 + x^{2} + x^{4})^{n} mod 2 |_{x=1}
= (1 + x^{3}) · (1 + x^{2} + x^{4})^{n} mod 2 |_{x=1}
= 1 · (1 + x^{2} + x^{4})^{n} mod 2 |_{x=1} + x^{3} · (1 + x^{2} + x^{4})^{n} mod 2 |_{x=1}
= (1 + x^{2} + x^{4})^{n} mod 2 |_{x=1} + (1 + x^{2} + x^{4})^{n} mod 2 |_{x=1}
= (1 + x + x^{2})^{n} mod 2 |_{x=1} + (1 + x + x^{2})^{n} mod 2 |_{x=1}= 2a_{1}(n) .

Hence

a_{2}(2n + 1) = 2a_{1}(n) .  (Recurrence2odd)

So the uninvited guest, a_{2}(n), did not invite further guests, and now we have a super-fast way to compute a_{1}(n) for large n, using the system

a_{1}(2n) = a_{1}(n) , a_{1}(2n + 1) = a_{1}(n) + a_{2}(n) ;
a_{2}(2n) = 2a_{1}(n) , a_{2}(2n + 1) = 2a_{1}(n) .  (System)

For certain “odd-rule” cellular automata, the sequence a_{1}(n), n ≥ 0 is completely determined by the subsequence b_{1}(k) := a_{1}(2^{k} − 1), k ≥ 0 [Sl], and the b_{1}(k), unlike the a_{1}(n), often have simple generating functions, which we can derive (rigorously) by these methods. With a_{1}(n) as defined above, let

f_{1}(t) := ∑_{k=0}^{∞} b_{1}(k)t^{k}

be the generating function for b_{1}(k), and similarly define b_{2}(k) := a_{2}(2^{k} − 1) and

f_{2}(t) := ∑_{k=0}^{∞} b_{2}(k)t^{k} .

From Eq. (System), we have

b_{1}(k) = b_{1}(k − 1) + b_{2}(k − 1) , b_{2}(k) = 2b_{1}(k − 1) ,

and since by direct computation, b_{1}(0) = 1, b_{2}(0) = 2, we arrive at a system of two linear equations for the unknowns f_{1}(t) and f_{2}(t):

{ f_{1}(t) = 1 + tf_{1}(t) + tf_{2}(t) , f_{2}(t) = 2 + 2tf_{1}(t) } ,

whose solution is

f_{1}(t) = \frac{1 + 2t}{(1 + t)(1 − 2t)} , f_{2}(t) = \frac{2}{(1 + t)(1 − 2t)}

(A001045, A014113 in [OEIS]). But we really don’t care about f_{2}(t), we just needed it in order to find f_{1}(t), so now we can safely discard it, and get the

Theorem:

f_{1}(t) = \frac{1 + 2t}{(1 + t)(1 − 2t)} .

The general case

Fix once and for all a prime p and a polynomial P = P (x_{1}, …, x_{k}) ∈ Z[x_{1}, …, x_{k}]. If A(x_{1}, …, x_{k}) is any element of Z[x_{1}, …, x_{k}], we define the functional

A(x_{1}, …, x_{k}) → A(x_{1}, …, x_{k}) mod p|_{x_{1}=1,...,x_{k}=1}  (Reduce)

to mean “expand A(x_{1}, …, x_{k}) as a sum of monomials, reduce the coefficients modulo p to one of the numbers {0, 1, …, p − 1} ∈ Z, and finally set all the variables x_{i} equal to 1”.

For any polynomial Q = Q(x_{1}, …, x_{k}) ∈ Z[x_{1}, …, x_{k}] whose degree in each of the variables is less than p, define

a_{Q}(n) := QP^{n} mod p|_{x_{1}=1,...,x_{k}=1} .

For 0 ≤ i < p, we have

a_{Q}(pn + i) = Q(x_{1}, …, x_{k})P (x_{1}, …, x_{k})^{pn+i} mod p |_{x_{1}=1,...,x_{k}=1}
= [ Q(x_{1}, …, x_{k})P (x_{1}, …, x_{k})^{i}]P (x_{1}, …, x_{k})^{np} mod p |_{x_{1}=1,...,x_{k}=1}
= [ Q(x_{1}, …, x_{k})P (x_{1}, …, x_{k})^{i}](P (x_{1}, …, x_{k})^{p})^{n} mod p |_{x_{1}=1,...,x_{k}=1}
= [ Q(x_{1}, …, x_{k})P (x_{1}, …, x_{k})^{i}]P (x_{1}^{p}, …, x_{k}^{p})^{n} mod p |_{x_{1}=1,...,x_{k}=1} .

Now write

Q(x_{1}, …, x_{k})P (x_{1}, …, x_{k})^{i} mod p = ∑_{(α_{1},...,α_{k})∈{0,...,p−1}^{k}} x_{1}^{α_{1}} ⋯ x_{k}^{α_{k}}R_{(α_{1},...,α_{k})}(x_{1}^{p}, …, x_{k}^{p}) .

(Here again “mod p” applies just to the coefficients, not the variables.) Hence

a_{Q}(np+i) = ∑_{(α_{1},...,α_{k})∈{0,...,p−1}^{k}} x_{1}^{α_{1}} ⋯ x_{k}^{α_{k}}R_{(α_{1},...,α_{k})}(x_{1}^{p}, …, x_{k}^{p})P (x_{1}^{p}, …, x_{k}^{p})^{n} mod p |_{x_{1}=1,...,x_{k}=1}

= ∑_{(α_{1},...,α_{k})∈{0,...,p−1}^{k}} R_{(α_{1},...,α_{k})}(x_{1}^{p}, …, x_{k}^{p})P (x_{1}^{p}, …, x_{k}^{p})^{n} mod p |_{x_{1}=1,...,x_{k}=1} ,

= ∑_{(α_{1},...,α_{k})∈{0,...,p−1}^{k}} R_{(α_{1},...,α_{k})}(x_{1}, …, x_{k})P (x_{1}, …, x_{k})^{n} mod p |_{x_{1}=1,...,x_{k}=1} ,

= ∑_{(α_{1},...,α_{k})∈{0,...,p−1}^{k}} a_{R_{(α_{1},...,α_{k})}}(n) .

In other words for any Q(x_{1}, …, x_{k}) and each of the residue classes i, 0 ≤ i ≤ p − 1, we can find a multiset of polynomials, let’s call it S_{i}(Q), such that

a_{Q}(np + i) = ∑_{R∈S_{i}(Q)} a_{R}(n) .

We really only care about the case Q = 1, but the algebra forces us to consider other Q’s, and they in turn force us to treat still other Q’s, and so on. However, by the pigeon-hole principle, this process must terminate, and we obtain a finite recurrence scheme, containing say m equations. Placing all the Q’s that appear into some arbitrary order, with Q_{1} = 1, we get a (logarithmic-time) recurrence scheme:

a_{j}(np + i) = ∑_{l∈S_{i}(j)} a_{l}(n) ,

for 1 ≤ j ≤ m, that enables the fast calculation of a_{1}(n) for any n.

Furthermore, by focusing only on i = p−1, and defining c_{j}(k) := a_{j}(p^{k} −1), we have, for 1 ≤ j ≤ m,

c_{j}(k) = ∑_{l∈S_{p−1}(j)} c_{l}(k − 1) .

Define the generating functions

f_{j}(t) := ∑_{k=0}^{∞} c_{j}(k)t^{k} (1 ≤ j ≤ m) .

Standard manipulations of generating functions convert the above recurrences into a system of m linear equations for the m unknowns f_{1}(t), …, f_{m}(t):

f_{j}(t) = c_{j}(0) + t ∑_{l∈S_{p−1}(j)} f_{l}(t) , 1 ≤ j ≤ m ,

that can be solved, at least in principle, yielding rigorous explicit expressions for all the f_{j}(t), and in particular for f_{1}(t), the one in which we are most interested. Note that this proves that the generating function, f_{1}(t), is always a rational function. If m is too large, and the system of equations cannot be solved, then one may try to use the recurrences to generate sufficiently many terms of the sequence c_{1}(k), and then guess the rational function f_{1}(t), using for example the Maple packgage gfun [SaZ]. It may then be possible to justify that guess, a posteriori, by finding upper bounds on the degree of the generating function.

Keeping track of the individual coefficients

If instead of the functional Eq. (Reduce), one uses, for some formal variables s_{1}, …, s_{p−1},

∑_{α} c_{α}x^{α} → ∑_{α} s_{c_{α}} ,

one can modify the above arguments and keep track of the number of occurrences of each i (i = 1, …, p − 1) as coefficients in the expansion of P (x_{1}, …, x_{k})^{n} mod p. The Maple package CAcount Everything discussed above is implemented in the Maple package CAcount, which can be downloaded from the web page for this article:

http://www.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/CAcount.html  , where there

are also many samples of input and output files that readers can use as templates for further computations. To see the list of the main procedures, type

ezra(); ,

or to see the list of procedures that handle the more refined case, where one keeps track of the individual coefficients (only useful for p > 2), type

ezraG(); .

To get instructions on using a particular procedure, type

ezra(ProcedureName); .

For example. procedure CAaut finds the recurrence ‘automaton’, and to get help with it, type

ezra(CAaut); .

For our toy example, type

CAaut([1+x+x**2,1],[x],2,2);

which produces as output the pair

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

where the first component,

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

is Maple’s way of encoding the recurrence a_{1}(2n) = a_{1}(n) , a_{1}(2n+1) = a_{2}(n)+a_{1}(n) ; a_{2}(2n) = a_{1}(n)+a_{1}(n) , a_{2}(2n+1) = a_{1}(n)+a_{1}(n) .

The second component

[1, 2]

is Maple’s way of encoding the initial conditions

a_{1}(1) = 1 , a_{2}(1) = 2 .

Procedure SeqF uses the scheme, once found, to compute as many terms as desired, while procedure ARLT (for anti-run-length-transform, see [Sl]) computes the sparse subsequence in the places p^{i} − 1. Procedure GFsP finds the proved generating function for that subsequence, and if the size of the system is too big, GFsG guesses it faster, and as we mentioned above, the guess can be justified a posteriori.

References

[3by3] Shalosh B. Ekhad, N. J. A. Sloane, and Doron Zeilberger, “Odd-Rule” Cellular Automata on the Square Grid, in preparation, March 2015.

[OEIS] The OEIS Foundation Inc., The On-Line Encyclopedia of Integer Sequences, https://oeis.org.

[RZ] Eric Rowland and Doron Zeilberger, A Case Study in Meta-AUTOMATION: AUTOMATIC Generation of Congruence AUTOMATA For Combinatorial Sequences, J. Difference Equations and Applications 20 (2014), 973–988;

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

[SaZ] Bruno Salvy and Paul Zimmermann, GFUN: a Maple package for the Manipulation of Generating and Holonomic Functions in One Variable, ACM Trans. Math. Software 20 (1994), 163–177.

[Sl] N. J. A. Sloane On the Number of ON Cells in Cellular Automata, 2015;

http://arxiv.org/abs/1503.01168.

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

N. J. A. Sloane, The OEIS Foundation Inc, 11 South Adelaide Ave, Highland Park, NJ 08904, USA, and Department of Mathematics, Rutgers University (New Brunswick); njasloane at gmail dot com ; http://neilsloane.com/ .

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

zeilberg at math dot rutgers dot edu ; http://www.math.rutgers.edu/~zeilberg/ .

Published in The Personal Journal of Shalosh B. Ekhad and Doron Zeilberger

(http://www.math.rutgers.edu/~zeilberg/pj.html), N. J. A. Sloane’s home page
(http://neilsloane.com/), and arxiv.org .

Mar 05, 2015