Search of linear recurrence correlation with constant integer coefficients in pre-set sequence by means of Eucludean algorithm

Recurrences pointed in the title are usual for the enumerative combinatoric problems that can be solved by using the transfer matrix method [2], [3], [4], [5]. There are some software tools for finding such recurrences – the so called Computer Algebra Systems (CAS). However, the testing has shown that the ability of Mathematica 7.0/8.0 is not enough to get repeated well-known results [8]. CAS Maple 14 is more effective but the increase of the recurrence order dramatically raises the amount of RAM to be used. For example, the attempt of finding the recurrence with the order of 2086 for the number of Hamiltonian circuits in P₁₂ × P_n [4] has failed the system.

The paper [1] was dedicated to the Dixon’s algorithm application [7] to the problem of finding linear recurrence. It was extraordinary faster than the known software tools because of the p-adic expansions’ usage. This paper shows that the Euclidean algorithm [6] modified in such a way that it uses modular arithmetic allows us to reduce the running time comparing to Dixon’s algorithm in some cases.

ALGORITHM

According to the Cayley – Hamilton theorem any sequence a₁, a₂, a₃… which has been obtained by use of transfer matrix method satisfies a linear homogeneous recurrence with the constant coefficients

« к = Е у=1 C • a_k _j , k > N, (1)

where N is known order. The coefficients of the recurrence c_j can be found out of the system of linear equations A _N ∙ c _N = b _N where

	( a .	a . -1    '	•    ai ^		( c .'		(a aN +1
A N =	a N +1	aN	a 2	, c N	c 1	, b N	a N +2
	V a 2 N -1	a 2 N -2   '	' aN J		v CN J		v a 2 N j

Note that the uniqueness of the solution is ensured by the inequality in (1). The system requires first 2N terms of the sequence to be known. So our case is significantly different from the other known algorithms for recovering the recurrence for terms in the given sequences with unknown properties.

If N is unknown, then we have to consider its maximum value, assuming N to be equal to the order of the transfer matrix which can significantly exceed the real order of our relation. In the Dixon’s algorithm we have to find the rank of A _N which shows the real order of the recurrence. The Euclidian algorithm works rather differently – it doesn’t need to determine the order before calculations. Let’s begin with its classical description [6].

Iterative scheme of the algorithm can be represented as follows. Let us assume s^x) = X2N, t^x) = Ef=o-1ai • %l, and

-*">«= Й °], then we repeat the following procedure until get degt(r)(x) < N – 1

9 % = R?J

Д(г)(%) = [1 -qw^J ^(r 1)(%)

[ s⁽’(%)      0        1         s ⁽ ’ ^-1) (%)l.

[t (-) cx)J   [° -e ^(r)to l к^-1) (%)]

After that, coefficients of the polynomial f (%) = —Д^(-1) • Л^ (%) where A= ^ 2’2 (0) correspond to the coefficients c_j of our recurrence (1).

However, the usage of rational numbers for coefficients of all involved polynomials brings extreme growth of numerator and denominator at the intermediate steps of the algorithm: the size of fractions have been grown up to one million of bits on hard tests. At the last step when calculating f (%) = —Д^(-1) • Д ’) (%) , the coefficients should become integer and rather small. For example, the number of Hamiltonian cycles in P₁₀×P_n are subjected to the recurrence of order 346 [8] (sequence

A180504). The required initial numbers in the sequence of this numbers can exceed 1000 decimal digits, while the biggest coefficient c_j for the recurrence has only 70 decimal digits. Therefore, it was decided to use modular arithmetic to avoid hard fraction calculations.

Let us consider that all calculations are performed by modulo, where p is a prime integer. First, we consider some small prime number, for example p = 2³¹ – 1 . Consider the modified Euclidean algorithm, where all polynomial operations are performed by modulo p . So, all coefficients of all polynomials at the intermediate steps of the algorithm are less than p . The modified algorithm will give us some solution c j^p) by modulo p. Then we have to restore the right coefficients (if possible).

Consider that the numbers c_j lie in the range (- [p/2] , [p/2] ). So if a number c® is greater or equal to [p/2] , then we have to decide C j = с ( ^р) - p , otherwise C j = cj^p\ On having done this, we should check the solution by substitution c_j into (1). If all the numbers a_k for k = N, N+1, …, 2N are satisfying the recurrence, we have found the right solution. Otherwise we have to greaten p . Good strategy is to find the smallest prime that is bigger than p² and to repeat all calculations again. Another strategy is to input the maximum number of bits b to the program. It finds the nearest to 2 ^b prime p and after having done all calculation shous if the answer is found or it isn’t. On the one hand, using the second strategy the user have to select the upper bound on b by hands, while the first strategy will make all selections automatically. On the other hand, the second strategy allows to select more accurate bound on b , then the first one.

EXAMPLE

Let us consider an example: we have a sequence of integers, 1, 14, 154, 1696, 18684, 205832, 2267544, 24980352, 275195536, 3031685984, … This sequence corresponds to the number of Hamiltonian cycles in P₅×P_2n for n = 1, 2, …, 10 [9] (sequence A006865). Let us apply the algorithm performing all operation by modulo p=13 . Let us assume

