Problem Summary

For integers \(n \ge 2\) and \(m \ge 1\), let \(T(n,m)\) be the number of length-\(m\) sequences \((a_1,\dots,a_m)\) such that

$$1 \le a_i \le n-1 \quad (1 \le i \le m),\qquad a_i+a_{i+1} \le n \quad (1 \le i \le m-1).$$

The task is to evaluate \(T(5000,10^{12})\) modulo \(10^9+7\). A direct dynamic program over all lengths up to \(10^{12}\) is impossible, so the solution first studies the finite transition system behind the neighbour rule, then reconstructs a short linear recurrence and jumps directly to the required term.

Mathematical Approach

The admissible sequences form a finite-state process on the values \(1,2,\dots,n-1\). Once that process is written in linear-algebra form, the distant term can be recovered from a much shorter recurrence.

Step 1: Encode the Neighbour Rule as a Transition Matrix

Use the last entry of the sequence as the state. If the current value is \(x\) and the next value is \(y\), the pair is valid exactly when

$$x+y \le n.$$

This gives a matrix \(A\in\{0,1\}^{(n-1)\times(n-1)}\) defined by

$$A_{y,x}=\begin{cases} 1,& x+y\le n,\\ 0,& x+y\gt n. \end{cases}$$

Each row is a consecutive block of ones followed by zeros. The first row contains \(n-1\) ones, the next row contains \(n-2\), and so on down to a single one. This triangular pattern is what makes the warm-up computation efficient.

Step 2: Write the Dynamic Programming Recurrence

Let \(f_t(y)\) be the number of valid sequences of length \(t\) ending at \(y\). For length \(1\), every state is allowed, so

$$f_1(y)=1 \qquad (1\le y\le n-1).$$

For \(t\ge 1\), the previous value can be any \(x\) with \(x\le n-y\), hence

$$f_{t+1}(y)=\sum_{x=1}^{n-y} f_t(x).$$

If we define prefix sums

$$P_t(k)=\sum_{x=1}^{k} f_t(x),$$

then the transition becomes

$$f_{t+1}(y)=P_t(n-y).$$

That identity explains the reverse-prefix update used by the implementation: one full layer can be produced in linear time. The total number of valid sequences of length \(t\) is

$$T(n,t)=\sum_{y=1}^{n-1} f_t(y).$$

In vector form, with \(F_t=(f_t(1),\dots,f_t(n-1))^T\) and \(\mathbf{1}\) the all-ones vector,

$$F_t=A^{t-1}\mathbf{1},\qquad T(n,t)=\mathbf{1}^T A^{t-1}\mathbf{1}.$$

Step 3: Why a Linear Recurrence Must Exist

The matrix \(A\) has size \(n-1\), so its minimal polynomial has degree at most \(n-1\). By the Cayley-Hamilton principle, every scalar sequence extracted from powers of \(A\), including \(T(n,t)\), must satisfy a linear recurrence of some order \(r\le n-1\):

$$T(n,k+r)=c_1T(n,k+r-1)+c_2T(n,k+r-2)+\cdots+c_rT(n,k)\pmod{10^9+7}.$$

This is the key structural fact. Instead of iterating the dynamic program all the way to \(m\), we only need enough initial terms to recover the recurrence and then evaluate the distant term quickly.

Step 4: Recover the Shortest Recurrence from Initial Terms

The implementation generates an initial segment

$$T(n,1),T(n,2),T(n,3),\dots$$

modulo \(10^9+7\), and then applies the Berlekamp-Massey algorithm to find the shortest linear recurrence that matches those values. If the recovered order is \(r\), its characteristic relation can be written as

$$x^r-c_1x^{r-1}-c_2x^{r-2}-\cdots-c_r=0.$$

Because \(r\le n-1\), any sample length comfortably above \(2(n-1)\) is enough for recovery; the implementations use \(2(n-1)+8\) warm-up terms.

Step 5: Compute the Distant Term by Polynomial Reduction

Once the recurrence is known, the problem becomes a standard fast nth-term computation. Reduce \(x^{m-1}\) modulo the recurrence polynomial:

$$x^{m-1}\equiv q_0+q_1x+\cdots+q_{r-1}x^{r-1}\pmod{x^r-c_1x^{r-1}-\cdots-c_r}.$$

Then the desired term is

$$T(n,m)\equiv q_0T(n,1)+q_1T(n,2)+\cdots+q_{r-1}T(n,r)\pmod{10^9+7}.$$

The coefficients \(q_i\) are obtained by binary exponentiation inside the quotient ring defined by the recurrence, so the huge exponent \(m=10^{12}\) contributes only a \(\log m\) factor.

Worked Example: \(n=3\)

Now the allowed values are \(1\) and \(2\). The neighbour rule forbids only the pair \((2,2)\), so

$$A=\begin{pmatrix} 1 & 1\\ 1 & 0 \end{pmatrix}.$$

The totals by length are

$$T(3,1)=2,\qquad T(3,2)=3,\qquad T(3,3)=5,\qquad T(3,4)=8.$$

In this case the recurrence is simply

$$T(3,t)=T(3,t-1)+T(3,t-2),$$

so the checkpoint \(T(3,4)=8\) follows immediately.

How the Code Works

The C++, Python, and Java implementations all follow the same mathematical pipeline. First they generate an initial block of totals using the end-state dynamic program. Instead of recomputing every state transition naively, they build prefix sums of the current layer and read those sums in reverse order, which turns one transition step into an \(O(n)\) operation.

Next they run Berlekamp-Massey modulo \(10^9+7\) on the warm-up totals and obtain the shortest recurrence for \(T(n,m)\). If the requested length already lies inside the precomputed prefix, the answer is returned directly. Otherwise the implementation performs binary exponentiation on reduced polynomials, producing the coefficients \(q_0,\dots,q_{r-1}\) and therefore the required distant term.

Complexity Analysis

Each warm-up dynamic-programming step costs \(O(n)\) time and \(O(n)\) memory. Since only \(O(n)\) initial terms are generated, the warm-up phase costs \(O(n^2)\) time overall. Berlekamp-Massey on a recurrence of order \(r\le n-1\) costs \(O(r^2)\), and the final nth-term stage costs \(O(r^2\log m)\). Therefore the total complexity is

$$O(n^2+r^2+r^2\log m)=O(n^2\log m),$$

with \(O(n)\) memory.

Footnotes and References

  1. Project Euler problem page: Problem 654: Neighbourly Constraints
  2. Berlekamp-Massey algorithm: Wikipedia - Berlekamp-Massey algorithm
  3. Fast computation of linear recurrences: cp-algorithms - Linear Recurrence
  4. Cayley-Hamilton theorem: Wikipedia - Cayley-Hamilton theorem
  5. Prefix sums: Wikipedia - Prefix sum

Problemzusammenfassung

Für ganze Zahlen \(n \ge 2\) und \(m \ge 1\) bezeichne \(T(n,m)\) die Anzahl der Folgen \((a_1,\dots,a_m)\) der Länge \(m\) mit

$$1 \le a_i \le n-1 \quad (1 \le i \le m),\qquad a_i+a_{i+1} \le n \quad (1 \le i \le m-1).$$

