Problem Summary

For a fixed \(n\), we count sequences \((a_1,a_2,\dots,a_n)\) with each \(a_i \in \{1,2,\dots,n\}\) and with the iterated-function condition

$$a_{i+1}=a_{a_i}\qquad (1 \le i < n).$$

Project Euler 977 asks for this count when \(n=10^6\), reported modulo \(10^9+7\). The brute-force search space has size \(n^n\), so the only viable route is to understand the structure forced by the recurrence and then sum the resulting combinatorial classes in closed form.

Mathematical Approach

The key idea is that the recurrence does not describe arbitrary sequences: it describes the forward orbit of a single self-map, and every valid sequence eventually becomes periodic. The solution classifies sequences by the first place where that periodic regime begins.

From the recurrence to an orbit

Define a function \(f:\{1,\dots,n\}\to\{1,\dots,n\}\) by \(f(i)=a_i\). Then

$$a_{i+1}=a_{a_i}=f(a_i)=f(f(i)).$$

Inductively this gives

$$a_i=f^i(1)\qquad (i \ge 1),$$

so the sequence is exactly the orbit of 1 under repeated application of \(f\). More generally, the orbit of index \(m\) is just a shifted suffix:

$$f^t(m)=a_{m+t-1}\qquad (t \ge 1).$$

This has an important consequence: if two positions carry the same value, then their entire future tails coincide. In particular, once a suffix starts repeating with some period, everything after that point is rigid.

Classify by the first periodic suffix

Let \(s\) be the first index such that the suffix

$$a_s,a_{s+1},\dots,a_n$$

is periodic from its first term. Write

$$N=n-s+1$$

for the suffix length, and let its exact period be \(l\). Then the suffix positions split into residue classes modulo \(l\):

$$S_u=\{\,s+u-1+ml : m \ge 0,\ s+u-1+ml \le n\,\}\qquad (1 \le u \le l).$$

If \(N=ql+r\) with \(0 \le r < l\), then the first \(r\) classes have size \(q+1\) and the remaining \(l-r\) classes have size \(q\).

Count one periodic suffix

Because the suffix has period \(l\), every position in the same class \(S_u\) carries the same value; call it \(c_u\). The recurrence forces a simple rule:

$$c_u \in S_{u+1}\quad (1 \le u < l),\qquad c_l \in S_1.$$

So choosing the suffix means choosing one element from the next residue class for each \(u\). The number of choices is therefore

$$P(N,l)=\prod_{u=1}^{l}|S_u|=q^{\,l-r}(q+1)^r,$$

where

$$q=\left\lfloor \frac{N}{l}\right\rfloor,\qquad r=N \bmod l.$$

This is the basic factor that appears everywhere in the implementations.

Attach the non-periodic prefix

If the periodic part starts at the very beginning, so \(s=1\) and \(N=n\), there is no prefix to attach. Those sequences contribute

$$A(n)=\sum_{l=1}^{n} P(n,l).$$

Now assume \(s>1\). Then \(a_{s-1}\) must be chosen so that

$$a_s=a_{a_{s-1}}=c_1.$$

The positions whose value is \(c_1\) are exactly the elements of \(S_1\), so \(a_{s-1}\) must lie in \(S_1\). But one choice is forbidden: the special element \(c_l \in S_1\). If we took \(a_{s-1}=c_l\), then the same \(l\)-periodic pattern would already start at position \(s-1\), contradicting the minimality of \(s\).

Therefore the number of admissible attachments is

$$|S_1|-1=\left\lceil \frac{N}{l}\right\rceil - 1,$$

which is \(q-1\) when \(r=0\) and \(q\) when \(r>0\). Once \(a_{s-1}\) is fixed, every earlier position must be the tautological forward link

$$a_i=i+1\qquad (1 \le i \le s-2),$$

because any earlier nontrivial copy would make the periodic regime begin even sooner.

So the full count is

$$F(n)=\sum_{l=1}^{n} P(n,l)+\sum_{N=1}^{n-1}\sum_{l=1}^{N}\left(\left\lceil \frac{N}{l}\right\rceil-1\right)P(N,l).$$

Worked Example: \(n=7\), suffix length \(N=4\), period \(l=2\)

Take \(s=4\), so the periodic suffix is \(a_4,a_5,a_6,a_7\). The residue classes are

$$S_1=\{4,6\},\qquad S_2=\{5,7\}.$$

To build a 2-periodic suffix, choose \(c_1 \in S_2\) and \(c_2 \in S_1\). There are

$$P(4,2)=2\cdot 2=4$$

choices. For example, \(c_1=5\) and \(c_2=6\) give the suffix

$$5,6,5,6.$$

The previous term \(a_3\) must lie in \(S_1\), but it cannot equal \(c_2=6\), otherwise the 2-periodic pattern would already begin at position 3. So \(a_3=4\) is forced, and then the earlier prefix is the rigid chain \(a_1=2\), \(a_2=3\). One valid sequence is therefore

$$ (2,3,4,5,6,5,6). $$

All four suffix choices work in the same way, so this class contributes

$$\left(\left\lceil \frac{4}{2}\right\rceil-1\right)P(4,2)=1 \cdot 4=4$$

sequences, exactly as the formula predicts.

Regroup the double sum by quotient blocks

The direct formula above is already correct, and the slower validation routines evaluate it exactly. The fast solver reorganizes the second double sum. For fixed \(l\), write

$$N=ql+r,\qquad 0 \le r < l,$$

so \(q=\lfloor N/l \rfloor\). Then

$$P(N,l)=q^{\,l-r}(q+1)^r.$$

When \(N\) runs through the block \(ql,ql+1,\dots,(q+1)l-1\), the attachment factor is \(q-1\) at \(r=0\) and \(q\) for \(r \ge 1\). If

$$R=\min(l-1,n-1-ql),$$

the whole block contributes

$$B_{l,q}=(q-1)q^l+\sum_{r=1}^{R} q^{\,l+1-r}(q+1)^r.$$

The inner sum is telescoping:

$$\sum_{r=1}^{R} q^{\,l+1-r}(q+1)^r=q^{\,l+1-R}(q+1)^{R+1}-q^{\,l+1}(q+1).$$

This is exactly the closed form used by the production implementations.

How the Code Works

Power-table precomputation

The C++, Python, and Java implementations precompute all powers that can appear later. For a base \(b\), the largest exponent ever needed is at most

$$\left\lfloor \frac{n-1}{b-1}\right\rfloor + 1,$$

so every required value \(b^e \bmod (10^9+7)\) can be read in \(O(1)\) time from a packed table. This avoids an enormous number of repeated modular exponentiations inside the main summation.

Two layers of counting

The implementations first evaluate

$$A(n)=\sum_{l=1}^{n}P(n,l),$$

which counts sequences whose periodic regime begins at position 1. They then add the correction term for later starts, but not as a naive triple loop over \((N,l,r)\). Instead, for each \(l\) they iterate over quotient plateaus \(q=\lfloor N/l \rfloor\), compute the block tail length \(R\), and add the closed block sum \(B_{l,q}\). That is why the code matches the mathematical formula but runs far faster.

Validation and parallel execution

Each implementation also contains small checkpoints: exhaustive enumeration confirms that the count for \(n=7\) is 174, and a direct evaluation of the ungathered double sum confirms that the count for \(n=100\) is 305741269. After that, the large computation for \(n=10^6\) uses the precomputed power table and the regrouped block formulas. The C++ and Java implementations split the \(l\)-range across worker threads, while the Python implementation performs the same arithmetic serially.

Complexity Analysis

The dominant work is not exponential anymore. Building the packed power tables costs

$$\sum_{b=2}^{n+1} O\!\left(\frac{n}{b-1}\right)=O(n \log n),$$

and the block summation for the correction term has the same order because

$$\sum_{l=1}^{n-1}\left\lfloor\frac{n-1}{l}\right\rfloor = O(n \log n).$$

So the fast method runs in \(O(n \log n)\) time and uses \(O(n \log n)\) memory for the power tables. In practice it is efficient because every block contribution is reduced to a handful of modular multiplications and table lookups.

Footnotes and References

  1. Problem page: https://projecteuler.net/problem=977
  2. Iterated function: Wikipedia - Iterated function
  3. Functional graph: Wikipedia - Functional graph
  4. Eventually periodic sequence: Wikipedia - Eventually periodic points
  5. Floor and ceiling functions: Wikipedia - Floor and ceiling functions
  6. Geometric series: Wikipedia - Geometric series
  7. Modular arithmetic: Wikipedia - Modular arithmetic

