Problem Summary

Let

$$T(x)=x-\frac{1}{x},\qquad x\neq 0.$$

For each exact period \(m\) with \(2\le m\le p\), the problem considers real periodic sequences generated by repeated application of \(T\). If one period of the sequence is \((x_0,x_1,\dots,x_{m-1})\), its range is

$$\max_{0\le i\lt m} x_i-\min_{0\le i\lt m} x_i.$$

The target quantity \(S(p)\) is the sum of these ranges over all primitive periodic sequences of lengths up to \(p\), counting different cyclic starting positions separately. Directly searching real periodic sequences is not practical, so the implementation re-encodes each orbit by inverse branches and solves one fixed-point problem per primitive orbit.

Mathematical Approach

The solver converts the periodic-sequence problem into a combinatorial enumeration of primitive binary necklaces, followed by a numerical fixed-point solve for each necklace.

Step 1: Write the Two Inverse Branches

If \(y=T(x)\), then \(x\) satisfies

$$x^2-yx-1=0.$$

So every value \(y\) has two real preimages, namely

$$g_{\pm}(y)=\frac{y\pm\sqrt{y^2+4}}{2}.$$

A periodic orbit can therefore be reconstructed backward by choosing, at each step, either the \(+\) branch or the \(-\) branch. A binary word \(\sigma=(s_0,s_1,\dots,s_{m-1})\) with \(s_i\in\{-1,+1\}\) determines the composition

$$G_\sigma=g_{s_{m-1}}\circ g_{s_{m-2}}\circ \cdots \circ g_{s_0}.$$

A period-\(m\) orbit appears when some starting value \(x_0\) satisfies the fixed-point equation

$$x_0=G_\sigma(x_0).$$

Step 2: Why Primitive Necklaces Are the Right Objects

Two binary words that differ only by a cyclic rotation describe the same orbit, just started at different points of the cycle. So the natural combinatorial object is not a raw binary word but a binary necklace.

If a word is a repetition of a shorter block, then the corresponding orbit has a smaller fundamental period. Therefore exact period \(m\) is represented by primitive binary necklaces of length \(m\). Their number is

$$P_m=\frac{1}{m}\sum_{d\mid m}\mu(d)\,2^{m/d},$$

where \(\mu\) is the Möbius function. This count is also used as a structural check before the numerical stage begins.

Step 3: Show That Each Branch Word Has a Unique Real Fixed Point

For \(s\in\{-1,+1\}\), differentiating one inverse branch gives

$$g_s'(y)=\frac{1}{2}\left(1+s\frac{y}{\sqrt{y^2+4}}\right).$$

Because

$$-1\lt\frac{y}{\sqrt{y^2+4}}\lt 1,$$

we obtain

$$0\lt g_s'(y)\lt 1.$$

So each inverse branch is an increasing contraction on \(\mathbb{R}\), and any finite composition \(G_\sigma\) is again a contraction. By the contraction mapping principle, every branch word \(\sigma\) has a unique real fixed point. That is why a fixed-point iteration converges reliably and why each primitive necklace corresponds to exactly one real primitive orbit.

Step 4: Recover the Orbit and Its Range

Once the fixed point \(x_0\) of \(G_\sigma\) is known, the orbit itself is obtained by applying the inverse branches in the chosen order:

$$x_{i+1}=g_{s_i}(x_i)\qquad (0\le i\le m-2).$$

The final branch returns from \(x_{m-1}\) to \(x_0\), closing the cycle. The range attached to \(\sigma\) is

$$R(\sigma)=\max_{0\le i\lt m} x_i-\min_{0\le i\lt m} x_i.$$

A primitive necklace represents one orbit up to rotation, but the original sum counts all \(m\) cyclic starting positions as distinct periodic sequences. Therefore one primitive necklace contributes

$$m\,R(\sigma),$$

and the full target becomes

$$S(p)=\sum_{m=2}^{p}\ \sum_{\sigma\in\mathcal{N}_m^{\mathrm{prim}}} m\,R(\sigma),$$

where \(\mathcal{N}_m^{\mathrm{prim}}\) denotes the primitive binary necklaces of length \(m\).

Step 5: Worked Example for \(m=2\)

There is exactly one primitive binary necklace of length \(2\): one \(+\) branch and one \(-\) branch, up to rotation. Its orbit is a 2-cycle \(\{x,-x\}\). Indeed, the condition \(T(x)=-x\) gives

$$x-\frac{1}{x}=-x,$$

so

$$2x=\frac{1}{x},\qquad x^2=\frac{1}{2}.$$

Thus the orbit points are

$$x=\pm\frac{1}{\sqrt{2}}.$$

The range is therefore

$$\frac{1}{\sqrt{2}}-\left(-\frac{1}{\sqrt{2}}\right)=\sqrt{2},$$

and because there are two cyclic starting positions, the total contribution is

$$2\sqrt{2}=2.8284271247\ldots$$

This matches the smallest validation value used by the implementation.

How the Code Works

The C++, Python, and Java implementations follow the same mathematical pipeline. For each \(m=2,3,\dots,p\), they generate all primitive binary necklaces of length \(m\) using a recursive necklace generator. Before any real-valued computation is performed, the generated count is compared with the Möbius-inversion formula above to confirm that only primitive representatives are present.

For each representative, the implementation starts from \(0\) and repeatedly applies the full inverse-branch composition \(G_\sigma\). This provides a stable initial approximation to the fixed point. It then refines that approximation with Newton's method applied to

$$\Phi(x)=G_\sigma(x)-x.$$

During this refinement, the derivative of the composition is accumulated by the chain rule while traversing the chosen branches, so the Newton correction uses the exact derivative of the composed inverse map.

After the fixed point is obtained, the implementation unfolds the orbit one branch at a time, tracks the minimum and maximum values encountered, computes \(R(\sigma)\), and adds \(mR(\sigma)\) to the running total. Compensated summation is used to reduce floating-point cancellation. The C++ and Python implementations distribute independent necklace blocks across workers, while the Java implementation delegates the numerical evaluation to the same underlying computation and returns the resulting decimal value.

Complexity Analysis

Let

$$P_m=\frac{1}{m}\sum_{d\mid m}\mu(d)\,2^{m/d}$$

be the number of primitive binary necklaces of length \(m\). Processing one necklace requires \(O(mI)\) arithmetic operations, where \(I\) is the small number of fixed-point and Newton iterations. Hence the total running time is

$$O\!\left(\sum_{m=2}^{p} P_m\,m\,I\right).$$

Since \(P_m\sim 2^m/m\), the growth is essentially exponential in \(p\), with the largest periods dominating the cost. Memory usage for a single length is \(O(P_m)\) to store the necklace representatives, plus \(O(1)\) orbit data and a small amount of worker-local accumulation state.

Footnotes and References

  1. Problem page: Project Euler 729
  2. Necklace counting: Wikipedia - Necklace (combinatorics)
  3. Möbius inversion: Wikipedia - Möbius inversion formula
  4. Contraction mapping principle: Wikipedia - Banach fixed-point theorem
  5. Newton's method: Wikipedia - Newton's method
  6. Compensated summation: Wikipedia - Kahan summation algorithm

Problemzusammenfassung

Sei