Gesucht ist \(T(5000,10^{12})\) modulo \(10^9+7\). Ein direktes dynamisches Programm bis zur Länge \(10^{12}\) ist ausgeschlossen. Daher nutzt die Lösung die endliche Übergangsstruktur der Nachbarschaftsbedingung, rekonstruiert daraus eine kurze lineare Rekurrenz und springt dann direkt zum gesuchten Term.

Mathematischer Ansatz

Die zulässigen Folgen bilden einen endlichen Zustandsprozess auf den Werten \(1,2,\dots,n-1\). Sobald dieser Prozess in linearer Form vorliegt, lässt sich der ferne Term über eine deutlich kürzere Rekurrenz bestimmen.

Schritt 1: Die Nachbarschaftsregel als Übergangsmatrix kodieren

Als Zustand dient jeweils der letzte Wert der Folge. Wenn der aktuelle Wert \(x\) und der nächste Wert \(y\) ist, dann ist das Paar genau dann erlaubt, wenn

$$x+y \le n.$$

Damit ergibt sich eine Matrix \(A\in\{0,1\}^{(n-1)\times(n-1)}\) mit

$$A_{y,x}=\begin{cases} 1,& x+y\le n,\\ 0,& x+y\gt n. \end{cases}$$

Jede Zeile besteht aus einem zusammenhängenden Block von Einsen und danach nur noch Nullen. In der ersten Zeile stehen \(n-1\) Einsen, in der nächsten \(n-2\), und so weiter bis hinunter zu einer einzigen Eins. Genau diese Dreiecksform macht die Aufwärmphase effizient.

Schritt 2: Die DP-Rekurrenz aufschreiben

Sei \(f_t(y)\) die Anzahl gültiger Folgen der Länge \(t\), die mit \(y\) enden. Für Länge \(1\) ist jeder Zustand erlaubt, also

$$f_1(y)=1 \qquad (1\le y\le n-1).$$

Für \(t\ge 1\) darf der vorherige Wert jedes \(x\) mit \(x\le n-y\) sein, also

$$f_{t+1}(y)=\sum_{x=1}^{n-y} f_t(x).$$

Definiert man Präfixsummen durch

$$P_t(k)=\sum_{x=1}^{k} f_t(x),$$

dann wird der Übergang zu

$$f_{t+1}(y)=P_t(n-y).$$

Genau das erklärt das im Programm verwendete Verfahren aus Präfixsummen und umgekehrtem Auslesen: Eine komplette Schicht lässt sich in linearer Zeit erzeugen. Die Gesamtzahl gültiger Folgen der Länge \(t\) ist

$$T(n,t)=\sum_{y=1}^{n-1} f_t(y).$$

In Vektorschreibweise, mit \(F_t=(f_t(1),\dots,f_t(n-1))^T\) und \(\mathbf{1}\) als Einsvektor, gilt

$$F_t=A^{t-1}\mathbf{1},\qquad T(n,t)=\mathbf{1}^T A^{t-1}\mathbf{1}.$$

Schritt 3: Warum zwingend eine lineare Rekurrenz existiert

Die Matrix \(A\) hat Größe \(n-1\), also besitzt ihr Minimalpolynom Grad höchstens \(n-1\). Nach dem Cayley-Hamilton-Prinzip erfüllt jede aus Potenzen von \(A\) gewonnene skalare Folge, insbesondere \(T(n,t)\), eine lineare Rekurrenz einer Ordnung \(r\le n-1\):

$$T(n,k+r)=c_1T(n,k+r-1)+c_2T(n,k+r-2)+\cdots+c_rT(n,k)\pmod{10^9+7}.$$

Das ist der strukturelle Kern der Lösung. Statt die DP bis \(m\) fortzusetzen, genügt es, genügend viele Anfangsterme zu berechnen, die Rekurrenz zu bestimmen und anschließend den fernen Term schnell auszuwerten.

Schritt 4: Die kürzeste Rekurrenz aus Anfangstermen rekonstruieren

Die Implementierung erzeugt zunächst ein Anfangsstück

$$T(n,1),T(n,2),T(n,3),\dots$$

modulo \(10^9+7\) und wendet darauf den Berlekamp-Massey-Algorithmus an. Dieser liefert die kürzeste lineare Rekurrenz, die zu allen berechneten Werten passt. Hat die gefundene Rekurrenz Ordnung \(r\), so lautet ihre charakteristische Beziehung

$$x^r-c_1x^{r-1}-c_2x^{r-2}-\cdots-c_r=0.$$

Weil \(r\le n-1\) gilt, reicht jede Stichprobe deutlich oberhalb von \(2(n-1)\); die Implementierungen verwenden \(2(n-1)+8\) Aufwärmterme.

Schritt 5: Den fernen Term per Polynomreduktion berechnen

Ist die Rekurrenz bekannt, wird das Problem zu einer Standardaufgabe der schnellen Berechnung des \(m\)-ten Terms. Man reduziert \(x^{m-1}\) modulo dem Rekurrenzpolynom:

$$x^{m-1}\equiv q_0+q_1x+\cdots+q_{r-1}x^{r-1}\pmod{x^r-c_1x^{r-1}-\cdots-c_r}.$$

Dann gilt

$$T(n,m)\equiv q_0T(n,1)+q_1T(n,2)+\cdots+q_{r-1}T(n,r)\pmod{10^9+7}.$$

Die Koeffizienten \(q_i\) entstehen durch binäre Exponentiation im Quotientenring der Rekurrenz. Deshalb verursacht der riesige Index \(m=10^{12}\) nur einen \(\log m\)-Faktor.

Durchgerechnetes Beispiel: \(n=3\)

Hier sind nur die Werte \(1\) und \(2\) erlaubt. Die Nachbarschaftsregel verbietet nur das Paar \((2,2)\), also ist

$$A=\begin{pmatrix} 1 & 1\\ 1 & 0 \end{pmatrix}.$$

Die Gesamtzahlen nach Länge sind

$$T(3,1)=2,\qquad T(3,2)=3,\qquad T(3,3)=5,\qquad T(3,4)=8.$$

In diesem Fall ist die Rekurrenz einfach

$$T(3,t)=T(3,t-1)+T(3,t-2),$$

also folgt der Kontrollwert \(T(3,4)=8\) sofort.

Wie der Code arbeitet

Die C++-, Python- und Java-Implementierungen folgen alle derselben mathematischen Pipeline. Zuerst erzeugen sie einen Anfangsblock der Totalsummen mit dem endzustandsbasierten dynamischen Programm. Statt jeden Zustandsübergang naiv neu auszurechnen, bilden sie Präfixsummen der aktuellen Schicht und lesen diese rückwärts aus. Dadurch kostet ein Übergangsschritt nur \(O(n)\) Zeit.

Anschließend wird Berlekamp-Massey modulo \(10^9+7\) auf diese Aufwärmwerte angewendet, um die kürzeste Rekurrenz für \(T(n,m)\) zu erhalten. Liegt die gesuchte Länge bereits im vorberechneten Präfix, wird der Wert direkt zurückgegeben. Andernfalls führt die Implementierung eine binäre Exponentiation auf reduzierten Polynomen aus, gewinnt so die Koeffizienten \(q_0,\dots,q_{r-1}\) und damit den gesuchten fernen Term.

