Problem Summary

The problem asks for the total potency of all magic bracelets whose bead labels do not exceed \(N\). A bracelet is treated as a cycle in which neighboring labels \(u\) and \(v\) are compatible exactly when

$$uv=x^2+1$$

for some integer \(x\). The implementations turn this into a graph problem: vertices are the admissible labels that can actually participate in this relation, and a bracelet is a cycle in that compatibility graph. Its potency is the sum of the labels on the cycle, counted once up to rotation and reflection.

The arithmetic filter used by the implementations keeps exactly \(1\), \(2\), and the families \(p^e\) and \(2p^e\) with \(p\equiv1\pmod4\). The whole solution is about building this graph canonically and then summing cycle weights without double counting.

Mathematical Approach

The key move is to replace each admissible label by a primitive lattice vector, so that the condition \(uv=x^2+1\) becomes a geometric determinant condition.

Step 1: Represent each admissible label by a primitive sum of two squares

For every prime \(p\equiv1\pmod4\), choose one decomposition

$$p=a^2+b^2.$$

In Gaussian integers this means

$$p=(a+bi)(a-bi).$$

Raising \(a+bi\) to the \(e\)-th power gives a primitive representative of \(p^e\):

$$z_{p^e}=(a+bi)^e,\qquad N(z_{p^e})=p^e.$$

Multiplying by \(1+i\) produces the doubled family:

$$N\bigl((1+i)z_{p^e}\bigr)=2p^e.$$

After taking absolute values and ordering coordinates so that \(x\ge y\ge0\), every stored label has one canonical primitive vector

$$n=x^2+y^2,\qquad \gcd(x,y)=1.$$

This is why the implementations can work with pairs \((x,y)\) instead of directly with the original square-plus-one condition.

Step 2: Convert compatibility into a determinant condition

Suppose two labels are represented by primitive vectors \((u,v)\) and \((s,t)\), so their labels are

$$m=u^2+v^2,\qquad n=s^2+t^2.$$

The Brahmagupta-Fibonacci identity gives

$$mn=(us+vt)^2+(ut-vs)^2.$$

Therefore, if the two vectors satisfy

$$|ut-vs|=1,$$

then

$$mn=(us+vt)^2+1,$$

so the two labels are compatible. Thus adjacency is encoded by a unimodular condition: neighboring beads correspond to primitive vectors whose determinant has absolute value \(1\).

This reformulation is the structural heart of the algorithm. Instead of searching for squares \(x^2+1\), the implementations search for determinant-\(1\) relationships among canonical primitive vectors.

Step 3: Recover the two smaller parent labels from Bézout coefficients

Take a canonical vector \((x,y)\) for some label \(n>1\). Because \(\gcd(x,y)=1\), the extended Euclidean algorithm provides integers \(\alpha,\beta\) such that

$$\alpha x+\beta y=1.$$

From this relation, the implementations build two smaller vectors

$$q_1=\bigl(|\beta|,|\alpha|\bigr),\qquad q_2=\bigl(x-|\beta|,\ y-|\alpha|\bigr).$$

These satisfy

$$q_1+q_2=(x,y),\qquad |\det(q_1,q_2)|=1.$$

Let their norms be

$$n_1=\|q_1\|^2,\qquad n_2=\|q_2\|^2.$$

Then \(n_1\) and \(n_2\) are the two parent labels of \(n\). Because \(q_1\) and \(q_2\) are unimodular, the parent labels are compatible. Since \((x,y)=q_1+q_2\), the current label is also compatible with each parent, because

$$\det(q_1,q_1+q_2)=\det(q_1,q_2),\qquad \det(q_1+q_2,q_2)=\det(q_1,q_2).$$

So every stored label sits naturally above an edge joining two smaller compatible labels. The special case \(n=2\) collapses to the single base edge from \(2\) down to \(1\).

Step 4: View the compatibility graph as a sparse DAG when oriented downward

Now orient every compatibility edge from the larger label to the smaller one. The parent construction above shows that every admissible label has at most two outgoing edges, and both targets are smaller. Hence the oriented graph is acyclic and very sparse.

A bracelet is an undirected cycle in the original compatibility graph. Once the edges are oriented downward, every such cycle has a unique largest label \(M\). Cutting the cycle at \(M\) leaves two descending chains from \(M\) to the endpoints of some lower compatibility edge. That observation gives a canonical place where each bracelet should be counted: exactly once, at its unique maximum.

Step 5: Count open descending paths and close them into bracelets

For every oriented edge \(i\to j\), the dynamic program stores two quantities:

$$C_{i\to j}=\text{number of open descending path fragments whose boundary edge is }(i,j),$$

$$T_{i\to j}=\text{sum of the interior label totals over those fragments}.$$

When the algorithm later reaches the edge \((i,j)\), every stored fragment can be closed by that edge into one bracelet. The added potency is

$$C_{i\to j}(i+j)+T_{i\to j}.$$

After closing the old fragments, the implementation inserts the trivial one-vertex fragment above the same edge, which seeds future longer paths.

If a label \(i\) has two parents \(a\) and \(b\), then every fragment above \(i\to a\) can be paired with every fragment above \(i\to b\). This creates an open path between \(a\) and \(b\) passing through \(i\). If the two incoming edge states are \((C_a,T_a)\) and \((C_b,T_b)\), the merged contribution is

$$C'=C_aC_b,$$

$$T'=C_aC_b\,i+C_aT_b+C_bT_a.$$

The count multiplies because the left and right choices are independent, and the interior sum adds the new vertex \(i\) plus the previously stored interior sums from both sides.

Worked Example: the triangle \(\{2,5,13\}\)

The label \(13\) has canonical vector

$$13=3^2+2^2,\qquad (x,y)=(3,2).$$

A Bézout relation is

$$1\cdot3+(-1)\cdot2=1.$$

So the two parent vectors are

$$q_1=(1,1),\qquad q_2=(2,1),$$

with norms

$$\|q_1\|^2=2,\qquad \|q_2\|^2=5.$$

Because

$$|\det(q_1,q_2)|=|1\cdot1-1\cdot2|=1,$$

the parent labels satisfy

$$2\cdot5=3^2+1,$$

so they are compatible. Also \(13\) is compatible with both \(2\) and \(5\), so \(\{2,5,13\}\) is a bracelet with potency

$$2+5+13=20.$$

The same mechanism gives \(5\) the parents \(1\) and \(2\), which yields the smaller triangle \(\{1,2,5\}\) with potency \(8\). These tiny examples are exactly what the edge-based dynamic program is reconstructing at scale.

How the Code Works

The C++, Python, and Java implementations follow the same plan. First they build a prime table up to \(N\) and a direct lookup table for exact squares. Then they generate the admissible labels together with one canonical primitive vector for each of them. Prime powers are produced by repeated multiplication inside \(\mathbb Z[i]\), and the doubled family is produced by one extra multiplication by \(1+i\).

Next, every admissible label greater than \(1\) is decomposed into its two parent labels via the extended Euclidean algorithm. Only parents that are themselves admissible are kept, so the result is a sparse directed graph from larger labels to smaller labels.

Finally, the implementations sweep the labels in descending order. For each downward edge they first close any previously accumulated open fragments, which adds finished bracelet potencies to the answer. Then they seed the trivial fragment on that edge. When a label has two parents, the two edge states are merged and forwarded to the lower edge joining those parents. Because every bracelet has a unique maximum label, this procedure counts each bracelet exactly once.

Complexity Analysis

Let \(B(N)\) be the number of admissible labels materialized by the algorithm. The prime sieve and the lookup tables use \(O(N)\) memory, and the sieve itself runs in \(O(N)\) time in these implementations. For each prime \(p\equiv1\pmod4\), the code searches for one representation \(p=a^2+b^2\) using exact-square tests; in the current implementation this costs at most \(O(\sqrt p)\) trials for that prime. After that, generating all stored labels, creating parent edges, and running the downward dynamic program are all \(O(B(N))\).

So the dominant graph phase is linear in the number of admissible labels, while the preprocessing cost is governed by the sieve, the square table, and the search for prime representations as sums of two squares.

Footnotes and References

  1. Problem page: https://projecteuler.net/problem=846
  2. Gaussian integers: Wikipedia - Gaussian integer
  3. Fermat's theorem on sums of two squares: Wikipedia - Fermat's theorem on sums of two squares
  4. Brahmagupta-Fibonacci identity: Wikipedia - Brahmagupta-Fibonacci identity
  5. Bézout's identity: Wikipedia - Bézout's identity

Problemzusammenfassung

Gesucht ist die gesamte Potenz aller magischen Armbänder, deren Beschriftungen höchstens \(N\) sind. Ein Armband wird als Zyklus betrachtet, in dem zwei benachbarte Werte \(u\) und \(v\) genau dann kompatibel sind, wenn

$$uv=x^2+1$$

