Computing the Cycle Structure of Finite Linear Systems

(1)

IFAC PapersOnLine 53-2 (2020) 4316–4321

ScienceDirect

Peer review under responsibility of International Federation of Automatic Control.

10.1016/j.ifacol.2020.12.2487

10.1016/j.ifacol.2020.12.2487 2405-8963

Computing the Cycle Structure of Finite Linear Systems

Eva Zerz^∗ and Hermann Giese^∗

∗Lehrstuhl D f¨ur Mathematik, RWTH Aachen University, 52062 Aachen, Germany (e-mail: eva.zerz@math.rwth-aachen.de)

Abstract:Consider a linear difference equation with constant coefficients in the ring of integers modulom. If the leading coefficient and the constant term are both units, then all trajectories are (purely) periodic. Moreover, the finite state set can be decomposed into disjoint cycles of various lengths. The following problems will be addressed: computing the cycle partition and determining the period w.r.t. a specific initial state. The latter question can often be reduced to calculating the order of an invertible matrix. If the prime factorization ofm is known, then it suffices to consider prime powers, by the Chinese remainder theorem. For primes, an efficient algorithm due to Leedham-Green may be used, which is available in group-theoretic computer algebra systems such as Magma or GAP. This approach will be extended to prime powers.

Finally, we will discuss how to relax the assumptions guaranteeing periodicity.

Keywords: Linear systems, algebraic systems theory, cycle length, periodicity, finite rings.

1. INTRODUCTION The Fibonacci equation

y(t+ 2) =y(t+ 1) +y(t),

where t is a nonnegative integer, is usually considered in characteristic zero. Choosing y(0) = 0 and y(1) = 1 as initial values, one obtains the famous sequence

y= (0,1,1,2,3,5,8,13,21,34, . . .).

In positive characteristic however, the sequence becomes (purely) periodic, that is, it returns to its initial values and thus runs into a loop.

Example: In characteristic 2, we obtain y= (0,1,1,0,1,1,0,1,1,0,1, . . .) and the length of the period is 3. In characteristic 3,

y= (0,1,1,2,0,2,2,1,0,1, . . .)

and the length of the period is 8. In characteristic 4, y= (0,1,1,2,3,1,0,1, . . .)

and the length of the period is 6.

The sequence of period lengths is known as the Pisano sequence

π= (1,3,8,6,20,24,16,12,24,60, . . .).

Its properties have been studied by several authors, see e.g.

Wall (1960), and it is popular in recreational mathematics.

However, it is also related to “serious” mathematical problems. A Wall-Sun-Sun prime is a prime p such that π(p) = π(p²). It has been checked experimentally that such a prime (if it exists) must be greater than 10¹⁴. On the other hand, we have π(n) =π(n²) forn∈ {6,12}. A proof of the nonexistence of Wall-Sun-Sun primes would yield an alternative proof of the first case of Fermat’s last theorem, see Sun and Sun (1992).

This work was supported by DFG-SFB/TRR 195.

In the present paper, we will study general monic difference equations of ordern, that is,

y(t+n) +a_n−1y(t+n−1) +. . .+a1y(t+ 1) +a0y(t) = 0 with ai ∈ Z. The case where all solutions modulo m are (purely) periodic will be characterized by the condition gcd(a0, m) = 1. We will then turn to computing the period of the solution evolving from specific initial values y(0), . . . , y(n−1). Next, we study when this period coin- cides with the (group-theoretic) order of the companion matrixAof the polynomialsⁿ+n−1

i=0 aisⁱ, ifAis considered as an element ofZⁿm^×ⁿ, whereZ^m:=Z/mZ. Note that A∈Zⁿm^×ⁿ is invertible if and only if gcd(det(A), m) = 1.

Due to the companion matrix structure, det(A)∈ {±a0}. The state setX := Zⁿm can be decomposed into disjoint cycles and we will show how to compute this partition.

Finally, the assumption of monicity (an = 1) will be relaxed.

2. MONIC EQUATIONS

