Detailed mathematical approach

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