Problemzusammenfassung

Gesucht ist die Anzahl der Folgen \((a_1,a_2,\dots,a_n)\) mit \(a_i \in \{1,2,\dots,n\}\), die die Bedingung

$$a_{i+1}=a_{a_i}\qquad (1 \le i < n)$$

erfüllen. In Problem 977 soll dieser Wert für \(n=10^6\) modulo \(10^9+7\) berechnet werden. Ein naiver Ansatz hätte Suchraumgröße \(n^n\) und ist damit völlig ausgeschlossen. Man muss die durch die Rekursion erzwungene Struktur direkt ausnutzen.

Mathematischer Ansatz

Die entscheidende Beobachtung ist, dass die Folge nicht frei ist, sondern als Iteration einer einzigen Abbildung entsteht. Dadurch wird sie nach endlich vielen Schritten periodisch, und genau diese erste periodische Phase kann man sauber zählen.

Von der Rekursion zur Bahn einer Abbildung

Setze \(f(i)=a_i\). Dann gilt

$$a_{i+1}=a_{a_i}=f(a_i)=f(f(i)).$$

Induktiv folgt daraus

$$a_i=f^i(1)\qquad (i \ge 1).$$

Die Folge ist also genau die Vorwärtsbahn von 1 unter \(f\). Sogar allgemeiner gilt für jeden Startindex \(m\)

$$f^t(m)=a_{m+t-1}\qquad (t \ge 1).$$

Wenn zwei Positionen denselben Wert tragen, dann stimmen daher ihre gesamten Fortsetzungen überein. Sobald also ein periodischer Abschnitt beginnt, ist der restliche Verlauf vollständig festgelegt.

Klassifikation nach dem ersten periodischen Suffix

Sei \(s\) der erste Index, ab dem das Suffix

$$a_s,a_{s+1},\dots,a_n$$

sofort periodisch ist. Schreibe

$$N=n-s+1$$

für seine Länge und \(l\) für seine exakte Periode. Dann zerfallen die Positionen dieses Suffixes in Restklassen modulo \(l\):

$$S_u=\{\,s+u-1+ml : m \ge 0,\ s+u-1+ml \le n\,\}\qquad (1 \le u \le l).$$

Ist \(N=ql+r\) mit \(0 \le r < l\), dann haben die ersten \(r\) Klassen Größe \(q+1\), die übrigen \(l-r\) Klassen Größe \(q\).

Wie viele periodische Suffixe gibt es?

Innerhalb jeder Klasse \(S_u\) ist der Wert konstant; nenne ihn \(c_u\). Die Rekursion erzwingt

$$c_u \in S_{u+1}\quad (1 \le u < l),\qquad c_l \in S_1.$$

Man wählt also für jede Klasse genau ein Element aus der jeweils nächsten Klasse. Daraus folgt unmittelbar

$$P(N,l)=\prod_{u=1}^{l}|S_u|=q^{\,l-r}(q+1)^r,$$

wobei

$$q=\left\lfloor \frac{N}{l}\right\rfloor,\qquad r=N \bmod l.$$

Dieses Produkt ist das Grundobjekt der ganzen Lösung.

Wie der nichtperiodische Präfix angehängt wird

Beginnt die periodische Phase sofort, also \(s=1\) und \(N=n\), gibt es keinen Präfix. Dann ist der Beitrag

$$A(n)=\sum_{l=1}^{n}P(n,l).$$

Nun sei \(s>1\). Dann muss \(a_{s-1}\) so gewählt werden, dass

$$a_s=a_{a_{s-1}}=c_1$$

gilt. Genau die Elemente aus \(S_1\) besitzen den Wert \(c_1\), also muss \(a_{s-1}\in S_1\) liegen. Eine Wahl ist jedoch verboten: das Element \(c_l \in S_1\). Würde man \(a_{s-1}=c_l\) setzen, dann begänne dieselbe \(l\)-Periodizität schon an Position \(s-1\).

Damit bleiben

$$|S_1|-1=\left\lceil \frac{N}{l}\right\rceil - 1$$

Möglichkeiten. Das ist \(q-1\), falls \(r=0\), und sonst \(q\). Sind diese festgelegt, dann müssen alle früheren Stellen die triviale Vorwärtskette bilden:

$$a_i=i+1\qquad (1 \le i \le s-2),$$

denn jede frühere nichttriviale Kopie würde die periodische Phase noch weiter nach vorne verschieben.

Die Gesamtformel lautet also

$$F(n)=\sum_{l=1}^{n} P(n,l)+\sum_{N=1}^{n-1}\sum_{l=1}^{N}\left(\left\lceil \frac{N}{l}\right\rceil-1\right)P(N,l).$$

Durchgerechnetes Beispiel: \(n=7\), \(N=4\), \(l=2\)

Setze \(s=4\). Das periodische Suffix ist also \(a_4,a_5,a_6,a_7\), und die Klassen sind

$$S_1=\{4,6\},\qquad S_2=\{5,7\}.$$

Für ein 2-periodisches Suffix wählt man \(c_1 \in S_2\) und \(c_2 \in S_1\). Das liefert

$$P(4,2)=2\cdot 2=4$$

Möglichkeiten. Wählt man etwa \(c_1=5\) und \(c_2=6\), dann lautet das Suffix

$$5,6,5,6.$$

Nun muss \(a_3 \in S_1\) liegen, darf aber nicht \(c_2=6\) sein, sonst würde die Periodizität schon bei Position 3 beginnen. Also ist \(a_3=4\), und davor bleibt nur die starre Kette \(a_1=2\), \(a_2=3\). Man erhält zum Beispiel

$$ (2,3,4,5,6,5,6). $$

Alle vier Suffixwahlen funktionieren genauso, also trägt diese Klasse genau

$$\left(\left\lceil \frac{4}{2}\right\rceil-1\right)P(4,2)=1\cdot 4=4$$

Folgen bei.

Umgruppierung nach Quotientenblöcken

Die Formel oben ist bereits korrekt und wird von den langsameren Prüfungen direkt ausgewertet. Für die schnelle Version wird der zweite Summand umgeordnet. Für festes \(l\) schreibe

$$N=ql+r,\qquad 0 \le r < l,$$

also \(q=\lfloor N/l \rfloor\). Dann ist

$$P(N,l)=q^{\,l-r}(q+1)^r.$$

Läuft \(N\) durch den Block \(ql,ql+1,\dots,(q+1)l-1\), so ist der Anheftungsfaktor bei \(r=0\) gleich \(q-1\) und für \(r \ge 1\) gleich \(q\). Mit

$$R=\min(l-1,n-1-ql)$$

wird daraus der Blockbeitrag

$$B_{l,q}=(q-1)q^l+\sum_{r=1}^{R} q^{\,l+1-r}(q+1)^r.$$

Die innere Summe teleskopiert zu

$$\sum_{r=1}^{R} q^{\,l+1-r}(q+1)^r=q^{\,l+1-R}(q+1)^{R+1}-q^{\,l+1}(q+1).$$

Genau diese Form wird im schnellen Teil der Implementierungen verwendet.

Wie der Code arbeitet

Vorberechnung aller benötigten Potenzen

Die C++-, Python- und Java-Implementierungen legen zunächst eine große Tabelle aller Potenzen an, die später vorkommen können. Für eine Basis \(b\) wird höchstens der Exponent

$$\left\lfloor \frac{n-1}{b-1}\right\rfloor + 1$$

benötigt. Dadurch kann jede vorkommende Potenz \(b^e \bmod (10^9+7)\) später in \(O(1)\) nachgeschlagen werden.

Zwei Zählebenen

Zuerst wird

$$A(n)=\sum_{l=1}^{n}P(n,l)$$

berechnet, also der Anteil der Folgen, deren periodische Phase sofort beginnt. Danach folgt der Korrekturterm für spätere Startpositionen. Dieser wird aber nicht naiv über alle \((N,l,r)\) iteriert, sondern nach Quotientenplateaus \(q=\lfloor N/l \rfloor\) gruppiert. Pro Block genügt dann eine konstante Anzahl modularer Operationen.

Prüfungen und Parallelisierung