$$T(x)=x-\frac{1}{x},\qquad x\neq 0.$$

Für jede exakte Periode \(m\) mit \(2\le m\le p\) betrachtet das Problem reelle periodische Folgen unter wiederholter Anwendung von \(T\). Ist eine Periode der Folge \((x_0,x_1,\dots,x_{m-1})\), dann ist ihre Spannweite

$$\max_{0\le i\lt m} x_i-\min_{0\le i\lt m} x_i.$$

Gesucht ist \(S(p)\), also die Summe dieser Spannweiten über alle primitiven periodischen Folgen der Längen bis \(p\), wobei verschiedene zyklische Startpositionen getrennt gezählt werden. Eine direkte Suche im Reellen wäre unpraktisch; deshalb kodiert die Implementierung jede Bahn über inverse Zweige und löst pro primitiver Bahn genau ein Fixpunktproblem.

Mathematischer Ansatz

Die Lösung ersetzt das direkte Durchsuchen periodischer Folgen durch eine Aufzählung primitiver binärer Halsketten und eine numerische Fixpunktberechnung für jede solche Halskette.

Schritt 1: Die beiden inversen Zweige aufschreiben

Aus \(y=T(x)\) folgt

$$x^2-yx-1=0.$$

Damit besitzt jeder Wert \(y\) genau zwei reelle Urbilder:

$$g_{\pm}(y)=\frac{y\pm\sqrt{y^2+4}}{2}.$$

Eine periodische Bahn kann also rückwärts rekonstruiert werden, indem man in jedem Schritt entweder den \(+\)- oder den \(-\)-Zweig wählt. Ein binäres Wort \(\sigma=(s_0,s_1,\dots,s_{m-1})\) mit \(s_i\in\{-1,+1\}\) definiert die Komposition

$$G_\sigma=g_{s_{m-1}}\circ g_{s_{m-2}}\circ \cdots \circ g_{s_0}.$$

Eine Bahn der Länge \(m\) entsteht genau dann, wenn ein Startwert \(x_0\) die Fixpunktgleichung

$$x_0=G_\sigma(x_0)$$

erfüllt.

Schritt 2: Warum primitive Halsketten die richtigen Objekte sind

Zwei binäre Wörter, die sich nur durch zyklische Rotation unterscheiden, beschreiben dieselbe Bahn, nur mit anderem Startpunkt. Deshalb sind binäre Halsketten die natürlichen Repräsentanten.

Ist ein Wort eine Wiederholung eines kürzeren Blocks, dann besitzt die zugehörige Bahn eine kleinere Grundperiode. Exakte Periode \(m\) wird daher durch primitive binäre Halsketten der Länge \(m\) beschrieben. Ihre Anzahl ist

$$P_m=\frac{1}{m}\sum_{d\mid m}\mu(d)\,2^{m/d},$$

wobei \(\mu\) die Möbius-Funktion bezeichnet. Genau diese Zahl dient auch als strukturelle Kontrolle vor der numerischen Berechnung.

Schritt 3: Jede Zweigfolge hat genau einen reellen Fixpunkt

Für \(s\in\{-1,+1\}\) gilt nach Differentiation

$$g_s'(y)=\frac{1}{2}\left(1+s\frac{y}{\sqrt{y^2+4}}\right).$$

Wegen

$$-1\lt\frac{y}{\sqrt{y^2+4}}\lt 1$$

folgt

$$0\lt g_s'(y)\lt 1.$$

Jeder inverse Zweig ist also eine wachsende Kontraktion auf \(\mathbb{R}\), und damit ist auch jede endliche Komposition \(G_\sigma\) eine Kontraktion. Nach dem Banachschen Fixpunktsatz besitzt daher jedes Zweigwort \(\sigma\) genau einen reellen Fixpunkt. Das erklärt, warum die Fixpunktiteration zuverlässig konvergiert und warum jede primitive Halskette genau einer primitiven reellen Bahn entspricht.

Schritt 4: Die Bahn und ihre Spannweite rekonstruieren

Ist der Fixpunkt \(x_0\) von \(G_\sigma\) bekannt, erhält man die Bahn durch sukzessives Anwenden der inversen Zweige:

$$x_{i+1}=g_{s_i}(x_i)\qquad (0\le i\le m-2).$$

Der letzte Zweig führt von \(x_{m-1}\) zurück nach \(x_0\). Die zu \(\sigma\) gehörige Spannweite ist

$$R(\sigma)=\max_{0\le i\lt m} x_i-\min_{0\le i\lt m} x_i.$$

Eine primitive Halskette repräsentiert eine Bahn nur bis auf Rotation, aber in der ursprünglichen Summe zählen alle \(m\) zyklischen Startpositionen separat. Deshalb trägt eine primitive Halskette

$$m\,R(\sigma)$$

bei, und insgesamt gilt

$$S(p)=\sum_{m=2}^{p}\ \sum_{\sigma\in\mathcal{N}_m^{\mathrm{prim}}} m\,R(\sigma).$$

Schritt 5: Durchgerechnetes Beispiel für \(m=2\)

Es gibt genau eine primitive binäre Halskette der Länge \(2\): einen \(+\)- und einen \(-\)-Zweig, bis auf Rotation. Die zugehörige Bahn ist ein 2-Zyklus \(\{x,-x\}\). Aus \(T(x)=-x\) folgt

$$x-\frac{1}{x}=-x,$$

also

$$2x=\frac{1}{x},\qquad x^2=\frac{1}{2}.$$

Damit sind die Bahnpunkte

$$x=\pm\frac{1}{\sqrt{2}}.$$

Die Spannweite ist also

$$\frac{1}{\sqrt{2}}-\left(-\frac{1}{\sqrt{2}}\right)=\sqrt{2},$$

und wegen der zwei zyklischen Startpositionen beträgt der Gesamtbeitrag

$$2\sqrt{2}=2.8284271247\ldots$$

Genau dieser Wert erscheint als kleinster Validierungspunkt der Implementierung.

Wie der Code arbeitet

Die C++-, Python- und Java-Implementierungen folgen demselben mathematischen Ablauf. Für jedes \(m=2,3,\dots,p\) erzeugen sie zunächst alle primitiven binären Halsketten der Länge \(m\) mit einem rekursiven Generator. Bevor mit reellen Zahlen gearbeitet wird, wird die erzeugte Anzahl mit der obigen Möbius-Formel verglichen, damit nur primitive Repräsentanten weiterverarbeitet werden.

Für jeden Repräsentanten startet die Berechnung bei \(0\) und wendet die vollständige inverse Zweigkomposition \(G_\sigma\) wiederholt an. Das liefert eine stabile erste Näherung des Fixpunkts. Danach verfeinert die Implementierung diese Näherung mit dem Newton-Verfahren für

$$\Phi(x)=G_\sigma(x)-x.$$

Dabei wird die Ableitung der gesamten Komposition über die Kettenregel während des Durchlaufs durch die Zweige mitgeführt, sodass die Newton-Korrektur die exakte Ableitung der zusammengesetzten inversen Abbildung verwendet.