Komplexitätsanalyse

Jeder Aufwärmschritt der DP kostet \(O(n)\) Zeit und \(O(n)\) Speicher. Da nur \(O(n)\) Anfangsterme erzeugt werden, benötigt diese Phase insgesamt \(O(n^2)\) Zeit. Berlekamp-Massey für eine Rekurrenz der Ordnung \(r\le n-1\) kostet \(O(r^2)\), und die abschließende Berechnung des \(m\)-ten Terms kostet \(O(r^2\log m)\). Insgesamt ergibt sich also

$$O(n^2+r^2+r^2\log m)=O(n^2\log m),$$

bei \(O(n)\) Speicherverbrauch.

Fußnoten und Referenzen

  1. Project-Euler-Problemseite: Problem 654: Neighbourly Constraints
  2. Berlekamp-Massey-Algorithmus: Wikipedia - Berlekamp-Massey algorithm
  3. Schnelle Berechnung linearer Rekurrenzen: cp-algorithms - Linear Recurrence
  4. Cayley-Hamilton-Theorem: Wikipedia - Cayley-Hamilton theorem
  5. Präfixsummen: Wikipedia - Prefix sum

Problem Özeti

\(n \ge 2\) ve \(m \ge 1\) için \(T(n,m)\),

$$1 \le a_i \le n-1 \quad (1 \le i \le m),\qquad a_i+a_{i+1} \le n \quad (1 \le i \le m-1)$$

koşullarını sağlayan uzunluğu \(m\) olan \((a_1,\dots,a_m)\) dizilerinin sayısı olsun. İstenen değer \(T(5000,10^{12}) \bmod (10^9+7)\) olduğundan, uzunluk boyunca sıradan DP yapmak mümkün değildir. Çözüm bunun yerine komşuluk kuralının oluşturduğu sonlu geçiş sistemini inceler, buradan kısa bir lineer rekürans çıkarır ve doğrudan hedef terime sıçrar.

Matematiksel Yaklaşım

Geçerli diziler, \(1,2,\dots,n-1\) değerleri üzerinde çalışan sonlu durumlu bir süreç oluşturur. Bu süreç doğrusal cebir diliyle yazılınca çok uzaktaki terim kısa bir rekürans üzerinden hesaplanabilir.

Adım 1: Komşuluk kuralını bir geçiş matrisi olarak yaz

Durum olarak dizinin son elemanını kullanalım. Mevcut değer \(x\), bir sonraki değer \(y\) ise geçerli olma koşulu tam olarak

$$x+y \le n$$

eşitsizliğidir. Böylece

$$A_{y,x}=\begin{cases} 1,& x+y\le n,\\ 0,& x+y\gt n \end{cases}$$

şeklinde tanımlanan \(A\in\{0,1\}^{(n-1)\times(n-1)}\) matrisi ortaya çıkar. Her satır solda art arda gelen birler ve onların ardından sıfırlardan oluşur. İlk satırda \(n-1\), sonraki satırda \(n-2\), en altta ise yalnızca \(1\) tane bir vardır. Uygulamadaki verimli ön hesaplama tam olarak bu üçgensel yapıdan yararlanır.

Adım 2: Dinamik programlama bağıntısını kur

\(f_t(y)\), uzunluğu \(t\) olan ve \(y\) ile biten geçerli dizi sayısı olsun. Uzunluk \(1\) için her durum mümkündür:

$$f_1(y)=1 \qquad (1\le y\le n-1).$$

\(t\ge 1\) iken önceki değer, \(x\le n-y\) koşulunu sağlayan herhangi bir sayı olabilir. Dolayısıyla

$$f_{t+1}(y)=\sum_{x=1}^{n-y} f_t(x).$$

Eğer

$$P_t(k)=\sum_{x=1}^{k} f_t(x)$$

önek toplamını tanımlarsak geçiş

$$f_{t+1}(y)=P_t(n-y)$$

haline gelir. Kodun kullandığı ters yönde okunan önek toplam güncellemesi tam olarak budur; böylece bir katman \(O(n)\) sürede üretilir. Uzunluğu \(t\) olan tüm geçerli dizilerin toplam sayısı ise

$$T(n,t)=\sum_{y=1}^{n-1} f_t(y)$$

olur. Vektör biçiminde, \(F_t=(f_t(1),\dots,f_t(n-1))^T\) ve \(\mathbf{1}\) birler vektörü olmak üzere

$$F_t=A^{t-1}\mathbf{1},\qquad T(n,t)=\mathbf{1}^T A^{t-1}\mathbf{1}.$$

Adım 3: Neden mutlaka lineer rekürans vardır

\(A\) matrisi \(n-1\) boyutludur; dolayısıyla minimal polinomunun derecesi en fazla \(n-1\)'dir. Cayley-Hamilton bakış açısından, \(A\)'nın kuvvetlerinden elde edilen her skaler dizi, özellikle de \(T(n,t)\), mertebesi \(r\le n-1\) olan bir lineer rekürans sağlar:

$$T(n,k+r)=c_1T(n,k+r-1)+c_2T(n,k+r-2)+\cdots+c_rT(n,k)\pmod{10^9+7}.$$

Asıl fikir budur. Yani DP'yi \(m\)'e kadar uzatmak yerine yeterli sayıda başlangıç terimi üretip reküransı bulmak ve ardından uzak terimi hızlı hesaplamak yeterlidir.

Adım 4: Başlangıç terimlerinden en kısa reküransı çıkar

Uygulama önce

$$T(n,1),T(n,2),T(n,3),\dots$$

değerlerinin bir başlangıç bölümünü \(10^9+7\) modunda üretir. Ardından Berlekamp-Massey algoritması bu değerlerle uyumlu en kısa lineer reküransı bulur. Bulunan mertebe \(r\) ise karakteristik ilişki

$$x^r-c_1x^{r-1}-c_2x^{r-2}-\cdots-c_r=0$$

şeklindedir. \(r\le n-1\) olduğundan \(2(n-1)\)'in biraz üzerindeki bir örnekleme yeterlidir; uygulamalar \(2(n-1)+8\) başlangıç terimi kullanır.

Adım 5: Uzak terimi polinom indirgeme ile hesapla

Rekürans elde edildikten sonra problem hızlı \(m\). terim hesabına dönüşür. \(x^{m-1}\) kuvvetini rekürans polinomu modunda indirgeriz:

$$x^{m-1}\equiv q_0+q_1x+\cdots+q_{r-1}x^{r-1}\pmod{x^r-c_1x^{r-1}-\cdots-c_r}.$$

Böylece

$$T(n,m)\equiv q_0T(n,1)+q_1T(n,2)+\cdots+q_{r-1}T(n,r)\pmod{10^9+7}.$$

\(q_i\) katsayıları reküransın tanımladığı bölüm halkasında ikili üs alma ile bulunduğundan, devasa \(m=10^{12}\) değeri sadece \(\log m\) çarpanı kadar maliyet getirir.