Alle drei Implementierungen enthalten kleine Selbsttests: vollständiges Durchzählen bestätigt \(174\) für \(n=7\), und die direkte ungruppierte Summation bestätigt \(305741269\) für \(n=100\). Erst danach wird die große Rechnung für \(n=10^6\) gestartet. C++ und Java teilen dabei den \(l\)-Bereich auf mehrere Threads auf; Python führt dieselbe Mathematik seriell aus.

Komplexitätsanalyse

Die Exponentialität ist vollständig verschwunden. Die Potenztabellen kosten

$$\sum_{b=2}^{n+1} O\!\left(\frac{n}{b-1}\right)=O(n \log n),$$

und für die Quotientenblöcke gilt analog

$$\sum_{l=1}^{n-1}\left\lfloor\frac{n-1}{l}\right\rfloor=O(n \log n).$$

Damit läuft die schnelle Methode in \(O(n \log n)\) Zeit und benötigt \(O(n \log n)\) Speicher für die vorberechneten Potenzen. Praktisch ist sie schnell, weil jede Blocksumme auf wenige Tabellenzugriffe und modulare Multiplikationen reduziert wird.

Fußnoten und Referenzen

  1. Problemseite: https://projecteuler.net/problem=977
  2. Iterierte Funktion: Wikipedia - Iterated function
  3. Funktionaler Graph: Wikipedia - Functional graph
  4. Schließlich periodische Folgen: Wikipedia - Eventually periodic points
  5. Gaußklammer und Aufrundung: Wikipedia - Floor and ceiling functions
  6. Geometrische Reihe: Wikipedia - Geometric series
  7. Modulare Arithmetik: Wikipedia - Modular arithmetic

Problem Özeti

Aranan şey, her terimi \(\{1,2,\dots,n\}\) içinde olan ve

$$a_{i+1}=a_{a_i}\qquad (1 \le i < n)$$

koşulunu sağlayan \((a_1,a_2,\dots,a_n)\) dizilerinin sayısıdır. Problem 977 bu sayıyı \(n=10^6\) için, \(10^9+7\) modunda istiyor. Düz kuvvet denemesi \(n^n\) olasılık gerektirdiği için çözüm tamamen yapısal bir sınıflandırmaya dayanmak zorundadır.

Matematiksel Yaklaşım

Temel fikir şu: Bu dizi koşulu rastgele bir yerel bağıntı değil, tek bir fonksiyonun 1 noktasından başlayan iterasyonunu tarif ediyor. Bu yüzden her geçerli dizi sonunda periyodik hale gelir ve doğru sınıflandırma, bu periyodik rejimin ilk başladığı yeri izlemektir.

Bağıntıyı bir iterasyon problemi olarak yazmak

\(f(i)=a_i\) tanımını yapalım. O zaman

$$a_{i+1}=a_{a_i}=f(a_i)=f(f(i)).$$

Bundan tümevarımla

$$a_i=f^i(1)\qquad (i \ge 1)$$

elde edilir. Yani dizi, 1'in \(f\) altındaki ileri yörüngesidir. Hatta her \(m\) için

$$f^t(m)=a_{m+t-1}\qquad (t \ge 1)$$

olur. Bu çok güçlü bir değişmezdir: iki konum aynı değeri taşıyorsa, o iki konumdan sonraki tüm kuyruk da aynıdır. Dolayısıyla periyodiklik başladıktan sonra geri kalan her şey katı biçimde belirlenir.

İlk periyodik son eke göre sınıflandırma

\(s\),

$$a_s,a_{s+1},\dots,a_n$$

son ekinin kendi ilk teriminden itibaren periyodik olduğu ilk indeks olsun. Son ekin uzunluğunu

$$N=n-s+1$$

ve tam periyodunu \(l\) ile gösterelim. O zaman bu son ekin konumları \(l\)'ye göre artık sınıflarına ayrılır:

$$S_u=\{\,s+u-1+ml : m \ge 0,\ s+u-1+ml \le n\,\}\qquad (1 \le u \le l).$$

Eğer \(N=ql+r\) ve \(0 \le r < l\) ise, ilk \(r\) sınıfın boyu \(q+1\), geri kalan \(l-r\) sınıfın boyu \(q\) olur.

Tek bir periyodik son ek nasıl sayılır?

Son ek \(l\)-periyodik olduğundan, aynı sınıftaki bütün konumların değeri aynıdır; buna \(c_u\) diyelim. Bağıntı şu zorunlu kuralı verir:

$$c_u \in S_{u+1}\quad (1 \le u < l),\qquad c_l \in S_1.$$

Yani her sınıf için bir sonraki sınıftan bir indeks seçiyoruz. Bu yüzden seçim sayısı doğrudan

$$P(N,l)=\prod_{u=1}^{l}|S_u|=q^{\,l-r}(q+1)^r$$

olur; burada

$$q=\left\lfloor \frac{N}{l}\right\rfloor,\qquad r=N \bmod l.$$

Koddaki bütün kuvvet terimleri bu ürünün farklı görünümleridir.

Periyodik olmayan önekin eklenmesi

Eğer periyodik bölüm en başta başlıyorsa, yani \(s=1\) ve \(N=n\), ekstra bir önek yoktur. Bu sınıfın toplamı

$$A(n)=\sum_{l=1}^{n}P(n,l)$$

şeklindedir.

Şimdi \(s>1\) olsun. O zaman

$$a_s=a_{a_{s-1}}=c_1$$

olmalı. Değeri \(c_1\) olan konumlar tam olarak \(S_1\)'deki elemanlardır; dolayısıyla \(a_{s-1}\) mutlaka \(S_1\) içinden seçilir. Fakat bu seçimlerden biri yasaktır: \(c_l \in S_1\). Eğer \(a_{s-1}=c_l\) seçilseydi aynı \(l\)-periyodik düzen bir adım daha erken, yani \(s-1\)'de başlamış olurdu.

Bu nedenle geçerli bağlanma sayısı

$$|S_1|-1=\left\lceil \frac{N}{l}\right\rceil - 1$$

olur. Bu da \(r=0\) için \(q-1\), \(r>0\) için \(q\)'dur. Bu seçim yapıldıktan sonra daha önceki bütün konumlar zorunlu olarak

$$a_i=i+1\qquad (1 \le i \le s-2)$$

şeklindeki ileri zincirdir; çünkü daha erken bir kopyalama olsaydı periyodik rejim daha da öne kayardı.

Sonuç olarak toplam sayı

$$F(n)=\sum_{l=1}^{n} P(n,l)+\sum_{N=1}^{n-1}\sum_{l=1}^{N}\left(\left\lceil \frac{N}{l}\right\rceil-1\right)P(N,l)$$

olur.

Çalışılmış örnek: \(n=7\), \(N=4\), \(l=2\)

\(s=4\) alalım. Periyodik son ek \(a_4,a_5,a_6,a_7\) olur ve sınıflar

$$S_1=\{4,6\},\qquad S_2=\{5,7\}$$

şeklindedir. 2-periyodik bir son ek kurmak için \(c_1 \in S_2\) ve \(c_2 \in S_1\) seçilir. Toplam

$$P(4,2)=2\cdot 2=4$$

olası son ek vardır. Örneğin \(c_1=5\), \(c_2=6\) seçimi

$$5,6,5,6$$

son ekini üretir. Önceki terim \(a_3\) mutlaka \(S_1\) içinde olmalıdır, ama \(c_2=6\) olamaz; aksi halde periyot 3. konumda başlamış olur. Dolayısıyla \(a_3=4\) zorunludur. Sonra önek zinciri \(a_1=2\), \(a_2=3\) gelir ve

$$ (2,3,4,5,6,5,6) $$

dizisini elde ederiz. Dört son ek seçiminin hepsi aynı mantıkla çalıştığından bu sınıfın katkısı

$$\left(\left\lceil \frac{4}{2}\right\rceil-1\right)P(4,2)=1\cdot 4=4$$

olur.

Bölüm bloklarına göre yeniden düzenleme

Yukarıdaki formül doğrudur ve yavaş doğrulama kodları onu aynen hesaplar. Hızlı çözüm ise ikinci toplamı bloklara ayırır. Sabit \(l\) için

$$N=ql+r,\qquad 0 \le r < l$$

yazalım; burada \(q=\lfloor N/l \rfloor\). O zaman