für ein ganzzahliges \(x\) gilt. Die Implementierungen übersetzen das in ein Graphenproblem: Knoten sind die zulässigen Werte, die in dieser Relation überhaupt auftreten können, und ein Armband ist ein Zyklus im Kompatibilitätsgraphen. Seine Potenz ist die Summe der Knotenwerte, wobei Drehungen und Spiegelungen nur einmal gezählt werden.

Der arithmetische Filter der Implementierungen behält genau \(1\), \(2\) sowie die Familien \(p^e\) und \(2p^e\) mit \(p\equiv1\pmod4\). Die Lösung besteht darin, diesen Graphen kanonisch aufzubauen und danach die Gewichte aller Zyklen ohne Doppelzählung zu summieren.

Mathematischer Ansatz

Der entscheidende Schritt ist, jeden zulässigen Wert durch einen primitiven Gittervektor zu ersetzen. Dann wird die Bedingung \(uv=x^2+1\) zu einer Determinantenbedingung.

Schritt 1: Jeder zulässige Wert als primitive Summe zweier Quadrate

Für jede Primzahl \(p\equiv1\pmod4\) wählen wir eine Darstellung

$$p=a^2+b^2.$$

In den gaußschen Zahlen bedeutet das

$$p=(a+bi)(a-bi).$$

Die \(e\)-te Potenz von \(a+bi\) liefert einen primitiven Repräsentanten von \(p^e\):

$$z_{p^e}=(a+bi)^e,\qquad N(z_{p^e})=p^e.$$

Durch Multiplikation mit \(1+i\) erhält man die verdoppelte Familie:

$$N\bigl((1+i)z_{p^e}\bigr)=2p^e.$$

Nach Absolutbetrag und Sortierung der Koordinaten mit \(x\ge y\ge0\) besitzt jeder gespeicherte Wert genau einen kanonischen primitiven Vektor

$$n=x^2+y^2,\qquad \gcd(x,y)=1.$$

Darum können die Implementierungen mit Paaren \((x,y)\) arbeiten, statt direkt nach Gleichungen der Form Quadrat plus \(1\) zu suchen.

Schritt 2: Kompatibilität wird zu \(|\det|=1\)

Seien zwei Werte durch primitive Vektoren \((u,v)\) und \((s,t)\) dargestellt. Dann sind ihre Werte

$$m=u^2+v^2,\qquad n=s^2+t^2.$$

Aus der Brahmagupta-Fibonacci-Identität folgt

$$mn=(us+vt)^2+(ut-vs)^2.$$

Gilt also

$$|ut-vs|=1,$$

dann erhält man

$$mn=(us+vt)^2+1,$$

also sind die beiden Werte kompatibel. Nachbarschaft im Armband entspricht daher einer unimodularen Bedingung: benachbarte Perlen haben primitive Vektoren mit Determinante \(\pm1\).

Das ist die eigentliche strukturelle Vereinfachung. Anstelle einer Suche nach Quadraten \(x^2+1\) suchen die Implementierungen nach Determinantenbeziehungen zwischen kanonischen primitiven Vektoren.

Schritt 3: Die zwei kleineren Eltern aus Bézout-Koeffizienten gewinnen

Betrachte den kanonischen Vektor \((x,y)\) eines Wertes \(n>1\). Wegen \(\gcd(x,y)=1\) liefert der erweiterte Euklid ganze Zahlen \(\alpha,\beta\) mit

$$\alpha x+\beta y=1.$$

Daraus konstruieren die Implementierungen zwei kleinere Vektoren

$$q_1=\bigl(|\beta|,|\alpha|\bigr),\qquad q_2=\bigl(x-|\beta|,\ y-|\alpha|\bigr).$$

Diese erfüllen

$$q_1+q_2=(x,y),\qquad |\det(q_1,q_2)|=1.$$

Mit ihren Normen

$$n_1=\|q_1\|^2,\qquad n_2=\|q_2\|^2$$

erhält man die beiden Elternwerte von \(n\). Da \(q_1\) und \(q_2\) unimodular sind, sind die Eltern kompatibel. Weil außerdem \((x,y)=q_1+q_2\) gilt, ist der aktuelle Wert auch mit jedem Elternwert kompatibel, denn

$$\det(q_1,q_1+q_2)=\det(q_1,q_2),\qquad \det(q_1+q_2,q_2)=\det(q_1,q_2).$$

Jeder gespeicherte Wert sitzt also über einer Kante, die zwei kleinere kompatible Werte verbindet. Der Sonderfall \(n=2\) fällt auf die eine Basiskante von \(2\) nach \(1\) zusammen.

Schritt 4: Der Kompatibilitätsgraph wird orientiert zu einem dünnen DAG

Nun orientieren wir jede Kompatibilitätskante vom größeren zum kleineren Wert. Die Elternkonstruktion zeigt, dass jeder zulässige Wert höchstens zwei ausgehende Kanten besitzt und beide Ziele kleiner sind. Damit ist der orientierte Graph azyklisch und sehr dünn.

Ein Armband ist ein ungerichteter Zyklus im ursprünglichen Kompatibilitätsgraphen. Nach der Orientierung besitzt jeder solcher Zyklus einen eindeutig größten Wert \(M\). Schneidet man dort auf, bleiben zwei absteigende Ketten von \(M\) zu den Endpunkten einer kleineren Kompatibilitätskante. Genau an diesem Maximum sollte also gezählt werden, und zwar genau einmal.

Schritt 5: Offene absteigende Pfade zählen und zu Zyklen schließen

Für jede orientierte Kante \(i\to j\) speichert die dynamische Programmierung zwei Größen:

$$C_{i\to j}=\text{Anzahl offener absteigender Pfadfragmente mit Randkante }(i,j),$$

$$T_{i\to j}=\text{Summe der inneren Knotenwerte über diese Fragmente}.$$

Erreicht der Algorithmus später die Kante \((i,j)\), so kann jedes gespeicherte Fragment durch diese Kante zu genau einem Armband geschlossen werden. Der Potenzbeitrag ist

$$C_{i\to j}(i+j)+T_{i\to j}.$$

Danach fügt die Implementierung auf derselben Kante das triviale Ein-Knoten-Fragment ein und erzeugt damit Saatgut für spätere längere Pfade.

Hat ein Wert \(i\) zwei Eltern \(a\) und \(b\), dann kann jedes Fragment über \(i\to a\) mit jedem Fragment über \(i\to b\) gepaart werden. So entsteht ein offener Pfad zwischen \(a\) und \(b\), der durch \(i\) verläuft. Sind die beiden Zustände \((C_a,T_a)\) und \((C_b,T_b)\), dann lautet der zusammengeführte Beitrag

$$C'=C_aC_b,$$

$$T'=C_aC_b\,i+C_aT_b+C_bT_a.$$

Die Anzahl multipliziert sich wegen der unabhängigen linken und rechten Wahlmöglichkeiten, und bei der inneren Summe kommen der neue Knoten \(i\) sowie die bisherigen inneren Summen beider Seiten hinzu.

Durchgerechnetes Beispiel: das Dreieck \(\{2,5,13\}\)

Der Wert \(13\) hat den kanonischen Vektor

$$13=3^2+2^2,\qquad (x,y)=(3,2).$$

Eine Bézout-Beziehung ist

$$1\cdot3+(-1)\cdot2=1.$$

Daraus folgen die beiden Elternvektoren

$$q_1=(1,1),\qquad q_2=(2,1),$$

mit Normen

$$\|q_1\|^2=2,\qquad \|q_2\|^2=5.$$

Weil

$$|\det(q_1,q_2)|=|1\cdot1-1\cdot2|=1,$$

gilt

$$2\cdot5=3^2+1,$$

also sind die Eltern kompatibel. Zugleich ist \(13\) mit \(2\) und \(5\) kompatibel, daher bildet \(\{2,5,13\}\) ein Armband mit Potenz

$$2+5+13=20.$$

Ganz analog hat \(5\) die Eltern \(1\) und \(2\), woraus das kleinere Dreieck \(\{1,2,5\}\) mit Potenz \(8\) entsteht. Genau diese kleinen Muster rekonstruiert die kantenbasierte DP im großen Maßstab.

Wie der Code arbeitet

Die C++-, Python- und Java-Implementierungen folgen demselben Plan. Zuerst bauen sie eine Primtabelle bis \(N\) und eine direkte Nachschlagetabelle für exakte Quadrate auf. Danach erzeugen sie die zulässigen Werte mitsamt je einem kanonischen primitiven Vektor. Primzahlpotenzen entstehen durch wiederholte Multiplikation in \(\mathbb Z[i]\), und die verdoppelte Familie durch eine zusätzliche Multiplikation mit \(1+i\).