Sobald der Fixpunkt feststeht, wird die Bahn Zweig für Zweig entfaltet, das Minimum und Maximum werden verfolgt, \(R(\sigma)\) wird berechnet, und \(mR(\sigma)\) wird zur laufenden Summe addiert. Kompensierte Summation verringert Rundungsfehler. Die C++- und Python-Implementierungen verteilen unabhängige Halskettenblöcke auf mehrere Arbeiter, während die Java-Implementierung dieselbe zugrunde liegende numerische Berechnung aufruft und den resultierenden Dezimalwert ausgibt.

Komplexitätsanalyse

Sei

$$P_m=\frac{1}{m}\sum_{d\mid m}\mu(d)\,2^{m/d}$$

die Anzahl primitiver binärer Halsketten der Länge \(m\). Die Verarbeitung einer einzelnen Halskette benötigt \(O(mI)\) arithmetische Operationen, wobei \(I\) die kleine Anzahl von Fixpunkt- und Newton-Iterationen ist. Somit ergibt sich insgesamt

$$O\!\left(\sum_{m=2}^{p} P_m\,m\,I\right).$$

Da \(P_m\sim 2^m/m\) gilt, wächst der Aufwand im Wesentlichen exponentiell in \(p\), und die größten Perioden dominieren die Laufzeit. Der Speicherbedarf für eine feste Länge ist \(O(P_m)\) zum Speichern der Halskettenrepräsentanten, zuzüglich \(O(1)\) für Bahndaten und etwas arbeiterlokalen Akkumulatorzustand.

Fußnoten und Referenzen

  1. Problemseite: Project Euler 729
  2. Halskettenzählung: Wikipedia - Necklace (combinatorics)
  3. Möbius-Inversion: Wikipedia - Möbius inversion formula
  4. Kontraktionsprinzip: Wikipedia - Banach fixed-point theorem
  5. Newton-Verfahren: Wikipedia - Newton's method
  6. Kompensierte Summation: Wikipedia - Kahan summation algorithm

Problem Özeti

Şunu tanımlayalım:

$$T(x)=x-\frac{1}{x},\qquad x\neq 0.$$

Problem, \(2\le m\le p\) aralığındaki her tam periyot için \(T\) dönüşümü altında oluşan reel periyodik dizileri inceler. Bir periyot \((x_0,x_1,\dots,x_{m-1})\) ise aralık

$$\max_{0\le i\lt m} x_i-\min_{0\le i\lt m} x_i$$

şeklinde tanımlanır. Amaç, uzunluğu en fazla \(p\) olan tüm ilkel periyodik diziler için bu aralıkları toplamak ve farklı döngüsel başlangıç noktalarını ayrı diziler olarak saymaktır. Reel sayılar üzerinde doğrudan arama yapmak verimli değildir; bu yüzden uygulama her yörüngeyi ters dallarla kodlayıp her ilkel yörünge için tek bir sabit nokta problemi çözer.

Matematiksel Yaklaşım

Çözüm, periyodik dizileri doğrudan üretmek yerine ilkel ikili boyunlukları sayar ve her boyunluk için sayısal bir sabit nokta hesabı yapar.

Adım 1: İki Ters Dalı Elde Et

\(y=T(x)\) ise

$$x^2-yx-1=0$$

olur. Dolayısıyla her \(y\) değeri için iki reel öncül vardır:

$$g_{\pm}(y)=\frac{y\pm\sqrt{y^2+4}}{2}.$$

Böylece bir periyodik yörünge geriye doğru kurulurken her adımda \(+\) veya \(-\) dalı seçilir. \(s_i\in\{-1,+1\}\) olacak şekilde \(\sigma=(s_0,s_1,\dots,s_{m-1})\) ikili kelimesi

$$G_\sigma=g_{s_{m-1}}\circ g_{s_{m-2}}\circ \cdots \circ g_{s_0}$$

bileşimini belirler. Periyot \(m\) olan bir yörünge,

$$x_0=G_\sigma(x_0)$$

sabit nokta denklemini sağlayan bir \(x_0\) ile elde edilir.

Adım 2: Neden İlkel Boyunluklar Kullanılır?

İki ikili kelime yalnızca döngüsel kaydırma ile birbirine dönüşüyorsa, aslında aynı yörüngeyi farklı başlangıç noktalarından yazıyorlardır. Bu nedenle ham kelimeler değil, döngüsel sınıflar yani boyunluklar doğal nesnedir.

Eğer bir kelime daha kısa bir bloğun tekrarıysa, ortaya çıkan yörüngenin gerçek temel periyodu daha küçüktür. O halde tam periyot \(m\) için doğru temsilciler, uzunluğu \(m\) olan ilkel ikili boyunluklardır. Bunların sayısı

$$P_m=\frac{1}{m}\sum_{d\mid m}\mu(d)\,2^{m/d}$$

şeklindedir; burada \(\mu\) Möbius fonksiyonudur. Uygulama sayısal hesaplamaya geçmeden önce üretilen temsilci sayısını bu formülle karşılaştırır.

Adım 3: Her Dal Kelimesinin Tek Bir Reel Sabit Noktası Vardır

\(s\in\{-1,+1\}\) için türev

$$g_s'(y)=\frac{1}{2}\left(1+s\frac{y}{\sqrt{y^2+4}}\right)$$

olur. Ayrıca

$$-1\lt\frac{y}{\sqrt{y^2+4}}\lt 1$$

olduğundan

$$0\lt g_s'(y)\lt 1$$

elde edilir. Yani her ters dal \(\mathbb{R}\) üzerinde artan bir daralmadır; dolayısıyla her sonlu bileşim \(G_\sigma\) da daralmadır. Banach sabit nokta teoremine göre her dal kelimesinin tek bir reel sabit noktası vardır. Sabit nokta iterasyonunun güvenilir çalışmasının ve her ilkel boyunluğun tek bir reel ilkel yörüngeye karşılık gelmesinin nedeni budur.

Adım 4: Yörüngeyi ve Aralığını Yeniden Kur

\(G_\sigma\)'nın sabit noktası \(x_0\) bulunduğunda, yörüngenin diğer noktaları seçilen ters dallar sırayla uygulanarak elde edilir:

$$x_{i+1}=g_{s_i}(x_i)\qquad (0\le i\le m-2).$$

Son dal tekrar \(x_0\)'a dönerek çevrimi kapatır. \(\sigma\)'ya ait aralık

$$R(\sigma)=\max_{0\le i\lt m} x_i-\min_{0\le i\lt m} x_i$$

şeklindedir. İlkel bir boyunluk yörüngeyi yalnızca dönme eşdeğerliği altında temsil eder; fakat problemde aynı döngünün \(m\) farklı başlangıç noktası ayrı periyodik diziler olarak sayılır. Bu yüzden tek bir ilkel boyunluğun katkısı

$$m\,R(\sigma)$$

olur ve toplam

$$S(p)=\sum_{m=2}^{p}\ \sum_{\sigma\in\mathcal{N}_m^{\mathrm{prim}}} m\,R(\sigma)$$

biçimine dönüşür.

Adım 5: \(m=2\) İçin Çalışılmış Örnek

Uzunluğu \(2\) olan yalnızca bir ilkel ikili boyunluk vardır: dönmeye kadar bir \(+\) ve bir \(-\) dalı. Buna karşılık gelen yörünge bir 2-çevrim \(\{x,-x\}\) olur. Gerçekten de \(T(x)=-x\) yazarsak