$$P(N,l)=q^{\,l-r}(q+1)^r$$

olur. \(N\), \(ql,ql+1,\dots,(q+1)l-1\) aralığında gezerken bağlanma katsayısı \(r=0\) için \(q-1\), diğer durumlarda \(q\) olur. Eğer

$$R=\min(l-1,n-1-ql)$$

ise, blok katkısı

$$B_{l,q}=(q-1)q^l+\sum_{r=1}^{R} q^{\,l+1-r}(q+1)^r$$

şeklindedir. İç toplam teleskopik olarak

$$\sum_{r=1}^{R} q^{\,l+1-r}(q+1)^r=q^{\,l+1-R}(q+1)^{R+1}-q^{\,l+1}(q+1)$$

olur. Uygulamaların hızlı kısmı tam olarak bu kapalı biçimi kullanır.

Kod Nasıl Çalışır

Önceden kuvvet tablosu kurma

C++, Python ve Java uygulamaları önce ileride ihtiyaç duyulacak bütün taban-kuvvet değerlerini paketlenmiş bir tabloda hazırlar. Bir \(b\) tabanı için gerekli en büyük üs

$$\left\lfloor \frac{n-1}{b-1}\right\rfloor + 1$$

kadardır. Böylece daha sonra gereken her \(b^e \bmod (10^9+7)\) değeri sabit zamanda okunur.

İki ayrı toplama katmanı

İlk olarak

$$A(n)=\sum_{l=1}^{n}P(n,l)$$

hesaplanır; bu, periyodik kısmı en başta başlayan dizileri sayar. Sonra daha geç başlayan sınıfların düzeltme terimi eklenir. Burada naif bir \((N,l,r)\) üçlü döngüsü kurulmaz; bunun yerine her \(l\) için sabit bölüm değeri veren \(q=\lfloor N/l \rfloor\) platoları dolaşılır ve her blok için tek kapalı form uygulanır.

Doğrulama ve paralellik

Her üç uygulama da küçük denetimler içerir: tam tarama ile \(n=7\) için sonucun 174 olduğu doğrulanır, sonra da gruplanmamış çift toplam kullanılarak \(n=100\) için 305741269 doğrulanır. Büyük hesaplama bundan sonra \(n=10^6\) için yapılır. C++ ve Java sürümleri \(l\) aralığını iş parçacıklarına böler; Python sürümü aynı aritmetiği tek iş parçacığında yürütür.

Karmaşıklık Analizi

Üstel maliyet tamamen ortadan kalkmıştır. Kuvvet tablolarının maliyeti

$$\sum_{b=2}^{n+1} O\!\left(\frac{n}{b-1}\right)=O(n \log n)$$

ve bölüm bloklarının toplam sayısı da

$$\sum_{l=1}^{n-1}\left\lfloor\frac{n-1}{l}\right\rfloor=O(n \log n)$$

mertebesindedir. Dolayısıyla hızlı yöntem \(O(n \log n)\) zamanda çalışır ve kuvvet tabloları için \(O(n \log n)\) bellek kullanır. Pratikte hızın kaynağı, her blok katkısının birkaç modüler çarpım ve tablo erişimiyle hesaplanabilmesidir.

Dipnotlar ve Kaynaklar

  1. Problem sayfası: https://projecteuler.net/problem=977
  2. İterasyonlu fonksiyon: Wikipedia - Iterated function
  3. Fonksiyonel grafik: Wikipedia - Functional graph
  4. Sonunda periyodik noktalar: Wikipedia - Eventually periodic points
  5. Taban ve tavan fonksiyonları: Wikipedia - Floor and ceiling functions
  6. Geometrik seri: Wikipedia - Geometric series
  7. Modüler aritmetik: Wikipedia - Modular arithmetic

Resumen del Problema

Se quiere contar cuantas secuencias \((a_1,a_2,\dots,a_n)\) con \(a_i \in \{1,2,\dots,n\}\) satisfacen

$$a_{i+1}=a_{a_i}\qquad (1 \le i < n).$$

En el problema 977 la evaluacion final es para \(n=10^6\), modulo \(10^9+7\). Un barrido directo sobre las \(n^n\) secuencias posibles es imposible; la solucion real tiene que extraer la estructura combinatoria escondida en la condicion funcional.

Enfoque Matemático

La recurrencia describe la iteracion de una sola funcion sobre el punto 1. Eso fuerza a que toda secuencia valida termine entrando en un regimen periodico. La forma mas natural de contar es clasificar por la primera posicion donde empieza ese sufijo periodico.

La recurrencia como orbita de una funcion

Definamos \(f(i)=a_i\). Entonces

$$a_{i+1}=a_{a_i}=f(a_i)=f(f(i)).$$

Por induccion se obtiene

$$a_i=f^i(1)\qquad (i \ge 1).$$

Es decir, la secuencia es exactamente la orbita hacia delante de 1 bajo \(f\). Mas aun, para cualquier indice \(m\),

$$f^t(m)=a_{m+t-1}\qquad (t \ge 1).$$

Esto implica que si dos posiciones tienen el mismo valor, entonces sus colas futuras coinciden por completo. Una vez que aparece la periodicidad, el resto de la secuencia queda rigidamente determinado.

Clasificar por el primer sufijo periodico

Sea \(s\) el primer indice tal que el sufijo

$$a_s,a_{s+1},\dots,a_n$$

es periodico desde su primer termino. Escribimos

$$N=n-s+1$$

para la longitud de ese sufijo y \(l\) para su periodo exacto. Las posiciones del sufijo se separan en clases de residuo modulo \(l\):

$$S_u=\{\,s+u-1+ml : m \ge 0,\ s+u-1+ml \le n\,\}\qquad (1 \le u \le l).$$

Si \(N=ql+r\) con \(0 \le r < l\), las primeras \(r\) clases tienen tamano \(q+1\) y las otras \(l-r\) tienen tamano \(q\).

Contar un sufijo periodico fijo

Dentro de cada clase \(S_u\) el valor es constante; llamemos \(c_u\) a ese valor. La condicion funcional obliga a que

$$c_u \in S_{u+1}\quad (1 \le u < l),\qquad c_l \in S_1.$$

Asi, construir el sufijo equivale a elegir para cada clase un elemento de la clase siguiente. El numero total de elecciones es

$$P(N,l)=\prod_{u=1}^{l}|S_u|=q^{\,l-r}(q+1)^r,$$

donde

$$q=\left\lfloor \frac{N}{l}\right\rfloor,\qquad r=N \bmod l.$$

Ese producto es el bloque elemental de toda la derivacion.

Como se engancha el prefijo no periodico

Si el regimen periodico empieza desde el inicio, es decir \(s=1\) y \(N=n\), no hay prefijo adicional. Esa parte aporta

$$A(n)=\sum_{l=1}^{n}P(n,l).$$

Supongamos ahora \(s>1\). Entonces el termino anterior debe cumplir

$$a_s=a_{a_{s-1}}=c_1.$$

Los indices cuyo valor es \(c_1\) son exactamente los elementos de \(S_1\), asi que \(a_{s-1}\) debe estar en \(S_1\). Pero una eleccion esta prohibida: \(c_l \in S_1\). Si se tomara \(a_{s-1}=c_l\), el mismo patron \(l\)-periodico empezaria ya en la posicion \(s-1\).

Por tanto, el numero de acoplamientos validos es

$$|S_1|-1=\left\lceil \frac{N}{l}\right\rceil - 1,$$

que vale \(q-1\) cuando \(r=0\) y \(q\) cuando \(r>0\). Fijado ese enlace, todas las posiciones anteriores quedan forzadas a la cadena trivial

$$a_i=i+1\qquad (1 \le i \le s-2),$$

porque cualquier copia no trivial antes de eso haria que la parte periodica empezara mas pronto.

En consecuencia,

$$F(n)=\sum_{l=1}^{n} P(n,l)+\sum_{N=1}^{n-1}\sum_{l=1}^{N}\left(\left\lceil \frac{N}{l}\right\rceil-1\right)P(N,l).$$

Ejemplo trabajado: \(n=7\), \(N=4\), \(l=2\)

Tomemos \(s=4\). Entonces el sufijo periodico es \(a_4,a_5,a_6,a_7\), con clases

$$S_1=\{4,6\},\qquad S_2=\{5,7\}.$$