Anschließend wird jeder zulässige Wert größer als \(1\) mithilfe des erweiterten Euklid in seine beiden Elternwerte zerlegt. Es bleiben nur solche Eltern erhalten, die selbst zulässig sind; dadurch entsteht ein dünner gerichteter Graph von großen zu kleinen Werten.

Zum Schluss laufen die Implementierungen die Werte absteigend durch. Auf jeder abwärts gerichteten Kante schließen sie zunächst alle zuvor angesammelten offenen Fragmente und addieren deren Potenzen zum Ergebnis. Danach setzen sie das triviale Fragment auf dieser Kante. Wenn ein Wert zwei Eltern hat, werden die beiden Kantenzustände verschmolzen und auf die tiefere Kante zwischen diesen Eltern weitergereicht. Weil jedes Armband einen eindeutigen Maximalwert besitzt, wird jedes davon genau einmal gezählt.

Komplexitätsanalyse

Sei \(B(N)\) die Anzahl der vom Algorithmus materialisierten zulässigen Werte. Primtabelle und Nachschlagetabellen benötigen \(O(N)\) Speicher, und das Sieb selbst läuft in diesen Implementierungen in \(O(N)\) Zeit. Für jede Primzahl \(p\equiv1\pmod4\) sucht der Code mit Exaktquadrat-Tests eine Darstellung \(p=a^2+b^2\); in der vorliegenden Implementierung kostet das im schlimmsten Fall \(O(\sqrt p)\) Versuche für diese Primzahl. Danach sind das Erzeugen aller gespeicherten Werte, das Aufbauen der Elternkanten und die absteigende DP jeweils \(O(B(N))\).

Die eigentliche Graphphase ist also linear in der Zahl der zulässigen Werte, während die Vorverarbeitung von Sieb, Quadrattabelle und der Suche nach Primdarstellungen als Summe zweier Quadrate bestimmt wird.

Fußnoten und Referenzen

  1. Problemseite: https://projecteuler.net/problem=846
  2. Gaußsche Zahlen: Wikipedia - Gaußsche Zahlen
  3. Satz von Fermat über Summen von zwei Quadraten: Wikipedia - Summe von zwei Quadraten
  4. Brahmagupta-Fibonacci-Identität: Wikipedia - Brahmagupta-Fibonacci-Identität
  5. Bézout-Identität: Wikipedia - Satz von Bézout

Problem Özeti

İstenen şey, boncuk etiketleri \(N\)'i aşmayan tüm sihirli bileziklerin toplam potency değeridir. Bir bilezik, komşu iki etiket \(u\) ve \(v\) için ancak

$$uv=x^2+1$$

olduğunda geçerli sayılan bir çevrim olarak ele alınır. C++, Python ve Java uygulamaları bu durumu bir grafik problemine çevirir: düğümler bu bağıntıya gerçekten girebilen uygun etiketlerdir, bilezik ise uyumluluk grafiğinde bir çevrimdir. Potency, çevrim üzerindeki etiketlerin toplamıdır; dönme ve yansıma aynı bilezik kabul edilir.

Uygulamaların kullandığı aritmetik filtre tam olarak \(1\), \(2\) ve \(p\equiv1\pmod4\) için \(p^e\) ile \(2p^e\) ailelerini bırakır. Çözümün özü, bu grafiği kanonik biçimde kurup sonra çevrim ağırlıklarını çakışmasız toplamaktır.

Matematiksel Yaklaşım

Ana fikir, her uygun etiketi ilkel bir kafes vektörü ile temsil etmektir. Böylece \(uv=x^2+1\) koşulu determinantı \(1\) olan bir geometrik koşula dönüşür.

Adım 1: Her uygun etiketi iki karenin ilkel toplamı olarak temsil et

Her \(p\equiv1\pmod4\) asalı için

$$p=a^2+b^2$$

şeklinde bir ayrışım seçilir. Gauss tamsayılarında bu,

$$p=(a+bi)(a-bi)$$

demektir. \(a+bi\)'nin \(e\). kuvveti \(p^e\) için ilkel bir temsil üretir:

$$z_{p^e}=(a+bi)^e,\qquad N(z_{p^e})=p^e.$$

\(1+i\) ile çarpmak ise iki katlı aileyi verir:

$$N\bigl((1+i)z_{p^e}\bigr)=2p^e.$$

Koordinatların mutlak değerleri alınıp \(x\ge y\ge0\) olacak şekilde sıralandığında, depolanan her etiket için tek bir kanonik ilkel vektör elde edilir:

$$n=x^2+y^2,\qquad \gcd(x,y)=1.$$

Bu yüzden uygulamalar, doğrudan kare artı \(1\) denklemleriyle uğraşmak yerine \((x,y)\) çiftleri üzerinde çalışır.

Adım 2: Uyumluluğu determinant koşuluna dönüştür

İki etiketin ilkel vektörleri \((u,v)\) ve \((s,t)\) olsun. Etiketlerin kendileri

$$m=u^2+v^2,\qquad n=s^2+t^2$$

şeklindedir. Brahmagupta-Fibonacci özdeşliği şunu verir:

$$mn=(us+vt)^2+(ut-vs)^2.$$

Dolayısıyla

$$|ut-vs|=1$$

ise

$$mn=(us+vt)^2+1$$

olur ve iki etiket uyumludur. Yani komşuluk, determinantın mutlak değerinin \(1\) olmasıyla kodlanır.

Bu dönüşüm algoritmanın kalbidir. Uygulamalar, \(x^2+1\) biçimindeki çarpımları aramak yerine kanonik ilkel vektörler arasındaki unimodüler ilişkileri arar.

Adım 3: Bézout katsayılarından iki küçük ebeveyni çıkar

\(n>1\) etiketi için kanonik vektör \((x,y)\) olsun. \(\gcd(x,y)=1\) olduğundan genişletilmiş Öklid algoritması

$$\alpha x+\beta y=1$$

eşitliğini sağlayan \(\alpha,\beta\) tamsayılarını verir. Uygulamalar buradan iki küçük vektör kurar:

$$q_1=\bigl(|\beta|,|\alpha|\bigr),\qquad q_2=\bigl(x-|\beta|,\ y-|\alpha|\bigr).$$

Bunlar

$$q_1+q_2=(x,y),\qquad |\det(q_1,q_2)|=1$$

özelliklerini sağlar. Normlarını

$$n_1=\|q_1\|^2,\qquad n_2=\|q_2\|^2$$

dersek, bunlar \(n\)'nin iki ebeveyn etiketidir. \(q_1\) ile \(q_2\) unimodüler olduğundan ebeveynler birbirleriyle uyumludur. Ayrıca \((x,y)=q_1+q_2\) olduğu için mevcut etiket her iki ebeveynle de uyumludur; çünkü

$$\det(q_1,q_1+q_2)=\det(q_1,q_2),\qquad \det(q_1+q_2,q_2)=\det(q_1,q_2).$$

Böylece depolanan her etiket, iki küçük uyumlu etiketi bağlayan bir kenarın üzerinde yer alır. \(n=2\) özel durumunda bu yapı, \(2\)'den \(1\)'e giden tek temel kenara çöker.

Adım 4: Uyumluluk grafiğini aşağı yönlü seyrek bir DAG olarak gör

Artık her uyumluluk kenarını büyük etiketten küçük etikete yönlendiriyoruz. Yukarıdaki ebeveyn yapısı, her uygun etiketin en fazla iki çıkış kenarı olduğunu ve her iki hedefin de daha küçük olduğunu gösterir. Bu yüzden yönlendirilmiş grafik döngüsüzdür ve çok seyrektir.

Bir bilezik, özgün uyumluluk grafiğinde yönsüz bir çevrimdir. Kenarlar aşağı yönlendirildiğinde her çevrimin tek bir maksimum etiketi \(M\) vardır. Çevrimi bu maksimumda kestiğimizde, \(M\)'den daha aşağıdaki bir uyumluluk kenarının iki ucuna inen iki azalan zincir kalır. Bu da her bileziğin tam olarak nerede sayılması gerektiğini gösterir: kendi benzersiz maksimumunda.

Adım 5: Açık azalan yolları say ve onları bileziğe kapat

Her yönlü \(i\to j\) kenarı için dinamik programlama iki büyüklük tutar:

$$C_{i\to j}=\text{sınır kenarı }(i,j)\text{ olan açık azalan yol parçası sayısı},$$

$$T_{i\to j}=\text{bu parçalardaki iç etiket toplamlarının toplamı}.$$

Algoritma daha sonra \((i,j)\) kenarına geldiğinde, saklanan her açık parça bu kenarla kapatılarak tek bir bilezik oluşturur. Eklenen potency miktarı

$$C_{i\to j}(i+j)+T_{i\to j}$$

olur. Eski parçalar kapatıldıktan sonra, uygulama aynı kenar üzerinde gelecekte büyüyecek olan trivial tek-düğüm parçasını ekler.