Çalışılmış Örnek: \(n=3\)

Bu durumda izin verilen değerler yalnızca \(1\) ve \(2\)'dir. Komşuluk kuralı sadece \((2,2)\) çiftini yasaklar. O halde

$$A=\begin{pmatrix} 1 & 1\\ 1 & 0 \end{pmatrix}.$$

Uzunluğa göre toplam sayılar

$$T(3,1)=2,\qquad T(3,2)=3,\qquad T(3,3)=5,\qquad T(3,4)=8$$

olur. Bu özel durumda rekürans

$$T(3,t)=T(3,t-1)+T(3,t-2)$$

şeklindedir; dolayısıyla kontrol değeri \(T(3,4)=8\) hemen elde edilir.

Kod Nasıl Çalışır

C++, Python ve Java uygulamalarının hepsi aynı matematiksel hattı izler. Önce son değere göre tutulan DP ile başlangıç toplamları üretilir. Her geçişi naif biçimde yeniden toplamak yerine mevcut katmanın önek toplamları alınır ve bunlar ters sırada okunur; bu sayede tek bir geçiş adımı \(O(n)\) zamanda tamamlanır.

Daha sonra Berlekamp-Massey algoritması bu başlangıç terimlerine \(10^9+7\) modunda uygulanır ve \(T(n,m)\) için en kısa rekürans bulunur. İstenen uzunluk önceden hesaplanan aralığın içindeyse sonuç doğrudan döndürülür. Aksi halde uygulama indirgenmiş polinomlar üzerinde ikili üs alma yapar, \(q_0,\dots,q_{r-1}\) katsayılarını elde eder ve hedef terimi hesaplar.

Karmaşıklık Analizi

Her DP ön hesaplama adımı \(O(n)\) zaman ve \(O(n)\) bellek kullanır. Yalnızca \(O(n)\) tane başlangıç terimi üretildiği için bu faz toplamda \(O(n^2)\) zamana mal olur. Mertebesi \(r\le n-1\) olan bir rekürans için Berlekamp-Massey maliyeti \(O(r^2)\), son \(m\). terimi hızlı hesaplama maliyeti ise \(O(r^2\log m)\)'dir. Toplam karmaşıklık böylece

$$O(n^2+r^2+r^2\log m)=O(n^2\log m)$$

ve bellek kullanımı \(O(n)\) olur.

Dipnotlar ve Kaynaklar

  1. Project Euler problem sayfası: Problem 654: Neighbourly Constraints
  2. Berlekamp-Massey algoritması: Wikipedia - Berlekamp-Massey algorithm
  3. Lineer reküransların hızlı hesaplanması: cp-algorithms - Linear Recurrence
  4. Cayley-Hamilton teoremi: Wikipedia - Cayley-Hamilton theorem
  5. Önek toplamlar: Wikipedia - Prefix sum

Resumen del Problema

Para enteros \(n \ge 2\) y \(m \ge 1\), sea \(T(n,m)\) el número de secuencias \((a_1,\dots,a_m)\) de longitud \(m\) tales que

$$1 \le a_i \le n-1 \quad (1 \le i \le m),\qquad a_i+a_{i+1} \le n \quad (1 \le i \le m-1).$$

La meta es calcular \(T(5000,10^{12})\) módulo \(10^9+7\). Un programa dinámico directo hasta longitud \(10^{12}\) es imposible, así que la solución estudia primero el sistema finito de transiciones inducido por la restricción entre vecinos, reconstruye una recurrencia lineal corta y luego salta directamente al término pedido.

Enfoque Matemático

Las secuencias válidas forman un proceso de estados finitos sobre los valores \(1,2,\dots,n-1\). Una vez escrito ese proceso en forma lineal, el término lejano puede recuperarse a partir de una recurrencia mucho más corta.

Paso 1: Codificar la regla de vecindad como una matriz de transición

Tomamos como estado el último valor de la secuencia. Si el valor actual es \(x\) y el siguiente es \(y\), el par es válido exactamente cuando

$$x+y \le n.$$

Esto define una matriz \(A\in\{0,1\}^{(n-1)\times(n-1)}\) dada por

$$A_{y,x}=\begin{cases} 1,& x+y\le n,\\ 0,& x+y\gt n. \end{cases}$$

Cada fila es un bloque consecutivo de unos seguido por ceros. La primera fila tiene \(n-1\) unos, la siguiente \(n-2\), y así sucesivamente hasta llegar a una sola unidad. Esa forma triangular es la razón de que la fase inicial sea eficiente.

Paso 2: Escribir la recurrencia de programación dinámica

Sea \(f_t(y)\) el número de secuencias válidas de longitud \(t\) que terminan en \(y\). Para longitud \(1\), cualquier estado sirve, así que

$$f_1(y)=1 \qquad (1\le y\le n-1).$$

Para \(t\ge 1\), el valor anterior puede ser cualquier \(x\) con \(x\le n-y\), por lo que

$$f_{t+1}(y)=\sum_{x=1}^{n-y} f_t(x).$$

Si definimos sumas prefijo por

$$P_t(k)=\sum_{x=1}^{k} f_t(x),$$

entonces la transición queda

$$f_{t+1}(y)=P_t(n-y).$$

Eso explica la actualización con prefijos leídos en orden inverso que usa la implementación: una capa completa se genera en tiempo lineal. El número total de secuencias válidas de longitud \(t\) es

$$T(n,t)=\sum_{y=1}^{n-1} f_t(y).$$

En forma vectorial, con \(F_t=(f_t(1),\dots,f_t(n-1))^T\) y \(\mathbf{1}\) el vector de unos,

$$F_t=A^{t-1}\mathbf{1},\qquad T(n,t)=\mathbf{1}^T A^{t-1}\mathbf{1}.$$

Paso 3: Por qué debe existir una recurrencia lineal

La matriz \(A\) tiene tamaño \(n-1\), así que su polinomio minimal tiene grado a lo sumo \(n-1\). Por la idea de Cayley-Hamilton, toda sucesión escalar extraída de potencias de \(A\), incluida \(T(n,t)\), satisface una recurrencia lineal de algún orden \(r\le n-1\):

$$T(n,k+r)=c_1T(n,k+r-1)+c_2T(n,k+r-2)+\cdots+c_rT(n,k)\pmod{10^9+7}.$$

Ese es el hecho estructural central. En vez de extender la DP hasta \(m\), basta con producir suficientes términos iniciales, recuperar la recurrencia y evaluar después el término lejano.

Paso 4: Recuperar la recurrencia mínima a partir de términos iniciales

La implementación genera primero un prefijo

$$T(n,1),T(n,2),T(n,3),\dots$$

módulo \(10^9+7\), y luego aplica el algoritmo de Berlekamp-Massey para hallar la recurrencia lineal más corta compatible con esos valores. Si el orden recuperado es \(r\), la relación característica puede escribirse como

$$x^r-c_1x^{r-1}-c_2x^{r-2}-\cdots-c_r=0.$$