Para construir un sufijo 2-periodico elegimos \(c_1 \in S_2\) y \(c_2 \in S_1\). Hay

$$P(4,2)=2\cdot 2=4$$

opciones. Por ejemplo, \(c_1=5\) y \(c_2=6\) producen

$$5,6,5,6.$$

Ahora \(a_3\) tiene que pertenecer a \(S_1\), pero no puede ser \(c_2=6\), porque entonces la periodicidad empezaria ya en la posicion 3. Asi que \(a_3=4\), y antes de eso queda forzada la cadena \(a_1=2\), \(a_2=3\). Una secuencia valida es

$$ (2,3,4,5,6,5,6). $$

Las cuatro elecciones del sufijo funcionan del mismo modo, de manera que esta clase aporta

$$\left(\left\lceil \frac{4}{2}\right\rceil-1\right)P(4,2)=1\cdot 4=4$$

secuencias.

Reagrupar la suma por bloques de cociente

La formula anterior ya es correcta y las rutinas lentas de comprobacion la evaluan tal cual. La version rapida reorganiza el segundo sumando. Para un \(l\) fijo escribimos

$$N=ql+r,\qquad 0 \le r < l,$$

de modo que \(q=\lfloor N/l \rfloor\). Entonces

$$P(N,l)=q^{\,l-r}(q+1)^r.$$

Cuando \(N\) recorre el bloque \(ql,ql+1,\dots,(q+1)l-1\), el factor de enganche vale \(q-1\) si \(r=0\) y \(q\) si \(r \ge 1\). Si

$$R=\min(l-1,n-1-ql),$$

el bloque completo contribuye

$$B_{l,q}=(q-1)q^l+\sum_{r=1}^{R} q^{\,l+1-r}(q+1)^r.$$

La suma interior es telescopica:

$$\sum_{r=1}^{R} q^{\,l+1-r}(q+1)^r=q^{\,l+1-R}(q+1)^{R+1}-q^{\,l+1}(q+1).$$

Esa es exactamente la identidad que usan las implementaciones rapidas.

Cómo Funciona el Código

Precomputacion de potencias

Las implementaciones en C++, Python y Java precalculan todas las potencias modulares que pueden aparecer. Para una base \(b\), el mayor exponente necesario es

$$\left\lfloor \frac{n-1}{b-1}\right\rfloor + 1,$$

de modo que cualquier \(b^e \bmod (10^9+7)\) se recupera despues en tiempo \(O(1)\) desde una tabla compacta.

Dos niveles de suma

Primero se evalua

$$A(n)=\sum_{l=1}^{n}P(n,l),$$

que cuenta las secuencias cuyo comportamiento periodico empieza en la posicion 1. Luego se agrega la correccion para los comienzos mas tardios. En vez de recorrer ingenuamente todos los triples \((N,l,r)\), el codigo agrupa por mesetas del cociente \(q=\lfloor N/l \rfloor\) y suma cada bloque con la forma cerrada anterior.

Verificacion y paralelismo

Cada implementacion contiene verificaciones pequenas: un barrido exhaustivo confirma que para \(n=7\) el conteo es 174, y la suma directa no reagruipada confirma que para \(n=100\) el valor es 305741269. Despues de eso se calcula el caso grande \(n=10^6\). Las versiones en C++ y Java reparten el rango de \(l\) entre varios hilos; la version en Python hace la misma aritmetica en serie.

Análisis de Complejidad

El coste exponencial desaparece por completo. La construccion de las tablas de potencias requiere

$$\sum_{b=2}^{n+1} O\!\left(\frac{n}{b-1}\right)=O(n \log n),$$

y el barrido de bloques de cociente tiene el mismo orden porque

$$\sum_{l=1}^{n-1}\left\lfloor\frac{n-1}{l}\right\rfloor=O(n \log n).$$

Por tanto, el metodo rapido funciona en \(O(n \log n)\) tiempo y usa \(O(n \log n)\) memoria para las tablas de potencias. En la practica es eficiente porque cada bloque se reduce a unas pocas multiplicaciones modulares y consultas de tabla.

Notas y Referencias

  1. Pagina del problema: https://projecteuler.net/problem=977
  2. Funcion iterada: Wikipedia - Iterated function
  3. Grafo funcional: Wikipedia - Functional graph
  4. Puntos eventualmente periodicos: Wikipedia - Eventually periodic points
  5. Funciones piso y techo: Wikipedia - Floor and ceiling functions
  6. Serie geometrica: Wikipedia - Geometric series
  7. Aritmetica modular: Wikipedia - Modular arithmetic

问题概述

我们要统计长度为 \(n\) 的序列 \((a_1,a_2,\dots,a_n)\),其中每个 \(a_i\) 都属于 \(\{1,2,\dots,n\}\),并且满足

$$a_{i+1}=a_{a_i}\qquad (1 \le i < n).$$

Project Euler 977 要求在 \(n=10^6\) 时求出这样的序列个数,并对 \(10^9+7\) 取模。直接枚举的规模是 \(n^n\),完全不可行,所以必须先把这个递推对应的真实结构拆开,再把计数化成闭式求和。

数学方法

这道题最重要的观察是:它并不是在描述一串彼此独立的数,而是在描述同一个自映射反复迭代时,从点 1 出发得到的轨道。因此每个合法序列最终都会进入一个周期段,而正确的分类方式就是找出“从哪里开始进入周期”。

把递推改写成函数迭代

定义函数 \(f:\{1,\dots,n\}\to\{1,\dots,n\}\),令

$$f(i)=a_i.$$

那么原条件变成

$$a_{i+1}=a_{a_i}=f(a_i)=f(f(i)).$$

由此可以归纳得到

$$a_i=f^i(1)\qquad (i \ge 1).$$

也就是说,整个序列就是点 1 在函数 \(f\) 下的前向轨道。更进一步,对任何下标 \(m\) 都有

$$f^t(m)=a_{m+t-1}\qquad (t \ge 1).$$

这意味着如果两个位置取到相同的值,那么从这两个位置开始往后的整段尾部也完全相同。所以一旦周期出现,后面的行为就被彻底锁死了。

按“最早的周期后缀”分类

设 \(s\) 是最小的下标,使得后缀

$$a_s,a_{s+1},\dots,a_n$$

从自己的第一项开始就是周期的。记这个后缀的长度为

$$N=n-s+1,$$

记它的精确周期为 \(l\)。于是后缀中的位置按照模 \(l\) 的余类分成

$$S_u=\{\,s+u-1+ml : m \ge 0,\ s+u-1+ml \le n\,\}\qquad (1 \le u \le l).$$

若 \(N=ql+r\) 且 \(0 \le r < l\),那么前 \(r\) 个余类大小为 \(q+1\),其余 \(l-r\) 个余类大小为 \(q\)。这正是代码里反复出现的整除商与余数结构。

固定一个周期后缀以后,如何计数?

因为这个后缀具有周期 \(l\),同一个余类里的所有位置取值都相同。记第 \(u\) 类的公共值为 \(c_u\)。递推关系强迫它满足

$$c_u \in S_{u+1}\quad (1 \le u < l),\qquad c_l \in S_1.$$

换句话说,构造这个后缀时,需要对每个余类,从“下一个余类”里挑一个具体下标作为它的值。于是总的选择数就是

$$P(N,l)=\prod_{u=1}^{l}|S_u|=q^{\,l-r}(q+1)^r,$$

其中

$$q=\left\lfloor \frac{N}{l}\right\rfloor,\qquad r=N \bmod l.$$

这就是实现里所有幂次项的来源。

如何把前面的非周期部分接上去

如果周期从一开始就出现,也就是 \(s=1\)、\(N=n\),那么根本没有前缀需要处理。这一部分的贡献是

$$A(n)=\sum_{l=1}^{n}P(n,l).$$

现在设 \(s>1\)。为了让位置 \(s\) 的值满足递推,必须有

$$a_s=a_{a_{s-1}}=c_1.$$

哪些位置的值等于 \(c_1\) 呢?恰好就是集合 \(S_1\) 中的所有位置。因此 \(a_{s-1}\) 必须从 \(S_1\) 里选。但是其中有一个位置不能选,那就是 \(c_l \in S_1\)。如果选了它,那么同样的 \(l\)-周期结构就会从 \(s-1\) 开始,和 \(s\) 的最小性矛盾。

因此可行的连接方式数为