LetZdenote the ring of integers, and letNdenote the set of nonnegative integers. Letnbe a positive integer, and let a0, . . . , an−1 ∈Z be given. Consider the linear difference equation

y(t+n)+an−1y(t+n−1)+. . .+a1y(t+1)+a0y(t) = 0, (1) wheret∈N. Lety0, . . . , y_n−1∈Zbe given. Together with the initial condition

y(0) =y0, . . . , y(n−1) =y_n−1 (2) the initial value problem (1), (2) has a unique solution (y(t))_t∈N ∈ Z^N. For each integer m > 1, one obtains an induced sequence (y(t))t∈N ∈ Z^Nm, where Z^m = Z/mZ denotes the ring of integers modulom. In the following, we will make no notational distinction between an integer and its residue class modulom, or between a polynomial inZ[s] and its residue class inZ^m[s].

Computing the Cycle Structure of Finite Linear Systems

y(t+ 2) =y(t+ 1) +y(t),

y= (0,1,1,2,3,5,8,13,21,34, . . .).

y= (0,1,1,2,0,2,2,1,0,1, . . .)

π= (1,3,8,6,20,24,16,12,24,60, . . .).

y(t+n) +an−1y(t+n−1) +. . .+a1y(t+ 1) +a0y(t) = 0 with ai ∈ Z. The case where all solutions modulo m are (purely) periodic will be characterized by the condition gcd(a0, m) = 1. We will then turn to computing the period of the solution evolving from specific initial values y(0), . . . , y(n−1). Next, we study when this period coin- cides with the (group-theoretic) order of the companion matrixAof the polynomialsⁿ+n−1

i=0 aisⁱ, ifAis considered as an element ofZⁿm^×ⁿ, whereZ^m:=Z/mZ. Note that A∈Z^n×nm is invertible if and only if gcd(det(A), m) = 1.

2. MONIC EQUATIONS

y(t+n)+an−1y(t+n−1)+. . .+a1y(t+1)+a0y(t) = 0, (1) wheret∈N. Lety0, . . . , yn−1∈Zbe given. Together with the initial condition

y(0) =y0, . . . , y(n−1) =yn−1 (2) the initial value problem (1), (2) has a unique solution (y(t))t∈N ∈ Z^N. For each integer m > 1, one obtains an induced sequence (y(t))_t∈N ∈ Z^Nm, where Z^m = Z/mZ denotes the ring of integers modulom. In the following, we will make no notational distinction between an integer and its residue class modulom, or between a polynomial inZ[s] and its residue class inZm[s].

Computing the Cycle Structure of Finite Linear Systems

y(t+ 2) =y(t+ 1) +y(t),

y= (0,1,1,2,3,5,8,13,21,34, . . .).

y= (0,1,1,2,0,2,2,1,0,1, . . .)

π= (1,3,8,6,20,24,16,12,24,60, . . .).

2. MONIC EQUATIONS

Computing the Cycle Structure of Finite Linear Systems

y(t+ 2) =y(t+ 1) +y(t),

y= (0,1,1,2,3,5,8,13,21,34, . . .).

y= (0,1,1,2,0,2,2,1,0,1, . . .)

π= (1,3,8,6,20,24,16,12,24,60, . . .).

y(t+n) +an−1y(t+n−1) +. . .+a1y(t+ 1) +a0y(t) = 0 with ai ∈ Z. The case where all solutions modulo m are (purely) periodic will be characterized by the condition gcd(a0, m) = 1. We will then turn to computing the period of the solution evolving from specific initial values y(0), . . . , y(n−1). Next, we study when this period coin- cides with the (group-theoretic) order of the companion matrixAof the polynomialsⁿ+n−1

i=0 aisⁱ, ifAis considered as an element ofZⁿm^×ⁿ, whereZ^m:=Z/mZ. Note that A∈Z^n×nm is invertible if and only if gcd(det(A), m) = 1.

2. MONIC EQUATIONS