Como \(r\le n-1\), cualquier muestra claramente superior a \(2(n-1)\) es suficiente; las implementaciones usan \(2(n-1)+8\) términos de arranque.

Paso 5: Calcular el término lejano por reducción de polinomios

Una vez conocida la recurrencia, el problema pasa a ser el cálculo rápido del término \(m\)-ésimo. Reducimos \(x^{m-1}\) módulo el polinomio de la recurrencia:

$$x^{m-1}\equiv q_0+q_1x+\cdots+q_{r-1}x^{r-1}\pmod{x^r-c_1x^{r-1}-\cdots-c_r}.$$

Entonces

$$T(n,m)\equiv q_0T(n,1)+q_1T(n,2)+\cdots+q_{r-1}T(n,r)\pmod{10^9+7}.$$

Los coeficientes \(q_i\) se obtienen por exponenciación binaria dentro del anillo cociente definido por la recurrencia. Por eso el enorme índice \(m=10^{12}\) solo aporta un factor \(\log m\).

Ejemplo Desarrollado: \(n=3\)

Aquí los valores permitidos son \(1\) y \(2\). La regla de vecindad prohíbe únicamente el par \((2,2)\), así que

$$A=\begin{pmatrix} 1 & 1\\ 1 & 0 \end{pmatrix}.$$

Los totales por longitud son

$$T(3,1)=2,\qquad T(3,2)=3,\qquad T(3,3)=5,\qquad T(3,4)=8.$$

En este caso la recurrencia es simplemente

$$T(3,t)=T(3,t-1)+T(3,t-2),$$

de modo que el punto de control \(T(3,4)=8\) se obtiene de inmediato.

Cómo Funciona el Código

Las implementaciones en C++, Python y Java siguen la misma tubería matemática. Primero generan un bloque inicial de totales con la DP basada en el último estado. En lugar de recomputar ingenuamente cada transición, construyen sumas prefijo de la capa actual y las leen en orden inverso, con lo que cada paso de transición cuesta \(O(n)\).

Después ejecutan Berlekamp-Massey módulo \(10^9+7\) sobre esos totales iniciales y obtienen la recurrencia más corta para \(T(n,m)\). Si la longitud pedida ya está dentro del prefijo precalculado, la respuesta se devuelve directamente. En caso contrario, la implementación hace exponenciación binaria sobre polinomios reducidos, obtiene los coeficientes \(q_0,\dots,q_{r-1}\) y calcula el término lejano requerido.

Análisis de Complejidad

Cada paso de la DP de arranque cuesta \(O(n)\) tiempo y \(O(n)\) memoria. Como solo se generan \(O(n)\) términos iniciales, esa fase cuesta \(O(n^2)\) en total. Berlekamp-Massey sobre una recurrencia de orden \(r\le n-1\) cuesta \(O(r^2)\), y la etapa final para el término \(m\)-ésimo cuesta \(O(r^2\log m)\). Por lo tanto, la complejidad total es

$$O(n^2+r^2+r^2\log m)=O(n^2\log m),$$

con memoria \(O(n)\).

Notas y Referencias

  1. Página del problema en Project Euler: Problem 654: Neighbourly Constraints
  2. Algoritmo de Berlekamp-Massey: Wikipedia - Berlekamp-Massey algorithm
  3. Cálculo rápido de recurrencias lineales: cp-algorithms - Linear Recurrence
  4. Teorema de Cayley-Hamilton: Wikipedia - Cayley-Hamilton theorem
  5. Sumas prefijo: Wikipedia - Prefix sum

问题概述

对整数 \(n \ge 2\) 与 \(m \ge 1\),记 \(T(n,m)\) 为满足下列条件的长度为 \(m\) 的序列 \((a_1,\dots,a_m)\) 的个数:

$$1 \le a_i \le n-1 \quad (1 \le i \le m),\qquad a_i+a_{i+1} \le n \quad (1 \le i \le m-1).$$

题目要求计算 \(T(5000,10^{12})\bmod (10^9+7)\)。如果沿着长度从 \(1\) 一直动态规划到 \(10^{12}\),计算量显然无法接受。因此解法先把相邻约束写成一个有限状态转移系统,再从这个系统中恢复出一个很短的线性递推,最后直接跳到目标项。

数学方法

所有合法序列都对应于值域 \(1,2,\dots,n-1\) 上的一个有限状态过程。只要把这个过程写成线性代数形式,远处的项就可以通过更短的递推关系恢复出来,而不必逐项推进到 \(m\)。

步骤 1:把相邻限制写成转移矩阵

把序列的最后一个值作为状态。若当前值为 \(x\),下一项为 \(y\),那么这一对合法当且仅当

$$x+y \le n.$$

于是得到一个 \(A\in\{0,1\}^{(n-1)\times(n-1)}\) 矩阵:

$$A_{y,x}=\begin{cases} 1,& x+y\le n,\\ 0,& x+y\gt n. \end{cases}$$

每一行都是一段连续的 \(1\),后面跟着若干个 \(0\)。第一行有 \(n-1\) 个 \(1\),第二行有 \(n-2\) 个,一直到最后一行只剩下 \(1\) 个。这种三角形结构非常重要,因为它让初始动态规划可以用前缀和在线性时间内完成一层转移。

步骤 2:写出动态规划递推

设 \(f_t(y)\) 表示长度为 \(t\)、并且最后一项等于 \(y\) 的合法序列数。长度为 \(1\) 时,每个状态都可取,因此

$$f_1(y)=1 \qquad (1\le y\le n-1).$$

当 \(t\ge 1\) 时,前一项可以是任意满足 \(x\le n-y\) 的值,所以

$$f_{t+1}(y)=\sum_{x=1}^{n-y} f_t(x).$$

如果定义前缀和

$$P_t(k)=\sum_{x=1}^{k} f_t(x),$$

那么转移就能改写成

$$f_{t+1}(y)=P_t(n-y).$$

这正是实现中“先做前缀和,再反向读取”的数学原因。于是一次完整的状态更新只需 \(O(n)\) 时间。长度为 \(t\) 的合法序列总数是

$$T(n,t)=\sum_{y=1}^{n-1} f_t(y).$$

用向量记号表示,若 \(F_t=(f_t(1),\dots,f_t(n-1))^T\),\(\mathbf{1}\) 表示全 \(1\) 向量,则

$$F_t=A^{t-1}\mathbf{1},\qquad T(n,t)=\mathbf{1}^T A^{t-1}\mathbf{1}.$$

步骤 3:为什么一定存在线性递推

矩阵 \(A\) 的维数是 \(n-1\),所以它的极小多项式次数至多为 \(n-1\)。由 Cayley-Hamilton 思想可知,从 \(A\) 的幂中提取出的任意标量序列都会满足某个线性递推,\(T(n,t)\) 当然也不例外。也就是说,存在某个 \(r\le n-1\),使得

$$T(n,k+r)=c_1T(n,k+r-1)+c_2T(n,k+r-2)+\cdots+c_rT(n,k)\pmod{10^9+7}.$$

这一步给出了解法的核心结构:真正需要的不是把 DP 做到第 \(m\) 层,而是先算出足够多的起始项,再把隐藏的递推关系恢复出来。