$$|S_1|-1=\left\lceil \frac{N}{l}\right\rceil - 1.$$

当 \(r=0\) 时它等于 \(q-1\),当 \(r>0\) 时它等于 \(q\)。一旦 \(a_{s-1}\) 确定,前面所有更早的位置都被迫形成最简单的前向链

$$a_i=i+1\qquad (1 \le i \le s-2),$$

因为只要更早的位置出现一次非平凡复制,周期段就会更早开始。

于是总计数公式是

$$F(n)=\sum_{l=1}^{n} P(n,l)+\sum_{N=1}^{n-1}\sum_{l=1}^{N}\left(\left\lceil \frac{N}{l}\right\rceil-1\right)P(N,l).$$

具体例子:\(n=7\),后缀长度 \(N=4\),周期 \(l=2\)

取 \(s=4\),也就是把周期段放在 \(a_4,a_5,a_6,a_7\)。此时

$$S_1=\{4,6\},\qquad S_2=\{5,7\}.$$

为了构造一个 2 周期后缀,我们需要选 \(c_1 \in S_2\)、\(c_2 \in S_1\)。所以共有

$$P(4,2)=2\cdot 2=4$$

种后缀。比如取 \(c_1=5\)、\(c_2=6\),就得到后缀

$$5,6,5,6.$$

这时 \(a_3\) 必须属于 \(S_1\),但又不能等于 \(c_2=6\),否则 2 周期会提前到位置 3 开始。所以 \(a_3=4\) 被唯一确定。再往前的位置只能是 \(a_1=2\)、\(a_2=3\),于是得到一个合法序列

$$ (2,3,4,5,6,5,6). $$

四种后缀都按同样方式工作,因此这个类别的总贡献正是

$$\left(\left\lceil \frac{4}{2}\right\rceil-1\right)P(4,2)=1\cdot 4=4.$$

把双重求和重排成商分块

上面的公式已经完整正确,较慢的校验版本就是直接按它计算。真正高效的实现会把第二项按整除商重排。对固定的 \(l\),写成

$$N=ql+r,\qquad 0 \le r < l,$$

其中 \(q=\lfloor N/l \rfloor\)。于是

$$P(N,l)=q^{\,l-r}(q+1)^r.$$

当 \(N\) 在区间 \(ql,ql+1,\dots,(q+1)l-1\) 内变化时,连接因子在 \(r=0\) 时是 \(q-1\),在 \(r \ge 1\) 时是 \(q\)。若

$$R=\min(l-1,n-1-ql),$$

则整块贡献为

$$B_{l,q}=(q-1)q^l+\sum_{r=1}^{R} q^{\,l+1-r}(q+1)^r.$$

内层和是一个望远镜式求和:

$$\sum_{r=1}^{R} q^{\,l+1-r}(q+1)^r=q^{\,l+1-R}(q+1)^{R+1}-q^{\,l+1}(q+1).$$

生产级实现正是用这个闭式,把大量重复工作折叠掉的。

代码如何工作

预计算幂表

C++、Python 和 Java 实现都会先建立一个紧凑的幂表。对某个底数 \(b\),真正可能用到的最大指数不超过

$$\left\lfloor \frac{n-1}{b-1}\right\rfloor + 1,$$

所以之后任意需要的 \(b^e \bmod (10^9+7)\) 都能用 \(O(1)\) 时间直接取出,而不必在主循环里反复做快速幂。

两层计数流程

第一层先算

$$A(n)=\sum_{l=1}^{n}P(n,l),$$

也就是周期从第一个位置就开始的所有情况。第二层再补上更晚才进入周期的那些类。这里没有直接三重枚举 \((N,l,r)\),而是对每个 \(l\) 按 \(q=\lfloor N/l \rfloor\) 的平台分组,再用块公式一次性求完整段的贡献。

校验与并行

三个版本都包含小规模校验:暴力枚举确认 \(n=7\) 时答案是 174,直接按原始双重求和计算确认 \(n=100\) 时答案是 305741269。之后才开始 \(n=10^6\) 的正式计算。C++ 与 Java 会把 \(l\) 的范围切分给多个线程,Python 则顺序执行同样的算术。

复杂度分析

指数级复杂度已经完全消失。幂表预处理的代价是

$$\sum_{b=2}^{n+1} O\!\left(\frac{n}{b-1}\right)=O(n \log n),$$

而商分块的总块数也满足

$$\sum_{l=1}^{n-1}\left\lfloor\frac{n-1}{l}\right\rfloor=O(n \log n).$$

因此快速算法的时间复杂度是 \(O(n \log n)\),空间复杂度也是 \(O(n \log n)\),主要由幂表占据。实际运行中它足够快,因为每个块都只需要少量模乘和数组读取。

注释与参考资料

  1. 题目页面: https://projecteuler.net/problem=977
  2. 迭代函数: Wikipedia - Iterated function
  3. 函数图: Wikipedia - Functional graph
  4. 最终进入周期的点: Wikipedia - Eventually periodic points
  5. 下取整与上取整: Wikipedia - Floor and ceiling functions
  6. 几何级数: Wikipedia - Geometric series
  7. 模运算: Wikipedia - Modular arithmetic

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

Нужно посчитать количество последовательностей \((a_1,a_2,\dots,a_n)\), где каждый \(a_i \in \{1,2,\dots,n\}\) и выполняется условие

$$a_{i+1}=a_{a_i}\qquad (1 \le i < n).$$

В задаче Project Euler 977 требуется найти это число для \(n=10^6\) по модулю \(10^9+7\). Полный перебор размера \(n^n\) невозможен, поэтому решение должно опираться не на перебор, а на точную структурную классификацию допустимых последовательностей.

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

Главная мысль состоит в том, что рекуррентное условие описывает итерации одной функции. Отсюда следует, что любая допустимая последовательность рано или поздно входит в периодический режим. Значит, естественно классифицировать объекты по тому, где именно этот режим начинается впервые.

Рекурсия как орбита одной функции

Положим \(f(i)=a_i\). Тогда

$$a_{i+1}=a_{a_i}=f(a_i)=f(f(i)).$$

Индукцией получаем

$$a_i=f^i(1)\qquad (i \ge 1).$$

То есть последовательность совпадает с прямой орбитой точки 1 при итерации \(f\). Более того, для любого индекса \(m\)

$$f^t(m)=a_{m+t-1}\qquad (t \ge 1).$$

Отсюда видно важное свойство: если в двух позициях стоит одно и то же значение, то и все последующие хвосты у этих позиций одинаковы. Поэтому после появления периодичности все оставшиеся элементы уже жестко определены.

Классификация по первому периодическому суффиксу

Пусть \(s\) — наименьший индекс, начиная с которого суффикс

$$a_s,a_{s+1},\dots,a_n$$

является периодическим сразу с первого члена. Обозначим его длину через

$$N=n-s+1,$$

а точный период через \(l\). Тогда позиции этого суффикса распадаются на классы по модулю \(l\):

$$S_u=\{\,s+u-1+ml : m \ge 0,\ s+u-1+ml \le n\,\}\qquad (1 \le u \le l).$$

Если \(N=ql+r\) и \(0 \le r < l\), то первые \(r\) классов имеют размер \(q+1\), а остальные \(l-r\) — размер \(q\).

Подсчет одного периодического суффикса

Так как суффикс \(l\)-периодичен, во всех позициях одного и того же класса \(S_u\) стоит одно и то же значение; обозначим его \(c_u\). Рекурсия заставляет выполняться правилу

$$c_u \in S_{u+1}\quad (1 \le u < l),\qquad c_l \in S_1.$$

Значит, чтобы построить такой суффикс, нужно для каждого класса выбрать один индекс из следующего класса. Число вариантов равно

$$P(N,l)=\prod_{u=1}^{l}|S_u|=q^{\,l-r}(q+1)^r,$$

где

$$q=\left\lfloor \frac{N}{l}\right\rfloor,\qquad r=N \bmod l.$$

Это базовая формула, на которой держится вся программа.

Как присоединяется непериодический префикс

Если периодический режим начинается сразу, то есть \(s=1\) и \(N=n\), никакого префикса нет. Тогда вклад равен

$$A(n)=\sum_{l=1}^{n}P(n,l).$$

Пусть теперь \(s>1\). Тогда предыдущий элемент должен удовлетворять условию