y(0) =y0, . . . , y(n−1) =yn−1 (2) the initial value problem (1), (2) has a unique solution (y(t))t∈N ∈ Z^N. For each integer m > 1, one obtains an induced sequence (y(t))_t∈N ∈ Z^Nm, where Z^m = Z/mZ denotes the ring of integers modulom. In the following, we will make no notational distinction between an integer and its residue class modulom, or between a polynomial inZ[s] and its residue class inZm[s].

Computing the Cycle Structure of Finite Linear Systems

y(t+ 2) =y(t+ 1) +y(t),

y= (0,1,1,2,3,5,8,13,21,34, . . .).

y= (0,1,1,2,0,2,2,1,0,1, . . .)

π= (1,3,8,6,20,24,16,12,24,60, . . .).

y(t+n) +an−1y(t+n−1) +. . .+a1y(t+ 1) +a0y(t) = 0 with ai ∈ Z. The case where all solutions modulo m are (purely) periodic will be characterized by the condition gcd(a0, m) = 1. We will then turn to computing the period of the solution evolving from specific initial values y(0), . . . , y(n−1). Next, we study when this period coin- cides with the (group-theoretic) order of the companion matrixAof the polynomialsⁿ+n−1

i=0 aisⁱ, ifAis considered as an element ofZ^n×nm , whereZ^m:=Z/mZ. Note that A∈Zⁿm^×ⁿ is invertible if and only if gcd(det(A), m) = 1.

2. MONIC EQUATIONS

LetZdenote the ring of integers, and letNdenote the set of nonnegative integers. Letnbe a positive integer, and let a0, . . . , a_n−1 ∈Z be given. Consider the linear difference equation

y(0) =y0, . . . , y(n−1) =yn−1 (2) the initial value problem (1), (2) has a unique solution (y(t))t∈N ∈ Z^N. For each integer m > 1, one obtains an induced sequence (y(t))t∈N ∈ Z^Nm, where Z^m = Z/mZ denotes the ring of integers modulom. In the following, we will make no notational distinction between an integer and its residue class modulom, or between a polynomial inZ[s] and its residue class inZ^m[s].

Computing the Cycle Structure of Finite Linear Systems

y(t+ 2) =y(t+ 1) +y(t),

y= (0,1,1,2,3,5,8,13,21,34, . . .).

y= (0,1,1,2,0,2,2,1,0,1, . . .)

π= (1,3,8,6,20,24,16,12,24,60, . . .).

2. MONIC EQUATIONS

Computing the Cycle Structure of Finite Linear Systems

y(t+ 2) =y(t+ 1) +y(t),

y= (0,1,1,2,3,5,8,13,21,34, . . .).

y= (0,1,1,2,0,2,2,1,0,1, . . .)

π= (1,3,8,6,20,24,16,12,24,60, . . .).

i=0 aisⁱ, ifAis considered as an element ofZⁿm^×ⁿ, whereZm:=Z/mZ. Note that A∈Zⁿm^×ⁿ is invertible if and only if gcd(det(A), m) = 1.

2. MONIC EQUATIONS

LetZdenote the ring of integers, and letNdenote the set of nonnegative integers. Letnbe a positive integer, and let a0, . . . , a_n−1 ∈Z be given. Consider the linear difference equation

y(t+n)+a_n−1y(t+n−1)+. . .+a1y(t+1)+a0y(t) = 0, (1) wheret∈N. Lety0, . . . , yn−1∈Zbe given. Together with the initial condition

y(0) =y0, . . . , y(n−1) =yn−1 (2) the initial value problem (1), (2) has a unique solution (y(t))t∈N ∈ Z^N. For each integer m > 1, one obtains an induced sequence (y(t))t∈N ∈ Z^Nm, where Z^m = Z/mZ denotes the ring of integers modulom. In the following, we will make no notational distinction between an integer and its residue class modulom, or between a polynomial inZ[s] and its residue class inZ^m[s].

Computing the Cycle Structure of Finite Linear Systems

y(t+ 2) =y(t+ 1) +y(t),

y= (0,1,1,2,3,5,8,13,21,34, . . .).

y= (0,1,1,2,0,2,2,1,0,1, . . .)

