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.
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)\).
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}$$
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)\).
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.
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.
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.
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.
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)\).
Sei \(T(n)\) die Anzahl der Pflasterungen eines Streifens der Lange \(n\) mit \(1\times 2\)-Dominos und \(1\times 1\)-Feldern, die eine von zehn Ziffern tragen konnen. Fur den leeren Streifen gilt \(T(0)=1\), und insgesamt erhalt man
$$T(0)=1,\qquad T(1)=10,\qquad T(n)=10T(n-1)+T(n-2)\quad(n\ge 2).$$
Gesucht ist
$$S(L)=\sum_{1\le a,b,c\le L}\gcd\bigl(T(c^a),T(c^b)\bigr)\pmod{987898789}$$
fur \(L=2000\). Eine direkte Auswertung ist ausgeschlossen, weil sowohl die Anzahl der Tripel als auch die Indizes \(c^a\) sehr schnell explodieren.
Betrachte die Lucas-Folge erster Art
$$U_0=0,\qquad U_1=1,\qquad U_{n+2}=10U_{n+1}+U_n.$$
Dann gilt
$$U_1=1,\quad U_2=10,\quad U_3=101,\quad U_4=1020,\dots$$
also genau
$$T(n)=U_{n+1}.$$
Fur diese Lucas-Folge mit \(Q=-1\) gilt die starke Teilbarkeitseigenschaft
$$\gcd(U_m,U_n)=U_{\gcd(m,n)}.$$
Damit folgt fur beliebige \(x,y\ge 0\)
$$\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).$$
Das ganze Problem reduziert sich also auf \(\gcd(c^a+1,c^b+1)\).
Schreibe
$$a=2^r u,\qquad b=2^s v,$$
wobei \(u\) und \(v\) ungerade sind. Entscheidend ist nur, ob \(r=s\) gilt, also ob \(\nu_2(a)=\nu_2(b)\).
Im Fall \(r=s\) sei \(d=\gcd(a,b)\). Dann gibt es Darstellungen
$$a=d\alpha,\qquad b=d\beta,\qquad \gcd(\alpha,\beta)=1,$$
und \(\alpha,\beta\) sind beide ungerade. Mit \(y=c^d\) ergibt sich
$$\gcd(c^a+1,c^b+1)=\gcd(y^\alpha+1,y^\beta+1)=y+1=c^d+1.$$
Im Fall \(r\ne s\) faktorisiert man wieder \(d=\gcd(a,b)\) heraus; dann ist einer der Quotienten \(a/d\), \(b/d\) ungerade und der andere gerade. Sei \(g\) ein gemeinsamer Teiler von \(y^{a/d}+1\) und \(y^{b/d}+1\). Dann gilt
$$y^{a/d}\equiv -1 \pmod g,\qquad y^{b/d}\equiv -1 \pmod g.$$
Erhebt man die erste Kongruenz auf die gerade Potenz \(b/d\) und die zweite auf die ungerade Potenz \(a/d\), so folgt
$$y^{ab/d^2}\equiv 1 \pmod g,\qquad y^{ab/d^2}\equiv -1 \pmod g,$$
also \(2\equiv 0\pmod g\). Jeder gemeinsame Teiler teilt also \(2\). Daher
$$\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{ und }c\text{ ungerade ist},\\ 1,& \nu_2(a)\ne \nu_2(b)\text{ und }c\text{ gerade ist}. \end{cases}$$
Setzt man dieses Ergebnis in die Lucas-Formel ein, erhalt man mit \(d=\gcd(a,b)\)
$$\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{ und }c\text{ ungerade ist},\\ T(0)=1,& \nu_2(a)\ne \nu_2(b)\text{ und }c\text{ gerade ist}. \end{cases}$$
Genau diese Fallunterscheidung benutzt das Programm.
Fur jedes \(d\in\{1,\dots,L\}\) definiere
$$C_d=\#\{(a,b):1\le a,b\le L,\ \nu_2(a)=\nu_2(b),\ \gcd(a,b)=d\}.$$
Erganzend sei
$$C_{\mathrm{diff}}=L^2-\sum_{d=1}^{L} C_d,$$
also die Anzahl geordneter Paare mit verschiedener 2-adischer Bewertung. Fur festes \(c\) setze
$$B(c)=\begin{cases} 10,& c\text{ ungerade},\\ 1,& c\text{ gerade}. \end{cases}$$
Dann hat dieses \(c\) den Gesamtbeitrag
$$\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).$$
Durch Summation uber alle \(c\) folgt
$$\boxed{S(L)=\sum_{c=1}^{L}\left(C_{\mathrm{diff}}\,B(c)+\sum_{d=1}^{L} C_d\,T(c^d)\right).}$$
Die Statistik der Paare hangt damit nicht mehr von \(c\) ab und wird nur einmal vorab berechnet.
Sei \(M=987898789\). Die Rekursion besitzt die Begleitmatrix
$$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\) ist der Zustandsraum endlich, also ist die Folge periodisch. Die Implementierungen bestimmen die exakte Periode uber die Ordnung dieser Matrix. Weil \(M\) prim ist und die Diskriminante
$$\Delta=10^2+4=104$$
quadratischer Rest modulo \(M\) ist, also
$$\left(\frac{104}{M}\right)=1,$$
liegen die Eigenwerte in \(\mathbb{F}_M\), und die Matrixordnung teilt \(M-1\). Deshalb startet der Algorithmus mit dem Kandidaten \(M-1=987898788\), faktorisiert ihn und versucht Primfaktoren zu entfernen, solange fur den kleineren Exponenten schon \(A^k=I\) gilt. Fur diesen Modul bleibt der Kandidat unverandert, also
$$\pi=987898788.$$
Daher darf jeder riesige Index reduziert werden:
$$T(n)\equiv T(n\bmod \pi)\pmod M.$$
Fur festes \(c\) werden die reduzierten Potenzen iterativ erzeugt durch
$$p_1\equiv c,\qquad p_{d+1}\equiv p_d\,c \pmod \pi,$$
sodass kein gigantisches \(c^d\) explizit gebildet werden muss.
Fur \(L=2\) lauten die geordneten Paare \((1,1)\), \((1,2)\), \((2,1)\), \((2,2)\). Gleiche 2-adische Bewertung tritt nur bei \((1,1)\) und \((2,2)\) auf. Also
$$C_1=1,\qquad C_2=1,\qquad C_{\mathrm{diff}}=2.$$
Bei \(c=1\) ist \(B(1)=10\) und \(T(1^d)=T(1)=10\), also
$$2\cdot 10 + 1\cdot 10 + 1\cdot 10 = 40.$$
Bei \(c=2\) gilt \(B(2)=1\), \(T(2)=101\) und \(T(4)=10301\), also
$$2\cdot 1 + 1\cdot 101 + 1\cdot 10301 = 10404.$$
Damit ergibt sich
$$S(2)=40+10404=10444,$$
genau der Kontrollwert aus den Implementierungen.
Die C++-, Python- und Java-Implementierungen folgen demselben mathematischen Plan. Zuerst wird fur jede Zahl von \(1\) bis \(L\) die 2-adische Bewertung vorab berechnet. Danach werden alle geordneten Paare \((a,b)\) durchlaufen. Bei gleicher Bewertung wandert das Paar in den Eimer mit Index \(\gcd(a,b)\); bei verschiedener Bewertung erhoht es nur den Restzahler \(C_{\mathrm{diff}}\).
Anschliessend werden die Periode modulo \(987898789\) und die binaren Potenzen der Begleitmatrix vorbereitet. Fur jedes \(c=1,\dots,L\) erzeugt die Implementierung dann die reduzierten Werte \(c^1,c^2,\dots,c^L\pmod \pi\) durch wiederholte Multiplikation. Jeder reduzierte Index wird per Matrixpotenz in \(T(\cdot)\bmod M\) umgewandelt, mit seinem vorab bekannten Gewicht multipliziert und aufsummiert. Die C++-Implementierung parallelisiert zusatzlich die aussere Schleife uber \(c\), mathematisch bleibt das Verfahren jedoch identisch.
Die Vorberechnung der Paarzahlen betrachtet alle geordneten Paare \((a,b)\) und kostet daher \(O(L^2)\) Zeit. Die Hauptsumme besitzt ebenfalls \(L^2\) Terme. Die Auswertung eines reduzierten Indexes benotigt eine feste Tabelle von Matrixquadraten und damit formal \(O(\log \pi)\) Zeit pro Term; da \(\pi=987898788\) fur diese Aufgabe fest ist, ist das praktisch ein kleiner Konstantfaktor. Insgesamt ergibt sich also \(O(L^2\log \pi)\), effektiv \(O(L^2)\), bei \(O(L)\) Speicher.
\(T(n)\), uzunlugu \(n\) olan bir seridin \(1\times 2\) domino taslari ve on farkli rakamdan biriyle boyanabilen \(1\times 1\) taslarla doseme sayisi olsun. Bos serit icin bir doseme vardir; dolayisiyla
$$T(0)=1,\qquad T(1)=10,\qquad T(n)=10T(n-1)+T(n-2)\quad(n\ge 2).$$
Hesaplanmasi istenen toplam
$$S(L)=\sum_{1\le a,b,c\le L}\gcd\bigl(T(c^a),T(c^b)\bigr)\pmod{987898789}$$
olup burada \(L=2000\). Dogrudan hesap yapmak mumkun degildir; hem \(L^3\) tane ucgenli secim vardir hem de \(c^a\) indisleri cok hizli buyur.
Birinci tur Lucas dizisini
$$U_0=0,\qquad U_1=1,\qquad U_{n+2}=10U_{n+1}+U_n$$
seklinde tanimlayalim. Bu durumda
$$U_1=1,\quad U_2=10,\quad U_3=101,\quad U_4=1020,\dots$$
ve dolayisiyla
$$T(n)=U_{n+1}$$
olur. \(Q=-1\) parametreli bu Lucas dizisi icin standart kuvvetli bolunebilme ozelligi vardir:
$$\gcd(U_m,U_n)=U_{\gcd(m,n)}.$$
Buradan tum \(x,y\ge 0\) icin
$$\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)$$
elde edilir. Boylece ana problem \(\gcd(c^a+1,c^b+1)\) hesabina indirgenir.
Su yazilisi kullanalim:
$$a=2^r u,\qquad b=2^s v,$$
burada \(u\) ve \(v\) tek sayilardir. Sonucun belirleyicisi \(r=s\) olup olmamasidir; yani \(\nu_2(a)=\nu_2(b)\) esitligi.
Eger \(r=s\) ise \(d=\gcd(a,b)\) olsun. O zaman
$$a=d\alpha,\qquad b=d\beta,\qquad \gcd(\alpha,\beta)=1$$
ve \(\alpha,\beta\) yine tektir. \(y=c^d\) dersek
$$\gcd(c^a+1,c^b+1)=\gcd(y^\alpha+1,y^\beta+1)=y+1=c^d+1$$
olur.
Eger \(r\ne s\) ise yine \(d=\gcd(a,b)\) faktorunu ayiralim. Bu kez \(a/d\) ve \(b/d\) oranlarindan biri tek, digeri cifttir. \(g\), \(y^{a/d}+1\) ile \(y^{b/d}+1\) ifadelerinin ortak boleni olsun. O halde
$$y^{a/d}\equiv -1 \pmod g,\qquad y^{b/d}\equiv -1 \pmod g.$$
Birinci kongruensin ustunu cift olan \(b/d\) kuvvetine, ikinci kongruensin ustunu tek olan \(a/d\) kuvvetine alirsak
$$y^{ab/d^2}\equiv 1 \pmod g,\qquad y^{ab/d^2}\equiv -1 \pmod g$$
elde edilir; dolayisiyla \(2\equiv 0\pmod g\) olur. Yani ortak bolen ancak \(2\)'yi bolebilir. Sonuc olarak
$$\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{ ve }c\text{ tek},\\ 1,& \nu_2(a)\ne \nu_2(b)\text{ ve }c\text{ cift}. \end{cases}$$
Yukaridaki sonuca Lucas baglantisini uygularsak, \(d=\gcd(a,b)\) icin
$$\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{ ve }c\text{ tek},\\ T(0)=1,& \nu_2(a)\ne \nu_2(b)\text{ ve }c\text{ cift}. \end{cases}$$
Kodun dayandigi temel kimlik tam olarak budur.
Her \(d\in\{1,\dots,L\}\) icin su sayiyi tanimlayalim:
$$C_d=\#\{(a,b):1\le a,b\le L,\ \nu_2(a)=\nu_2(b),\ \gcd(a,b)=d\}.$$
Ayrica
$$C_{\mathrm{diff}}=L^2-\sum_{d=1}^{L} C_d$$
farkli 2-adik degerlemeye sahip sirali ciftlerin sayisi olsun. Sabit bir \(c\) icin
$$B(c)=\begin{cases} 10,& c\text{ tek},\\ 1,& c\text{ cift} \end{cases}$$
olarak yazarsak, bu tek \(c\)'nin tum katkisi
$$\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)$$
olur. Boylece
$$\boxed{S(L)=\sum_{c=1}^{L}\left(C_{\mathrm{diff}}\,B(c)+\sum_{d=1}^{L} C_d\,T(c^d)\right)}$$
sonucuna ulasiriz. Yani \((a,b)\) ciftlerine ait sayim bilgisi \(c\)'den bagimsizdir ve bir kez onceden hesaplanabilir.
\(M=987898789\) olsun. Rekurans su eslik matrisiyle yazilir:
$$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\) ortaminda durum uzayi sonludur; dolayisiyla dizi periyodiktir. Uygulamalar bu periyodu matrisin mertebesinden hesaplar. \(M\) asaldir ve diskriminant
$$\Delta=10^2+4=104$$
icin
$$\left(\frac{104}{M}\right)=1$$
oldugundan, \(A\)'nin ozdegerleri \(\mathbb{F}_M\) icindedir ve matris mertebesi \(M-1\)'i boler. Bu nedenle algoritma \(M-1=987898788\) adayindan baslar, bunu carpanlara ayirir ve daha kucuk bir kuvvette \(A^k=I\) elde edildigi surece asal carpanlari tek tek siler. Bu problemde hicbir indirgeme olmaz; dolayisiyla
$$\pi=987898788.$$
Boylece devasa indisler
$$T(n)\equiv T(n\bmod \pi)\pmod M$$
ile kucultulur. Sabit \(c\) icin indirgenmis kuvvetler de
$$p_1\equiv c,\qquad p_{d+1}\equiv p_d\,c \pmod \pi$$
seklinde uretilir; yani program hicbir zaman gercek \(c^d\) sayisini acik bicimde kurmaz.
\(L=2\) icin sirali ciftler \((1,1)\), \((1,2)\), \((2,1)\), \((2,2)\) olur. Esit 2-adik degerleme sadece \((1,1)\) ve \((2,2)\) ciftlerinde vardir. Bu nedenle
$$C_1=1,\qquad C_2=1,\qquad C_{\mathrm{diff}}=2.$$
\(c=1\) icin \(B(1)=10\) ve tum \(T(1^d)=T(1)=10\) oldugundan katki
$$2\cdot 10 + 1\cdot 10 + 1\cdot 10 = 40$$
olur. \(c=2\) icin \(B(2)=1\), \(T(2)=101\) ve \(T(4)=10301\), dolayisiyla
$$2\cdot 1 + 1\cdot 101 + 1\cdot 10301 = 10404.$$
Toplam
$$S(2)=40+10404=10444$$
cikar; bu deger implementasyonlardaki kontrol sonucu ile aynidir.
C++, Python ve Java implementasyonlari ayni matematiksel akisi izler. Once \(1\) ile \(L\) arasindaki tum sayilar icin 2-adik degerleme onceden bulunur. Sonra tum sirali \((a,b)\) ciftleri taranir. Eger degerlemeler ayniysa cift, \(\gcd(a,b)\) ile indekslenen kovaya eklenir; farkliysa yalnizca kalan sayacini arttirir.
Daha sonra \(987898789\) modundaki periyot ve eslik matrisinin ikili kuvvetleri hazirlanir. Her \(c=1,\dots,L\) icin indirgenmis indisler \(c,c^2,\dots,c^L\pmod \pi\) tekrar carpma ile uretilir. Her indirgenmis indis matrix us alma yardimiyla \(T(\cdot)\bmod M\) degerine cevrilir, onceden hesaplanan agirlik ile carpilarak toplama eklenir. C++ implementasyonu distaki \(c\) dongusunu paralellestirir; ancak aritmetik yontem uc dilde de aynidir.
Cift sayim on islemi tum sirali \((a,b)\) ciftlerini gezdigi icin \(O(L^2)\) zamandir. Ana toplamin kendisi de \(L^2\) terim icerir. Tek bir indirgenmis indisi hesaplamak sabit boyutlu bir matris-kareleri tablosu kullandigi icin resmen \(O(\log \pi)\) zaman alir; burada \(\pi=987898788\) sabit oldugundan bu pratikte kucuk bir katsayidir. Sonuc olarak toplam sure \(O(L^2\log \pi)\), bu problem ozelinde etkin olarak \(O(L^2)\), bellek kullanimi ise \(O(L)\) olur.
Sea \(T(n)\) el numero de recubrimientos de una tira de longitud \(n\) usando fichas \(1\times 2\) y fichas \(1\times 1\) coloreadas con uno de diez digitos. La tira vacia aporta una configuracion, de modo que
$$T(0)=1,\qquad T(1)=10,\qquad T(n)=10T(n-1)+T(n-2)\quad(n\ge 2).$$
Hay que calcular
$$S(L)=\sum_{1\le a,b,c\le L}\gcd\bigl(T(c^a),T(c^b)\bigr)\pmod{987898789}$$
para \(L=2000\). Un recorrido directo no sirve: hay \(L^3\) ternas y los indices \(c^a\) se hacen gigantescos enseguida.
Consideremos la sucesion de Lucas de primera especie
$$U_0=0,\qquad U_1=1,\qquad U_{n+2}=10U_{n+1}+U_n.$$
Entonces
$$U_1=1,\quad U_2=10,\quad U_3=101,\quad U_4=1020,\dots$$
y por tanto
$$T(n)=U_{n+1}.$$
Para esta sucesion con \(Q=-1\) vale la propiedad de divisibilidad fuerte
$$\gcd(U_m,U_n)=U_{\gcd(m,n)}.$$
Luego, para cualesquiera \(x,y\ge 0\),
$$\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).$$
El problema completo queda reducido a entender \(\gcd(c^a+1,c^b+1)\).
Escribimos
$$a=2^r u,\qquad b=2^s v,$$
con \(u\) y \(v\) impares. Todo depende de si \(r=s\), es decir, de si \(\nu_2(a)=\nu_2(b)\).
Si \(r=s\), sea \(d=\gcd(a,b)\). Entonces
$$a=d\alpha,\qquad b=d\beta,\qquad \gcd(\alpha,\beta)=1,$$
y tanto \(\alpha\) como \(\beta\) son impares. Si ponemos \(y=c^d\), obtenemos
$$\gcd(c^a+1,c^b+1)=\gcd(y^\alpha+1,y^\beta+1)=y+1=c^d+1.$$
Si \(r\ne s\), volvemos a factorizar \(d=\gcd(a,b)\), pero ahora uno de los cocientes \(a/d\), \(b/d\) es impar y el otro es par. Sea \(g\) un divisor comun de \(y^{a/d}+1\) y \(y^{b/d}+1\). Entonces
$$y^{a/d}\equiv -1 \pmod g,\qquad y^{b/d}\equiv -1 \pmod g.$$
Elevando la primera congruencia a la potencia par \(b/d\) y la segunda a la potencia impar \(a/d\), se llega a
$$y^{ab/d^2}\equiv 1 \pmod g,\qquad y^{ab/d^2}\equiv -1 \pmod g,$$
asi que \(2\equiv 0\pmod g\). Por tanto, todo divisor comun divide a \(2\). En consecuencia,
$$\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{ y }c\text{ es impar},\\ 1,& \nu_2(a)\ne \nu_2(b)\text{ y }c\text{ es par}. \end{cases}$$
Insertando la identidad anterior en la formula de Lucas, con \(d=\gcd(a,b)\), resulta
$$\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{ y }c\text{ es impar},\\ T(0)=1,& \nu_2(a)\ne \nu_2(b)\text{ y }c\text{ es par}. \end{cases}$$
Esta es exactamente la ley aritmetica que aprovechan las implementaciones.
Para cada \(d\in\{1,\dots,L\}\), definimos
$$C_d=\#\{(a,b):1\le a,b\le L,\ \nu_2(a)=\nu_2(b),\ \gcd(a,b)=d\}.$$
Tambien definimos
$$C_{\mathrm{diff}}=L^2-\sum_{d=1}^{L} C_d,$$
que cuenta los pares ordenados con distinta valuacion 2-adica. Para un \(c\) fijo, escribimos
$$B(c)=\begin{cases} 10,& c\text{ impar},\\ 1,& c\text{ par}. \end{cases}$$
Entonces la contribucion total de ese \(c\) es
$$\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).$$
Al sumar sobre \(c\), obtenemos
$$\boxed{S(L)=\sum_{c=1}^{L}\left(C_{\mathrm{diff}}\,B(c)+\sum_{d=1}^{L} C_d\,T(c^d)\right).}$$
Asi, la parte combinatoria dependiente de \((a,b)\) queda precomputada y ya no depende de \(c\).
Sea \(M=987898789\). La recurrencia se codifica con la matriz companera
$$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\), el espacio de estados es finito y la sucesion es periodica. Las implementaciones obtienen el periodo exacto a partir del orden de la matriz. Como \(M\) es primo y el discriminante
$$\Delta=10^2+4=104$$
satisface
$$\left(\frac{104}{M}\right)=1,$$
los autovalores de \(A\) viven en \(\mathbb{F}_M\), y por eso el orden de \(A\) divide a \(M-1\). El algoritmo toma como candidato inicial \(M-1=987898788\), lo factoriza y va quitando factores primos cuando una potencia menor ya produce \(A^k=I\). En este modulo no se elimina nada, de modo que
$$\pi=987898788.$$
Por tanto,
$$T(n)\equiv T(n\bmod \pi)\pmod M.$$
Para cada \(c\), las potencias reducidas se generan iterativamente:
$$p_1\equiv c,\qquad p_{d+1}\equiv p_d\,c \pmod \pi,$$
y asi nunca hace falta construir el entero gigante \(c^d\).
Si \(L=2\), los pares ordenados son \((1,1)\), \((1,2)\), \((2,1)\), \((2,2)\). Solo \((1,1)\) y \((2,2)\) tienen la misma valuacion 2-adica, luego
$$C_1=1,\qquad C_2=1,\qquad C_{\mathrm{diff}}=2.$$
Para \(c=1\), se tiene \(B(1)=10\) y \(T(1^d)=T(1)=10\), asi que la contribucion es
$$2\cdot 10 + 1\cdot 10 + 1\cdot 10 = 40.$$
Para \(c=2\), \(B(2)=1\), \(T(2)=101\) y \(T(4)=10301\), por lo que
$$2\cdot 1 + 1\cdot 101 + 1\cdot 10301 = 10404.$$
En total,
$$S(2)=40+10404=10444,$$
que coincide con el punto de control usado por las implementaciones.
Las implementaciones en C++, Python y Java siguen exactamente el mismo esquema matematico. Primero precalculan la valuacion 2-adica de todos los enteros entre \(1\) y \(L\), y luego recorren todos los pares ordenados \((a,b)\). Si las valuaciones coinciden, el par se acumula en el contenedor indexado por \(\gcd(a,b)\); si no coinciden, solo incrementa el contador residual.
Despues calculan el periodo de la recurrencia modulo \(987898789\), preparan las potencias binarias de la matriz companera y recorren \(c=1,\dots,L\). Para cada \(c\), los indices reducidos \(c,c^2,\dots,c^L\pmod \pi\) se generan por multiplicacion repetida. Cada indice reducido se evalua mediante exponenciacion matricial, se multiplica por su peso precomputado y se suma al acumulado final. La implementacion en C++ ademas reparte el bucle exterior entre varios hilos, pero la aritmetica es la misma en los tres lenguajes.
La precuenta de pares examina todos los pares ordenados \((a,b)\), luego cuesta \(O(L^2)\). La suma principal tambien contiene \(L^2\) terminos. Evaluar un indice reducido usa una tabla fija de cuadrados de matrices, asi que formalmente cada termino cuesta \(O(\log \pi)\); aqui \(\pi=987898788\) es fijo, por lo que en la practica es un factor constante pequeno. El tiempo total es \(O(L^2\log \pi)\), efectivamente \(O(L^2)\) para este problema, y la memoria utilizada es \(O(L)\).
设 \(T(n)\) 表示长度为 \(n\) 的带状区域的铺法数量，允许使用 \(1\times 2\) 骨牌以及带有 10 种数字颜色的 \(1\times 1\) 小砖。空带有 1 种铺法，因此
$$T(0)=1,\qquad T(1)=10,\qquad T(n)=10T(n-1)+T(n-2)\quad(n\ge 2).$$
题目要求计算
$$S(L)=\sum_{1\le a,b,c\le L}\gcd\bigl(T(c^a),T(c^b)\bigr)\pmod{987898789}$$
其中 \(L=2000\)。直接枚举完全不可行，因为一方面有 \(L^3\) 组 \((a,b,c)\)，另一方面指数 \(c^a\) 很快就会变得极大。
先定义第一类 Lucas 序列
$$U_0=0,\qquad U_1=1,\qquad U_{n+2}=10U_{n+1}+U_n.$$
于是
$$U_1=1,\quad U_2=10,\quad U_3=101,\quad U_4=1020,\dots$$
所以铺砖数正好只是平移一位：
$$T(n)=U_{n+1}.$$
对于这里的 Lucas 序列，标准强可整除性质成立：
$$\gcd(U_m,U_n)=U_{\gcd(m,n)}.$$
因此对任意 \(x,y\ge 0\)，都有
$$\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).$$
原问题于是转化为研究 \(\gcd(c^a+1,c^b+1)\)。
把指数写成
$$a=2^r u,\qquad b=2^s v,$$
其中 \(u,v\) 都是奇数。结论只取决于 \(r\) 和 \(s\) 是否相等，也就是 \(\nu_2(a)\) 与 \(\nu_2(b)\) 是否相同。
若 \(r=s\)，记 \(d=\gcd(a,b)\)。则可写成
$$a=d\alpha,\qquad b=d\beta,\qquad \gcd(\alpha,\beta)=1,$$
并且 \(\alpha,\beta\) 都是奇数。令 \(y=c^d\)，便有
$$\gcd(c^a+1,c^b+1)=\gcd(y^\alpha+1,y^\beta+1)=y+1=c^d+1.$$
若 \(r\ne s\)，同样提出 \(d=\gcd(a,b)\)，此时 \(a/d\) 与 \(b/d\) 中一个是奇数、另一个是偶数。设 \(g\) 是 \(y^{a/d}+1\) 与 \(y^{b/d}+1\) 的公因子，则
$$y^{a/d}\equiv -1 \pmod g,\qquad y^{b/d}\equiv -1 \pmod g.$$
把第一式提升到偶数次幂 \(b/d\)，把第二式提升到奇数次幂 \(a/d\)，得到
$$y^{ab/d^2}\equiv 1 \pmod g,\qquad y^{ab/d^2}\equiv -1 \pmod g.$$
于是 \(2\equiv 0\pmod g\)，所以任何公因子都只能整除 \(2\)。因此
$$\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{ 且 }c\text{ 为奇数},\\ 1,& \nu_2(a)\ne \nu_2(b)\text{ 且 }c\text{ 为偶数}. \end{cases}$$
把上面的结论代回 Lucas 公式。若 \(d=\gcd(a,b)\)，则
$$\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{ 且 }c\text{ 为奇数},\\ T(0)=1,& \nu_2(a)\ne \nu_2(b)\text{ 且 }c\text{ 为偶数}. \end{cases}$$
这就是实现所依赖的核心等式。真正困难的三重 gcd 问题，被压缩成了“同 2 进层级时取 \(T(c^d)\)，否则取常数 10 或 1”。
对每个 \(d\in\{1,\dots,L\}\)，定义
$$C_d=\#\{(a,b):1\le a,b\le L,\ \nu_2(a)=\nu_2(b),\ \gcd(a,b)=d\}.$$
再定义补集计数
$$C_{\mathrm{diff}}=L^2-\sum_{d=1}^{L} C_d,$$
它表示 2 进赋值不同的有序对数量。对固定的 \(c\)，记
$$B(c)=\begin{cases} 10,& c\text{ 为奇数},\\ 1,& c\text{ 为偶数}. \end{cases}$$
那么这个 \(c\) 的总贡献就是
$$\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).$$
再对 \(c\) 求和，可得
$$\boxed{S(L)=\sum_{c=1}^{L}\left(C_{\mathrm{diff}}\,B(c)+\sum_{d=1}^{L} C_d\,T(c^d)\right).}$$
这样一来，与 \((a,b)\) 有关的组合信息只需预处理一次，之后就能被所有 \(c\) 共享。
设 \(M=987898789\)。递推可以写成伴随矩阵形式：
$$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}.$$
在模 \(M\) 下，状态空间有限，所以序列必然周期性。实现通过矩阵阶来求精确周期。由于 \(M\) 是素数，且判别式
$$\Delta=10^2+4=104$$
满足
$$\left(\frac{104}{M}\right)=1,$$
所以 \(A\) 的特征值位于 \(\mathbb{F}_M\) 中，矩阵阶必整除 \(M-1\)。实现于是从候选值 \(M-1=987898788\) 出发，对其分解质因数，并不断尝试删去可删除的质因子，只要更小的指数已经满足 \(A^k=I\)。在这个模数下没有任何质因子可以去掉，因此
$$\pi=987898788.$$
于是所有巨大下标都可约简为
$$T(n)\equiv T(n\bmod \pi)\pmod M.$$
对固定 \(c\)，约简后的幂通过递推生成：
$$p_1\equiv c,\qquad p_{d+1}\equiv p_d\,c \pmod \pi,$$
这样程序从头到尾都不需要显式构造真正的 \(c^d\)。
当 \(L=2\) 时，有序对为 \((1,1)\)、\((1,2)\)、\((2,1)\)、\((2,2)\)。只有 \((1,1)\) 和 \((2,2)\) 的 2 进赋值相同，因此
$$C_1=1,\qquad C_2=1,\qquad C_{\mathrm{diff}}=2.$$
对 \(c=1\)，有 \(B(1)=10\)，且所有 \(T(1^d)=T(1)=10\)，于是贡献为
$$2\cdot 10 + 1\cdot 10 + 1\cdot 10 = 40.$$
对 \(c=2\)，有 \(B(2)=1\)，而 \(T(2)=101\)、\(T(4)=10301\)，所以贡献为
$$2\cdot 1 + 1\cdot 101 + 1\cdot 10301 = 10404.$$
因此
$$S(2)=40+10404=10444,$$
这与实现中的校验值完全一致。
C++、Python 和 Java 实现遵循完全相同的数学流程。首先预处理 \(1\) 到 \(L\) 每个整数的 2 进赋值，然后枚举所有有序对 \((a,b)\)。如果两者赋值相同，就把该对加入按 \(\gcd(a,b)\) 分类的桶；否则只增加“不同赋值”计数。
接着，程序求出模 \(987898789\) 下的周期，并预先构造伴随矩阵的二进制幂表。随后对每个 \(c=1,\dots,L\)，通过重复乘法生成 \(c,c^2,\dots,c^L\pmod \pi\)。每个约简后的指数都用矩阵快速幂转成 \(T(\cdot)\bmod M\)，再乘以预处理好的权重并累加到答案中。C++ 版本还把最外层的 \(c\) 循环分配到多个线程，但三种语言的数论与递推逻辑完全一致。
配对预处理要检查所有有序对 \((a,b)\)，因此时间是 \(O(L^2)\)。主求和本身也有 \(L^2\) 个项。对单个约简指数求 \(T(\cdot)\) 时，需要扫描固定长度的矩阵平方表，所以形式上每项代价为 \(O(\log \pi)\)；在这里 \(\pi=987898788\) 由模数固定，因此这是很小的常数因子。总时间复杂度为 \(O(L^2\log \pi)\)，在本题中可视为 \(O(L^2)\)，空间复杂度为 \(O(L)\)。
Пусть \(T(n)\) обозначает число разбиений полосы длины \(n\) с помощью домино \(1\times 2\) и плиток \(1\times 1\), окрашенных в один из десяти цифровых цветов. Для пустой полосы имеется одно разбиение, поэтому
$$T(0)=1,\qquad T(1)=10,\qquad T(n)=10T(n-1)+T(n-2)\quad(n\ge 2).$$
Требуется вычислить
$$S(L)=\sum_{1\le a,b,c\le L}\gcd\bigl(T(c^a),T(c^b)\bigr)\pmod{987898789}$$
при \(L=2000\). Прямой перебор не годится: троек слишком много, а показатели \(c^a\) почти сразу становятся астрономическими.
Рассмотрим последовательность Лукаса первого рода
$$U_0=0,\qquad U_1=1,\qquad U_{n+2}=10U_{n+1}+U_n.$$
Тогда
$$U_1=1,\quad U_2=10,\quad U_3=101,\quad U_4=1020,\dots$$
и, следовательно,
$$T(n)=U_{n+1}.$$
Для такой Lucas-последовательности с параметром \(Q=-1\) выполняется сильное свойство делимости:
$$\gcd(U_m,U_n)=U_{\gcd(m,n)}.$$
Значит, для любых \(x,y\ge 0\),
$$\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).$$
Тем самым задача сводится к вычислению \(\gcd(c^a+1,c^b+1)\).
Запишем
$$a=2^r u,\qquad b=2^s v,$$
где \(u\) и \(v\) нечетны. Ответ зависит только от того, совпадают ли \(r\) и \(s\), то есть равны ли \(\nu_2(a)\) и \(\nu_2(b)\).
Если \(r=s\), положим \(d=\gcd(a,b)\). Тогда
$$a=d\alpha,\qquad b=d\beta,\qquad \gcd(\alpha,\beta)=1,$$
причем \(\alpha\) и \(\beta\) нечетны. Обозначив \(y=c^d\), получаем
$$\gcd(c^a+1,c^b+1)=\gcd(y^\alpha+1,y^\beta+1)=y+1=c^d+1.$$
Если \(r\ne s\), снова выносим \(d=\gcd(a,b)\), но теперь одно из чисел \(a/d\), \(b/d\) нечетно, а другое четно. Пусть \(g\) является общим делителем чисел \(y^{a/d}+1\) и \(y^{b/d}+1\). Тогда
$$y^{a/d}\equiv -1 \pmod g,\qquad y^{b/d}\equiv -1 \pmod g.$$
Возведем первое сравнение в четную степень \(b/d\), а второе в нечетную степень \(a/d\). Получим
$$y^{ab/d^2}\equiv 1 \pmod g,\qquad y^{ab/d^2}\equiv -1 \pmod g,$$
откуда следует \(2\equiv 0\pmod g\). Значит, любой общий делитель обязан делить \(2\). Следовательно,
$$\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{ и }c\text{ нечетно},\\ 1,& \nu_2(a)\ne \nu_2(b)\text{ и }c\text{ четно}. \end{cases}$$
Подставляя найденную формулу в Lucas-идентичность, при \(d=\gcd(a,b)\) получаем
$$\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{ и }c\text{ нечетно},\\ T(0)=1,& \nu_2(a)\ne \nu_2(b)\text{ и }c\text{ четно}. \end{cases}$$
Именно эту развилку использует алгоритм: либо нужно значение \(T(c^d)\), либо появляется постоянный вклад \(10\) или \(1\).
Для каждого \(d\in\{1,\dots,L\}\) определим
$$C_d=\#\{(a,b):1\le a,b\le L,\ \nu_2(a)=\nu_2(b),\ \gcd(a,b)=d\}.$$
Также введем
$$C_{\mathrm{diff}}=L^2-\sum_{d=1}^{L} C_d,$$
то есть число упорядоченных пар с разными 2-адическими показателями. Для фиксированного \(c\) положим
$$B(c)=\begin{cases} 10,& c\text{ нечетно},\\ 1,& c\text{ четно}. \end{cases}$$
Тогда вклад данного \(c\) равен
$$\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).$$
Суммируя по всем \(c\), получаем
$$\boxed{S(L)=\sum_{c=1}^{L}\left(C_{\mathrm{diff}}\,B(c)+\sum_{d=1}^{L} C_d\,T(c^d)\right).}$$
Комбинаторная часть, связанная только с парами \((a,b)\), тем самым вычисляется один раз и затем используется для всех \(c\).
Пусть \(M=987898789\). Рекуррентное соотношение кодируется матрицей-компаньоном
$$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}.$$
По модулю \(M\) множество состояний конечно, следовательно, последовательность периодична. Реализации находят точный период через порядок матрицы. Так как \(M\) простое, а дискриминант
$$\Delta=10^2+4=104$$
удовлетворяет
$$\left(\frac{104}{M}\right)=1,$$
собственные значения матрицы лежат в \(\mathbb{F}_M\), и потому порядок \(A\) делит \(M-1\). Алгоритм начинает с кандидата \(M-1=987898788\), раскладывает его на простые множители и удаляет те множители, для которых меньшая степень уже дает \(A^k=I\). В данном модуле ничего сократить нельзя, поэтому
$$\pi=987898788.$$
Отсюда
$$T(n)\equiv T(n\bmod \pi)\pmod M.$$
Для фиксированного \(c\) редуцированные степени строятся итеративно:
$$p_1\equiv c,\qquad p_{d+1}\equiv p_d\,c \pmod \pi,$$
так что огромные числа \(c^d\) никогда не создаются явно.
При \(L=2\) имеем пары \((1,1)\), \((1,2)\), \((2,1)\), \((2,2)\). Одинаковая 2-адическая оценка есть только у \((1,1)\) и \((2,2)\), значит
$$C_1=1,\qquad C_2=1,\qquad C_{\mathrm{diff}}=2.$$
Для \(c=1\) имеем \(B(1)=10\) и \(T(1^d)=T(1)=10\), поэтому вклад равен
$$2\cdot 10 + 1\cdot 10 + 1\cdot 10 = 40.$$
Для \(c=2\) получаем \(B(2)=1\), \(T(2)=101\) и \(T(4)=10301\), значит вклад равен
$$2\cdot 1 + 1\cdot 101 + 1\cdot 10301 = 10404.$$
Итак,
$$S(2)=40+10404=10444,$$
что в точности совпадает с контрольным значением из реализаций.
Реализации на C++, Python и Java следуют одному и тому же математическому плану. Сначала для всех чисел от \(1\) до \(L\) вычисляется 2-адическая оценка. Затем перебираются все упорядоченные пары \((a,b)\). Если оценки равны, пара добавляется в корзину с номером \(\gcd(a,b)\); если различны, она увеличивает только остаточный счетчик.
После этого программа находит период рекурсии по модулю \(987898789\), строит таблицу двоичных степеней матрицы-компаньона и проходит по всем \(c=1,\dots,L\). Для каждого \(c\) значения \(c,c^2,\dots,c^L\pmod \pi\) генерируются повторным умножением. Каждый редуцированный индекс преобразуется в \(T(\cdot)\bmod M\) матричным возведением в степень, умножается на заранее найденный вес и добавляется к ответу. В версии на C++ внешний цикл по \(c\) дополнительно распараллелен, но математическая схема во всех трех языках одинакова.
Предварительный подсчет пар рассматривает все упорядоченные пары \((a,b)\), то есть занимает \(O(L^2)\) времени. Основная сумма тоже содержит \(L^2\) слагаемых. Вычисление одного редуцированного индекса использует фиксированную таблицу квадратов матрицы, поэтому формально стоимость одного слагаемого равна \(O(\log \pi)\); здесь \(\pi=987898788\) фиксировано модулем, так что это небольшой постоянный множитель. Итого получаем \(O(L^2\log \pi)\), что для данной задачи практически эквивалентно \(O(L^2)\), при памяти \(O(L)\).
لتكن \(T(n)\) عدد طرق تبليط شريط طوله \(n\) باستخدام قطع دومينو \(1\times 2\) وقطع \(1\times 1\) يمكن تلوينها بعشرة رموز رقمية مختلفة. ولان الشريط الفارغ يملك طريقة واحدة، فان المتتالية تحقق
$$T(0)=1,\qquad T(1)=10,\qquad T(n)=10T(n-1)+T(n-2)\quad(n\ge 2).$$
المطلوب هو حساب
$$S(L)=\sum_{1\le a,b,c\le L}\gcd\bigl(T(c^a),T(c^b)\bigr)\pmod{987898789}$$
عند \(L=2000\). الحساب المباشر غير عملي لان عدد الثلاثيات هو \(L^3\)، كما ان الفهارس \(c^a\) تصبح هائلة بسرعة كبيرة.
نعرف متتالية لوكاس من النوع الاول كما يلي:
$$U_0=0,\qquad U_1=1,\qquad U_{n+2}=10U_{n+1}+U_n.$$
وعندئذ نحصل على
$$U_1=1,\quad U_2=10,\quad U_3=101,\quad U_4=1020,\dots$$
ومن ثم
$$T(n)=U_{n+1}.$$
لهذه المتتالية مع \(Q=-1\) توجد خاصية القسمة القوية المعروفة:
$$\gcd(U_m,U_n)=U_{\gcd(m,n)}.$$
وبالتالي، لاي \(x,y\ge 0\)،
$$\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).$$
وهكذا تختزل المسألة كلها الى فهم \(\gcd(c^a+1,c^b+1)\).
نكتب
$$a=2^r u,\qquad b=2^s v,$$
حيث \(u\) و\(v\) فرديان. النتيجة تعتمد فقط على ما اذا كان \(r=s\)، اي على ما اذا كان \(\nu_2(a)=\nu_2(b)\).
اذا كان \(r=s\)، فلتكن \(d=\gcd(a,b)\). عندئذ يمكن كتابة
$$a=d\alpha,\qquad b=d\beta,\qquad \gcd(\alpha,\beta)=1,$$
وتكون \(\alpha,\beta\) فرديتين ايضا. بوضع \(y=c^d\) نحصل على
$$\gcd(c^a+1,c^b+1)=\gcd(y^\alpha+1,y^\beta+1)=y+1=c^d+1.$$
اما اذا كان \(r\ne s\)، فنستخرج \(d=\gcd(a,b)\) مرة اخرى، لكن يصبح احد العددين \(a/d\) و\(b/d\) فرديا والاخر زوجيا. ولنفترض ان \(g\) قاسم مشترك للعددين \(y^{a/d}+1\) و\(y^{b/d}+1\). عندئذ
$$y^{a/d}\equiv -1 \pmod g,\qquad y^{b/d}\equiv -1 \pmod g.$$
برفع العلاقة الاولى الى القوة الزوجية \(b/d\)، ورفع الثانية الى القوة الفردية \(a/d\)، نحصل على
$$y^{ab/d^2}\equiv 1 \pmod g,\qquad y^{ab/d^2}\equiv -1 \pmod g,$$
ومن ثم \(2\equiv 0\pmod g\). وهذا يعني ان كل قاسم مشترك يجب ان يقسم \(2\). لذلك
$$\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),\ c\text{ odd},\\ 1,& \nu_2(a)\ne \nu_2(b),\ c\text{ even}. \end{cases}$$
بالتعويض في صيغة لوكاس السابقة، ومع \(d=\gcd(a,b)\)، نحصل على
$$\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),\ c\text{ odd},\\ T(0)=1,& \nu_2(a)\ne \nu_2(b),\ c\text{ even}. \end{cases}$$
وهذه هي القاعدة الحسابية المركزية التي تستعملها التنفيذات الثلاثة.
لكل \(d\in\{1,\dots,L\}\) نعرف
$$C_d=\#\{(a,b):1\le a,b\le L,\ \nu_2(a)=\nu_2(b),\ \gcd(a,b)=d\}.$$
ونعرف ايضا
$$C_{\mathrm{diff}}=L^2-\sum_{d=1}^{L} C_d,$$
وهو عدد الازواج المرتبة التي تختلف فيها التقييمات 2-الادية. ولكل \(c\) ثابت نضع
$$B(c)=\begin{cases} 10,& c\text{ odd},\\ 1,& c\text{ even}. \end{cases}$$
عندئذ يكون اسهام هذا \(c\) هو
$$\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).$$
وبالجمع على جميع قيم \(c\) نحصل على
$$\boxed{S(L)=\sum_{c=1}^{L}\left(C_{\mathrm{diff}}\,B(c)+\sum_{d=1}^{L} C_d\,T(c^d)\right).}$$
وبذلك تصبح احصاءات الازواج \((a,b)\) مستقلة عن \(c\)، ويمكن حسابها مرة واحدة فقط.
لنكتب \(M=987898789\). يمكن تمثيل العلاقة التراجعية بالمصفوفة المرافقة
$$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}.$$
بترديد \(M\) يكون فضاء الحالات منتهيا، ولذلك تصبح المتتالية دورية. التنفيذات تحسب الدورة الدقيقة من رتبة هذه المصفوفة. وبما ان \(M\) اولي، والمميز
$$\Delta=10^2+4=104$$
يحقق
$$\left(\frac{104}{M}\right)=1,$$
فان القيم الذاتية للمصفوفة تنتمي الى \(\mathbb{F}_M\)، ومن ثم رتبة \(A\) تقسم \(M-1\). لهذا السبب يبدا التنفيذ من المرشح \(M-1=987898788\)، ثم يفككه الى عوامله الاولية ويحاول حذف كل عامل اولي ما دام الاس الاصغر ما زال يعطي \(A^k=I\). في هذا الموديل لا يحدث اي اختزال، ولذلك
$$\pi=987898788.$$
وبالتالي يمكن اختزال كل فهرس ضخم وفق
$$T(n)\equiv T(n\bmod \pi)\pmod M.$$
ولكل \(c\) ثابت تتولد القوى المختزلة تكراريا عبر
$$p_1\equiv c,\qquad p_{d+1}\equiv p_d\,c \pmod \pi,$$
ولذلك لا تبني الخوارزمية العدد الضخم \(c^d\) صراحة.
عندما \(L=2\)، تكون الازواج المرتبة هي \((1,1)\)، \((1,2)\)، \((2,1)\)، \((2,2)\). التقييم 2-الادي المتساوي يظهر فقط في \((1,1)\) و\((2,2)\)، ولذلك
$$C_1=1,\qquad C_2=1,\qquad C_{\mathrm{diff}}=2.$$
عند \(c=1\) لدينا \(B(1)=10\) وكل قيمة \(T(1^d)=T(1)=10\)، ولذلك يكون الاسهام
$$2\cdot 10 + 1\cdot 10 + 1\cdot 10 = 40.$$
وعند \(c=2\) يكون \(B(2)=1\)، و\(T(2)=101\)، و\(T(4)=10301\)، ومن ثم
$$2\cdot 1 + 1\cdot 101 + 1\cdot 10301 = 10404.$$
وعليه
$$S(2)=40+10404=10444,$$
وهو تماما مقدار التحقق الذي تستعمله التنفيذات.
تنفيذات C++ وPython وJava تتبع الخطة الرياضية نفسها. اولا تحسب التقييم 2-الادي لكل عدد من \(1\) الى \(L\)، ثم تمر على جميع الازواج المرتبة \((a,b)\). اذا كان التقييمان متساويين يضاف الزوج الى الحاوية المفهرسة بواسطة \(\gcd(a,b)\)، واذا اختلفا يزداد فقط العداد المتبقي.
بعد ذلك تحسب الدورة بترديد \(987898789\)، وتبني جدول القوى الثنائية للمصفوفة المرافقة، ثم تمر على \(c=1,\dots,L\). لكل \(c\)، تتولد الفهارس المختزلة \(c,c^2,\dots,c^L\pmod \pi\) بضرب متكرر. ثم يحول كل فهرس مختزل الى \(T(\cdot)\bmod M\) بواسطة الاس المصفوفي السريع، ويضرب في وزنه المحسوب مسبقا ويضاف الى الجواب. تنفيذ C++ يوزع الحلقة الخارجية على عدة خيوط، لكن البنية العددية نفسها مشتركة بين اللغات الثلاث.
مرحلة عد الازواج تفحص جميع الازواج المرتبة \((a,b)\)، ولذلك تكلف \(O(L^2)\) زمنا. والمجموع الرئيسي يحتوي هو ايضا على \(L^2\) حد. حساب قيمة واحدة بعد اختزال الفهرس يستخدم جدولا ثابتا من مربعات المصفوفة، ولذلك تكون الكلفة النظرية لكل حد \(O(\log \pi)\)؛ وهنا \(\pi=987898788\) ثابتة بحكم الموديل، فيتحول ذلك عمليا الى عامل صغير ثابت. وبناء على ذلك فالتعقيد الكلي هو \(O(L^2\log \pi)\)، وهو فعليا \(O(L^2)\) في هذه المسألة، مع استهلاك ذاكرة \(O(L)\).