Eğer \(i\)'nin iki ebeveyni \(a\) ve \(b\) ise, \(i\to a\) üzerindeki her parça \(i\to b\) üzerindeki her parça ile eşleştirilebilir. Böylece \(a\) ile \(b\) arasında, \(i\)'den geçen yeni bir açık yol oluşur. İki giriş durumu \((C_a,T_a)\) ve \((C_b,T_b)\) ise birleşim formülü

$$C'=C_aC_b,$$

$$T'=C_aC_b\,i+C_aT_b+C_bT_a$$

şeklindedir. Sayıların çarpılması sol ve sağ seçimlerin bağımsızlığından gelir; iç toplam ise yeni düğüm \(i\) ile her iki taraftan gelen önceki iç toplamları birleştirir.

Çalışılmış Örnek: \(\{2,5,13\}\) üçgeni

\(13\) etiketinin kanonik vektörü

$$13=3^2+2^2,\qquad (x,y)=(3,2)$$

şeklindedir. Bir Bézout bağıntısı

$$1\cdot3+(-1)\cdot2=1$$

olduğundan ebeveyn vektörleri

$$q_1=(1,1),\qquad q_2=(2,1)$$

ve normları

$$\|q_1\|^2=2,\qquad \|q_2\|^2=5$$

olur. Çünkü

$$|\det(q_1,q_2)|=|1\cdot1-1\cdot2|=1,$$

dolayısıyla

$$2\cdot5=3^2+1$$

ve ebeveynler uyumludur. Ayrıca \(13\), hem \(2\) hem de \(5\) ile uyumludur. Böylece \(\{2,5,13\}\) potency değeri

$$2+5+13=20$$

olan bir bileziktir. Aynı mekanizma \(5\) için ebeveynleri \(1\) ve \(2\) yapar; oradan da potency değeri \(8\) olan \(\{1,2,5\}\) üçgeni çıkar. Kenar tabanlı DP tam olarak bu küçük yapıları büyük ölçekte yeniden kurar.

Kod Nasıl Çalışır

C++, Python ve Java uygulamalarının planı aynıdır. Önce \(N\)'e kadar asal tablosu ve tam kareler için doğrudan erişim tablosu kurulur. Sonra her uygun etiket için bir kanonik ilkel vektör üretilir. Asal kuvvetler \(\mathbb Z[i]\) içinde tekrarlı çarpımla, iki katlı aile ise bir kez daha \(1+i\) ile çarpılarak elde edilir.

Daha sonra \(1\)'den büyük her uygun etiket, genişletilmiş Öklid yardımıyla iki ebeveyn etiketine ayrılır. Sadece kendileri de uygun olan ebeveynler tutulur; böylece büyükten küçüğe giden seyrek bir yönlü grafik oluşur.

Son aşamada etiketler büyükten küçüğe taranır. Her aşağı yönlü kenarda önce birikmiş açık parçalar kapatılır ve bunların bilezik potency değerleri sonuca eklenir. Ardından aynı kenara trivial parça yerleştirilir. Bir etikette iki ebeveyn varsa, iki kenar durumu birleştirilir ve ortaya çıkan açık yol daha alttaki ebeveyn kenarına aktarılır. Her bileziğin tek bir maksimum etiketi olduğundan, her biri tam bir kez sayılır.

Karmaşıklık Analizi

\(B(N)\), algoritmanın gerçekten oluşturduğu uygun etiket sayısı olsun. Asal tablosu ve yardımcı tablolar \(O(N)\) bellek kullanır; elek aşaması da bu uygulamalarda \(O(N)\) zamanda çalışır. Her \(p\equiv1\pmod4\) asalı için kod, \(p=a^2+b^2\) ayrışımını tam kare testleriyle arar; mevcut uygulamada bunun maliyeti o asal için en fazla \(O(\sqrt p)\) denemedir. Bundan sonraki etiket üretimi, ebeveyn kenarlarının kurulması ve aşağı yönlü DP aşaması ise toplamda \(O(B(N))\) olur.

Özetle grafik aşaması uygun etiket sayısına göre lineerdir; ön işleme maliyeti ise eleğe, kare tablosuna ve asalların iki kare toplamı ayrışımlarını bulma işine bağlıdır.

Dipnotlar ve Kaynaklar

  1. Problem sayfası: https://projecteuler.net/problem=846
  2. Gauss tamsayıları: Wikipedia - Gauss tamsayıları
  3. İki kare toplamı teoremi: Wikipedia - Fermat's theorem on sums of two squares
  4. Brahmagupta-Fibonacci özdeşliği: Wikipedia - Brahmagupta-Fibonacci identity
  5. Bézout özdeşliği: Wikipedia - Bézout özdeşliği

Resumen del Problema

Se pide la potency total de todas las pulseras mágicas cuyas etiquetas no superan \(N\). Una pulsera se trata como un ciclo en el que dos etiquetas vecinas \(u\) y \(v\) son compatibles exactamente cuando

$$uv=x^2+1$$

para algún entero \(x\). Las implementaciones convierten esto en un problema de grafos: los vértices son las etiquetas admisibles que realmente pueden participar en esa relación, y una pulsera es un ciclo en el grafo de compatibilidad. Su potency es la suma de las etiquetas del ciclo, contando una sola vez hasta rotación y reflexión.

El filtro aritmético usado por las implementaciones conserva exactamente \(1\), \(2\) y las familias \(p^e\) y \(2p^e\) con \(p\equiv1\pmod4\). Toda la solución consiste en construir ese grafo de forma canónica y después sumar los pesos de sus ciclos sin contar ninguno dos veces.

Enfoque Matemático

La idea central es reemplazar cada etiqueta admisible por un vector primitivo de la retícula. Entonces la condición \(uv=x^2+1\) pasa a ser una condición sobre determinantes.

Paso 1: Representar cada etiqueta admisible como suma primitiva de dos cuadrados

Para cada primo \(p\equiv1\pmod4\), elegimos una descomposición

$$p=a^2+b^2.$$

En los enteros gaussianos esto significa

$$p=(a+bi)(a-bi).$$

Al elevar \(a+bi\) a la potencia \(e\) obtenemos un representante primitivo de \(p^e\):

$$z_{p^e}=(a+bi)^e,\qquad N(z_{p^e})=p^e.$$

Multiplicar por \(1+i\) produce la familia duplicada:

$$N\bigl((1+i)z_{p^e}\bigr)=2p^e.$$

Tras tomar valores absolutos y ordenar las coordenadas para que \(x\ge y\ge0\), cada etiqueta almacenada queda asociada a un único vector primitivo canónico

$$n=x^2+y^2,\qquad \gcd(x,y)=1.$$

Por eso las implementaciones trabajan con pares \((x,y)\) en lugar de buscar directamente ecuaciones del tipo cuadrado más \(1\).

Paso 2: Convertir la compatibilidad en una condición de determinante

Supongamos que dos etiquetas están representadas por los vectores primitivos \((u,v)\) y \((s,t)\). Entonces sus valores son

$$m=u^2+v^2,\qquad n=s^2+t^2.$$

La identidad de Brahmagupta-Fibonacci da

$$mn=(us+vt)^2+(ut-vs)^2.$$

Por tanto, si

$$|ut-vs|=1,$$

entonces

$$mn=(us+vt)^2+1,$$

y las dos etiquetas son compatibles. Así, la adyacencia queda codificada por una condición unimodular: dos cuentas vecinas corresponden a vectores primitivos cuyo determinante tiene valor absoluto \(1\).

Esta reformulación es el núcleo estructural del algoritmo. En vez de buscar cuadrados \(x^2+1\), las implementaciones buscan relaciones de determinante \(1\) entre vectores primitivos canónicos.

Paso 3: Recuperar los dos padres menores a partir de coeficientes de Bézout

Tomemos el vector canónico \((x,y)\) de una etiqueta \(n>1\). Como \(\gcd(x,y)=1\), el algoritmo extendido de Euclides produce enteros \(\alpha,\beta\) tales que

$$\alpha x+\beta y=1.$$

A partir de ellos, las implementaciones construyen dos vectores menores

$$q_1=\bigl(|\beta|,|\alpha|\bigr),\qquad q_2=\bigl(x-|\beta|,\ y-|\alpha|\bigr).$$

Estos cumplen

$$q_1+q_2=(x,y),\qquad |\det(q_1,q_2)|=1.$$

Si definimos sus normas por

$$n_1=\|q_1\|^2,\qquad n_2=\|q_2\|^2,$$

obtenemos las dos etiquetas padre de \(n\). Como \(q_1\) y \(q_2\) son unimodulares, los padres son compatibles entre sí. Además, como \((x,y)=q_1+q_2\), la etiqueta actual también es compatible con cada padre, ya que

$$\det(q_1,q_1+q_2)=\det(q_1,q_2),\qquad \det(q_1+q_2,q_2)=\det(q_1,q_2).$$