π= (1,3,8,6,20,24,16,12,24,60, . . .).

2. MONIC EQUATIONS

y(0) =y0, . . . , y(n−1) =yn−1 (2) the initial value problem (1), (2) has a unique solution (y(t))_t∈N ∈ Z^N. For each integer m > 1, one obtains an induced sequence (y(t))t∈N ∈ Z^Nm, where Z^m = Z/mZ denotes the ring of integers modulom. In the following, we will make no notational distinction between an integer and its residue class modulom, or between a polynomial inZ[s] and its residue class inZ^m[s].

Computing the Cycle Structure of Finite Linear Systems

y(t+ 2) =y(t+ 1) +y(t),

y= (0,1,1,2,3,5,8,13,21,34, . . .).

y= (0,1,1,2,0,2,2,1,0,1, . . .)

π= (1,3,8,6,20,24,16,12,24,60, . . .).

2. MONIC EQUATIONS

y(0) =y0, . . . , y(n−1) =y_n−1 (2) the initial value problem (1), (2) has a unique solution (y(t))t∈N ∈ Z^N. For each integer m > 1, one obtains an induced sequence (y(t))t∈N ∈ Z^Nm, where Z^m = Z/mZ denotes the ring of integers modulom. In the following, we will make no notational distinction between an integer and its residue class modulom, or between a polynomial inZ[s] and its residue class inZ^m[s].

Theorem 1. For each m > 0, there exist integers k ≥ 0 andp >0 such that

y(t+p) =y(t) for allt≥k. If gcd(a0, m) = 1, thenkcan be taken to be zero.

Proof: Define

x(t) := [y(t), . . . , y(t+n−1)]^T. Thenx(t+ 1) =Ax(t), where

A=







0 1

... . ..

0 1

−a0 −a1 . . . −an−1







is the companion matrix of the polynomial χA = sⁿ + n−1

i=0 aisⁱ. Since the ringZ^mis finite, the sequencey, and hence alsox, can only take finitely many values. Thus there exists k ≥ 0 andp > 0 such that x(k+p) = x(k). This implies thatx(t+p) =x(t) for all t≥k.

If gcd(a0, m) = 1, then a0 is a unit in Z^m. This implies that the polynomial indeterminatesis a unit in the finite ring R := Zm[s]/χA. Hence it has finite order, that is, there existsp >0 such thats^p= 1 holds inR. This means that the polynomial s^p−1 is a multiple of χA in Z^m[s]

and thusy(t+p) =y(t) holds for allt∈N. Givenmwith gcd(a0, m) = 1, we will address the compu- tation of the smallest numberp >0 such that

(1) x(p) =x(0). Since this implies thatx(t+p) =x(t) for allt ∈N, this number will be called the period of x with respect to the initial valuex0:= [y0, . . . , yn−1]^T. It is the smallest p >0 with A^px0 =x0 and will be denoted byπm(x0).

(2) A^p = I. This p is the order of the invertible matrix A∈Zⁿm^×ⁿ, denoted by ordm(A).

We observe that πm(x0) must be a divisor of ordm(A).

Moreover,

ordm(A) = lcm{πm(x0)| x0∈Zⁿm}.

By the Chinese remainder theorem, it suffices to consider the case wheremis a prime power.

Example: Consider

y(t+3) =y(t+2)+y(t+1)+y(t), y(0) =y(1) =y(2) = 1.

• Modulo 2, we obtain the constant sequence 1, and henceπ2([1,1,1]^T) = 1, but ord2(A) = 4. Indeed, we have

x^T₀ π2(x0) [0,0,0] 1 [0,0,1] 4 [0,1,0] 2 [0,1,1] 4 [1,0,0] 4 [1,0,1] 2 [1,1,0] 4 [1,1,1] 1.

Adopting the notation used by Deng (2015), the cycle structure isX = Z³2 = 2C1∪C2∪C4. This means that there are two fixed points (i.e., cycles of length one), one cycle of length two, and one cycle of length four, corresponding to a partition of the state set into cycles according to 8=1+1+2+4.