$$x-\frac{1}{x}=-x$$

ve buradan

$$2x=\frac{1}{x},\qquad x^2=\frac{1}{2}$$

çıkar. Demek ki yörünge noktaları

$$x=\pm\frac{1}{\sqrt{2}}$$

olur. Aralık

$$\frac{1}{\sqrt{2}}-\left(-\frac{1}{\sqrt{2}}\right)=\sqrt{2}$$

ve iki döngüsel başlangıç noktası olduğu için toplam katkı

$$2\sqrt{2}=2.8284271247\ldots$$

şeklindedir. Bu değer, uygulamanın en küçük doğrulama sonucuyla birebir örtüşür.

Kod Nasıl Çalışır?

C++, Python ve Java uygulamaları aynı matematiksel akışı izler. Her \(m=2,3,\dots,p\) için önce uzunluğu \(m\) olan tüm ilkel ikili boyunluklar özyinelemeli bir üreteçle oluşturulur. Reel sayı hesabına başlamadan önce üretilen temsilci sayısı yukarıdaki Möbius formülüyle karşılaştırılarak yalnızca ilkel temsilcilerin kaldığı doğrulanır.

Her temsilci için hesap \(0\)'dan başlar ve tam ters-dal bileşimi \(G_\sigma\) art arda uygulanır. Bu işlem sabit nokta için kararlı bir ilk yaklaşım verir. Ardından uygulama

$$\Phi(x)=G_\sigma(x)-x$$

fonksiyonuna Newton yöntemi uygular. Bu sırada bileşimin türevi, seçilen dallar boyunca zincir kuralı kullanılarak biriktirilir; böylece Newton düzeltmesi birleşik ters dönüşümün gerçek türeviyle yapılır.

Sabit nokta elde edildikten sonra yörünge dal dal açılır, görülen minimum ve maksimum değerler izlenir, \(R(\sigma)\) hesaplanır ve \(mR(\sigma)\) toplam değere eklenir. Kayan nokta hatasını azaltmak için telafili toplama kullanılır. C++ ve Python sürümleri bağımsız boyunluk bloklarını çalışanlara dağıtır; Java sürümü ise aynı temel sayısal hesabı çağırıp son ondalık çıktıyı döndürür.

Karmaşıklık Analizi

Uzunluğu \(m\) olan ilkel ikili boyunluk sayısı

$$P_m=\frac{1}{m}\sum_{d\mid m}\mu(d)\,2^{m/d}$$

olsun. Tek bir boyunluğu işlemek \(O(mI)\) aritmetik işlem gerektirir; burada \(I\) sabit nokta ve Newton iterasyonlarının küçük sayısıdır. Toplam çalışma süresi

$$O\!\left(\sum_{m=2}^{p} P_m\,m\,I\right)$$

şeklindedir. \(P_m\sim 2^m/m\) olduğundan büyüme özünde \(p\)'ye göre üsseldir ve maliyetin büyük kısmı en uzun periyotlardan gelir. Tek bir uzunluk için bellek kullanımı, boyunluk temsilcilerini saklamak üzere \(O(P_m)\), ayrıca \(O(1)\) yörünge verisi ve küçük miktarda çalışan-yerel birikim durumudur.

Dipnotlar ve Kaynaklar

  1. Problem sayfası: Project Euler 729
  2. Boyunluk sayımı: Wikipedia - Necklace (combinatorics)
  3. Möbius tersleme: Wikipedia - Möbius inversion formula
  4. Daralma dönüşümü ilkesi: Wikipedia - Banach fixed-point theorem
  5. Newton yöntemi: Wikipedia - Newton's method
  6. Telafili toplama: Wikipedia - Kahan summation algorithm

Resumen del Problema

Sea

$$T(x)=x-\frac{1}{x},\qquad x\neq 0.$$

Para cada periodo exacto \(m\) con \(2\le m\le p\), el problema considera las sucesiones periódicas reales generadas por iterar \(T\). Si un periodo es \((x_0,x_1,\dots,x_{m-1})\), su rango es

$$\max_{0\le i\lt m} x_i-\min_{0\le i\lt m} x_i.$$

La cantidad buscada \(S(p)\) suma esos rangos sobre todas las sucesiones periódicas primitivas de longitudes hasta \(p\), contando por separado las distintas posiciones cíclicas de inicio. Buscar directamente entre todas las sucesiones reales no es viable, así que la implementación recodifica cada órbita mediante ramas inversas y resuelve un problema de punto fijo por cada órbita primitiva.

Enfoque Matemático

La idea central es sustituir la enumeración directa de sucesiones periódicas por una enumeración de collares binarios primitivos y un cálculo numérico de punto fijo para cada collar.

Paso 1: Escribir las Dos Ramas Inversas

Si \(y=T(x)\), entonces \(x\) satisface

$$x^2-yx-1=0.$$

Por lo tanto, cada valor \(y\) tiene dos preimágenes reales:

$$g_{\pm}(y)=\frac{y\pm\sqrt{y^2+4}}{2}.$$

Una órbita periódica puede reconstruirse hacia atrás eligiendo en cada paso la rama \(+\) o la rama \(-\). Una palabra binaria \(\sigma=(s_0,s_1,\dots,s_{m-1})\) con \(s_i\in\{-1,+1\}\) define la composición

$$G_\sigma=g_{s_{m-1}}\circ g_{s_{m-2}}\circ \cdots \circ g_{s_0}.$$

Existe una órbita de periodo \(m\) cuando un valor inicial \(x_0\) verifica

$$x_0=G_\sigma(x_0).$$

Paso 2: Por Qué los Collares Primitivos Son el Objeto Correcto

Dos palabras binarias que solo difieren por una rotación cíclica describen la misma órbita, pero con distinto punto de arranque. Así, el objeto combinatorio correcto es un collar binario y no una palabra lineal.

Si una palabra es repetición de un bloque más corto, la órbita resultante tiene un periodo fundamental menor. En consecuencia, el periodo exacto \(m\) queda representado por los collares binarios primitivos de longitud \(m\). Su número viene dado por

$$P_m=\frac{1}{m}\sum_{d\mid m}\mu(d)\,2^{m/d},$$

donde \(\mu\) es la función de Möbius. Esa misma fórmula sirve también como verificación estructural antes de empezar la parte numérica.

Paso 3: Cada Palabra de Ramas Tiene un Único Punto Fijo Real

Para \(s\in\{-1,+1\}\), la derivada de una rama inversa es

$$g_s'(y)=\frac{1}{2}\left(1+s\frac{y}{\sqrt{y^2+4}}\right).$$

Como

$$-1\lt\frac{y}{\sqrt{y^2+4}}\lt 1,$$

se obtiene

$$0\lt g_s'(y)\lt 1.$$

Esto significa que cada rama inversa es una contracción creciente en \(\mathbb{R}\), y por tanto toda composición finita \(G_\sigma\) también lo es. Por el teorema del punto fijo de Banach, cada palabra de ramas posee un único punto fijo real. Esa es la razón matemática por la cual la iteración de punto fijo converge bien y por la que cada collar primitivo produce exactamente una órbita real primitiva.

Paso 4: Reconstruir la Órbita y su Rango

