Detailed mathematical approach

Problem Summary

Let \(T(n)\) denote the number of tilings of a strip of length \(n\) using \(1\times 2\) dominoes and \(1\times 1\) tiles colored with one of ten digits. The empty strip contributes one tiling, so

$T(0)=1,\qquad T(1)=10,\qquad T(n)=10T(n-1)+T(n-2)\quad(n\ge 2).$

The task is to evaluate

$S(L)=\sum_{1\le a,b,c\le L}\gcd\bigl(T(c^a),T(c^b)\bigr)\pmod{987898789}$

for \(L=2000\). A direct computation is impossible because there are \(L^3\) triples and the indices \(c^a\) become enormous almost immediately.

Mathematical Approach

Step 1: Rewrite the tiling sequence as a Lucas sequence

Define a Lucas sequence of the first kind by

$U_0=0,\qquad U_1=1,\qquad U_{n+2}=10U_{n+1}+U_n.$

Then

$U_1=1,\quad U_2=10,\quad U_3=101,\quad U_4=1020,\dots$

so the tiling numbers are simply shifted by one position:

$T(n)=U_{n+1}.$

This is the key structural observation, because Lucas sequences with parameter \(Q=-1\) satisfy the strong divisibility law

$\gcd(U_m,U_n)=U_{\gcd(m,n)}.$

Therefore, for arbitrary nonnegative integers \(x\) and \(y\),

$\gcd(T(x),T(y))=\gcd(U_{x+1},U_{y+1})=U_{\gcd(x+1,y+1)}=T\bigl(\gcd(x+1,y+1)-1\bigr).$

So the original problem is reduced to the arithmetic of \(\gcd(c^a+1,c^b+1)\).

Step 2: Evaluate \(\gcd(c^a+1,c^b+1)\) with 2-adic valuations

Write

$a=2^r u,\qquad b=2^s v,$

where \(u\) and \(v\) are odd. The answer depends only on whether \(r\) and \(s\) are equal, i.e. on whether \(\nu_2(a)=\nu_2(b)\).

If \(r=s\), let \(d=\gcd(a,b)\). Then we may write

$a=d\alpha,\qquad b=d\beta,\qquad \gcd(\alpha,\beta)=1,$

and both \(\alpha\) and \(\beta\) are odd. Setting \(y=c^d\), we get

$\gcd(c^a+1,c^b+1)=\gcd(y^\alpha+1,y^\beta+1)=y+1=c^d+1.$

If \(r\ne s\), factor out \(d=\gcd(a,b)\) again, but now one of the quotients \(a/d\), \(b/d\) is odd and the other is even. Put \(y=c^d\), and let \(g\) be a common divisor of \(y^{a/d}+1\) and \(y^{b/d}+1\). Then

$y^{a/d}\equiv -1 \pmod g,\qquad y^{b/d}\equiv -1 \pmod g.$

Raise the first congruence to the even exponent \(b/d\), and the second to the odd exponent \(a/d\). This gives

$y^{ab/d^2}\equiv 1 \pmod g,\qquad y^{ab/d^2}\equiv -1 \pmod g,$

hence \(2\equiv 0\pmod g\). So every common divisor must divide \(2\). Therefore

$\gcd(c^a+1,c^b+1)=\begin{cases} c^{\gcd(a,b)}+1,& \nu_2(a)=\nu_2(b),\\ 2,& \nu_2(a)\ne \nu_2(b)\text{ and }c\text{ is odd},\\ 1,& \nu_2(a)\ne \nu_2(b)\text{ and }c\text{ is even}. \end{cases}$

Step 3: Translate that gcd back to tiling numbers

Substitute the previous identity into the Lucas-sequence formula. If \(d=\gcd(a,b)\), then

$\gcd\bigl(T(c^a),T(c^b)\bigr)=\begin{cases} T(c^d),& \nu_2(a)=\nu_2(b),\\ T(1)=10,& \nu_2(a)\ne \nu_2(b)\text{ and }c\text{ is odd},\\ T(0)=1,& \nu_2(a)\ne \nu_2(b)\text{ and }c\text{ is even}. \end{cases}$

This is exactly the dichotomy exploited by the implementations. The huge indices disappear from the gcd itself; only the equal-valuation pairs still need values of the form \(T(c^d)\).

Step 4: Rearrange the triple sum by counting exponent pairs once

For each \(d\in\{1,\dots,L\}\), define the number of ordered pairs

$C_d=\#\{(a,b):1\le a,b\le L,\ \nu_2(a)=\nu_2(b),\ \gcd(a,b)=d\}.$

Also define the complementary count

$C_{\mathrm{diff}}=L^2-\sum_{d=1}^{L} C_d,$

which counts the ordered pairs with different 2-adic valuations. For a fixed base \(c\), let