步骤 4:从初始项中恢复最短递推

实现会先生成一段前缀

$$T(n,1),T(n,2),T(n,3),\dots$$

并全部在模 \(10^9+7\) 下处理,然后对这串值应用 Berlekamp-Massey 算法,求出与它们完全吻合的最短线性递推。若恢复出的阶数为 \(r\),则相应的特征关系可以写成

$$x^r-c_1x^{r-1}-c_2x^{r-2}-\cdots-c_r=0.$$

由于 \(r\le n-1\),所以只要初始样本数明显超过 \(2(n-1)\) 就足够稳定恢复;实现中使用的是 \(2(n-1)+8\) 个预热项。

步骤 5:通过多项式约化直接求远处的第 \(m\) 项

一旦递推已经知道,问题就变成标准的“快速求线性递推第 \(m\) 项”。把 \(x^{m-1}\) 对递推多项式取模:

$$x^{m-1}\equiv q_0+q_1x+\cdots+q_{r-1}x^{r-1}\pmod{x^r-c_1x^{r-1}-\cdots-c_r}.$$

于是目标值就是

$$T(n,m)\equiv q_0T(n,1)+q_1T(n,2)+\cdots+q_{r-1}T(n,r)\pmod{10^9+7}.$$

这里的 \(q_i\) 系数可以通过递推多项式所定义的商环中的二进制快速幂求出,因此即使 \(m=10^{12}\) 很大,额外代价也只是一个 \(\log m\) 因子。

例子:\(n=3\)

此时允许的值只有 \(1\) 和 \(2\)。相邻约束只禁止 \((2,2)\) 这一对,因此

$$A=\begin{pmatrix} 1 & 1\\ 1 & 0 \end{pmatrix}.$$

按长度统计的总数为

$$T(3,1)=2,\qquad T(3,2)=3,\qquad T(3,3)=5,\qquad T(3,4)=8.$$

在这个小例子里,递推恰好就是

$$T(3,t)=T(3,t-1)+T(3,t-2),$$

所以检查点 \(T(3,4)=8\) 会立刻得到。

代码如何工作

C++、Python 和 Java 实现都遵循同一条数学路线。第一阶段先用“按最后一个值分类”的动态规划生成一批起始总数。为了避免朴素地为每个状态重新求和,它们先对当前层做前缀和,再把这些前缀和按反向顺序读出,于是一次状态转移只需要 \(O(n)\) 时间。

第二阶段对这些预热项在模 \(10^9+7\) 下运行 Berlekamp-Massey,得到 \(T(n,m)\) 的最短线性递推。如果所求长度已经落在预计算前缀之内,就直接返回对应值。否则实现会在约化多项式空间中做二进制快速幂,求出 \(q_0,\dots,q_{r-1}\) 这些系数,并据此合成最终答案。

复杂度分析

每一步预热动态规划需要 \(O(n)\) 时间和 \(O(n)\) 空间。由于只会生成 \(O(n)\) 个初始项,预热阶段总时间是 \(O(n^2)\)。对阶数 \(r\le n-1\) 的递推做 Berlekamp-Massey 需要 \(O(r^2)\),最后的快速求第 \(m\) 项需要 \(O(r^2\log m)\)。因此总复杂度为

$$O(n^2+r^2+r^2\log m)=O(n^2\log m),$$

空间复杂度为 \(O(n)\)。

注释与参考资料

  1. Project Euler 题目页面:Problem 654: Neighbourly Constraints
  2. Berlekamp-Massey 算法:Wikipedia - Berlekamp-Massey algorithm
  3. 线性递推的快速计算:cp-algorithms - Linear Recurrence
  4. Cayley-Hamilton 定理:Wikipedia - Cayley-Hamilton theorem
  5. 前缀和:Wikipedia - Prefix sum

Краткое описание задачи

Для целых \(n \ge 2\) и \(m \ge 1\) обозначим через \(T(n,m)\) число последовательностей \((a_1,\dots,a_m)\) длины \(m\), удовлетворяющих условиям

$$1 \le a_i \le n-1 \quad (1 \le i \le m),\qquad a_i+a_{i+1} \le n \quad (1 \le i \le m-1).$$

Нужно вычислить \(T(5000,10^{12})\) по модулю \(10^9+7\). Прямое динамическое программирование по длине до \(10^{12}\) невозможно, поэтому решение сначала анализирует конечную систему переходов, заданную ограничением на соседние элементы, затем восстанавливает короткую линейную рекуррентность и уже по ней быстро получает нужный член.

Математический подход

Все допустимые последовательности образуют конечный процесс состояний на значениях \(1,2,\dots,n-1\). Как только этот процесс записан в линейной форме, далёкий член можно выразить через намного более короткую рекуррентность.

Шаг 1: Записать соседнее ограничение как матрицу переходов

Состоянием считаем последний элемент последовательности. Если текущий элемент равен \(x\), а следующий равен \(y\), то пара допустима тогда и только тогда, когда

$$x+y \le n.$$

Отсюда возникает матрица \(A\in\{0,1\}^{(n-1)\times(n-1)}\):

$$A_{y,x}=\begin{cases} 1,& x+y\le n,\\ 0,& x+y\gt n. \end{cases}$$

Каждая строка состоит из непрерывного блока единиц, за которым следуют нули. В первой строке стоит \(n-1\) единиц, во второй \(n-2\), и так далее, пока в последней не останется одна единица. Именно эта треугольная форма позволяет эффективно выполнять стартовые переходы.

Шаг 2: Выписать рекуррентность динамического программирования

Пусть \(f_t(y)\) обозначает число допустимых последовательностей длины \(t\), оканчивающихся значением \(y\). Для длины \(1\) каждый допустимый конечный элемент возможен, поэтому

$$f_1(y)=1 \qquad (1\le y\le n-1).$$

При \(t\ge 1\) предыдущим значением может быть любой \(x\), удовлетворяющий \(x\le n-y\). Следовательно,

$$f_{t+1}(y)=\sum_{x=1}^{n-y} f_t(x).$$

Если ввести префиксные суммы

$$P_t(k)=\sum_{x=1}^{k} f_t(x),$$

то переход переписывается как

$$f_{t+1}(y)=P_t(n-y).$$

Именно поэтому реализация использует префиксные суммы с обратным чтением: один полный шаг получается за линейное время. Общее число допустимых последовательностей длины \(t\) равно

$$T(n,t)=\sum_{y=1}^{n-1} f_t(y).$$

В векторной форме, если \(F_t=(f_t(1),\dots,f_t(n-1))^T\), а \(\mathbf{1}\) есть вектор из единиц, то

$$F_t=A^{t-1}\mathbf{1},\qquad T(n,t)=\mathbf{1}^T A^{t-1}\mathbf{1}.$$

Шаг 3: Почему линейная рекуррентность обязана существовать

Матрица \(A\) имеет размер \(n-1\), значит её минимальный многочлен имеет степень не выше \(n-1\). По идее теоремы Кэли-Гамильтона любая скалярная последовательность, извлечённая из степеней \(A\), и в частности \(T(n,t)\), обязана удовлетворять линейной рекуррентности некоторого порядка \(r\le n-1\):