Una vez conocido el punto fijo \(x_0\) de \(G_\sigma\), los demás puntos de la órbita se obtienen aplicando las ramas inversas en el orden elegido:

$$x_{i+1}=g_{s_i}(x_i)\qquad (0\le i\le m-2).$$

La última rama devuelve de \(x_{m-1}\) a \(x_0\), cerrando el ciclo. El rango asociado a \(\sigma\) es

$$R(\sigma)=\max_{0\le i\lt m} x_i-\min_{0\le i\lt m} x_i.$$

Un collar primitivo representa una órbita solo hasta rotación, pero la suma original cuenta las \(m\) posiciones cíclicas de inicio como sucesiones distintas. Por eso la contribución de un collar es

$$m\,R(\sigma),$$

y el objetivo completo queda

$$S(p)=\sum_{m=2}^{p}\ \sum_{\sigma\in\mathcal{N}_m^{\mathrm{prim}}} m\,R(\sigma).$$

Paso 5: Ejemplo Trabajado para \(m=2\)

Existe exactamente un collar binario primitivo de longitud \(2\): una rama \(+\) y una rama \(-\), salvo rotación. La órbita correspondiente es un 2-ciclo \(\{x,-x\}\). En efecto, si \(T(x)=-x\), entonces

$$x-\frac{1}{x}=-x,$$

de donde

$$2x=\frac{1}{x},\qquad x^2=\frac{1}{2}.$$

Así, los puntos de la órbita son

$$x=\pm\frac{1}{\sqrt{2}}.$$

El rango vale

$$\frac{1}{\sqrt{2}}-\left(-\frac{1}{\sqrt{2}}\right)=\sqrt{2},$$

y como hay dos posiciones cíclicas de inicio, la contribución total es

$$2\sqrt{2}=2.8284271247\ldots$$

Ese valor coincide con la comprobación más pequeña utilizada por la implementación.

Cómo Funciona el Código

Las implementaciones en C++, Python y Java siguen el mismo flujo matemático. Para cada \(m=2,3,\dots,p\), generan primero todos los collares binarios primitivos de longitud \(m\) mediante un generador recursivo. Antes de hacer cálculos en coma flotante, comparan la cantidad generada con la fórmula de Möbius para confirmar que solo hay representantes primitivos.

Para cada representante, la implementación comienza en \(0\) y aplica repetidamente la composición completa de ramas inversas \(G_\sigma\). Eso proporciona una aproximación inicial estable del punto fijo. Luego refina esa aproximación con el método de Newton aplicado a

$$\Phi(x)=G_\sigma(x)-x.$$

Durante ese refinamiento, la derivada de toda la composición se acumula con la regla de la cadena al recorrer las ramas elegidas, de modo que la corrección de Newton usa la derivada exacta del mapa inverso compuesto.

Una vez encontrado el punto fijo, la implementación despliega la órbita paso a paso, mantiene el mínimo y el máximo observados, calcula \(R(\sigma)\) y añade \(mR(\sigma)\) al total acumulado. Se usa suma compensada para reducir la cancelación numérica. Las versiones de C++ y Python reparten bloques independientes de collares entre trabajadores, mientras que la versión de Java delega la evaluación numérica al mismo núcleo subyacente y devuelve el valor decimal final.

Análisis de Complejidad

Sea

$$P_m=\frac{1}{m}\sum_{d\mid m}\mu(d)\,2^{m/d}$$

el número de collares binarios primitivos de longitud \(m\). Procesar un collar requiere \(O(mI)\) operaciones aritméticas, donde \(I\) es el pequeño número de iteraciones de punto fijo y de Newton. Por tanto, el tiempo total es

$$O\!\left(\sum_{m=2}^{p} P_m\,m\,I\right).$$

Como \(P_m\sim 2^m/m\), el crecimiento es esencialmente exponencial en \(p\), dominado por los periodos más grandes. La memoria necesaria para una longitud fija es \(O(P_m)\) para guardar los representantes de collares, más \(O(1)\) para los datos de la órbita y una pequeña cantidad de estado local por trabajador.

Notas y Referencias

  1. Página del problema: Project Euler 729
  2. Conteo de collares: Wikipedia - Necklace (combinatorics)
  3. Inversión de Möbius: Wikipedia - Möbius inversion formula
  4. Principio de contracción: Wikipedia - Banach fixed-point theorem
  5. Método de Newton: Wikipedia - Newton's method
  6. Suma compensada: Wikipedia - Kahan summation algorithm

问题概述

$$T(x)=x-\frac{1}{x},\qquad x\neq 0.$$

对于每个满足 \(2\le m\le p\) 的精确周期 \(m\),题目研究在反复作用 \(T\) 后得到的实数周期序列。如果一个周期写成 \((x_0,x_1,\dots,x_{m-1})\),那么它的区间长度定义为

$$\max_{0\le i\lt m} x_i-\min_{0\le i\lt m} x_i.$$

目标量 \(S(p)\) 就是把所有长度不超过 \(p\) 的原始周期序列的区间长度加总起来,而且同一个周期轨道的不同循环起点要分别计数。直接在实数范围里枚举所有周期序列并不现实,所以实现采用另一种表示:先用逆分支把每条轨道编码成一个二进制模式,再对每个原始模式求一个固定点,从而恢复整条周期轨道。

数学方法

核心思路是把“找所有周期序列”转换成“枚举所有原始二元项链”,然后对每个项链做一次数值固定点求解。

步骤 1:写出两个逆分支

若 \(y=T(x)\),则 \(x\) 满足二次方程

$$x^2-yx-1=0.$$

因此每个 \(y\) 都有两个实数原像:

$$g_{\pm}(y)=\frac{y\pm\sqrt{y^2+4}}{2}.$$

这说明如果我们沿着一条周期轨道“逆着走”,每一步都要在 \(+\) 分支和 \(-\) 分支之间做一次选择。于是,一个长度为 \(m\) 的二进制词 \(\sigma=(s_0,s_1,\dots,s_{m-1})\),其中 \(s_i\in\{-1,+1\}\),就决定了一个复合映射

$$G_\sigma=g_{s_{m-1}}\circ g_{s_{m-2}}\circ \cdots \circ g_{s_0}.$$

如果存在某个起点 \(x_0\) 满足

$$x_0=G_\sigma(x_0),$$

那么沿着这组逆分支选择走一圈之后又回到原点,就得到一个周期为 \(m\) 的候选轨道。

步骤 2:为什么要用原始项链而不是普通二进制串

同一条周期轨道如果只是把起点往后平移一位,得到的分支序列只是循环位移,轨道本身并没有变化。因此真正应该区分的是“循环等价类”,也就是二元项链,而不是线性排列的二进制串。

另外,如果一个二进制串是某个更短模式的重复,那么对应轨道的真正基本周期也更短,不应算作长度 \(m\) 的原始对象。所以,精确周期 \(m\) 对应的正是长度为 \(m\) 的原始二元项链。它们的个数由 Möbius 反演给出:

$$P_m=\frac{1}{m}\sum_{d\mid m}\mu(d)\,2^{m/d},$$

其中 \(\mu\) 是 Möbius 函数。实现会先用这条公式验证生成出的代表数目是否正确,再进入后续的实数计算。