De este modo, cada etiqueta almacenada queda situada sobre una arista que une dos etiquetas compatibles más pequeñas. El caso especial \(n=2\) se reduce a la arista base única de \(2\) hacia \(1\).

Paso 4: Ver el grafo de compatibilidad como un DAG disperso orientado hacia abajo

Ahora orientamos cada arista de compatibilidad desde la etiqueta mayor hacia la menor. La construcción anterior muestra que cada etiqueta admisible tiene como mucho dos aristas salientes, y ambas apuntan a valores más pequeños. Por eso el grafo orientado es acíclico y muy disperso.

Una pulsera es un ciclo no dirigido en el grafo original de compatibilidad. Una vez orientado, cada ciclo tiene una etiqueta máxima única \(M\). Si se corta el ciclo en \(M\), quedan dos cadenas descendentes desde \(M\) hasta los extremos de una arista compatible más baja. Eso da un lugar canónico donde contar cada pulsera: exactamente una vez, en su máximo único.

Paso 5: Contar caminos abiertos descendentes y cerrarlos en pulseras

Para cada arista orientada \(i\to j\), la programación dinámica guarda dos cantidades:

$$C_{i\to j}=\text{número de fragmentos de camino abierto descendente con arista frontera }(i,j),$$

$$T_{i\to j}=\text{suma de las sumas interiores de etiquetas en esos fragmentos}.$$

Cuando el algoritmo alcanza más tarde la arista \((i,j)\), cada fragmento almacenado puede cerrarse con esa arista y formar una pulsera. La potency añadida es

$$C_{i\to j}(i+j)+T_{i\to j}.$$

Después de cerrar los fragmentos antiguos, la implementación inserta el fragmento trivial de un solo vértice sobre la misma arista para sembrar fragmentos más largos en pasos futuros.

Si una etiqueta \(i\) tiene dos padres \(a\) y \(b\), entonces cada fragmento sobre \(i\to a\) puede combinarse con cada fragmento sobre \(i\to b\). Eso crea un nuevo camino abierto entre \(a\) y \(b\) que pasa por \(i\). Si los dos estados son \((C_a,T_a)\) y \((C_b,T_b)\), la contribución fusionada es

$$C'=C_aC_b,$$

$$T'=C_aC_b\,i+C_aT_b+C_bT_a.$$

El recuento se multiplica porque las elecciones izquierda y derecha son independientes, y la suma interior añade el nuevo vértice \(i\) junto con las sumas interiores ya acumuladas en ambos lados.

Ejemplo Trabajado: el triángulo \(\{2,5,13\}\)

La etiqueta \(13\) tiene como vector canónico

$$13=3^2+2^2,\qquad (x,y)=(3,2).$$

Una relación de Bézout es

$$1\cdot3+(-1)\cdot2=1.$$

Por tanto, los vectores padre son

$$q_1=(1,1),\qquad q_2=(2,1),$$

con normas

$$\|q_1\|^2=2,\qquad \|q_2\|^2=5.$$

Como

$$|\det(q_1,q_2)|=|1\cdot1-1\cdot2|=1,$$

se tiene

$$2\cdot5=3^2+1,$$

de modo que los padres son compatibles. Además, \(13\) es compatible tanto con \(2\) como con \(5\), así que \(\{2,5,13\}\) es una pulsera con potency

$$2+5+13=20.$$

El mismo mecanismo da a \(5\) los padres \(1\) y \(2\), produciendo el triángulo menor \(\{1,2,5\}\) con potency \(8\). Estos ejemplos pequeños son exactamente la estructura que la DP basada en aristas reconstruye a gran escala.

Cómo Funciona el Código

Las implementaciones en C++, Python y Java siguen el mismo esquema. Primero construyen una tabla de primos hasta \(N\) y una tabla de consulta directa para cuadrados perfectos. Después generan las etiquetas admisibles junto con un vector primitivo canónico para cada una. Las potencias primas se obtienen mediante multiplicación repetida en \(\mathbb Z[i]\), y la familia duplicada con una multiplicación adicional por \(1+i\).

Luego, cada etiqueta admisible mayor que \(1\) se descompone en sus dos padres mediante el algoritmo extendido de Euclides. Solo se conservan los padres que también son admisibles, lo que produce un grafo dirigido muy disperso desde etiquetas grandes hacia etiquetas pequeñas.

Por último, las implementaciones recorren las etiquetas en orden descendente. En cada arista hacia abajo primero cierran todos los fragmentos abiertos acumulados y suman sus potencies al resultado. Después siembran el fragmento trivial sobre esa misma arista. Cuando una etiqueta tiene dos padres, los dos estados de arista se fusionan y el nuevo fragmento abierto se envía a la arista inferior que une a esos padres. Como cada pulsera tiene una etiqueta máxima única, cada una se cuenta exactamente una vez.

Análisis de Complejidad

Sea \(B(N)\) el número de etiquetas admisibles materializadas por el algoritmo. La criba de primos y las tablas auxiliares usan \(O(N)\) memoria, y la propia criba corre en \(O(N)\) tiempo en estas implementaciones. Para cada primo \(p\equiv1\pmod4\), el código busca una representación \(p=a^2+b^2\) mediante pruebas de cuadrados exactos; en la implementación actual eso cuesta a lo sumo \(O(\sqrt p)\) intentos para ese primo. Después, generar todas las etiquetas almacenadas, crear las aristas a padres y ejecutar la programación dinámica descendente cuesta \(O(B(N))\).

En consecuencia, la fase de grafo es lineal en el número de etiquetas admisibles, mientras que el preprocesamiento está dominado por la criba, la tabla de cuadrados y la búsqueda de representaciones primas como suma de dos cuadrados.

Notas y Referencias

  1. Página del problema: https://projecteuler.net/problem=846
  2. Enteros gaussianos: Wikipedia - Entero gaussiano
  3. Teorema de Fermat sobre suma de dos cuadrados: Wikipedia - Teorema de Fermat sobre suma de dos cuadrados
  4. Identidad de Brahmagupta-Fibonacci: Wikipedia - Identidad de Brahmagupta-Fibonacci
  5. Identidad de Bézout: Wikipedia - Identidad de Bézout

问题概述

题目要求计算所有标签不超过 \(N\) 的魔法手链的 potency 总和。手链被视为一个环;若相邻两个标签 \(u\) 和 \(v\) 满足

$$uv=x^2+1$$

对某个整数 \(x\) 成立,则这两个标签可以相邻。实现代码先把问题改写成图论问题:顶点是那些真正可能出现在这种关系中的可行标签,手链则是这个兼容图中的一个环。手链的 potency 是环上所有标签之和,并且旋转与镜像只计一次。

实现中保留下来的标签恰好是 \(1\)、\(2\),以及满足 \(p\equiv1\pmod4\) 时的 \(p^e\) 和 \(2p^e\) 两族。整个求解过程就是先以规范方式构造这个图,再在不重复计数的前提下把所有环的权值总和累加出来。

数学方法

核心思想是把每个可行标签替换成一个原始格点向量。这样,原来的 \(uv=x^2+1\) 会变成一个关于行列式的几何条件。

步骤 1:把每个可行标签写成原始的两平方和表示

对每个满足 \(p\equiv1\pmod4\) 的素数,选取一个分解

$$p=a^2+b^2.$$

在高斯整数中,这等价于

$$p=(a+bi)(a-bi).$$

将 \(a+bi\) 提到 \(e\) 次幂,就得到 \(p^e\) 的一个原始表示:

$$z_{p^e}=(a+bi)^e,\qquad N(z_{p^e})=p^e.$$

再乘上 \(1+i\) 就得到翻倍后的那一族:

$$N\bigl((1+i)z_{p^e}\bigr)=2p^e.$$

随后把坐标取绝对值,并调整为 \(x\ge y\ge0\) 的顺序,于是每个被保存的标签都对应一个规范的原始向量

$$n=x^2+y^2,\qquad \gcd(x,y)=1.$$

这就是为什么实现代码不直接在“某个乘积是否等于平方加一”上做搜索,而是改成在 \((x,y)\) 这样的向量上工作。

步骤 2:把兼容条件变成行列式等于 \(\pm1\)

设两个标签的原始向量分别为 \((u,v)\) 和 \((s,t)\),则它们对应的标签值是

$$m=u^2+v^2,\qquad n=s^2+t^2.$$

Brahmagupta-Fibonacci 恒等式给出

$$mn=(us+vt)^2+(ut-vs)^2.$$

因此,如果

$$|ut-vs|=1,$$

那么就有

$$mn=(us+vt)^2+1,$$

这正是题目要求的兼容关系。于是,“两个标签可以相邻”就等价于“它们的原始向量的行列式绝对值等于 \(1\)”这一 unimodular 条件。