N = 3, s ⁽⁰⁾( x ) = x ⁶, t ⁽⁰⁾( x ) = 1 + x + 11 x ² + 6 x ³ + 3 x ⁴ + 5 x ⁵, and №[0 01

Begin calculating

Q(1)(x) = 9x + 4

^H? 4 + J

s ⁽¹⁾( x ) = 1 + x + 11 x ² + 6 x ³ + 3 x ⁴ + 3 x ⁵ t ⁽¹⁾( x ) = 9 + 12 x ² + 7 x ³ + 12 x ⁴

We will continue until the degree of the polynomial t ( x ) is less than N – 1 = 3 – 1 = 2.

Q ⁽²⁾( x ) = 10 x +2

A (2)           1         4 + 9% J

⁽ %)   111 + 3% 9 + 6% + 12%²J

s ⁽²⁾( x ) = 9 + 12 x + 7 x ³ +12 x ⁴ t ⁽²⁾( x ) = 9 + 2 x + 2 x ³

Q ⁽³⁾( x ) = 6 x + 10

Л^{(3) (} %) = [,

11 + 3%    9 + 6% + 12%2

. 8 + 8% + 8%2 10 + 7% + 6%3

s ⁽³⁾( x ) = 9 + 2 x + 2 x ³ t ⁽³⁾( x ) = 10 + 4 x .

We see that the degree of polynomial t ⁽³⁾ ( x ) is equal to 1 and it is less than 2. So we have / ⁽ %) = _л^(_о) m ^odp = 11% ⁺ 2%³ and respectively we have got the recurrence a_n = 11 a_{n –1} + 2 a_{n –3}. Having done this, we should check the solution by substitution c_j into (1). In our case all the numbers a_n for n = deg Q +1,…,2N are satisfying the recurrence, we can argue that the solution is right. If we chose p=7, for example, we would get a_n = 4a_n-1+2a_n-3. Obviously, this is the wrong answer. That is why we have chosen another p .

TEST RESULTS

The testing was conducted on four tasks, the first, the second and the third problems are enumeration the number of Hamiltonian circles in rectangular lattices Pm × Pn [4], the number Hamiltonian circles in thin cylinder Cm × Pn [4] and the number Hamiltonian circles in torus Cm × Cn [4] and the fourth one is the placement of kings in a rectangular chessboard of size 2m × 2n [2].

Table 1

Lattice size, the corresponding order of the relation and the spent time (s) for the problem of finding the number of Hamiltonian circles in rectangular lattices

m	N	Mathematica 8	Maple 14	Dixon	Euclid
7	18	0,3	0,5	0,0	0,0
8	66	25,9	0,5	0,0	0,1
9	104	–	2,0	0,2	0,2
10	346	–	14,2	1,7	1,6
11	671	–	555,9	32,0	32,3
12	2086	–	–	821,8	650,0

Table 2

Cylinder size, the corresponding order of the relation and the spent time (s) for the problem of finding the number of Hamiltonian circles in thin cylinder

m	N	Mathematica 8	Maple 14	Dixon	Euclid
8	20	0,4	0,0	0,0	0,0
9	51	8,6	0,0	0,0	0,1
10	74	92,5	1,0	0,0	0,1
11	246	–	6,8	0,7	0,6
12	303	–	26,4	1,3	1,2
13	1320	–	–	197,1	160,5
14	1514	–	–	346,8	267,5

Table 3

Table 4

Torus size, the corresponding order of the relation and the spent time(s) for the problem of finding the number of Hamiltonian circles in torus

m	N	Mathematica 8	Maple 14	Dixon	Euclid
4	28	0,4	0,0	0,0	0,088
5	84	74,4	0,6	0,0	0,113
6	257	–	26,3	0,7	0,492
7	856	–	430,8	29,5	19,612
8	2785	–	–	1849,4	>2 GiB RAM

Board size, the corresponding order of the relation and the spent time (s) for the problem of placement of kings in a rectangular chessboard

m	N	Mathematica 8	Maple 14	Dixon	Euclid
4	17	0,4	0,0	0,0	0,049
5	31	0,5	0,0	0,0	0,045
6	75	38,9	0,0	0,2	0,194
7	124	812,6	1,5	0,1	0,255
8	307	–	15,4	1,5	2,699
9	548	–	774,7	10,5	17,413
10	1318	–	–	455,1	>2 GiB RAM

We have used the computer with Intel Intel Core 2 Duo E-8400 @ 3.00 GHz, 4 GiB RAM, 32-bit operating system. Tables 1, 2, 3 and 4 show the runtime of Mathematica 8, Maple 14, Dixon’s algorithm [1] and the described Euclidean algorithm. The runtime is given in seconds.

As we can see that a described Euclidean algorithm is faster than other known tools in our test cases. Of course, one could test this algorithm in other test cases and get the other results. We had to show that in case of the leak of RAM Euclidean algorithm works better with exponentially growing integer sequences.

* Работа выполнена при поддержке Программы стратегического развития (ПСР) ПетрГУ в рамках реализации комплекса мероприятий по развитию научно-исследовательской деятельности на 2012–2016 гг.