步骤 3:证明每个分支词都只有一个实固定点

对 \(s\in\{-1,+1\}\) 求导,可得

$$g_s'(y)=\frac{1}{2}\left(1+s\frac{y}{\sqrt{y^2+4}}\right).$$

$$-1\lt\frac{y}{\sqrt{y^2+4}}\lt 1,$$

所以总有

$$0\lt g_s'(y)\lt 1.$$

这说明每个逆分支在整个实轴上都是单调递增的压缩映射,任何有限次复合 \(G_\sigma\) 仍然是压缩映射。根据压缩映射原理,每个 \(\sigma\) 都存在唯一的实固定点。这一点非常关键:它解释了为什么固定点迭代会稳定收敛,也保证了每个原始项链恰好对应一条唯一的实原始周期轨道,而不会出现多解或漏解。

步骤 4:由固定点恢复整条轨道并计算区间长度

求出 \(G_\sigma\) 的固定点 \(x_0\) 之后,就可以按选定的逆分支顺序逐步恢复轨道上的其他点:

$$x_{i+1}=g_{s_i}(x_i)\qquad (0\le i\le m-2).$$

最后一个分支会把 \(x_{m-1}\) 再送回 \(x_0\),于是周期闭合。与项链 \(\sigma\) 对应的区间长度记为

$$R(\sigma)=\max_{0\le i\lt m} x_i-\min_{0\le i\lt m} x_i.$$

这里还要注意一个计数细节:一个原始项链只代表“忽略起点旋转后的同一条轨道”,但题目要求把同一条轨道的 \(m\) 个循环起点都当作不同的周期序列来统计。因此,一条原始项链对总和的贡献不是 \(R(\sigma)\),而是

$$m\,R(\sigma).$$

于是最终目标可以写成

$$S(p)=\sum_{m=2}^{p}\ \sum_{\sigma\in\mathcal{N}_m^{\mathrm{prim}}} m\,R(\sigma),$$

其中 \(\mathcal{N}_m^{\mathrm{prim}}\) 表示长度为 \(m\) 的原始二元项链集合。

步骤 5:\(m=2\) 的完整示例

长度为 \(2\) 的原始二元项链只有一个,也就是一个 \(+\) 分支和一个 \(-\) 分支,循环位移视为同一个对象。对应的轨道是一个 2-周期 \(\{x,-x\}\)。因为此时必须满足 \(T(x)=-x\),所以

$$x-\frac{1}{x}=-x.$$

整理得

$$2x=\frac{1}{x},\qquad x^2=\frac{1}{2}.$$

于是轨道上的两个点就是

$$x=\pm\frac{1}{\sqrt{2}}.$$

它的区间长度为

$$\frac{1}{\sqrt{2}}-\left(-\frac{1}{\sqrt{2}}\right)=\sqrt{2}.$$

由于同一条 2-周期轨道有两个不同的循环起点,因此它对总和的贡献是

$$2\sqrt{2}=2.8284271247\ldots$$

这正好对应实现里最小的校验值,也说明“项链去重,再乘以周期长度”的思路是正确的。

代码如何工作

C++、Python 和 Java 实现共享同一条数学主线。对每个 \(m=2,3,\dots,p\),程序先用递归项链生成方法列出所有长度为 \(m\) 的原始二元项链。随后,它把生成数量与上面的 Möbius 公式进行比较,确认没有把非原始模式混入计算。

对每个代表项链,程序从 \(0\) 出发,不断应用整段逆分支复合映射 \(G_\sigma\),先得到一个稳定的固定点近似值。之后,再对

$$\Phi(x)=G_\sigma(x)-x$$

做 Newton 修正,以更快地提高精度。Newton 步骤中需要整个复合映射的导数,程序在遍历各个逆分支时用链式法则把这些导数逐段乘起来,因此修正量使用的是完整复合映射的真实导数。

固定点求出后,程序按照项链编码的逆分支顺序把一整个周期展开出来,同时维护当前最小值和最大值,从而得到 \(R(\sigma)\)。然后把 \(mR(\sigma)\) 加入总和。为了抑制大量浮点数累加时的舍入误差,代码使用了补偿求和。C++ 和 Python 版本都会把互不相关的项链块分配给不同工作单元并行处理;Java 版本则调用相同的数值核心,并返回最终的小数结果。

复杂度分析

记长度为 \(m\) 的原始二元项链数量为

$$P_m=\frac{1}{m}\sum_{d\mid m}\mu(d)\,2^{m/d}.$$

处理一条项链需要 \(O(mI)\) 次算术操作,其中 \(I\) 是固定点迭代和 Newton 修正所需的少量步数。因此总时间复杂度为

$$O\!\left(\sum_{m=2}^{p} P_m\,m\,I\right).$$

由于 \(P_m\sim 2^m/m\),整体增长本质上是关于 \(p\) 的指数级,主要成本由最大周期贡献。对某个固定长度 \(m\),内存复杂度是 \(O(P_m)\),用于存储该长度的项链代表;除此之外只需要 \(O(1)\) 的轨道数据和少量工作单元局部累加状态。

注释与参考资料

  1. 题目页面:Project Euler 729
  2. 项链计数:Wikipedia - Necklace (combinatorics)
  3. Möbius 反演:Wikipedia - Möbius inversion formula
  4. 压缩映射定理:Wikipedia - Banach fixed-point theorem
  5. Newton 方法:Wikipedia - Newton's method
  6. 补偿求和:Wikipedia - Kahan summation algorithm

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

Пусть

$$T(x)=x-\frac{1}{x},\qquad x\neq 0.$$

Для каждого точного периода \(m\), где \(2\le m\le p\), рассматриваются вещественные периодические последовательности, возникающие при повторном применении \(T\). Если один период записывается как \((x_0,x_1,\dots,x_{m-1})\), то его размах равен

$$\max_{0\le i\lt m} x_i-\min_{0\le i\lt m} x_i.$$

Искомая величина \(S(p)\) складывает эти размахи по всем примитивным периодическим последовательностям длины не больше \(p\), причем разные циклические начальные позиции считаются отдельно. Прямой перебор вещественных периодических последовательностей неэффективен, поэтому реализация кодирует каждую орбиту через выбор обратных ветвей и решает одну задачу на фиксированную точку для каждой примитивной орбиты.

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

Вместо прямого перечисления всех периодических последовательностей решение перечисляет примитивные бинарные ожерелья и для каждого из них численно находит фиксированную точку.

Шаг 1: Выписать две обратные ветви

Если \(y=T(x)\), то число \(x\) удовлетворяет уравнению

$$x^2-yx-1=0.$$

Значит, у каждого \(y\) есть два вещественных прообраза:

$$g_{\pm}(y)=\frac{y\pm\sqrt{y^2+4}}{2}.$$

Следовательно, периодическую орбиту можно восстанавливать назад, выбирая на каждом шаге либо ветвь \(+\), либо ветвь \(-\). Бинарное слово \(\sigma=(s_0,s_1,\dots,s_{m-1})\), где \(s_i\in\{-1,+1\}\), задает композицию

$$G_\sigma=g_{s_{m-1}}\circ g_{s_{m-2}}\circ \cdots \circ g_{s_0}.$$