• Modulo 3, we obtain a sequence of period 13, that is, π3([1,1,1]^T) = 13, and we also have ord3(A) = 13. Indeed, all x0= 0 yieldπ3(x0) = 13.

• Modulo 4, we obtain a sequence of period 4, that is, π4([1,1,1]^T) = 4, but ord4(A) = 8.

Theorem 2. Given x0 = [y0, . . . , yn−1]^T ∈ Zⁿ, compute yn, . . . , y2n−2 ∈Zaccording to (1) and define the Hankel matrix

H =







y0 y1 . . . yn−1

y1 yn

... . .. ... yn−1 yn . . . y2n−2





∈Zⁿ^×ⁿ.

If gcd(det(H), m) = 1, then πm(x0) = ordm(A). Since x0 =en = [0, . . . ,0,1]^T yields a matrixH with det(H)∈ {±1}, we conclude that there always exists x0 with πm(x0) = ordm(A).

Proof: Suppose thatA^px0 =x0. Let H₋i denote the i-th column ofH. Noting thatx0=H₋₁andAHi =H₋(i+1), we may conclude that A^pH = H. By assumption, H is invertible as a matrix overZ^m, and henceA^p = I holds

overZ^m.

The prime divisors of det(H) are called bad primes of the initial value problem (1),(2).

Example: For y(t+ 3) = y(t+ 2) +y(t+ 1) +y(t) and x0 = [1,1,1]^T, the only bad prime is 2. For odd m, we haveπm(x0) = ordm(A).

3. COMPUTING THE ORDER

We will now discuss how to compute ordm(A) for a matrix A∈Z^n×nthat is invertible when considered as an element of Zⁿm^×ⁿ. Note that we do not restrict to companion matrices. If the prime factorization of m is known, the problem can be reduced to the case where m is prime power, by the Chinese remainder theorem. For primesp, an efficient algorithm by Celler and Leedham-Green (1997) can be used to compute ordp(A).

Theorem 3. LetA∈Zⁿ^×ⁿ and a primepwithpdet(A) be given. ThenAis invertible over Zp^k for allk≥1. We have

ordp^k(A)|ordp(A)p^k⁻¹.

Proof: SinceA^r=I implies that the order ofAdivides r, it suffices to show that

A^ord^p^(A)p^k⁻¹≡I (modp^k). SetB :=A^ord^p^(A). Then we have

B ≡I (mod p). We need to prove that

B^p^k−1≡I (modp^k)

for all k ≥ 1. This will be done by induction on k. The statement is clearly true for k = 1. Let’s assume that it holds fork. This means that

B^p^k⁻¹ =I+p^kC

for some matrixC∈Zⁿ^×ⁿ. By the binomial theorem,

(2)

Theorem 1. For each m > 0, there exist integers k ≥ 0 andp >0 such that

y(t+p) =y(t) for allt≥k.

If gcd(a0, m) = 1, thenkcan be taken to be zero.

Proof: Define

x(t) := [y(t), . . . , y(t+n−1)]^T. Thenx(t+ 1) =Ax(t), where

A=







0 1

... . ..

0 1

−a0 −a1 . . . −an−1







is the companion matrix of the polynomial χA = sⁿ + n−1

i=0 aisⁱ. Since the ringZ^mis finite, the sequencey, and hence alsox, can only take finitely many values. Thus there exists k ≥ 0 andp > 0 such that x(k+p) = x(k). This implies thatx(t+p) =x(t) for all t≥k.

If gcd(a0, m) = 1, then a0 is a unit in Z^m. This implies that the polynomial indeterminatesis a unit in the finite ring R := Zm[s]/χA. Hence it has finite order, that is, there existsp >0 such thats^p= 1 holds inR. This means that the polynomial s^p −1 is a multiple of χA in Z^m[s]

and thusy(t+p) =y(t) holds for allt∈N. Givenmwith gcd(a0, m) = 1, we will address the compu- tation of the smallest numberp >0 such that