$B(c)=\begin{cases} 10,& c\text{ odd},\\ 1,& c\text{ even}. \end{cases}$

Then the entire contribution of this single \(c\) is

$\sum_{1\le a,b\le L}\gcd\bigl(T(c^a),T(c^b)\bigr)=C_{\mathrm{diff}}\,B(c)+\sum_{d=1}^{L} C_d\,T(c^d).$

Summing over \(c\) yields the final rearrangement

$\boxed{S(L)=\sum_{c=1}^{L}\left(C_{\mathrm{diff}}\,B(c)+\sum_{d=1}^{L} C_d\,T(c^d)\right).}$

Now the pair statistics \(\{C_d\}\) are independent of \(c\), so they can be precomputed once.

Step 5: Compute \(T(c^d)\bmod M\) by period reduction

Let \(M=987898789\). The recurrence can be encoded with the companion matrix

$A=\begin{bmatrix}10&1\\ 1&0\end{bmatrix},\qquad \begin{bmatrix}T(n+1)\\ T(n)\end{bmatrix}=A^n\begin{bmatrix}10\\ 1\end{bmatrix}.$

Modulo \(M\), the state space is finite, so the sequence is periodic. The implementations compute the exact period from the matrix order. Since \(M\) is prime and the discriminant

$\Delta=10^2+4=104$

is a quadratic residue modulo \(M\), namely

$\left(\frac{104}{M}\right)=1,$

the eigenvalues of \(A\) lie in \(\mathbb{F}_M\), so the order of \(A\) divides \(M-1\). The algorithm starts from the candidate \(M-1=987898788\), factors it, and removes prime factors whenever the smaller exponent already gives \(A^k=I\). For this modulus no reduction occurs, so the period is

$\pi=987898788.$

Hence every huge index is reduced as

$T(n)\equiv T(n\bmod \pi)\pmod M.$

For each fixed \(c\), the reduced powers are generated iteratively:

$p_1\equiv c,\qquad p_{d+1}\equiv p_d\,c \pmod \pi,$

so the program never forms \(c^d\) as a gigantic integer.

Worked example: \(L=2\)

When \(L=2\), the ordered pairs \((a,b)\) are \((1,1)\), \((1,2)\), \((2,1)\), \((2,2)\). Equal 2-adic valuation occurs only for \((1,1)\) and \((2,2)\), so

$C_1=1,\qquad C_2=1,\qquad C_{\mathrm{diff}}=2.$

For \(c=1\), every term \(T(1^d)\) equals \(T(1)=10\), and \(B(1)=10\). Thus

$2\cdot 10 + 1\cdot 10 + 1\cdot 10 = 40.$

For \(c=2\), we have \(B(2)=1\), \(T(2)=101\), and \(T(4)=10301\), so

$2\cdot 1 + 1\cdot 101 + 1\cdot 10301 = 10404.$

Therefore

$S(2)=40+10404=10444,$

which matches the checkpoint used by the implementations.

How the Code Works

The C++, Python, and Java implementations all follow the same mathematical plan. They first precompute the 2-adic valuation of every integer from \(1\) to \(L\), then scan all ordered pairs \((a,b)\). If the valuations are equal, the pair is added to the bucket indexed by \(\gcd(a,b)\); otherwise it contributes only to the residual count \(C_{\mathrm{diff}}\).

Next they determine the recurrence period modulo \(987898789\), build a table of binary powers of the companion matrix, and loop over \(c=1,\dots,L\). For each fixed \(c\), the reduced indices \(c^1,c^2,\dots,c^L\pmod \pi\) are generated by repeated multiplication. Each reduced index is converted to \(T(\cdot)\bmod M\) by matrix exponentiation, multiplied by its precomputed pair count, and added to the running sum. The C++ implementation additionally parallelizes the outer loop over \(c\), but the arithmetic is identical in all three languages.

Complexity Analysis

The pair-count precomputation examines all ordered pairs \((a,b)\), so it costs \(O(L^2)\) time. The main accumulation also has \(L^2\) terms. Evaluating one reduced index uses a fixed table of matrix squares, so formally the cost is \(O(\log \pi)\) per term; here \(\pi=987898788\) is fixed by the modulus, so in practice this is a small constant. Thus the total running time is \(O(L^2\log \pi)\), which is effectively \(O(L^2)\) for this problem, and the memory usage is \(O(L)\).

References

1. Problem page: https://projecteuler.net/problem=440
2. Lucas sequences: Wikipedia — Lucas sequence
3. Linear recurrences and companion matrices: Wikipedia — Linear recurrence with constant coefficients
4. Exponentiation by squaring: Wikipedia — Exponentiation by squaring
5. Legendre symbol: Wikipedia — Legendre symbol