Орбита периода \(m\) возникает тогда, когда некоторый стартовый элемент \(x_0\) удовлетворяет уравнению

$$x_0=G_\sigma(x_0).$$

Шаг 2: Почему нужны именно примитивные ожерелья

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

Если слово является повторением более короткого блока, то соответствующая орбита имеет меньший фундаментальный период. Значит, точному периоду \(m\) соответствуют именно примитивные бинарные ожерелья длины \(m\). Их число равно

$$P_m=\frac{1}{m}\sum_{d\mid m}\mu(d)\,2^{m/d},$$

где \(\mu\) обозначает функцию Мёбиуса. Эта формула используется и как структурная проверка перед началом вещественных вычислений.

Шаг 3: У каждого слова ветвей есть единственная вещественная фиксированная точка

Для \(s\in\{-1,+1\}\) производная одной обратной ветви равна

$$g_s'(y)=\frac{1}{2}\left(1+s\frac{y}{\sqrt{y^2+4}}\right).$$

Так как

$$-1\lt\frac{y}{\sqrt{y^2+4}}\lt 1,$$

получаем

$$0\lt g_s'(y)\lt 1.$$

Значит, каждая обратная ветвь является возрастающим сжимающим отображением на \(\mathbb{R}\), а любая конечная композиция \(G_\sigma\) тоже является сжатием. По теореме Банаха о неподвижной точке каждое слово \(\sigma\) имеет единственную вещественную фиксированную точку. Именно это обосновывает устойчивую сходимость итерации по фиксированной точке и гарантирует, что каждому примитивному ожерелью соответствует ровно одна вещественная примитивная орбита.

Шаг 4: Восстановить орбиту и ее размах

Когда фиксированная точка \(x_0\) для \(G_\sigma\) найдена, остальные точки орбиты получаются последовательным применением выбранных обратных ветвей:

$$x_{i+1}=g_{s_i}(x_i)\qquad (0\le i\le m-2).$$

Последняя ветвь переводит \(x_{m-1}\) обратно в \(x_0\), замыкая цикл. Размах, связанный с \(\sigma\), равен

$$R(\sigma)=\max_{0\le i\lt m} x_i-\min_{0\le i\lt m} x_i.$$

Примитивное ожерелье представляет орбиту лишь с точностью до циклического сдвига, но в исходной сумме все \(m\) циклических стартовых позиции считаются различными периодическими последовательностями. Поэтому вклад одного примитивного ожерелья равен

$$m\,R(\sigma),$$

а вся целевая сумма записывается как

$$S(p)=\sum_{m=2}^{p}\ \sum_{\sigma\in\mathcal{N}_m^{\mathrm{prim}}} m\,R(\sigma).$$

Шаг 5: Подробный пример для \(m=2\)

Для длины \(2\) существует ровно одно примитивное бинарное ожерелье: одна ветвь \(+\) и одна ветвь \(-\), с точностью до циклического сдвига. Соответствующая орбита является 2-циклом \(\{x,-x\}\). Если \(T(x)=-x\), то

$$x-\frac{1}{x}=-x,$$

откуда

$$2x=\frac{1}{x},\qquad x^2=\frac{1}{2}.$$

Следовательно, точки орбиты равны

$$x=\pm\frac{1}{\sqrt{2}}.$$

Ее размах равен

$$\frac{1}{\sqrt{2}}-\left(-\frac{1}{\sqrt{2}}\right)=\sqrt{2},$$

а поскольку существует две циклические стартовые позиции, полный вклад составляет

$$2\sqrt{2}=2.8284271247\ldots$$

Именно это значение совпадает с самым маленьким контрольным результатом, используемым в реализации.

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

Реализации на C++, Python и Java следуют одной и той же математической схеме. Для каждого \(m=2,3,\dots,p\) сначала рекурсивно генерируются все примитивные бинарные ожерелья длины \(m\). Затем число полученных представителей сравнивается с формулой Мёбиуса, чтобы убедиться, что в вычисление не попали непервичные шаблоны.

Для каждого представителя вычисление начинается из точки \(0\), после чего многократно применяется полная композиция обратных ветвей \(G_\sigma\). Это дает устойчивое начальное приближение к фиксированной точке. Далее используется метод Ньютона для функции

$$\Phi(x)=G_\sigma(x)-x.$$

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

После нахождения фиксированной точки код разворачивает всю орбиту шаг за шагом, отслеживает минимум и максимум, вычисляет \(R(\sigma)\) и прибавляет \(mR(\sigma)\) к общей сумме. Чтобы уменьшить потерю точности при суммировании большого числа вещественных вкладов, используется компенсированное суммирование. Версии на C++ и Python распределяют независимые блоки ожерелий между рабочими задачами, а версия на Java обращается к тому же численному ядру и возвращает итоговое десятичное значение.

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

Пусть

$$P_m=\frac{1}{m}\sum_{d\mid m}\mu(d)\,2^{m/d}$$

обозначает число примитивных бинарных ожерелий длины \(m\). Обработка одного ожерелья требует \(O(mI)\) арифметических операций, где \(I\) - небольшое число итераций по фиксированной точке и шагов Ньютона. Поэтому суммарное время работы равно

$$O\!\left(\sum_{m=2}^{p} P_m\,m\,I\right).$$

Поскольку \(P_m\sim 2^m/m\), рост по сути экспоненциален по \(p\), и наибольший вклад дают самые длинные периоды. Память для одной фиксированной длины составляет \(O(P_m)\) на хранение представителей ожерелий, плюс \(O(1)\) на данные текущей орбиты и небольшой объем локального состояния для рабочих накопителей.

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

  1. Страница задачи: Project Euler 729
  2. Подсчет ожерелий: Wikipedia - Necklace (combinatorics)
  3. Формула обращения Мёбиуса: Wikipedia - Möbius inversion formula
  4. Теорема Банаха о неподвижной точке: Wikipedia - Banach fixed-point theorem
  5. Метод Ньютона: Wikipedia - Newton's method
  6. Компенсированное суммирование: Wikipedia - Kahan summation algorithm

ملخص المشكلة

لنعرّف

$$T(x)=x-\frac{1}{x},\qquad x\neq 0.$$

لكل فترة دقيقة \(m\) حيث \(2\le m\le p\)، تبحث المسألة في المتتاليات الدورية الحقيقية الناتجة عن تكرار تطبيق \(T\). إذا كُتبت دورة واحدة على الصورة \((x_0,x_1,\dots,x_{m-1})\)، فإن مداها يساوي

$$\max_{0\le i\lt m} x_i-\min_{0\le i\lt m} x_i.$$

القيمة المطلوبة \(S(p)\) هي مجموع هذه المديات على جميع المتتاليات الدورية الأولية ذات الأطوال حتى \(p\)، مع اعتبار نقاط البداية الدورية المختلفة متتاليات متميزة. البحث المباشر بين جميع المتتاليات الحقيقية غير عملي، ولذلك تعيد المنهجية ترميز كل مدار عبر فروع عكسية، ثم تحل مسألة نقطة تثبيت واحدة لكل مدار أولي.

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

الفكرة الأساسية هي استبدال التعداد المباشر للمتتاليات الدورية بتعداد القلائد الثنائية الأولية، ثم حل نقطة تثبيت عدديا لكل قلادة.