(1) x(p) =x(0). Since this implies thatx(t+p) =x(t) for allt ∈N, this number will be called the period of x with respect to the initial valuex0:= [y0, . . . , yn−1]^T. It is the smallest p >0 with A^px0 =x0 and will be denoted byπm(x0).

(2) A^p = I. This p is the order of the invertible matrix A∈Zⁿm^×ⁿ, denoted by ordm(A).

We observe that πm(x0) must be a divisor of ordm(A).

Moreover,

ordm(A) = lcm{πm(x0)|x0∈Zⁿm}.

By the Chinese remainder theorem, it suffices to consider the case wheremis a prime power.

Example: Consider

y(t+3) =y(t+2)+y(t+1)+y(t), y(0) =y(1) =y(2) = 1.

• Modulo 2, we obtain the constant sequence 1, and henceπ2([1,1,1]^T) = 1, but ord2(A) = 4. Indeed, we have

x^T₀ π2(x0) [0,0,0] 1 [0,0,1] 4 [0,1,0] 2 [0,1,1] 4 [1,0,0] 4 [1,0,1] 2 [1,1,0] 4 [1,1,1] 1.

Adopting the notation used by Deng (2015), the cycle structure is X = Z³2 = 2C1∪C2∪C4. This means that there are two fixed points (i.e., cycles of length one), one cycle of length two, and one cycle of length four, corresponding to a partition of the state set into cycles according to 8=1+1+2+4.

• Modulo 3, we obtain a sequence of period 13, that is, π3([1,1,1]^T) = 13, and we also have ord3(A) = 13.

Indeed, all x0= 0 yieldπ3(x0) = 13.

• Modulo 4, we obtain a sequence of period 4, that is, π4([1,1,1]^T) = 4, but ord4(A) = 8.

Theorem 2. Given x0 = [y0, . . . , yn−1]^T ∈ Zⁿ, compute yn, . . . , y2n−2 ∈Zaccording to (1) and define the Hankel matrix

H =







y0 y1 . . . yn−1

y1 yn

... . .. ... yn−1 yn . . . y2n−2





∈Zⁿ^×ⁿ.

If gcd(det(H), m) = 1, then πm(x0) = ordm(A). Since x0 =en = [0, . . . ,0,1]^T yields a matrixH with det(H)∈ {±1}, we conclude that there always exists x0 with πm(x0) = ordm(A).

Proof: Suppose thatA^px0 =x0. Let H₋i denote the i-th column ofH. Noting thatx0=H₋₁andAHi =H₋(i+1), we may conclude that A^pH = H. By assumption, H is invertible as a matrix overZ^m, and hence A^p =I holds

overZ^m.

The prime divisors of det(H) are called bad primes of the initial value problem (1),(2).

Example: For y(t+ 3) = y(t+ 2) +y(t+ 1) +y(t) and x0 = [1,1,1]^T, the only bad prime is 2. For odd m, we haveπm(x0) = ordm(A).

3. COMPUTING THE ORDER

We will now discuss how to compute ordm(A) for a matrix A∈Z^n×nthat is invertible when considered as an element of Zⁿm^×ⁿ. Note that we do not restrict to companion matrices. If the prime factorization of m is known, the problem can be reduced to the case where m is prime power, by the Chinese remainder theorem. For primesp, an efficient algorithm by Celler and Leedham-Green (1997) can be used to compute ordp(A).

Theorem 3. Let A∈Zⁿ^×ⁿ and a primepwithpdet(A) be given. ThenAis invertible over Zp^k for allk ≥1. We have

ordp^k(A)|ordp(A)p^k⁻¹.

Proof: SinceA^r=I implies that the order ofAdivides r, it suffices to show that

A^ord^p^(A)p^k⁻¹≡I (modp^k).

SetB :=A^ord^p^(A). Then we have B≡I (mod p).

We need to prove that

B^p^k−1≡I (modp^k)

for all k ≥ 1. This will be done by induction on k. The statement is clearly true for k = 1. Let’s assume that it holds fork. This means that

B^p^k⁻¹ =I+p^kC

for some matrixC∈Zⁿ^×ⁿ. By the binomial theorem,