$$T(n,k+r)=c_1T(n,k+r-1)+c_2T(n,k+r-2)+\cdots+c_rT(n,k)\pmod{10^9+7}.$$

Это главный структурный факт. Вместо того чтобы продолжать DP до самого \(m\), достаточно вычислить достаточно длинный начальный префикс, восстановить скрытую рекуррентность и затем быстро получить далёкий член.

Шаг 4: Восстановить кратчайшую рекуррентность по начальным членам

Сначала реализация строит префикс

$$T(n,1),T(n,2),T(n,3),\dots$$

по модулю \(10^9+7\), после чего применяет алгоритм Берлекэмпа-Мэсси и находит кратчайшую линейную рекуррентность, согласованную с этими значениями. Если найденный порядок равен \(r\), то характеристическое соотношение можно записать как

$$x^r-c_1x^{r-1}-c_2x^{r-2}-\cdots-c_r=0.$$

Так как \(r\le n-1\), достаточно запаса чуть больше чем \(2(n-1)\) членов; в реализациях берётся \(2(n-1)+8\) стартовых значений.

Шаг 5: Вычислить далёкий член через редукцию многочлена

Когда рекуррентность уже известна, задача превращается в стандартное быстрое вычисление \(m\)-го члена линейной рекуррентности. Нужно сократить \(x^{m-1}\) по модулю характеристического многочлена:

$$x^{m-1}\equiv q_0+q_1x+\cdots+q_{r-1}x^{r-1}\pmod{x^r-c_1x^{r-1}-\cdots-c_r}.$$

Тогда

$$T(n,m)\equiv q_0T(n,1)+q_1T(n,2)+\cdots+q_{r-1}T(n,r)\pmod{10^9+7}.$$

Коэффициенты \(q_i\) получаются бинарным возведением в степень в факторкольце, заданном рекуррентностью. Поэтому огромный индекс \(m=10^{12}\) даёт лишь множитель \(\log m\) во времени работы.

Разобранный пример: \(n=3\)

Здесь разрешены только значения \(1\) и \(2\). Ограничение на соседние элементы запрещает лишь пару \((2,2)\), значит

$$A=\begin{pmatrix} 1 & 1\\ 1 & 0 \end{pmatrix}.$$

Общие количества по длинам равны

$$T(3,1)=2,\qquad T(3,2)=3,\qquad T(3,3)=5,\qquad T(3,4)=8.$$

В этом случае рекуррентность просто совпадает с формулой Фибоначчи:

$$T(3,t)=T(3,t-1)+T(3,t-2),$$

поэтому контрольное значение \(T(3,4)=8\) получается сразу.

Как работает код

Реализации на C++, Python и Java следуют одной и той же математической схеме. Сначала они строят стартовый блок суммарных значений с помощью DP по последнему элементу. Вместо наивного пересчёта каждого перехода они вычисляют префиксные суммы текущего слоя и читают их в обратном порядке, так что один переходный шаг занимает \(O(n)\) времени.

Затем к этим стартовым значениям по модулю \(10^9+7\) применяется алгоритм Берлекэмпа-Мэсси, и получается кратчайшая рекуррентность для \(T(n,m)\). Если требуемая длина уже лежит внутри предвычисленного префикса, ответ возвращается сразу. Иначе реализация выполняет бинарное возведение в степень на редуцированных многочленах, находит коэффициенты \(q_0,\dots,q_{r-1}\) и собирает итоговый далёкий член.

Анализ сложности

Каждый шаг стартового динамического программирования требует \(O(n)\) времени и \(O(n)\) памяти. Поскольку генерируется только \(O(n)\) начальных членов, эта фаза стоит \(O(n^2)\) времени. Алгоритм Берлекэмпа-Мэсси для рекуррентности порядка \(r\le n-1\) работает за \(O(r^2)\), а финальное быстрое вычисление \(m\)-го члена занимает \(O(r^2\log m)\). Следовательно, общая сложность равна

$$O(n^2+r^2+r^2\log m)=O(n^2\log m),$$

при использовании \(O(n)\) памяти.

Примечания и ссылки

  1. Страница задачи Project Euler: Problem 654: Neighbourly Constraints
  2. Алгоритм Берлекэмпа-Мэсси: Wikipedia - Berlekamp-Massey algorithm
  3. Быстрое вычисление линейных рекуррентностей: cp-algorithms - Linear Recurrence
  4. Теорема Кэли-Гамильтона: Wikipedia - Cayley-Hamilton theorem
  5. Префиксные суммы: Wikipedia - Prefix sum

ملخص المشكلة

لعددين صحيحين \(n \ge 2\) و\(m \ge 1\)، لنعرّف \(T(n,m)\) بأنه عدد المتتاليات \((a_1,\dots,a_m)\) ذات الطول \(m\) التي تحقق

$$1 \le a_i \le n-1 \quad (1 \le i \le m),\qquad a_i+a_{i+1} \le n \quad (1 \le i \le m-1).$$

المطلوب هو حساب \(T(5000,10^{12})\) بترديد \(10^9+7\). من المستحيل تنفيذ برمجة ديناميكية مباشرة حتى طول \(10^{12}\)، لذلك يعتمد الحل على فهم نظام الانتقالات المحدود الذي تفرضه علاقة الجوار، ثم استعادة علاقة عودية خطية قصيرة، ثم القفز مباشرة إلى الحد المطلوب.

المنهج الرياضي

كل المتتاليات المسموح بها تشكل عملية ذات حالات منتهية على القيم \(1,2,\dots,n-1\). وبمجرد كتابة هذه العملية بصيغة جبر خطي، يمكن استخراج الحد البعيد من علاقة عودية أقصر بكثير بدل متابعة الحساب خطوة خطوة حتى \(m\).

الخطوة 1: تمثيل شرط الجوار بمصفوفة انتقال

نأخذ آخر عنصر في المتتالية على أنه الحالة الحالية. إذا كانت القيمة الحالية \(x\) والقيمة التالية \(y\)، فإن الزوج يكون صالحًا بالضبط عندما

$$x+y \le n.$$

وهذا يعطينا مصفوفة \(A\in\{0,1\}^{(n-1)\times(n-1)}\) معرفة كما يلي:

$$A_{y,x}=\begin{cases} 1,& x+y\le n,\\ 0,& x+y\gt n. \end{cases}$$

كل صف فيها عبارة عن كتلة متصلة من الواحدات ثم أصفار بعدها. الصف الأول يحتوي \(n-1\) واحدًا، ثم \(n-2\) في الصف التالي، وهكذا حتى نصل إلى صف فيه واحد واحد فقط. هذا الشكل المثلثي هو السبب في أن مرحلة التمهيد الحسابية تصبح فعالة.

الخطوة 2: كتابة علاقة البرمجة الديناميكية

لتكن \(f_t(y)\) عدد المتتاليات الصحيحة ذات الطول \(t\) والمنتهية بالقيمة \(y\). عندما يكون الطول \(1\)، فجميع الحالات ممكنة، ولذلك