这一步是整个算法的结构核心。实现代码不再显式搜索平方数,而是在规范原始向量之间寻找行列式为 \(1\) 的关系。

步骤 3:用 Bézout 系数把较大的标签唯一拆成两个更小的父节点

考虑某个标签 \(n>1\) 的规范向量 \((x,y)\)。由于 \(\gcd(x,y)=1\),扩展欧几里得算法可以找到整数 \(\alpha,\beta\),使得

$$\alpha x+\beta y=1.$$

实现代码利用这组系数构造两个更小的向量

$$q_1=\bigl(|\beta|,|\alpha|\bigr),\qquad q_2=\bigl(x-|\beta|,\ y-|\alpha|\bigr).$$

它们满足

$$q_1+q_2=(x,y),\qquad |\det(q_1,q_2)|=1.$$

再记它们的范数为

$$n_1=\|q_1\|^2,\qquad n_2=\|q_2\|^2,$$

这两个数就是当前标签 \(n\) 的两个父标签。因为 \(q_1\) 与 \(q_2\) 的行列式绝对值是 \(1\),所以 \(n_1\) 与 \(n_2\) 彼此兼容。又因为 \((x,y)=q_1+q_2\),当前标签也会分别与两个父标签兼容,因为

$$\det(q_1,q_1+q_2)=\det(q_1,q_2),\qquad \det(q_1+q_2,q_2)=\det(q_1,q_2).$$

因此,每个被保存的标签都自然地“坐落”在一条连接两个更小兼容标签的边之上。唯一的退化情况是 \(n=2\),此时两个父节点重合,只剩下一条从 \(2\) 指向 \(1\) 的基础边。

步骤 4:把兼容图向下定向后得到一个稀疏 DAG

现在把每条兼容边都从较大的标签指向较小的标签。上面的父节点分解说明:每个可行标签最多只有两条出边,而且都指向更小的值。因此这个有向图没有有向环,并且非常稀疏。

原问题中的手链就是原始兼容图中的一个无向环。把边定向以后,每个这样的环都有唯一的最大标签 \(M\)。如果在 \(M\) 处把这个环“切开”,剩下的就是两条从 \(M\) 向下走到某条更低兼容边两个端点的下降链。于是每个手链都有一个规范的计数位置:就在它的唯一最大标签处计一次即可。

步骤 5:统计开放下降路径,再用边把它们闭合成手链

对每条有向边 \(i\to j\),动态规划保存两个量:

$$C_{i\to j}=\text{边界边为 }(i,j)\text{ 的开放下降路径片段个数},$$

$$T_{i\to j}=\text{这些路径片段内部标签和的总和}.$$

当算法之后真正处理到边 \((i,j)\) 时,之前存下来的每个开放片段都可以通过这条边闭合成一条手链。新增的 potency 为

$$C_{i\to j}(i+j)+T_{i\to j}.$$

关闭旧片段之后,程序再在同一条边上放入一个最简单的单顶点片段,为后续更长路径的拼接提供种子。

如果某个标签 \(i\) 有两个父节点 \(a\) 和 \(b\),那么 \(i\to a\) 上的每个片段都可以与 \(i\to b\) 上的每个片段配对,从而形成一条经过 \(i\)、连接 \(a\) 与 \(b\) 的新开放路径。若两边的状态分别是 \((C_a,T_a)\) 和 \((C_b,T_b)\),则合并后的贡献为

$$C'=C_aC_b,$$

$$T'=C_aC_b\,i+C_aT_b+C_bT_a.$$

这里路径数相乘,是因为左右两侧的选择彼此独立;内部标签和则是把新加入的顶点 \(i\) 与原先两边已经累计的内部和一起算进去。

示例:三角形 \(\{2,5,13\}\)

标签 \(13\) 的规范向量是

$$13=3^2+2^2,\qquad (x,y)=(3,2).$$

一组 Bézout 关系为

$$1\cdot3+(-1)\cdot2=1.$$

于是两个父向量为

$$q_1=(1,1),\qquad q_2=(2,1),$$

其范数分别是

$$\|q_1\|^2=2,\qquad \|q_2\|^2=5.$$

因为

$$|\det(q_1,q_2)|=|1\cdot1-1\cdot2|=1,$$

所以

$$2\cdot5=3^2+1,$$

说明 \(2\) 与 \(5\) 兼容。同时,\(13\) 也分别与 \(2\) 和 \(5\) 兼容,因此 \(\{2,5,13\}\) 构成一条手链,其 potency 是

$$2+5+13=20.$$

同样的方法还会把 \(5\) 拆成父节点 \(1\) 和 \(2\),于是得到更小的三角形 \(\{1,2,5\}\),其 potency 为 \(8\)。这正是边状态动态规划在大规模情形下不断重建的基本模式。

代码如何工作

C++、Python 和 Java 实现遵循同一套流程。首先建立到 \(N\) 为止的素数表,以及一个用于判断完全平方数的直接查表结构。然后为每个可行标签构造一个规范的原始向量。素数幂通过在 \(\mathbb Z[i]\) 中反复相乘得到,而翻倍后的那一族则通过额外再乘一次 \(1+i\) 得到。

接下来,每个大于 \(1\) 的可行标签都会借助扩展欧几里得算法分解成两个父标签。只有那些本身也可行的父标签才会被保留,因此最终得到的是一个从大标签指向小标签的稀疏有向图。

最后,程序按标签从大到小扫描。对每条向下的边,先把之前累计的开放片段闭合掉,并把对应的手链 potency 加到答案中;然后在同一条边上放入最简单的种子片段。当某个标签有两个父节点时,就把两侧边状态合并,再把新形成的开放路径传给更低的一条父边。由于每条手链都有唯一的最大标签,所以每条手链都会且只会被统计一次。

复杂度分析

设 \(B(N)\) 为算法实际构造出的可行标签个数。素数筛和辅助查表需要 \(O(N)\) 内存,而在这些实现中筛法本身运行时间也是 \(O(N)\)。对每个满足 \(p\equiv1\pmod4\) 的素数,代码需要通过完全平方测试去寻找一个表示 \(p=a^2+b^2\);在当前实现里,这一步对该素数最多需要 \(O(\sqrt p)\) 次尝试。之后,所有标签的生成、父边的建立以及向下的动态规划都只需要 \(O(B(N))\) 时间。

因此,真正的图阶段对可行标签数是线性的,而预处理成本主要来自筛法、平方查表以及为相关素数寻找两平方和表示。

注释与参考资料

  1. 题目页面:https://projecteuler.net/problem=846
  2. 高斯整数:Wikipedia - 高斯整数
  3. 费马两平方和定理:Wikipedia - 两个平方和定理
  4. Brahmagupta-Fibonacci 恒等式:Wikipedia - Brahmagupta-Fibonacci identity
  5. Bézout 等式:Wikipedia - 贝祖等式

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

Нужно найти суммарную potency всех магических браслетов, чьи метки не превосходят \(N\). Браслет рассматривается как цикл, в котором две соседние метки \(u\) и \(v\) совместимы тогда и только тогда, когда

$$uv=x^2+1$$

для некоторого целого \(x\). Реализации переводят это в задачу на графе: вершины - это допустимые метки, которые вообще могут участвовать в таком соотношении, а браслет - цикл в графе совместимости. Его potency равна сумме меток на цикле; повороты и отражения считаются одним и тем же браслетом.

Арифметический фильтр в реализациях оставляет ровно \(1\), \(2\), а также семейства \(p^e\) и \(2p^e\) для простых \(p\equiv1\pmod4\). Всё решение состоит в каноническом построении этого графа и последующем суммировании весов всех циклов без двойного счёта.

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

Главная идея - заменить каждую допустимую метку примитивным вектором решётки. Тогда условие \(uv=x^2+1\) превращается в условие на определитель.

Шаг 1: Представить каждую допустимую метку как примитивную сумму двух квадратов

Для каждого простого \(p\equiv1\pmod4\) выбираем разложение

$$p=a^2+b^2.$$

В гауссовых целых это означает

$$p=(a+bi)(a-bi).$$

Степень \(e\) элемента \(a+bi\) даёт примитивного представителя числа \(p^e\):

$$z_{p^e}=(a+bi)^e,\qquad N(z_{p^e})=p^e.$$

Умножение на \(1+i\) порождает удвоенное семейство:

$$N\bigl((1+i)z_{p^e}\bigr)=2p^e.$$

После взятия модулей координат и упорядочивания так, чтобы \(x\ge y\ge0\), каждая сохранённая метка получает единственный канонический примитивный вектор

$$n=x^2+y^2,\qquad \gcd(x,y)=1.$$

Именно поэтому реализации работают с парами \((x,y)\), а не пытаются напрямую перебирать равенства вида «квадрат плюс единица».

Шаг 2: Совместимость становится условием \(|\det|=1\)