الخطوة 1: كتابة الفرعين العكسيين

إذا كان \(y=T(x)\)، فإن \(x\) يحقق

$$x^2-yx-1=0.$$

ومن ثم يملك كل \(y\) صورتين سابقتين حقيقيتين:

$$g_{\pm}(y)=\frac{y\pm\sqrt{y^2+4}}{2}.$$

هذا يعني أن إعادة بناء المدار إلى الخلف تتطلب اختيار الفرع \(+\) أو الفرع \(-\) في كل خطوة. إذا أخذنا كلمة ثنائية \(\sigma=(s_0,s_1,\dots,s_{m-1})\) حيث \(s_i\in\{-1,+1\}\)، فإنها تحدد التركيب

$$G_\sigma=g_{s_{m-1}}\circ g_{s_{m-2}}\circ \cdots \circ g_{s_0}.$$

ويظهر مدار ذو فترة \(m\) عندما توجد قيمة ابتدائية \(x_0\) تحقق

$$x_0=G_\sigma(x_0).$$

الخطوة 2: لماذا القلائد الأولية هي التمثيل الصحيح

إذا اختلفت كلمتان ثنائيتان فقط بإزاحة دورية، فهما تصفان المدار نفسه لكن مع اختيار نقطة بداية مختلفة. لذلك فالكائن التوافقي الطبيعي هنا هو القلادة الثنائية لا الكلمة الخطية الخام.

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

$$P_m=\frac{1}{m}\sum_{d\mid m}\mu(d)\,2^{m/d},$$

حيث \(\mu\) هي دالة موبيوس. وتُستخدم هذه الصيغة أيضا كاختبار بنيوي قبل بدء الحسابات العددية.

الخطوة 3: لكل كلمة فروع نقطة تثبيت حقيقية وحيدة

بالنسبة إلى \(s\in\{-1,+1\}\)، فإن مشتقة الفرع العكسي هي

$$g_s'(y)=\frac{1}{2}\left(1+s\frac{y}{\sqrt{y^2+4}}\right).$$

وبما أن

$$-1\lt\frac{y}{\sqrt{y^2+4}}\lt 1,$$

فنحصل على

$$0\lt g_s'(y)\lt 1.$$

أي أن كل فرع عكسي هو تقليص متزايد على \(\mathbb{R}\)، وبالتالي فإن أي تركيب منته \(G_\sigma\) يبقى تقليصا. وبحسب مبدأ نقطة التثبيت للتطبيقات المتقلصة، تمتلك كل كلمة فروع \(\sigma\) نقطة تثبيت حقيقية وحيدة. هذا يبرر رياضيا نجاح تكرار نقطة التثبيت، كما يضمن أن كل قلادة أولية تقابل مدارا حقيقيا أوليا واحدا فقط.

الخطوة 4: إعادة بناء المدار وحساب المدى

بعد إيجاد نقطة التثبيت \(x_0\) للتركيب \(G_\sigma\)، نحصل على بقية نقاط المدار بتطبيق الفروع العكسية بالترتيب المختار:

$$x_{i+1}=g_{s_i}(x_i)\qquad (0\le i\le m-2).$$

الفرع الأخير يعيد \(x_{m-1}\) إلى \(x_0\) فيغلق الدورة. والمدى المرتبط بـ \(\sigma\) هو

$$R(\sigma)=\max_{0\le i\lt m} x_i-\min_{0\le i\lt m} x_i.$$

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

$$m\,R(\sigma),$$

وبالتالي يصبح الهدف الكلي

$$S(p)=\sum_{m=2}^{p}\ \sum_{\sigma\in\mathcal{N}_m^{\mathrm{prim}}} m\,R(\sigma).$$

الخطوة 5: مثال محلول للحالة \(m=2\)

توجد قلادة ثنائية أولية واحدة فقط بطول \(2\): فرع \(+\) وفرع \(-\) مع اعتبار الدوران متماثلا. والمدار الموافق هو دورة بطول \(2\) من الشكل \(\{x,-x\}\). فعلا إذا كان \(T(x)=-x\)، نحصل على

$$x-\frac{1}{x}=-x,$$

ومنها

$$2x=\frac{1}{x},\qquad x^2=\frac{1}{2}.$$

إذن نقطتا المدار هما

$$x=\pm\frac{1}{\sqrt{2}}.$$

ومن ثم يكون المدى

$$\frac{1}{\sqrt{2}}-\left(-\frac{1}{\sqrt{2}}\right)=\sqrt{2},$$

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

$$2\sqrt{2}=2.8284271247\ldots$$

وهذه هي بالضبط أصغر قيمة تحقق داخلية تستخدمها المنهجية.

كيف يعمل الكود

تسير تطبيقات C++ وPython وJava على الخط الرياضي نفسه. لكل \(m=2,3,\dots,p\)، يجري أولا توليد جميع القلائد الثنائية الأولية ذات الطول \(m\) بواسطة مولد تراجعي. وقبل إجراء أي حسابات عائمة، تتم مقارنة العدد الناتج مع صيغة موبيوس السابقة للتأكد من أن التوليد لم يحتو إلا على ممثلين أوليين.

بالنسبة إلى كل ممثل، يبدأ التنفيذ من \(0\) ويطبق التركيب الكامل للفروع العكسية \(G_\sigma\) مرارا للحصول على تقريب أولي ثابت لنقطة التثبيت. بعد ذلك تُستخدم طريقة نيوتن على الدالة

$$\Phi(x)=G_\sigma(x)-x.$$

وأثناء هذا التحسين، تُجمّع مشتقة التركيب كله باستعمال قاعدة السلسلة خلال المرور على الفروع المختارة، بحيث تعتمد خطوة نيوتن على المشتقة الحقيقية للتركيب العكسي الكامل.

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

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

إذا كان

$$P_m=\frac{1}{m}\sum_{d\mid m}\mu(d)\,2^{m/d}$$

هو عدد القلائد الثنائية الأولية ذات الطول \(m\)، فإن معالجة قلادة واحدة تتطلب \(O(mI)\) عملية حسابية، حيث إن \(I\) هو العدد الصغير من تكرارات نقطة التثبيت وخطوات نيوتن. لذلك فإن الزمن الكلي يساوي

$$O\!\left(\sum_{m=2}^{p} P_m\,m\,I\right).$$

وبما أن \(P_m\sim 2^m/m\)، فإن النمو في جوهره أسي بالنسبة إلى \(p\)، وتأتي الكلفة المسيطرة من أكبر الفترات. أما الذاكرة المطلوبة لطول واحد ثابت فهي \(O(P_m)\) لتخزين ممثلي القلائد، إضافة إلى \(O(1)\) لبيانات المدار وحالة تجميع محلية صغيرة لكل عامل.

هوامش ومراجع

  1. صفحة المسألة: Project Euler 729
  2. عدّ القلائد: Wikipedia - Necklace (combinatorics)
  3. عكس موبيوس: Wikipedia - Möbius inversion formula
  4. مبرهنة باناخ لنقطة التثبيت: Wikipedia - Banach fixed-point theorem
  5. طريقة نيوتن: Wikipedia - Newton's method
  6. الجمع التعويضي: Wikipedia - Kahan summation algorithm