$$f_1(y)=1 \qquad (1\le y\le n-1).$$

ولـ \(t\ge 1\)، يمكن أن تكون القيمة السابقة أي \(x\) يحقق \(x\le n-y\)، ولذلك نحصل على

$$f_{t+1}(y)=\sum_{x=1}^{n-y} f_t(x).$$

إذا عرّفنا المجاميع البادئة

$$P_t(k)=\sum_{x=1}^{k} f_t(x),$$

فإن الانتقال يصبح

$$f_{t+1}(y)=P_t(n-y).$$

وهذا يفسر مباشرة لماذا تستخدم التنفيذات مجاميع بادئة مع قراءة معكوسة: فطبقة انتقال كاملة تُبنى في زمن \(O(n)\). أما العدد الكلي للمتتاليات الصحيحة ذات الطول \(t\) فهو

$$T(n,t)=\sum_{y=1}^{n-1} f_t(y).$$

وبصيغة المتجهات، إذا كان \(F_t=(f_t(1),\dots,f_t(n-1))^T\) و\(\mathbf{1}\) هو متجه الواحدات، فإن

$$F_t=A^{t-1}\mathbf{1},\qquad T(n,t)=\mathbf{1}^T A^{t-1}\mathbf{1}.$$

الخطوة 3: لماذا لا بد أن توجد علاقة عودية خطية

حجم المصفوفة \(A\) هو \(n-1\)، ولذلك فدرجة كثير حدودها الأدنى لا تتجاوز \(n-1\). وبفكرة كايلي-هاملتون، فإن أي متتالية عددية مستخرجة من قوى \(A\)، ومنها \(T(n,t)\)، يجب أن تحقق علاقة عودية خطية من رتبة ما \(r\le n-1\):

$$T(n,k+r)=c_1T(n,k+r-1)+c_2T(n,k+r-2)+\cdots+c_rT(n,k)\pmod{10^9+7}.$$

هذه هي الفكرة البنيوية الأساسية. بدلًا من تمديد البرمجة الديناميكية حتى \(m\)، يكفي إنتاج عدد مناسب من الحدود الأولى، ثم استخراج العلاقة العودية المخفية، ثم استخدام تلك العلاقة للوصول سريعًا إلى الحد البعيد.

الخطوة 4: استعادة أقصر علاقة عودية من الحدود الابتدائية

تقوم التنفيذات أولًا بتوليد مقطع ابتدائي من

$$T(n,1),T(n,2),T(n,3),\dots$$

بترديد \(10^9+7\)، ثم تطبق خوارزمية Berlekamp-Massey لاستخراج أقصر علاقة عودية خطية توافق هذه القيم. إذا كانت الرتبة المستعادة هي \(r\)، فيمكن كتابة العلاقة المميزة على الصورة

$$x^r-c_1x^{r-1}-c_2x^{r-2}-\cdots-c_r=0.$$

وبما أن \(r\le n-1\)، فإن أي عدد من الحدود يزيد بوضوح على \(2(n-1)\) يكفي للاستعادة؛ والتنفيذات تستخدم \(2(n-1)+8\) حدود تمهيدية.

الخطوة 5: حساب الحد البعيد باختزال كثيرات الحدود

بعد معرفة العلاقة العودية، تتحول المسألة إلى حساب سريع للحد رقم \(m\). نختزل \(x^{m-1}\) بترديد كثير حدود العودية:

$$x^{m-1}\equiv q_0+q_1x+\cdots+q_{r-1}x^{r-1}\pmod{x^r-c_1x^{r-1}-\cdots-c_r}.$$

وعندئذ يكون

$$T(n,m)\equiv q_0T(n,1)+q_1T(n,2)+\cdots+q_{r-1}T(n,r)\pmod{10^9+7}.$$

وتُحسب المعاملات \(q_i\) بواسطة الأسّ الثنائي داخل الحلقة القسمة التي تحددها العودية، ولذلك فإن الفهرس الضخم \(m=10^{12}\) يضيف فقط عامل \(\log m\) في الزمن.

مثال محلول: \(n=3\)

في هذه الحالة القيم المسموح بها هي \(1\) و\(2\) فقط. وقيد الجوار يمنع الزوج \((2,2)\) وحده، لذا تكون

$$A=\begin{pmatrix} 1 & 1\\ 1 & 0 \end{pmatrix}.$$

وتكون المجاميع بحسب الطول هي

$$T(3,1)=2,\qquad T(3,2)=3,\qquad T(3,3)=5,\qquad T(3,4)=8.$$

وفي هذا المثال الصغير تكون العودية ببساطة

$$T(3,t)=T(3,t-1)+T(3,t-2),$$

ومن ثم نحصل مباشرة على قيمة التحقق \(T(3,4)=8\).

كيف تعمل الشيفرة

تنفيذات C++ وPython وJava تتبع جميعها الخطة الرياضية نفسها. في المرحلة الأولى تُولِّد كتلة من الحدود الابتدائية باستعمال برمجة ديناميكية مصنفة بحسب القيمة الأخيرة. وبدل إعادة حساب كل انتقال بشكل ساذج، تُحسب المجاميع البادئة للطبقة الحالية ثم تُقرأ بترتيب معكوس، وبذلك تصبح كلفة خطوة انتقال واحدة \(O(n)\).

بعد ذلك تُطبَّق خوارزمية Berlekamp-Massey بترديد \(10^9+7\) على هذه الحدود الابتدائية للحصول على أقصر عودية تحقق \(T(n,m)\). إذا كان الطول المطلوب موجودًا أصلًا داخل البادئة المحسوبة، تُعاد القيمة مباشرة. وإلا تُجري الشيفرة أسًّا ثنائيًا على كثيرات حدود مختزلة، فتستخرج المعاملات \(q_0,\dots,q_{r-1}\) ومن ثم تبني الحد النهائي المطلوب.

تحليل التعقيد

كل خطوة من خطوات البرمجة الديناميكية التمهيدية تكلف \(O(n)\) زمنًا و\(O(n)\) ذاكرة. وبما أن عدد الحدود الابتدائية المولدة هو \(O(n)\)، فإن هذه المرحلة تكلف \(O(n^2)\) زمنًا إجمالًا. خوارزمية Berlekamp-Massey لعلاقة من رتبة \(r\le n-1\) تكلف \(O(r^2)\)، ومرحلة الحساب السريع للحد رقم \(m\) تكلف \(O(r^2\log m)\). لذلك تكون الكلفة الكلية

$$O(n^2+r^2+r^2\log m)=O(n^2\log m),$$

مع ذاكرة \(O(n)\).

حواش ومراجع

  1. صفحة المسألة في Project Euler: Problem 654: Neighbourly Constraints
  2. خوارزمية Berlekamp-Massey: Wikipedia - Berlekamp-Massey algorithm
  3. الحساب السريع للعلاقات العودية الخطية: cp-algorithms - Linear Recurrence
  4. مبرهنة كايلي-هاملتون: Wikipedia - Cayley-Hamilton theorem
  5. المجاميع البادئة: Wikipedia - Prefix sum