Пусть две метки представлены примитивными векторами \((u,v)\) и \((s,t)\). Тогда сами метки равны

$$m=u^2+v^2,\qquad n=s^2+t^2.$$

Из тождества Брахмагупты-Фибоначчи следует

$$mn=(us+vt)^2+(ut-vs)^2.$$

Если выполнено

$$|ut-vs|=1,$$

то получаем

$$mn=(us+vt)^2+1,$$

и, значит, метки совместимы. Итак, соседство в браслете кодируется unimodular-условием: примитивные векторы соседних меток имеют определитель с модулем \(1\).

Это и есть структурное ядро алгоритма. Вместо поиска квадратов \(x^2+1\) реализации ищут отношения с определителем \(1\) между каноническими примитивными векторами.

Шаг 3: Восстановить двух меньших родителей из коэффициентов Бéзу

Возьмём канонический вектор \((x,y)\) для некоторой метки \(n>1\). Так как \(\gcd(x,y)=1\), расширенный алгоритм Евклида даёт целые \(\alpha,\beta\), удовлетворяющие

$$\alpha x+\beta y=1.$$

Из них реализации строят два меньших вектора

$$q_1=\bigl(|\beta|,|\alpha|\bigr),\qquad q_2=\bigl(x-|\beta|,\ y-|\alpha|\bigr).$$

Для них верно

$$q_1+q_2=(x,y),\qquad |\det(q_1,q_2)|=1.$$

Пусть их нормы равны

$$n_1=\|q_1\|^2,\qquad n_2=\|q_2\|^2.$$

Тогда это и есть две родительские метки числа \(n\). Поскольку \(q_1\) и \(q_2\) unimodular, родители совместимы между собой. А так как \((x,y)=q_1+q_2\), текущая метка тоже совместима с каждым из родителей, потому что

$$\det(q_1,q_1+q_2)=\det(q_1,q_2),\qquad \det(q_1+q_2,q_2)=\det(q_1,q_2).$$

Значит, каждая сохранённая метка естественным образом расположена над ребром, соединяющим две меньшие совместимые метки. Особый случай \(n=2\) вырождается в единственное базовое ребро от \(2\) к \(1\).

Шаг 4: После ориентации вниз граф совместимости становится разреженным DAG

Теперь ориентируем каждое ребро совместимости от большей метки к меньшей. Предыдущее построение показывает, что у каждой допустимой метки не более двух исходящих рёбер, и оба ведут к меньшим значениям. Следовательно, ориентированный граф ацикличен и очень разрежен.

Браслет - это неориентированный цикл в исходном графе совместимости. После ориентации у каждого такого цикла есть единственная максимальная метка \(M\). Если разрезать цикл в \(M\), останутся две нисходящие цепочки от \(M\) к концам некоторого более низкого совместимого ребра. Это даёт каноническое место подсчёта: каждый браслет следует учитывать ровно один раз, в его уникальном максимуме.

Шаг 5: Считать открытые нисходящие пути и замыкать их в браслеты

Для каждого ориентированного ребра \(i\to j\) динамика хранит две величины:

$$C_{i\to j}=\text{число открытых нисходящих фрагментов пути с граничным ребром }(i,j),$$

$$T_{i\to j}=\text{сумма внутренних сумм меток по этим фрагментам}.$$

Когда алгоритм позже доходит до ребра \((i,j)\), каждый сохранённый фрагмент можно замкнуть этим ребром в один браслет. Добавляемая potency равна

$$C_{i\to j}(i+j)+T_{i\to j}.$$

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

Если у метки \(i\) два родителя \(a\) и \(b\), то каждый фрагмент над \(i\to a\) можно сочетать с каждым фрагментом над \(i\to b\). Так получается новый открытый путь между \(a\) и \(b\), проходящий через \(i\). Если состояния двух сторон равны \((C_a,T_a)\) и \((C_b,T_b)\), то объединённый вклад имеет вид

$$C'=C_aC_b,$$

$$T'=C_aC_b\,i+C_aT_b+C_bT_a.$$

Число путей перемножается, потому что левая и правая части выбираются независимо, а внутренняя сумма складывается из нового узла \(i\) и уже накопленных внутренних сумм с обеих сторон.

Разобранный пример: треугольник \(\{2,5,13\}\)

Для метки \(13\) канонический вектор равен

$$13=3^2+2^2,\qquad (x,y)=(3,2).$$

Одно из соотношений Бéзу:

$$1\cdot3+(-1)\cdot2=1.$$

Отсюда получаем два родительских вектора

$$q_1=(1,1),\qquad q_2=(2,1),$$

с нормами

$$\|q_1\|^2=2,\qquad \|q_2\|^2=5.$$

Так как

$$|\det(q_1,q_2)|=|1\cdot1-1\cdot2|=1,$$

имеем

$$2\cdot5=3^2+1,$$

то есть родители совместимы. Кроме того, \(13\) совместима и с \(2\), и с \(5\), поэтому \(\{2,5,13\}\) образует браслет с potency

$$2+5+13=20.$$

Тот же механизм даёт для \(5\) родителей \(1\) и \(2\), а значит и меньший треугольник \(\{1,2,5\}\) с potency \(8\). Именно такие локальные структуры реберная динамика восстанавливает в общем случае.

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

Реализации на C++, Python и Java следуют одному и тому же плану. Сначала они строят таблицу простых чисел до \(N\) и таблицу прямого доступа для точных квадратов. Затем генерируют допустимые метки вместе с одним каноническим примитивным вектором для каждой из них. Простые степени получаются повторным умножением в \(\mathbb Z[i]\), а удвоенное семейство - ещё одним умножением на \(1+i\).

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

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

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

Пусть \(B(N)\) - число допустимых меток, которые реально материализует алгоритм. Решето простых и вспомогательные таблицы занимают \(O(N)\) памяти, а само решето работает за \(O(N)\) времени в данных реализациях. Для каждого простого \(p\equiv1\pmod4\) код ищет представление \(p=a^2+b^2\) с помощью проверок точных квадратов; в текущей реализации это требует не более \(O(\sqrt p)\) попыток для данного простого. После этого генерация всех меток, построение родительских рёбер и нисходящая динамика занимают \(O(B(N))\).

Следовательно, собственно графовая часть линейна по числу допустимых меток, а стоимость предобработки определяется решетом, таблицей квадратов и поиском представлений простых как суммы двух квадратов.

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

  1. Страница задачи: https://projecteuler.net/problem=846
  2. Гауссовы целые: Wikipedia - Гауссовы целые
  3. Теорема о сумме двух квадратов: Wikipedia - Сумма двух квадратов
  4. Тождество Брахмагупты-Фибоначчи: Wikipedia - Тождество Брахмагупты-Фибоначчи
  5. Тождество Бéзу: Wikipedia - Теорема Безу

ملخص المسألة

المطلوب هو حساب مجموع قيم potency لكل الأساور السحرية التي لا تتجاوز قيم خرزاتها \(N\). يُنظر إلى السوار على أنه دورة، وتكون القيمتان المتجاورتان \(u\) و\(v\) متوافقتين إذا وفقط إذا تحقق

$$uv=x^2+1$$

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

الترشيح الحسابي المستخدم في التطبيقات يُبقي بالضبط على \(1\) و\(2\) والعائلتين \(p^e\) و\(2p^e\) عندما يكون \(p\equiv1\pmod4\). وجوهر الحل هو بناء هذا الرسم بطريقة معيارية، ثم جمع أوزان جميع دوراته من دون أي ازدواج في العد.

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

الفكرة الأساسية هي استبدال كل قيمة صالحة بمتجه شبكي أولي. عندها تتحول العلاقة \(uv=x^2+1\) إلى شرط هندسي على المحدد.

الخطوة 1: تمثيل كل قيمة صالحة كمجموع أولي لمربعين

لكل أولي \(p\equiv1\pmod4\) نختار تفكيكاً من الشكل

$$p=a^2+b^2.$$

وفي أعداد غاوس هذا يعني

$$p=(a+bi)(a-bi).$$

رفع \(a+bi\) إلى القوة \(e\) يعطينا ممثلاً أولياً للقيمة \(p^e\):

$$z_{p^e}=(a+bi)^e,\qquad N(z_{p^e})=p^e.$$

وضرب هذا العدد في \(1+i\) يولد العائلة المضاعفة:

$$N\bigl((1+i)z_{p^e}\bigr)=2p^e.$$

بعد أخذ القيم المطلقة للإحداثيين وترتيبهما بحيث \(x\ge y\ge0\)، تمتلك كل قيمة مخزنة متجهاً معيارياً أولياً وحيداً

$$n=x^2+y^2,\qquad \gcd(x,y)=1.$$

ولهذا السبب تعمل التطبيقات على الأزواج \((x,y)\) بدلاً من البحث المباشر في معادلات من نوع مربع زائد \(1\).

