Let \(T(N)\) be the number of tours on a \(4\times N\) board that start at the upper-left square, end at the lower-left square, and visit every square exactly once. In graph-theoretic language, we are counting Hamiltonian paths on the \(4\times N\) grid graph with those two prescribed endpoints.
The Project Euler instance asks for \(T(10^{12}) \bmod 10^8\). A direct search is hopeless, and even stepping through a recurrence one column at a time would still take \(10^{12}\) iterations. The solution therefore needs two ingredients: a finite-width transfer argument to get a fixed recurrence, and fast matrix exponentiation to jump to the trillionth column length.
The decisive simplification is that the board is extremely long but only 4 cells tall. Fixed width means a left-to-right sweep can leave only finitely many possible boundary patterns, so the counting problem becomes a small linear dynamical system.
Sweep the board from left to right and imagine that the first \(n\) columns have already been processed. Every square strictly inside that processed region must already have its final degree in the Hamiltonian path, namely degree 2, except for the two prescribed endpoints in the first column. Any unfinished information can only live on the cut between columns \(n\) and \(n+1\).
That cut meets only 4 squares, so the unfinished path can expose only a tiny number of dangling connections. Because both global endpoints lie in the first column, the frontier can contain only \(0\), \(2\), or \(4\) dangling continuations to the right. In addition, any partial configuration that already forms a closed cycle, or traps a component so that it can never reconnect to the future, must be discarded immediately. These degree and connectivity rules leave only a small finite set of legal frontier profiles.
If we let \(S_n\) denote the vector counting those legal profiles after \(n\) columns, then adding one more column depends only on the current profile, not on the earlier history. Therefore there is a fixed transfer matrix \(B\) such that
$$S_{n+1}=B\,S_n.$$
One coordinate of \(S_n\) is the fully completed profile with no dangling edges; that coordinate is exactly \(T(n)\). The implementation does not store the whole frontier-state matrix explicitly. Instead it uses the smaller recurrence obtained after the auxiliary frontier counts are eliminated.
For this particular width and endpoint choice, the completed-tour count satisfies the fourth-order linear recurrence
$$T_n=2T_{n-1}+2T_{n-2}-2T_{n-3}+T_{n-4},\qquad n\ge 5,$$
with initial values
$$T_1=1,\qquad T_2=1,\qquad T_3=4,\qquad T_4=8.$$
The negative coefficient does not mean that tours are being counted negatively. It appears only after the auxiliary frontier states are algebraically removed from the transfer system. The underlying transfer counts remain nonnegative throughout; the sign change is a feature of the elimination step, not of the combinatorics itself.
Once the first four values are known, the recurrence generates the sequence immediately. For example,
$$T_5=2\cdot 8+2\cdot 4-2\cdot 1+1=23,$$
and then
$$T_6=2\cdot 23+2\cdot 8-2\cdot 4+1=55.$$
Continuing in the same way gives
$$T_7=144,\qquad T_8=360,\qquad T_9=921,\qquad T_{10}=2329.$$
That \(T_{10}=2329\) value is also used as a checkpoint in the implementations, so it serves as a convenient sanity test for both the mathematics and the code.
A fourth-order recurrence can be written as a first-order matrix recurrence by packaging four consecutive values together:
$$X_n=\begin{bmatrix}T_n\\T_{n-1}\\T_{n-2}\\T_{n-3}\end{bmatrix}.$$
Then
$$X_n= \begin{bmatrix} 2 & 2 & -2 & 1\\ 1 & 0 & 0 & 0\\ 0 & 1 & 0 & 0\\ 0 & 0 & 1 & 0 \end{bmatrix} X_{n-1}.$$
So for every \(N\ge 4\),
$$X_N=A^{N-4}X_4,\qquad X_4=\begin{bmatrix}8\\4\\1\\1\end{bmatrix}.$$
Now the huge exponent \(N-4\) can be handled by exponentiation by squaring, which needs only \(O(\log N)\) matrix multiplications. No closed form and no division are required; the method works directly modulo \(10^8\), even though that modulus is not prime.
The C++, Python, and Java implementations follow exactly the matrix formulation above. They return the four base values directly for \(N=1,2,3,4\). Otherwise they build the \(4\times 4\) companion matrix \(A\), initialize an identity matrix as the accumulator, and store the base vector \(X_4=[8,4,1,1]^T\).
They then raise \(A\) to the power \(N-4\) by binary exponentiation. Whenever the current binary digit of the exponent is 1, the accumulator is multiplied by the current matrix power. After that, the matrix is squared and the exponent is halved. Once the loop ends, the final matrix is applied to the base vector, and the first coordinate is the desired value \(T(N)\bmod M\).
All arithmetic is reduced modulo \(M\) after every multiply-add. That detail matters because the recurrence contains the coefficient \(-2\), so intermediate values can become negative before reduction. The implementations normalize such values back into the range \(0,\dots,M-1\). The C++ version also includes optional small-board checkpoint tests: exact brute-force counts for the tiniest boards and the known value \(T(10)=2329\).
The matrix size is constant, so both matrix-matrix multiplication and matrix-vector multiplication cost \(O(1)\) in asymptotic terms. Exponentiation by squaring performs \(O(\log N)\) such multiplications, giving a total running time of \(O(\log N)\).
The algorithm stores only a handful of \(4\times 4\) matrices and 4-entry vectors, so the extra memory usage is \(O(1)\). This is the key gain over a naive recurrence evaluation, which would still require \(O(N)\) steps and is therefore useless for \(N=10^{12}\).
Sei \(T(N)\) die Anzahl der Touren auf einem \(4\times N\)-Brett, die links oben beginnen, links unten enden und jedes Feld genau einmal besuchen. Graphentheoretisch zählen wir also Hamilton-Pfade im \(4\times N\)-Gittergraphen mit diesen beiden fest vorgegebenen Endpunkten.
Im Project-Euler-Fall ist \(T(10^{12}) \bmod 10^8\) gefragt. Eine direkte Suche ist ausgeschlossen, und selbst eine lineare Auswertung einer Rekurrenz über alle Spalten wäre bei \(10^{12}\) Schritten noch viel zu langsam. Die Lösung kombiniert deshalb ein Transferargument für feste Breite mit schneller Matrixpotenzierung.
Der entscheidende Punkt ist: Das Brett ist extrem lang, aber nur 4 Felder hoch. Bei fester Breite kann ein Links-nach-Rechts-Durchlauf nur endlich viele Randzustände hinterlassen. Damit wird aus dem kombinatorischen Problem ein kleines lineares Dynamiksystem.
Wir verarbeiten das Brett spaltenweise von links nach rechts und nehmen an, dass die ersten \(n\) Spalten bereits festgelegt sind. Jedes Feld im Inneren dieses bereits bearbeiteten Bereichs muss im endgültigen Hamilton-Pfad schon seinen endgültigen Grad 2 besitzen, abgesehen von den beiden vorgegebenen Endpunkten in der ersten Spalte. Alle noch offenen Informationen können also nur auf dem Schnitt zwischen Spalte \(n\) und Spalte \(n+1\) liegen.
Dieser Schnitt trifft nur 4 Felder. Deshalb kann der unfertige Pfad nur sehr wenige offene Verbindungen nach rechts besitzen. Da beide globalen Endpunkte schon in der ersten Spalte liegen, kann die Front nur \(0\), \(2\) oder \(4\) nach rechts fortzusetzende Enden enthalten. Außerdem müssen alle Teilkonfigurationen sofort verworfen werden, die bereits einen geschlossenen Zyklus bilden oder eine Komponente einschließen, die später nicht mehr sinnvoll verbunden werden kann. Aus diesen Grad- und Zusammenhangsbedingungen bleibt nur eine kleine endliche Menge zulässiger Randprofile übrig.
Bezeichnen wir mit \(S_n\) den Vektor der Anzahlen dieser zulässigen Profile nach \(n\) Spalten, dann hängt das Anhängen einer weiteren Spalte nur vom aktuellen Profil ab, nicht von der älteren Geschichte. Also existiert eine feste Transfermatrix \(B\) mit
$$S_{n+1}=B\,S_n.$$
Eine Komponente von \(S_n\) beschreibt den vollständig abgeschlossenen Zustand ohne offene Kanten; genau diese Komponente ist \(T(n)\). Die Implementierungen speichern jedoch nicht die volle Randzustandsmatrix, sondern die kleinere Rekurrenz, die nach Eliminierung der Hilfszustände entsteht.
Für genau diese Brettbreite und diese Endpunktwahl erfüllt die Zahl vollständiger Touren die Rekurrenz vierter Ordnung
$$T_n=2T_{n-1}+2T_{n-2}-2T_{n-3}+T_{n-4},\qquad n\ge 5,$$
mit den Anfangswerten
$$T_1=1,\qquad T_2=1,\qquad T_3=4,\qquad T_4=8.$$
Der negative Koeffizient bedeutet nicht, dass Touren negativ gezählt würden. Er entsteht erst dann, wenn man die Hilfsprofile algebraisch aus dem Transfersystem eliminiert. Die eigentlichen Transferzahlen bleiben nichtnegativ; das Minuszeichen ist nur ein Effekt der Eliminierung.
Sobald die ersten vier Werte bekannt sind, erzeugt die Rekurrenz die Folge sofort weiter. Zum Beispiel gilt
$$T_5=2\cdot 8+2\cdot 4-2\cdot 1+1=23,$$
und danach
$$T_6=2\cdot 23+2\cdot 8-2\cdot 4+1=55.$$
Setzt man das fort, erhält man
$$T_7=144,\qquad T_8=360,\qquad T_9=921,\qquad T_{10}=2329.$$
Der Wert \(T_{10}=2329\) wird auch als Kontrollwert in den Implementierungen verwendet und ist daher ein nützlicher Plausibilitätstest für Herleitung und Programm.
Eine Rekurrenz vierter Ordnung lässt sich als Rekurrenz erster Ordnung für einen 4-dimensionalen Zustandsvektor schreiben:
$$X_n=\begin{bmatrix}T_n\\T_{n-1}\\T_{n-2}\\T_{n-3}\end{bmatrix}.$$
Dann gilt
$$X_n= \begin{bmatrix} 2 & 2 & -2 & 1\\ 1 & 0 & 0 & 0\\ 0 & 1 & 0 & 0\\ 0 & 0 & 1 & 0 \end{bmatrix} X_{n-1}.$$
Also für jedes \(N\ge 4\):
$$X_N=A^{N-4}X_4,\qquad X_4=\begin{bmatrix}8\\4\\1\\1\end{bmatrix}.$$
Damit reduziert sich das riesige \(N-4\) auf Potenzieren per Exponentiation by Squaring, also auf nur \(O(\log N)\) Matrixmultiplikationen. Weder eine geschlossene Formel noch Divisionen sind nötig; das Verfahren arbeitet direkt modulo \(10^8\), obwohl dieser Modul nicht prim ist.
Die C++-, Python- und Java-Implementierungen setzen genau diese Matrixform um. Für \(N=1,2,3,4\) geben sie die vier Startwerte direkt zurück. Sonst bauen sie die \(4\times 4\)-Begleitmatrix \(A\), initialisieren eine Einheitsmatrix als Akkumulator und speichern den Basisvektor \(X_4=[8,4,1,1]^T\).
Danach wird \(A\) mit binärer Potenzierung auf die Potenz \(N-4\) gebracht. Immer wenn das aktuelle Bit des Exponenten gleich 1 ist, wird der Akkumulator mit der aktuellen Matrixpotenz multipliziert. Anschließend wird die Matrix quadriert und der Exponent halbiert. Nach Ende der Schleife wird die resultierende Matrix auf den Basisvektor angewendet; die erste Komponente ist genau \(T(N)\bmod M\).
Alle Rechenschritte werden nach jeder Multiplikation und Addition modulo \(M\) reduziert. Das ist wichtig, weil in der Rekurrenz der Koeffizient \(-2\) vorkommt und Zwischenwerte daher vorübergehend negativ werden können. Die Implementierungen normieren solche Werte zurück in den Bereich \(0,\dots,M-1\). Die C++-Version enthält zusätzlich kleine Kontrolltests: exakte Brute-Force-Zählungen auf den kleinsten Brettern sowie den bekannten Wert \(T(10)=2329\).
Die Matrixgröße ist konstant. Damit kosten sowohl Matrix-Matrix- als auch Matrix-Vektor-Multiplikationen asymptotisch \(O(1)\). Die binäre Potenzierung benötigt \(O(\log N)\) solcher Multiplikationen, also insgesamt \(O(\log N)\) Zeit.
Gespeichert werden nur wenige \(4\times 4\)-Matrizen und Vektoren der Länge 4, also \(O(1)\) zusätzlicher Speicher. Genau darin liegt der Vorteil gegenüber einer naiven Rekurrenzauswertung, die immer noch \(O(N)\) Schritte bräuchte und bei \(N=10^{12}\) unbrauchbar wäre.
\(T(N)\), \(4\times N\) tahtasında sol üst kareden başlayıp sol alt karede biten ve her kareyi tam bir kez ziyaret eden tur sayısı olsun. Graf kuramı dilinde bu, uç noktaları önceden belirlenmiş bir \(4\times N\) ızgara grafı üzerinde Hamilton yolu saymaktır.
Project Euler sürümünde istenen değer \(T(10^{12}) \bmod 10^8\)'dir. Doğrudan arama imkansızdır; hatta bir lineer yinelemeyi bile sütun sütun ilerleterek hesaplamak \(10^{12}\) adım gerektireceği için pratik değildir. Bu yüzden çözüm iki adıma dayanır: sabit genişlikten gelen sonlu bir geçiş modeli ve bu modelin logaritmik zamanda hesaplanmasını sağlayan matris üs alma tekniği.
Asıl sadeleştirme şuradan gelir: tahta çok uzun olabilir ama yüksekliği yalnızca 4'tür. Genişlik sabit olduğunda, soldan sağa süpürme yaparken sınırda görülebilecek durum sayısı sonludur. Böylece problem küçük boyutlu lineer bir dinamik sisteme iner.
Tahtayı soldan sağa işlediğimizi ve ilk \(n\) sütunun tamamlandığını düşünelim. Bu işlenmiş bölgenin tam içindeki her karenin, nihai Hamilton yolundaki son derecesi artık belirlenmiş olmalıdır; yani ilk sütundaki iki sabit uç dışında bütün bu karelerin derecesi 2'dir. O halde çözülmemiş tek bilgi, ancak \(n\). sütun ile \(n+1\). sütun arasındaki kesitte kalabilir.
Bu kesit yalnızca 4 kareyi kestiği için, sağa doğru devam eden açık bağlantı sayısı çok küçüktür. Üstelik iki küresel uç nokta da ilk sütunda bulunduğundan, sınırda yalnızca \(0\), \(2\) veya \(4\) tane sağa uzayan yarım bağlantı olabilir. Bunun yanında, daha bütün kareler kullanılmadan kapalı bir çevrim oluşturan ya da gelecekte tekrar bağlanması artık imkansız olan bir bileşen yaratan her kısmi düzen derhal elenir. Bu derece ve bağlantılılık kuralları sonunda elimizde yalnızca az sayıda geçerli sınır profili kalır.
İlk \(n\) sütundan sonra bu geçerli profillerin sayısını bir vektör olarak \(S_n\) ile gösterirsek, yeni bir sütun eklemek yalnızca mevcut profile bağlıdır; daha eski geçmişe bağlı değildir. Bu yüzden sabit bir geçiş matrisi \(B\) vardır ve
$$S_{n+1}=B\,S_n$$
yazabiliriz. \(S_n\)'nin bileşenlerinden biri, açık kenar bırakmayan tamamlanmış tur durumudur; bu bileşen tam olarak \(T(n)\)'dir. Uygulamalar bütün sınır-durumu matrisini saklamaz; onun yerine yardımcı durumlar elimine edildikten sonra elde edilen daha küçük yinelemeyi kullanır.
Bu tahta yüksekliği ve bu uç seçimi için tamamlanmış tur sayısı şu dördüncü dereceden lineer bağıntıyı sağlar:
$$T_n=2T_{n-1}+2T_{n-2}-2T_{n-3}+T_{n-4},\qquad n\ge 5,$$
başlangıç değerleri ise
$$T_1=1,\qquad T_2=1,\qquad T_3=4,\qquad T_4=8.$$
Buradaki negatif katsayı, turların negatif sayıldığı anlamına gelmez. Bu işaret yalnızca yardımcı sınır durumları cebirsel olarak sistemden çıkarıldığında ortaya çıkar. Alttaki geçiş sayıları baştan sona negatife düşmez; eksi işareti yalnızca eliminasyonun sonucudur.
İlk dört değer bilindikten sonra dizi kolayca ilerletilir. Örneğin
$$T_5=2\cdot 8+2\cdot 4-2\cdot 1+1=23,$$
ve ardından
$$T_6=2\cdot 23+2\cdot 8-2\cdot 4+1=55.$$
Aynı biçimde devam edilirse
$$T_7=144,\qquad T_8=360,\qquad T_9=921,\qquad T_{10}=2329$$
elde edilir. \(T_{10}=2329\) değeri uygulamalarda da kontrol noktası olarak kullanıldığı için, hem matematik hem de kod için güzel bir doğrulama örneğidir.
Dördüncü dereceden bir yineleme, ardışık dört değeri tek vektörde toplayarak birinci dereceden matris yinelemesine çevrilebilir:
$$X_n=\begin{bmatrix}T_n\\T_{n-1}\\T_{n-2}\\T_{n-3}\end{bmatrix}.$$
Böylece
$$X_n= \begin{bmatrix} 2 & 2 & -2 & 1\\ 1 & 0 & 0 & 0\\ 0 & 1 & 0 & 0\\ 0 & 0 & 1 & 0 \end{bmatrix} X_{n-1}$$
olur. Dolayısıyla \(N\ge 4\) için
$$X_N=A^{N-4}X_4,\qquad X_4=\begin{bmatrix}8\\4\\1\\1\end{bmatrix}.$$
Bundan sonra devasa \(N-4\) üssü ikili üs alma ile yalnızca \(O(\log N)\) matris çarpımıyla hesaplanabilir. Kapalı form gerekmez, bölme gerekmez; yöntem doğrudan \(10^8\) modunda çalışır ve modülün asal olmaması bir sorun yaratmaz.
C++, Python ve Java uygulamaları yukarıdaki matris formülünü doğrudan uygular. \(N=1,2,3,4\) için dört başlangıç değeri doğrudan döndürülür. Daha büyük \(N\) için \(4\times 4\) companion matrisi \(A\) kurulur, birim matris çarpan olarak hazırlanır ve taban vektörü \(X_4=[8,4,1,1]^T\) saklanır.
Sonra \(A\) matrisi \(N-4\) üssüne binary exponentiation ile yükseltilir. Üssün o andaki ikilik basamağı 1 ise biriken matris mevcut kuvvetle çarpılır. Ardından matris karesine çıkarılır ve üs ikiye bölünür. Döngü bittiğinde elde edilen matris taban vektörüne uygulanır; ilk bileşen aranan \(T(N)\bmod M\) değeridir.
Tüm işlemlerde her çarpma-toplama adımından sonra mod alınır. Bu nokta önemlidir, çünkü yinelemede \(-2\) katsayısı bulunduğundan ara değerler geçici olarak negatif olabilir. Uygulamalar bu değerleri \(0,\dots,M-1\) aralığına geri normalize eder. C++ sürümü ayrıca küçük doğrulama testleri de içerir: en küçük tahtalar için açık sayım ve bilinen \(T(10)=2329\) değeri.
Matris boyutu sabit olduğundan hem matris-matris hem de matris-vektör çarpımı asimptotik olarak \(O(1)\) maliyetlidir. İkili üs alma \(O(\log N)\) tane böyle çarpım yaptığı için toplam süre \(O(\log N)\) olur.
Algoritma yalnızca birkaç adet \(4\times 4\) matris ve 4 elemanlı vektör tuttuğu için ek bellek kullanımı \(O(1)\)'dir. Bu, sütun sütun ilerleyen saf yineleme hesabına göre büyük fark yaratır; o yöntem \(O(N)\) adım ister ve \(N=10^{12}\) için kullanışsızdır.
Sea \(T(N)\) el número de recorridos sobre un tablero \(4\times N\) que empiezan en la casilla superior izquierda, terminan en la inferior izquierda y visitan cada casilla exactamente una vez. En lenguaje de teoría de grafos, estamos contando caminos hamiltonianos del grafo de rejilla \(4\times N\) con esos dos extremos fijados.
La instancia de Project Euler pide \(T(10^{12}) \bmod 10^8\). Una búsqueda directa es imposible, e incluso avanzar una recurrencia columna por columna seguiría necesitando \(10^{12}\) pasos. Por eso la solución combina un argumento de transferencia para anchura fija con exponenciación matricial rápida.
La simplificación decisiva es que el tablero puede ser larguísimo, pero solo tiene altura 4. Cuando la anchura es fija, un barrido de izquierda a derecha solo puede dejar un número finito de patrones en la frontera. Eso convierte el problema en un pequeño sistema lineal.
Imaginemos que procesamos el tablero de izquierda a derecha y que las primeras \(n\) columnas ya están decididas. Toda casilla situada estrictamente dentro de esa región procesada ya debe tener su grado final dentro del camino hamiltoniano, es decir grado 2, salvo los dos extremos prescritos de la primera columna. Por tanto, la única información pendiente puede vivir en el corte entre las columnas \(n\) y \(n+1\).
Ese corte atraviesa solo 4 casillas, así que el camino parcial solo puede dejar muy pocas conexiones abiertas hacia la derecha. Como ambos extremos globales están en la primera columna, la frontera solo puede contener \(0\), \(2\) o \(4\) continuaciones colgantes hacia la derecha. Además, cualquier configuración parcial que ya forme un ciclo cerrado, o que encierre una componente que luego no pueda reconectarse, debe descartarse de inmediato. Estas restricciones de grado y conectividad dejan solo un pequeño conjunto finito de perfiles legales.
Si \(S_n\) denota el vector que cuenta esos perfiles legales tras \(n\) columnas, entonces añadir una columna nueva depende solo del perfil actual, no de toda la historia previa. Por tanto existe una matriz de transferencia fija \(B\) tal que
$$S_{n+1}=B\,S_n.$$
Una de las coordenadas de \(S_n\) corresponde al perfil completamente cerrado, sin aristas colgantes; esa coordenada es exactamente \(T(n)\). Las implementaciones no almacenan la matriz completa de estados de frontera, sino la recurrencia más pequeña que se obtiene al eliminar los estados auxiliares.
Para esta anchura y esta elección de extremos, el número de recorridos completos satisface la recurrencia lineal de orden cuatro
$$T_n=2T_{n-1}+2T_{n-2}-2T_{n-3}+T_{n-4},\qquad n\ge 5,$$
con valores iniciales
$$T_1=1,\qquad T_2=1,\qquad T_3=4,\qquad T_4=8.$$
El coeficiente negativo no significa que se estén contando recorridos de forma negativa. Aparece solo después de eliminar algebraicamente los estados auxiliares del sistema de transferencia. Los conteos subyacentes siguen siendo no negativos; el signo menos es un efecto de la eliminación, no de la combinatoria básica.
Una vez conocidos los primeros cuatro valores, la recurrencia genera la sucesión de inmediato. Por ejemplo,
$$T_5=2\cdot 8+2\cdot 4-2\cdot 1+1=23,$$
y después
$$T_6=2\cdot 23+2\cdot 8-2\cdot 4+1=55.$$
Continuando del mismo modo se obtiene
$$T_7=144,\qquad T_8=360,\qquad T_9=921,\qquad T_{10}=2329.$$
Ese valor \(T_{10}=2329\) también aparece como comprobación en las implementaciones, así que sirve como ejemplo concreto y como prueba de consistencia.
Una recurrencia de orden cuatro puede escribirse como una recurrencia matricial de primer orden si agrupamos cuatro valores consecutivos:
$$X_n=\begin{bmatrix}T_n\\T_{n-1}\\T_{n-2}\\T_{n-3}\end{bmatrix}.$$
Entonces
$$X_n= \begin{bmatrix} 2 & 2 & -2 & 1\\ 1 & 0 & 0 & 0\\ 0 & 1 & 0 & 0\\ 0 & 0 & 1 & 0 \end{bmatrix} X_{n-1}.$$
Por tanto, para todo \(N\ge 4\),
$$X_N=A^{N-4}X_4,\qquad X_4=\begin{bmatrix}8\\4\\1\\1\end{bmatrix}.$$
Ahora el enorme exponente \(N-4\) se maneja mediante exponenciación binaria, que solo necesita \(O(\log N)\) multiplicaciones de matrices. No hace falta una fórmula cerrada ni dividir por nada; el método funciona directamente módulo \(10^8\), aunque ese módulo no sea primo.
Las implementaciones en C++, Python y Java siguen exactamente la formulación matricial anterior. Para \(N=1,2,3,4\) devuelven directamente los cuatro valores base. En otro caso construyen la matriz compañera \(A\), inicializan una matriz identidad como acumulador y guardan el vector base \(X_4=[8,4,1,1]^T\).
Después elevan \(A\) a la potencia \(N-4\) mediante exponenciación binaria. Cada vez que el bit actual del exponente vale 1, multiplican el acumulador por la potencia actual de \(A\). Luego elevan la matriz al cuadrado y dividen el exponente entre 2. Al terminar el bucle, aplican la matriz resultante al vector base y la primera coordenada es el valor buscado \(T(N)\bmod M\).
Toda la aritmética se reduce módulo \(M\) después de cada multiplicación y suma. Esto importa porque la recurrencia contiene el coeficiente \(-2\), de modo que algunos valores intermedios pueden volverse negativos antes de la reducción. Las implementaciones normalizan esas cantidades al rango \(0,\dots,M-1\). La versión en C++ también incluye comprobaciones pequeñas: conteos exactos por fuerza bruta para los tableros más pequeños y el valor conocido \(T(10)=2329\).
El tamaño de la matriz es constante, así que tanto la multiplicación matriz-matriz como la multiplicación matriz-vector cuestan \(O(1)\) en términos asintóticos. La exponenciación binaria realiza \(O(\log N)\) multiplicaciones de ese tipo, por lo que el tiempo total es \(O(\log N)\).
El algoritmo solo guarda unas pocas matrices \(4\times 4\) y vectores de longitud 4, de modo que el uso extra de memoria es \(O(1)\). Esa es la gran mejora frente a evaluar la recurrencia paso a paso, lo que seguiría costando \(O(N)\) y no serviría para \(N=10^{12}\).
设 \(T(N)\) 表示在 \(4\times N\) 棋盘上的巡游条数：路径从左上角出发，在左下角结束，并且每个方格恰好经过一次。用图论语言来说，这就是在 \(4\times N\) 的网格图上，统计具有这两个指定端点的 Hamilton 路径数量。
Project Euler 这一题要求计算 \(T(10^{12}) \bmod 10^8\)。直接搜索显然不可能完成；即使已经知道一个线性递推，如果还按列从 \(1\) 一直推到 \(10^{12}\)，时间也完全不可接受。因此解法必须先利用固定宽度得到一个小规模状态转移，再用快速矩阵幂把指数级长度压缩成对数级计算。
最关键的观察是：棋盘可以非常长，但高度始终只有 4。高度固定意味着从左到右扫描时，切口上可能出现的边界形状是有限的，于是原问题可以压缩成一个很小的线性系统。
想象按列从左向右处理棋盘，并假设前 \(n\) 列已经完全处理完毕。对于这个已处理区域内部的每个格子，除了第一列中预先指定的两个端点之外，它在最终 Hamilton 路径中的度数都必须已经固定为 2。因此，所有尚未决定的信息只能出现在第 \(n\) 列与第 \(n+1\) 列之间的那条竖直切口上。
这条切口只穿过 4 个格子，所以未完成路径向右延伸的“悬空连接”数量非常有限。又因为全局的两个端点都位于第一列，所以边界上只能出现 \(0\)、\(2\) 或 \(4\) 个继续向右的悬空端。与此同时，任何部分构型如果已经形成了一个提前闭合的环，或者把某个连通块困死在左侧、以后再也无法与右侧重新连通，就必须立刻判为非法。度数约束和连通性约束一起，最终只留下很少几个合法的边界轮廓。
记 \(S_n\) 为处理完前 \(n\) 列之后，各种合法边界轮廓的计数向量。再加入一列时，新的计数只依赖当前边界轮廓，而与更早之前的细节无关，因此存在一个固定转移矩阵 \(B\)，满足
$$S_{n+1}=B\,S_n.$$
\(S_n\) 的某一个分量对应“没有任何悬空边、当前已经形成完整巡游”的状态，这个分量正是 \(T(n)\)。实现代码并没有显式保存完整的边界状态矩阵，而是使用把辅助状态消去之后得到的更小递推。
对于本题这种高度和端点安排，完整巡游数满足四阶线性递推
$$T_n=2T_{n-1}+2T_{n-2}-2T_{n-3}+T_{n-4},\qquad n\ge 5,$$
初值为
$$T_1=1,\qquad T_2=1,\qquad T_3=4,\qquad T_4=8.$$
这里出现的负系数并不表示“负数条路径”。真正的边界状态转移计数始终是非负的；\(-2T_{n-3}\) 只是把若干辅助状态从线性系统中代数消元以后出现的结果。因此负号属于公式整理的产物，而不是组合计数本身的异常现象。
一旦前四项已知，后面的值就可以直接由递推推出。例如
$$T_5=2\cdot 8+2\cdot 4-2\cdot 1+1=23,$$
接着
$$T_6=2\cdot 23+2\cdot 8-2\cdot 4+1=55.$$
继续向后可得
$$T_7=144,\qquad T_8=360,\qquad T_9=921,\qquad T_{10}=2329.$$
\(T_{10}=2329\) 这个数值也被实现代码用作检查点，因此它既是递推的直接示例，也是验证程序正确性的方便样例。
四阶递推可以改写成一阶矩阵递推，只要把连续四项打包成一个向量：
$$X_n=\begin{bmatrix}T_n\\T_{n-1}\\T_{n-2}\\T_{n-3}\end{bmatrix}.$$
于是有
$$X_n= \begin{bmatrix} 2 & 2 & -2 & 1\\ 1 & 0 & 0 & 0\\ 0 & 1 & 0 & 0\\ 0 & 0 & 1 & 0 \end{bmatrix} X_{n-1}.$$
因此对所有 \(N\ge 4\)，
$$X_N=A^{N-4}X_4,\qquad X_4=\begin{bmatrix}8\\4\\1\\1\end{bmatrix}.$$
这样一来，原本巨大的指数 \(N-4\) 就可以用二分幂在 \(O(\log N)\) 次矩阵乘法内完成。整个过程只需要加法和乘法，不需要求闭式，也不需要除法；即使模数 \(10^8\) 不是素数，方法仍然完全适用。
C++、Python 和 Java 实现都严格按照上面的矩阵模型来写。对于 \(N=1,2,3,4\)，程序直接返回四个初值。否则先建立 \(4\times 4\) 的伴随矩阵 \(A\)，再用单位矩阵作为累积器，并准备好基向量 \(X_4=[8,4,1,1]^T\)。
接下来用二进制快速幂计算 \(A^{N-4}\)。如果当前指数的二进制最低位是 1，就把累积器乘上当前矩阵；然后把矩阵自乘一次，同时把指数右移一位。循环结束后，再把得到的矩阵作用到基向量上，结果向量的第一项就是所求的 \(T(N)\bmod M\)。
所有乘加操作之后都会立刻取模 \(M\)。这一步非常重要，因为递推里含有 \(-2\) 这个系数，中间结果在取模前可能暂时变成负数。实现会把这些值重新规范到 \(0,\dots,M-1\) 的范围内。C++ 版本还额外加入了很小规模的正确性检查：对最小棋盘做精确暴力计数，并核对已知值 \(T(10)=2329\)。
矩阵维度固定为 4，所以矩阵与矩阵相乘、矩阵与向量相乘在渐近意义下都属于 \(O(1)\) 操作。快速幂需要做 \(O(\log N)\) 次这样的乘法，因此总时间复杂度是 \(O(\log N)\)。
程序只保存少量 \(4\times 4\) 矩阵和长度为 4 的向量，所以额外空间复杂度是 \(O(1)\)。这正是它相对于朴素递推的优势所在：朴素方法虽然已经比暴力搜索好得多，但仍要 \(O(N)\) 步，对于 \(N=10^{12}\) 依然无法使用。
Пусть \(T(N)\) обозначает число маршрутов на доске \(4\times N\), которые начинаются в верхней левой клетке, заканчиваются в нижней левой клетке и посещают каждую клетку ровно один раз. С точки зрения теории графов это число гамильтоновых путей в прямоугольной решетке \(4\times N\) с двумя заранее фиксированными концами.
В задаче Project Euler требуется найти \(T(10^{12}) \bmod 10^8\). Полный перебор невозможен, а даже пошаговое вычисление линейной рекурсии по одной колонке за раз все равно потребовало бы \(10^{12}\) шагов. Поэтому решение строится в два слоя: сначала из фиксированной высоты выводится конечная система переходов, а затем эта система ускоряется матричным возведением в степень.
Главное упрощение состоит в том, что доска может быть очень длинной, но ее высота всегда равна 4. При фиксированной высоте вертикальный разрез между обработанной и необработанной частью может иметь лишь конечное число допустимых конфигураций, и задача сводится к маленькой линейной динамической системе.
Будем мысленно обрабатывать доску слева направо и предположим, что первые \(n\) столбцов уже разобраны. Каждая клетка строго внутри этой обработанной области уже должна иметь окончательную степень в гамильтоновом пути, а именно степень 2, кроме двух заданных концов в первом столбце. Значит, вся еще не определенная информация может находиться только на разрезе между столбцами \(n\) и \(n+1\).
Этот разрез проходит всего через 4 клетки, поэтому число незавершенных продолжений пути вправо очень мало. Поскольку оба глобальных конца расположены в первом столбце, на границе может быть только \(0\), \(2\) или \(4\) висящих продолжения вправо. Кроме того, любую частичную конфигурацию, которая уже образует замкнутый цикл или запирает компоненту так, что ее невозможно будет снова соединить с будущими столбцами, нужно сразу отбросить. Из ограничений на степени и связность остается лишь небольшой конечный набор допустимых профилей границы.
Обозначим через \(S_n\) вектор количеств этих допустимых профилей после \(n\) столбцов. Тогда добавление еще одного столбца зависит только от текущего профиля, а не от всей предыстории, поэтому существует фиксированная матрица перехода \(B\), для которой
$$S_{n+1}=B\,S_n.$$
Одна из координат \(S_n\) соответствует полностью завершенному профилю без висящих ребер; именно эта координата и равна \(T(n)\). В реализации, однако, хранится не полная матрица граничных состояний, а более компактная рекурсия, полученная после исключения вспомогательных состояний.
Для данной высоты и данного выбора концов число завершенных маршрутов удовлетворяет линейной рекурсии четвертого порядка
$$T_n=2T_{n-1}+2T_{n-2}-2T_{n-3}+T_{n-4},\qquad n\ge 5,$$
с начальными значениями
$$T_1=1,\qquad T_2=1,\qquad T_3=4,\qquad T_4=8.$$
Отрицательный коэффициент не означает, что какие-то маршруты считаются с минусом. Он появляется только после алгебраического исключения вспомогательных профилей из системы переходов. Сами базовые количества переходов неотрицательны; знак минус является следствием устранения промежуточных состояний.
Как только известны первые четыре значения, рекурсия сразу продолжает последовательность. Например,
$$T_5=2\cdot 8+2\cdot 4-2\cdot 1+1=23,$$
а затем
$$T_6=2\cdot 23+2\cdot 8-2\cdot 4+1=55.$$
Если продолжить вычисление, получаем
$$T_7=144,\qquad T_8=360,\qquad T_9=921,\qquad T_{10}=2329.$$
Значение \(T_{10}=2329\) используется и в коде как контрольная точка, так что это одновременно и наглядный пример работы рекурсии, и удобная проверка корректности.
Рекурсию четвертого порядка удобно переписать как рекурсию первого порядка для вектора из четырех последовательных значений:
$$X_n=\begin{bmatrix}T_n\\T_{n-1}\\T_{n-2}\\T_{n-3}\end{bmatrix}.$$
Тогда
$$X_n= \begin{bmatrix} 2 & 2 & -2 & 1\\ 1 & 0 & 0 & 0\\ 0 & 1 & 0 & 0\\ 0 & 0 & 1 & 0 \end{bmatrix} X_{n-1}.$$
Следовательно, для любого \(N\ge 4\)
$$X_N=A^{N-4}X_4,\qquad X_4=\begin{bmatrix}8\\4\\1\\1\end{bmatrix}.$$
Теперь огромный показатель \(N-4\) вычисляется методом бинарного возведения в степень, который требует всего \(O(\log N)\) матричных умножений. Замкнутая формула не нужна, деления тоже нет; метод напрямую работает по модулю \(10^8\), хотя этот модуль не является простым.
Реализации на C++, Python и Java буквально следуют приведенной матричной форме. Для \(N=1,2,3,4\) они сразу возвращают четыре базовых значения. В остальных случаях строится матрица-компаньон \(A\) размера \(4\times 4\), матрица-аккумулятор инициализируется единичной матрицей, а в качестве стартового вектора берется \(X_4=[8,4,1,1]^T\).
Далее вычисляется степень \(A^{N-4}\) с помощью двоичного разложения показателя. Если текущий бит показателя равен 1, аккумулятор умножается на текущую степень матрицы. Затем матрица возводится в квадрат, а показатель делится пополам. После завершения цикла результирующая матрица применяется к стартовому вектору, и первая координата дает ответ \(T(N)\bmod M\).
После каждого умножения и сложения выполняется взятие по модулю \(M\). Это важно из-за коэффициента \(-2\) в рекурсии: промежуточные значения могут временно стать отрицательными, и их нужно вернуть в диапазон \(0,\dots,M-1\). В версии на C++ есть также небольшие контрольные проверки: точный перебор для самых маленьких досок и сверка со значением \(T(10)=2329\).
Размер матрицы постоянен, поэтому умножение матрицы на матрицу и матрицы на вектор имеет асимптотическую стоимость \(O(1)\). Метод быстрого возведения в степень делает \(O(\log N)\) таких операций, значит общее время равно \(O(\log N)\).
Алгоритм хранит лишь несколько матриц \(4\times 4\) и несколько векторов длины 4, так что дополнительная память равна \(O(1)\). В этом и заключается выигрыш по сравнению с наивным пошаговым вычислением рекурсии, которое по-прежнему стоило бы \(O(N)\) и было бы бесполезно при \(N=10^{12}\).
لتكن \(T(N)\) عدد المسارات على لوحة \(4\times N\) التي تبدأ من المربع العلوي الأيسر، وتنتهي في المربع السفلي الأيسر، وتزور كل مربع مرة واحدة بالضبط. بلغة نظرية الرسوم البيانية نحن نعد المسارات الهاملتونية في شبكة \(4\times N\) مع تثبيت نقطتي البداية والنهاية مسبقًا.
في نسخة Project Euler المطلوب هو حساب \(T(10^{12}) \bmod 10^8\). البحث المباشر مستحيل عمليًا، وحتى لو عرفنا علاقة تراجعية خطية فإن تطبيقها عمودًا بعد عمود سيظل يحتاج إلى \(10^{12}\) خطوة. لذلك يعتمد الحل على خطوتين مترابطتين: اختزال المسألة إلى عدد محدود من حالات الحد بسبب ثبات الارتفاع، ثم استخدام رفع المصفوفات السريع للوصول إلى طول هائل في زمن لوغاريتمي.
الفكرة الحاسمة هي أن اللوحة قد تكون طويلة جدًا، لكنها لا ترتفع إلا 4 خلايا فقط. وعندما يكون العرض ثابتًا بهذا الشكل، فإن المسح من اليسار إلى اليمين لا يمكن أن يترك إلا عددًا محدودًا من الأنماط على الحد الفاصل، وبذلك تتحول المسألة إلى نظام خطي صغير.
نتخيل أننا نعالج اللوحة عمودًا بعد عمود من اليسار إلى اليمين، وأن الأعمدة \(1\) حتى \(n\) قد حُسمت بالفعل. كل مربع يقع داخل الجزء المعالَج تمامًا يجب أن يكون قد أخذ درجته النهائية في المسار الهاملتوني، أي الدرجة 2، باستثناء النهايتين المحددتين في العمود الأول. لذلك فكل المعلومات غير المحسومة لا يمكن أن تظهر إلا على القطع العمودي بين العمود \(n\) والعمود \(n+1\).
هذا القطع يمر عبر 4 خلايا فقط، لذا فإن عدد الامتدادات المفتوحة نحو اليمين صغير جدًا. وبما أن طرفي المسار الكلي موجودان أصلًا في العمود الأول، فإن الحد لا يمكن أن يحمل إلا \(0\) أو \(2\) أو \(4\) نهايات معلقة تتجه إلى اليمين. كذلك يجب حذف أي تشكيل جزئي يُغلق دورة مبكرًا أو يعزل مكوّنًا لن يعود قادرًا على الاتصال بما تبقى من اللوحة. قيود الدرجات والاتصال هذه تترك لنا مجموعة صغيرة ومحدودة من بروفايلات الحدود الممكنة.
إذا رمزنا بـ \(S_n\) إلى متجه أعداد هذه البروفايلات القانونية بعد معالجة \(n\) أعمدة، فإن إضافة عمود جديد تعتمد فقط على بروفايل الحد الحالي، لا على التاريخ الكامل السابق. لذلك توجد مصفوفة انتقال ثابتة \(B\) تحقق
$$S_{n+1}=B\,S_n.$$
أحد مركبات \(S_n\) يمثل الحالة المكتملة تمامًا من دون أي حواف معلقة، وهذا المركب هو بالضبط \(T(n)\). التطبيق لا يخزن مصفوفة جميع حالات الحد صراحة، بل يستخدم العلاقة الأصغر الناتجة بعد حذف الحالات المساعدة جبريًا.
لهذا الارتفاع المحدد ولهذا الاختيار من النهايات، يحقق عدد الجولات المكتملة العلاقة الخطية من الرتبة الرابعة
$$T_n=2T_{n-1}+2T_{n-2}-2T_{n-3}+T_{n-4},\qquad n\ge 5,$$
مع القيم الابتدائية
$$T_1=1,\qquad T_2=1,\qquad T_3=4,\qquad T_4=8.$$
المعامل السالب هنا لا يعني أن بعض الجولات تُحسب بقيم سالبة. هذا الحد يظهر فقط بعد حذف الحالات المساعدة من نظام الانتقال على نحو جبري. أما العدود الأساسية نفسها فتبقى غير سالبة طوال الوقت؛ والإشارة السالبة مجرد أثر لعملية الاختزال.
ما إن تصبح القيم الأربع الأولى معروفة حتى تولد العلاقة بقية المتتالية مباشرة. على سبيل المثال
$$T_5=2\cdot 8+2\cdot 4-2\cdot 1+1=23,$$
ثم
$$T_6=2\cdot 23+2\cdot 8-2\cdot 4+1=55.$$
وبالاستمرار بالطريقة نفسها نحصل على
$$T_7=144,\qquad T_8=360,\qquad T_9=921,\qquad T_{10}=2329.$$
القيمة \(T_{10}=2329\) تُستخدم أيضًا كنقطة تحقق داخل التطبيقات، لذا فهي مثال عملي على صحة الصيغة وعلى سلامة التنفيذ معًا.
يمكن كتابة علاقة من الرتبة الرابعة على صورة علاقة مصفوفية من الرتبة الأولى إذا جمعنا أربع قيم متتالية في متجه واحد:
$$X_n=\begin{bmatrix}T_n\\T_{n-1}\\T_{n-2}\\T_{n-3}\end{bmatrix}.$$
وعندئذ يصبح
$$X_n= \begin{bmatrix} 2 & 2 & -2 & 1\\ 1 & 0 & 0 & 0\\ 0 & 1 & 0 & 0\\ 0 & 0 & 1 & 0 \end{bmatrix} X_{n-1}.$$
إذن لكل \(N\ge 4\) لدينا
$$X_N=A^{N-4}X_4,\qquad X_4=\begin{bmatrix}8\\4\\1\\1\end{bmatrix}.$$
وهكذا يمكن التعامل مع الأس الهائل \(N-4\) بواسطة الرفع بالتربيع، الذي يحتاج فقط إلى \(O(\log N)\) من ضرب المصفوفات. لا نحتاج إلى صيغة مغلقة ولا إلى أي قسمة؛ والطريقة تعمل مباشرة تحت الحساب modulo \(10^8\) رغم أن هذا الموديلو ليس أوليًا.
تنفذ نسخ C++ وPython وJava الصياغة المصفوفية نفسها حرفيًا. إذا كان \(N=1,2,3,4\) تُعاد القيم الابتدائية الأربع مباشرة. وإلا تُبنى مصفوفة المرافقة \(A\) ذات البعد \(4\times 4\)، وتُهيأ مصفوفة الوحدة كمُراكِم، ويُخزَّن المتجه الأساسي \(X_4=[8,4,1,1]^T\).
بعد ذلك تُحسب \(A^{N-4}\) باستخدام الرفع الثنائي. كلما كانت الخانة الثنائية الحالية للأس تساوي 1 ضُرب المُراكِم في القوة الحالية للمصفوفة. ثم تُربَّع المصفوفة ويُنصَّف الأس. وعند انتهاء الحلقة تُطبَّق المصفوفة الناتجة على المتجه الأساسي، ويكون المركب الأول هو القيمة المطلوبة \(T(N)\bmod M\).
كل عملية ضرب وجمع يتبعها أخذ الباقي modulo \(M\). هذه النقطة مهمة لأن العلاقة تحتوي على المعامل \(-2\)، ولذلك قد تصبح بعض القيم الوسيطة سالبة مؤقتًا قبل التطبيع. تقوم التطبيقات بإعادة هذه القيم إلى المجال \(0,\dots,M-1\). كما تتضمن نسخة C++ اختبارات تحقق صغيرة: عدًا صريحًا لأصغر الألواح، إضافة إلى مطابقة القيمة المعروفة \(T(10)=2329\).
بما أن أبعاد المصفوفة ثابتة، فإن ضرب المصفوفات وضرب المصفوفة في متجه كلاهما يكلّف \(O(1)\) من الناحية التقاربية. والرفع بالتربيع ينفذ \(O(\log N)\) من هذه العمليات، لذا يصبح الزمن الكلي \(O(\log N)\).
لا تخزن الخوارزمية إلا عددًا قليلًا من المصفوفات \(4\times 4\) والمتجهات ذات الطول 4، ولذلك فإن الذاكرة الإضافية هي \(O(1)\). وهذا هو المكسب الحقيقي مقارنة بالتقييم الساذج للعلاقة التراجعية خطوة خطوة، لأنه سيظل يحتاج إلى \(O(N)\) عملية ولن يكون مناسبًا عندما يكون \(N=10^{12}\).