$$a_s=a_{a_{s-1}}=c_1.$$

Значение \(c_1\) встречается ровно в позициях из \(S_1\), поэтому \(a_{s-1}\) обязан лежать в \(S_1\). Но один выбор запрещен: элемент \(c_l \in S_1\). Если взять \(a_{s-1}=c_l\), то тот же самый период уже начнется с позиции \(s-1\), что противоречит минимальности \(s\).

Следовательно, число допустимых присоединений равно

$$|S_1|-1=\left\lceil \frac{N}{l}\right\rceil - 1.$$

Это \(q-1\), когда \(r=0\), и \(q\), когда \(r>0\). После выбора \(a_{s-1}\) все более ранние позиции вынужденно образуют простую цепочку

$$a_i=i+1\qquad (1 \le i \le s-2),$$

потому что любая более ранняя нетривиальная копия сдвинула бы начало периодичности еще левее.

Итак, окончательная формула имеет вид

$$F(n)=\sum_{l=1}^{n} P(n,l)+\sum_{N=1}^{n-1}\sum_{l=1}^{N}\left(\left\lceil \frac{N}{l}\right\rceil-1\right)P(N,l).$$

Разобранный пример: \(n=7\), \(N=4\), \(l=2\)

Положим \(s=4\). Тогда периодический суффикс — это \(a_4,a_5,a_6,a_7\), а классы равны

$$S_1=\{4,6\},\qquad S_2=\{5,7\}.$$

Чтобы получить 2-периодический суффикс, нужно выбрать \(c_1 \in S_2\) и \(c_2 \in S_1\). Поэтому

$$P(4,2)=2\cdot 2=4.$$

Например, выбор \(c_1=5\), \(c_2=6\) дает суффикс

$$5,6,5,6.$$

Теперь \(a_3\) должен принадлежать \(S_1\), но не может равняться \(c_2=6\), иначе период начался бы уже с позиции 3. Значит, \(a_3=4\), а раньше остается только цепочка \(a_1=2\), \(a_2=3\). Получаем, например, последовательность

$$ (2,3,4,5,6,5,6). $$

Все четыре выбора суффикса работают так же, поэтому вклад этого класса равен

$$\left(\left\lceil \frac{4}{2}\right\rceil-1\right)P(4,2)=1\cdot 4=4.$$

Перегруппировка суммы по блокам одинакового частного

Формула выше уже точна, и медленные проверочные версии считают именно ее. Быстрый алгоритм просто иначе организует вторую сумму. Для фиксированного \(l\) пишем

$$N=ql+r,\qquad 0 \le r < l,$$

где \(q=\lfloor N/l \rfloor\). Тогда

$$P(N,l)=q^{\,l-r}(q+1)^r.$$

Когда \(N\) пробегает блок \(ql,ql+1,\dots,(q+1)l-1\), коэффициент присоединения равен \(q-1\) при \(r=0\) и \(q\) при \(r \ge 1\). Если

$$R=\min(l-1,n-1-ql),$$

то вклад целого блока равен

$$B_{l,q}=(q-1)q^l+\sum_{r=1}^{R} q^{\,l+1-r}(q+1)^r.$$

Внутренняя сумма телескопируется:

$$\sum_{r=1}^{R} q^{\,l+1-r}(q+1)^r=q^{\,l+1-R}(q+1)^{R+1}-q^{\,l+1}(q+1).$$

Именно эта закрытая форма и используется в производительном решении.

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

Предварительное вычисление степеней

Реализации на C++, Python и Java сначала строят большую таблицу всех степеней, которые потом могут понадобиться. Для основания \(b\) достаточно хранить показатели до

$$\left\lfloor \frac{n-1}{b-1}\right\rfloor + 1,$$

поэтому любая величина \(b^e \bmod (10^9+7)\) позже извлекается за \(O(1)\) без повторных вызовов быстрого возведения в степень.

Два слоя суммирования

Сначала вычисляется

$$A(n)=\sum_{l=1}^{n}P(n,l),$$

то есть количество последовательностей, у которых периодическая часть начинается сразу. Затем добавляется поправка для более позднего старта периода. Вместо наивного перебора всех троек \((N,l,r)\) код проходит по плато частного \(q=\lfloor N/l \rfloor\) и суммирует каждый блок по закрытой формуле.

Проверки и параллелизм

Во всех трех версиях есть небольшие проверки корректности: полный перебор подтверждает значение 174 при \(n=7\), а прямой расчет не перегруппированной формулы подтверждает 305741269 при \(n=100\). После этого уже выполняется большой запуск для \(n=10^6\). В версиях на C++ и Java диапазон по \(l\) делится между несколькими потоками; версия на Python выполняет ту же математику последовательно.

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

Экспоненциальный рост полностью устранен. Подготовка таблиц степеней стоит

$$\sum_{b=2}^{n+1} O\!\left(\frac{n}{b-1}\right)=O(n \log n),$$

а число блоков в основной сумме имеет тот же порядок, потому что

$$\sum_{l=1}^{n-1}\left\lfloor\frac{n-1}{l}\right\rfloor=O(n \log n).$$

Значит, быстрый метод работает за \(O(n \log n)\) времени и требует \(O(n \log n)\) памяти для таблиц степеней. На практике он эффективен, потому что вклад каждого блока сводится к нескольким модульным умножениям и обращениям к массивам.

Сноски и ссылки

  1. Страница задачи: https://projecteuler.net/problem=977
  2. Итерируемая функция: Wikipedia - Iterated function
  3. Функциональный граф: Wikipedia - Functional graph
  4. В конечном итоге периодические точки: Wikipedia - Eventually periodic points
  5. Функции пола и потолка: Wikipedia - Floor and ceiling functions
  6. Геометрическая прогрессия: Wikipedia - Geometric series
  7. Модульная арифметика: Wikipedia - Modular arithmetic

ملخص المسألة

نريد عد جميع المتتاليات \((a_1,a_2,\dots,a_n)\) التي تحقق \(a_i \in \{1,2,\dots,n\}\) ولكل \(i\) مع \(1 \le i < n\) يكون

$$a_{i+1}=a_{a_i}.$$

في المسألة 977 المطلوب هو هذا العدد عندما \(n=10^6\)، مع اخذ الناتج بترديد \(10^9+7\). البحث المباشر حجمه \(n^n\)، لذلك لا بد من فهم البنية التي تفرضها العلاقة نفسها ثم تحويل العد الى صيغ مغلقة قابلة للحساب السريع.

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

الفكرة الاساسية هي ان هذه العلاقة لا تصف اعدادا مستقلة، بل تصف مدار نقطة واحدة تحت تكرار دالة واحدة. ولهذا السبب كل متتالية صحيحة تدخل في النهاية في جزء دوري. التصنيف الصحيح هو بحسب اول موضع يبدأ عنده هذا الجزء الدوري.

تحويل العلاقة الى مدار دالة

نعرف دالة \(f\) على \(\{1,\dots,n\}\) بواسطة

$$f(i)=a_i.$$

عندئذ تصبح العلاقة

$$a_{i+1}=a_{a_i}=f(a_i)=f(f(i)).$$

وبالاستقراء نحصل على

$$a_i=f^i(1)\qquad (i \ge 1).$$

اي ان المتتالية كلها هي المدار الامامي للنقطة 1 تحت \(f\). بل يوجد ايضا

$$f^t(m)=a_{m+t-1}\qquad (t \ge 1)$$

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

التصنيف بحسب اول لاحقة دورية

لنفرض ان \(s\) هو اصغر فهرس تكون منه اللاحقة

$$a_s,a_{s+1},\dots,a_n$$

دورية ابتداء من حدها الاول. نكتب

$$N=n-s+1$$

لطول هذه اللاحقة، و\(l\) لدورتها الدقيقة. عندئذ تتجزا مواضع اللاحقة الى فئات توافقية بترديد \(l\):

$$S_u=\{\,s+u-1+ml : m \ge 0,\ s+u-1+ml \le n\,\}\qquad (1 \le u \le l).$$

اذا كتبنا \(N=ql+r\) مع \(0 \le r < l\)، فان اول \(r\) فئات حجمها \(q+1\)، وبقية \(l-r\) فئات حجمها \(q\).

عد لاحقة دورية واحدة