الخطوة 2: تحويل التوافق إلى شرط \(|\det|=1\)

لنفترض أن قيمتين تمثلان بالمتجهين الأوليين \((u,v)\) و\((s,t)\). عندئذ تكون القيمتان

$$m=u^2+v^2,\qquad n=s^2+t^2.$$

وتعطي هوية Brahmagupta-Fibonacci العلاقة

$$mn=(us+vt)^2+(ut-vs)^2.$$

لذلك إذا تحقق

$$|ut-vs|=1,$$

فإننا نحصل على

$$mn=(us+vt)^2+1,$$

وبالتالي تكون القيمتان متوافقتين. أي إن التجاور بين خرزتين يُشفَّر بشرط أحادي المعيار: المتجهان الأوليان لهما محدد قيمته المطلقة \(1\).

هذه إعادة الصياغة هي القلب البنيوي للخوارزمية. فبدلاً من البحث عن أعداد من شكل \(x^2+1\)، تبحث التطبيقات عن علاقات محددها \(1\) بين متجهات أولية معيارية.

الخطوة 3: استخراج الوالدين الأصغر من معاملات Bézout

خذ المتجه المعياري \((x,y)\) لقيمة ما \(n>1\). بما أن \(\gcd(x,y)=1\)، فإن خوارزمية إقليدس الموسعة تعطي عددين صحيحين \(\alpha,\beta\) بحيث

$$\alpha x+\beta y=1.$$

ومن هذه العلاقة تبني التطبيقات متجهين أصغر هما

$$q_1=\bigl(|\beta|,|\alpha|\bigr),\qquad q_2=\bigl(x-|\beta|,\ y-|\alpha|\bigr).$$

وهما يحققان

$$q_1+q_2=(x,y),\qquad |\det(q_1,q_2)|=1.$$

إذا كانت نورماتهما

$$n_1=\|q_1\|^2,\qquad n_2=\|q_2\|^2,$$

فإن \(n_1\) و\(n_2\) هما القيمتان الوالدتان للقيمة \(n\). وبما أن \(q_1\) و\(q_2\) يحققان شرط المحدد الأحادي، فالوالدان متوافقان فيما بينهما. كذلك بما أن \((x,y)=q_1+q_2\)، فإن القيمة الحالية متوافقة أيضاً مع كل واحد من والديها، لأن

$$\det(q_1,q_1+q_2)=\det(q_1,q_2),\qquad \det(q_1+q_2,q_2)=\det(q_1,q_2).$$

وبذلك تتموضع كل قيمة مخزنة فوق ضلع يصل بين قيمتين أصغر متوافقتين. والحالة الخاصة \(n=2\) تنكمش إلى الضلع الأساسي الوحيد من \(2\) إلى \(1\).

الخطوة 4: بعد توجيه الحواف إلى الأسفل نحصل على DAG متناثر

نوجّه الآن كل ضلع توافق من القيمة الأكبر إلى القيمة الأصغر. ويُظهر بناء الوالدين السابق أن كل قيمة صالحة لها على الأكثر ضلعان خارجان، وكلاهما يتجه إلى قيمة أصغر. لذلك يكون الرسم الموجّه خالياً من الدورات ومتناثراً جداً.

السوار نفسه هو دورة غير موجّهة في رسم التوافق الأصلي. وبعد التوجيه يصبح لكل دورة من هذا النوع قيمة عظمى وحيدة \(M\). وإذا قطعنا الدورة عند \(M\)، بقيت سلسلتان هابطتان من \(M\) إلى طرفي ضلع توافق أدنى. وهذا يعطينا موضعاً معيارياً للعد: كل سوار يُحسب مرة واحدة فقط عند قيمته العظمى الوحيدة.

الخطوة 5: عدّ المسارات المفتوحة الهابطة ثم إغلاقها إلى أساور

لكل ضلع موجّه \(i\to j\)، تخزّن البرمجة الديناميكية مقدارين:

$$C_{i\to j}=\text{count of open descending fragments on }(i,j),$$

$$T_{i\to j}=\text{sum of interior label totals for those fragments}.$$

وعندما تصل الخوارزمية لاحقاً إلى الضلع \((i,j)\)، يمكن إغلاق كل شظية مخزنة بهذا الضلع لتكوين سوار واحد. ومقدار potency المضاف هو

$$C_{i\to j}(i+j)+T_{i\to j}.$$

بعد إغلاق الشظايا القديمة، تضيف التطبيقات الشظية التافهة ذات الرأس الواحد فوق الضلع نفسه لتكون بذرة لمسارات أطول لاحقاً.

إذا كانت للقيمة \(i\) قيمتان والدتان \(a\) و\(b\)، فيمكن ضم كل شظية فوق \(i\to a\) مع كل شظية فوق \(i\to b\). وينتج عن ذلك مسار مفتوح جديد بين \(a\) و\(b\) يمر عبر \(i\). وإذا كانت حالتا الضلعين هما \((C_a,T_a)\) و\((C_b,T_b)\)، فإن مساهمة الدمج تكون

$$C'=C_aC_b,$$

$$T'=C_aC_b\,i+C_aT_b+C_bT_a.$$

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

مثال محلول: المثلث \(\{2,5,13\}\)

للقيمة \(13\) المتجه المعياري

$$13=3^2+2^2,\qquad (x,y)=(3,2).$$

وإحدى علاقات Bézout هي

$$1\cdot3+(-1)\cdot2=1.$$

ومنها نحصل على متجهي الوالدين

$$q_1=(1,1),\qquad q_2=(2,1),$$

ونورماتهما

$$\|q_1\|^2=2,\qquad \|q_2\|^2=5.$$

وبما أن

$$|\det(q_1,q_2)|=|1\cdot1-1\cdot2|=1,$$

فإن

$$2\cdot5=3^2+1,$$

أي إن الوالدين متوافقان. كما أن \(13\) متوافقة مع \(2\) ومع \(5\)، ولذلك فإن \(\{2,5,13\}\) تشكل سواراً بقيمة potency تساوي

$$2+5+13=20.$$

والآلية نفسها تعطي للقيمة \(5\) الوالدين \(1\) و\(2\)، ومن ثم المثلث الأصغر \(\{1,2,5\}\) ذي القيمة \(8\). وهذه الأنماط الصغيرة هي بالضبط ما تعيد البرمجة الديناميكية على الحواف بناءه على المدى الكبير.

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

تتبع تطبيقات C++ وPython وJava الخطة نفسها. في البداية تُبنى قائمة أولية حتى \(N\) مع جدول وصول مباشر للمربعات الكاملة. ثم تُنشأ القيم الصالحة مع متجه أولي معياري واحد لكل قيمة. وتنتج قوى الأعداد الأولية من الضرب المتكرر داخل \(\mathbb Z[i]\)، بينما تنتج العائلة المضاعفة من ضربة إضافية بعامل \(1+i\).

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

وفي المرحلة الأخيرة تُمسح القيم تنازلياً. عند كل ضلع هابط، تُغلق أولاً جميع الشظايا المفتوحة المتراكمة ويُضاف ما يقابلها من قيم potency إلى الناتج. ثم تُزرع الشظية التافهة على الضلع نفسه. وإذا كانت للقيمة والدان، تُدمج حالتا الضلعين ويُمرَّر المسار المفتوح الجديد إلى الضلع الأدنى بين الوالدين. ولأن لكل سوار قيمة عظمى وحيدة، فإن كل سوار يُحسب مرة واحدة فقط.

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

لتكن \(B(N)\) هي عدد القيم الصالحة التي تُنشئها الخوارزمية فعلياً. يستخدم غربال الأعداد الأولية والجداول المساعدة ذاكرة من رتبة \(O(N)\)، كما أن الغربال نفسه يعمل بزمن \(O(N)\) في هذه التطبيقات. ولكل أولي \(p\equiv1\pmod4\)، يبحث الكود عن تمثيل من الشكل \(p=a^2+b^2\) بواسطة اختبارات المربعات الكاملة؛ وفي التنفيذ الحالي تكلف هذه العملية في أسوأ الأحوال \(O(\sqrt p)\) محاولة لذلك الأولي. وبعد ذلك تكون عملية توليد جميع القيم المخزنة، وبناء حواف الوالدين، وتنفيذ البرمجة الديناميكية الهابطة كلها من رتبة \(O(B(N))\).

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

حواشٍ ومراجع

  1. صفحة المسألة: https://projecteuler.net/problem=846
  2. الأعداد الصحيحة الغاوسية: Wikipedia - عدد غاوسي
  3. مبرهنة مجموع مربعين: Wikipedia - Fermat's theorem on sums of two squares
  4. هوية Brahmagupta-Fibonacci: Wikipedia - Brahmagupta-Fibonacci identity
  5. هوية Bézout: Wikipedia - متطابقة بيزو