لان اللاحقة دورية بطول \(l\)، فان جميع المواضع داخل الفئة نفسها تحمل القيمة نفسها. نسمي قيمة الفئة \(S_u\) بالرمز \(c_u\). العلاقة تعطي الشرط

$$c_u \in S_{u+1}\quad (1 \le u < l),\qquad c_l \in S_1.$$

اي ان بناء اللاحقة يكافئ اختيار عنصر واحد من الفئة التالية لكل \(u\). ومن هنا يكون عدد الاختيارات

$$P(N,l)=\prod_{u=1}^{l}|S_u|=q^{\,l-r}(q+1)^r,$$

حيث

$$q=\left\lfloor \frac{N}{l}\right\rfloor,\qquad r=N \bmod l.$$

هذا هو الجسم الرياضي الحقيقي الذي تستعمله جميع الصيغ في الحل.

كيف نصل اللاحقة الدورية بالمقدمة غير الدورية

اذا بدا الجزء الدوري من الموضع الاول، اي عندما \(s=1\) و\(N=n\)، فلا توجد مقدمة اضافية. مساهمة هذه الحالة هي

$$A(n)=\sum_{l=1}^{n}P(n,l).$$

اما اذا كان \(s>1\)، فيجب ان يحقق الحد السابق

$$a_s=a_{a_{s-1}}=c_1.$$

المواضع التي تساوي قيمتها \(c_1\) هي بالضبط عناصر \(S_1\)، ولذلك يجب ان يكون \(a_{s-1}\) عنصرا من \(S_1\). لكن هناك اختيارا واحدا ممنوعا، وهو \(c_l \in S_1\). لو اخذنا \(a_{s-1}=c_l\)، فان النمط الدوري نفسه سيبدا فعليا من الموضع \(s-1\)، وهذا يناقض اصغرية \(s\).

اذن عدد طرق الوصل الصحيحة هو

$$|S_1|-1=\left\lceil \frac{N}{l}\right\rceil - 1.$$

وهذا يساوي \(q-1\) اذا كان \(r=0\)، ويساوي \(q\) اذا كان \(r>0\). وبعد تثبيت \(a_{s-1}\)، تصبح جميع المواضع السابقة مجبرة على السلسلة البسيطة

$$a_i=i+1\qquad (1 \le i \le s-2),$$

لان اي نسخ غير بدهي في موضع ابكر سيجعل الجزء الدوري يبدا قبل ذلك.

وبالتالي تكون الصيغة النهائية

$$F(n)=\sum_{l=1}^{n} P(n,l)+\sum_{N=1}^{n-1}\sum_{l=1}^{N}\left(\left\lceil \frac{N}{l}\right\rceil-1\right)P(N,l).$$

مثال عملي: \(n=7\)، وطول اللاحقة \(N=4\)، والدورة \(l=2\)

لنأخذ \(s=4\). عندئذ تكون اللاحقة الدورية هي \(a_4,a_5,a_6,a_7\)، والفئتان هما

$$S_1=\{4,6\},\qquad S_2=\{5,7\}.$$

لبناء لاحقة دورية بطول 2 نختار \(c_1 \in S_2\) و\(c_2 \in S_1\)، ولذلك يوجد

$$P(4,2)=2\cdot 2=4$$

اختيارات. فمثلا اذا اخترنا \(c_1=5\) و\(c_2=6\) حصلنا على اللاحقة

$$5,6,5,6.$$

الان يجب ان يكون \(a_3\) عنصرا من \(S_1\)، لكنه لا يجوز ان يساوي \(c_2=6\)، لان ذلك يجعل الدورية تبدا من الموضع 3. لذلك يثبت \(a_3=4\)، ثم تصبح المقدمة مجبرة على \(a_1=2\) و\(a_2=3\). فنحصل مثلا على المتتالية

$$ (2,3,4,5,6,5,6). $$

وكل اختيار من الاختيارات الاربع يعمل بالطريقة نفسها، لذا تكون مساهمة هذا الصنف

$$\left(\left\lceil \frac{4}{2}\right\rceil-1\right)P(4,2)=1\cdot 4=4.$$

اعادة ترتيب المجموع بحسب كتل القسمة الصحيحة

الصيغة السابقة صحيحة بالكامل، والنسخ البطيئة تستعملها مباشرة. اما النسخة السريعة فتعيد تنظيم الحد الثاني. مع \(l\) ثابت نكتب

$$N=ql+r,\qquad 0 \le r < l,$$

حيث \(q=\lfloor N/l \rfloor\). عندئذ

$$P(N,l)=q^{\,l-r}(q+1)^r.$$

عندما يتحرك \(N\) داخل الكتلة \(ql,ql+1,\dots,(q+1)l-1\)، يكون معامل الوصل \(q-1\) عند \(r=0\)، و\(q\) عند \(r \ge 1\). اذا عرفنا

$$R=\min(l-1,n-1-ql),$$

فان مساهمة الكتلة كلها تصبح

$$B_{l,q}=(q-1)q^l+\sum_{r=1}^{R} q^{\,l+1-r}(q+1)^r.$$

والمجموع الداخلي يتبسط على شكل متسلسل متلاش:

$$\sum_{r=1}^{R} q^{\,l+1-r}(q+1)^r=q^{\,l+1-R}(q+1)^{R+1}-q^{\,l+1}(q+1).$$

وهذه هي الصيغة المغلقة التي تعتمد عليها التطبيقات السريعة.

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

التحضير المسبق لجدول القوى

تقوم تطبيقات C++ وPython وJava ببناء جدول مضغوط لكل القوى المطلوبة مسبقا. بالنسبة للاساس \(b\)، اكبر اس مطلوب لا يتجاوز

$$\left\lfloor \frac{n-1}{b-1}\right\rfloor + 1,$$

ولذلك يمكن قراءة اي قيمة \(b^e \bmod (10^9+7)\) في زمن \(O(1)\) من الجدول، من غير اعادة حسابها داخل الحلقات الرئيسية.

طبقتان من العد

اولا تحسب التطبيقات

$$A(n)=\sum_{l=1}^{n}P(n,l),$$

وهو عدد المتتاليات التي يبدا فيها الجزء الدوري من الموضع الاول. ثم تضيف التصحيح الخاص بالحالات التي تبدا فيها الدورية لاحقا. وبدلا من الدوران الساذج على جميع الثلاثيات \((N,l,r)\)، يجري التجميع بحسب مستويات القسمة \(q=\lfloor N/l \rfloor\)، ثم تجمع كل كتلة بصيغتها المغلقة.

التحقق والتنفيذ المتوازي

تحتوي جميع النسخ على اختبارات صغيرة: العد الشامل يؤكد ان القيمة عند \(n=7\) هي 174، كما ان التقييم المباشر للصيغة غير المعاد تجميعها يؤكد ان القيمة عند \(n=100\) هي 305741269. بعد ذلك فقط تبدأ العملية الكبيرة عند \(n=10^6\). نسختا C++ وJava تقسمان مجال \(l\) على عدة خيوط، بينما تنفذ نسخة Python الحساب نفسه بصورة تسلسلية.

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

اختفى التعقيد الاسّي تماما. بناء جداول القوى يكلف

$$\sum_{b=2}^{n+1} O\!\left(\frac{n}{b-1}\right)=O(n \log n),$$

كما ان عدد كتل القسمة في المجموع الرئيسي له المرتبة نفسها، لان

$$\sum_{l=1}^{n-1}\left\lfloor\frac{n-1}{l}\right\rfloor=O(n \log n).$$

اذن تعمل الطريقة السريعة في زمن \(O(n \log n)\) وتستهلك ذاكرة \(O(n \log n)\)، ويذهب معظمها الى جداول القوى. وعمليا تكون الطريقة فعالة لان مساهمة كل كتلة تختزل الى عدد صغير من الضربات المعيارية وعمليات القراءة من المصفوفات.

حواش ومراجع

  1. صفحة المسألة: https://projecteuler.net/problem=977
  2. الدالة التكرارية: Wikipedia - Iterated function
  3. الرسم الوظيفي: Wikipedia - Functional graph
  4. النقاط التي تصبح دورية في النهاية: Wikipedia - Eventually periodic points
  5. دالتي الجزء الصحيح والسقف: Wikipedia - Floor and ceiling functions
  6. المتسلسلة الهندسية: Wikipedia - Geometric series
  7. الحساب المعياري: Wikipedia - Modular arithmetic