Problem 261: Pivotal Square Sums

Project Euler Problem 261 Solution

EulerSolve provides an optimized solution for Project Euler Problem 261, Pivotal Square Sums, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary A positive integer $k$ is called a square-pivot if there exist integers $m>0$ and $n\ge k$ such that $$ (k-m)^2+(k-m+1)^2+\cdots+k^2=(n+1)^2+(n+2)^2+\cdots+(n+m)^2. $$ So the $(m+1)$ consecutive squares ending at $k$ must equal the $m$ consecutive squares starting just after $n$. The task is to sum all distinct square-pivots $$ k\le L, $$ with $L=10^{10}$ in the full problem. The code does not brute-force all triples $(k,m,n)$; it converts the condition into a Pell-equation pipeline. Mathematical Approach 1. Expanding the Square-Sum Identity Start from $$ \sum_{i=0}^{m}(k-i)^2=\sum_{j=1}^{m}(n+j)^2. $$ Using the standard sum-of-squares formula, the left side becomes $$ (m+1)k^2-km(m+1)+\frac{m(m+1)(2m+1)}{6}, $$ while the right side becomes $$ mn^2+m(m+1)n+\frac{m(m+1)(2m+1)}{6}. $$ The cubic tail cancels immediately, leaving the key Diophantine relation $$ (m+1)k(k-m)=m\,n(n+m+1). $$ This is the real starting point of the implementation. 2. Quadratic Viewpoint and the Brute-Force Check For fixed $k$ and $m$, we can treat the previous identity as a quadratic equation in $n$: $$ n^2+(m+1)n-\frac{(m+1)k(k-m)}{m}=0. $$ Therefore $n$ is integral if and only if the discriminant $$ \Delta=(m+1)^2+4\frac{(m+1)k(k-m)}{m} $$ is a perfect square and the parity matches....

Detailed mathematical approach

Problem Summary

A positive integer $k$ is called a square-pivot if there are integers $m\gt0$ and $n\ge k$ such that the $m+1$ consecutive squares ending with $k^2$ have the same sum as the $m$ consecutive squares starting with $(n+1)^2$:

$$(k-m)^2+\cdots+k^2=(n+1)^2+\cdots+(n+m)^2.$$

The statement lists $3^2+4^2=5^2$, $20^2+21^2=29^2$, $21^2+22^2+23^2+24^2$${}=25^2+26^2+27^2$ and $108^2+109^2+110^2$${}=133^2+134^2$, so $4$, $21$, $24$ and $110$ are square-pivots. We must find the sum of all distinct square-pivots $k\le10^{10}$; a pivot that works with several values of $m$ is counted once. ^[1]

The whole solution in one sentence: for each of the $70710$ possible values of $m$, equating the two sums and changing variables turns the condition into a Pell equation $x^2-sy^2=1$, whose solutions can be listed one after another, and the solutions that pass a simple parity and size test give the pivots, whose distinct values up to $10^{10}$ we add up.

Mathematical Approach

We collect small examples, rewrite the equality of the two sums as one quadratic equation and change variables until it becomes Pell's equation, whose solutions are known completely. A check of parity and size then shows exactly which solutions belong to pivots.

1. Small cases

For $m=1$ the condition reads $(k-1)^2+k^2=(n+1)^2$: the legs of a right triangle with integer sides differ by $1$. Such triangles are $(3,4,5)$, $(20,21,29)$, $(119,120,169)$, $(696,697,985)$ and so on, so $4$, $21$, $120$ and $697$ are pivots. ^[2] A brute force over all $k\le1000$, straight from the definition, finds $27$ pivots.

Twenty-one of them are the numbers $2m(m+1)$ for $m=1,\dots,21$, that is $4$, $12$, $24$, $40$, …, $924$, each with $n=k$; for example $10^2+11^2+12^2=13^2+14^2$. This family is classical. ^[3] The other six, $21$, $110$, $120$, $315$, $697$ and $820$, need $n\gt k$, and $684=2\cdot18\cdot19$ has a second representation with $m=4$ and $n=760$. The method must produce both kinds of pivots and count $684$ only once.

2. One equation instead of two sums

We use $1+2+\cdots+m=m(m+1)/2$ and the square pyramidal number $P_m=1^2+2^2+\cdots+m^2$, which equals $m(m+1)(2m+1)/6$. ^[4] Expanding $(k-i)^2$ and $(n+j)^2$ gives

$$\begin{aligned}\sum_{i=0}^{m}(k-i)^2&=(m+1)k^2-m(m+1)k+P_m,\\ \sum_{j=1}^{m}(n+j)^2&=mn^2+m(m+1)n+P_m.\end{aligned}$$

The terms $P_m$ cancel, and the condition becomes

$$(m+1)\,k\,(k-m)=m\,n\,(n+m+1).\qquad(1)$$

For $k=n=4$ and $m=1$, both sides equal $24$.

3. Every pivot is at least $2m(m+1)$

The product $n(n+m+1)$ grows with $n$, so $n\ge k$ and (1) give $(m+1)k(k-m)\ge mk(k+m+1)$. Dividing by $k$ and collecting terms leaves

$$k\ge2m(m+1),$$

with equality exactly when $n=k$. Conversely, $k=n=2m(m+1)$ satisfies (1): then $k-m=m(2m+1)$ and $k+m+1=(m+1)(2m+1)$, and both sides equal $2m^2(m+1)^2(2m+1)$. So every $m\ge1$ has the pivot $2m(m+1)$, the family of Section 1. Two further consequences:

$k-m\ge2m^2+m\gt0$, so the left-hand side never contains the square of $0$ or of a negative number;
a pivot $k\le L$ needs $2m(m+1)\le L$, that is $m\le(\sqrt{2L+1}-1)/2$. For $L=10^{10}$ this means $m\le70710$, since $2\cdot70710\cdot70711=9\,999\,949\,620$ is at most $10^{10}$ and $2\cdot70711\cdot70712$ is not.

4. Completing the squares

Put $A=2n+m+1$ and $B=2k-m$. Then $4n(n+m+1)=A^2-(m+1)^2$ and $4k(k-m)=B^2-m^2$, so four times (1) reads $(m+1)(B^2-m^2)$${}=m(A^2-(m+1)^2)$, or

$$mA^2-(m+1)B^2=m(m+1).\qquad(2)$$

For $k=n=4$ and $m=1$ we get $A=10$, $B=7$ and $100-98=2$.

5. A change of variables

Put $x=A-B=2(n-k+m)+1$, so that $A=B+x$. Substituting this into the left-hand side of (2) and completing the square once more gives

$$\begin{aligned}mA^2-(m+1)B^2&=mx^2+2mxB-B^2\\ &=m(m+1)x^2-(B-mx)^2.\end{aligned}$$

This suggests the new variable $z=B-mx=(m+1)B-mA$. Conversely $B=mx+z$ and $A=(m+1)x+z$, so the integer pairs $(A,B)$ and $(x,z)$ correspond one-to-one, and (2) becomes

$$z^2=m(m+1)\,(x^2-1).\qquad(3)$$

In the example, $x=3$, $z=7-3=4$ and $16=2\cdot8$.

6. The squarefree part and Pell's equation

A number is squarefree if no square of a prime divides it, and every positive integer is, in exactly one way, a squarefree number times a square. ^[5] Write $m(m+1)=sq^2$ with $s$ squarefree; for example $1\cdot2=2\cdot1^2$, $2\cdot3=6\cdot1^2$, $3\cdot4=3\cdot2^2$, $4\cdot5=5\cdot2^2$ and $8\cdot9=2\cdot6^2$. Since $m^2\lt m(m+1)\lt(m+1)^2$, the number $m(m+1)$ is never a square, so $s\ge2$ and $s$ is not a square. And since $sq^2$ is even, $s$ or $q$ is even, so $qs$ is even.

Now (3) says $z^2=sq^2(x^2-1)$. As $q^2$ divides $z^2$, $q$ divides $z$; put $z=qw$, so that $w^2=s(x^2-1)$. Every prime factor of $s$ divides $w^2$, hence $w$, and since $s$ is squarefree, $s$ divides $w$; put $w=sy$. Then $s^2y^2=s(x^2-1)$, that is,

$$x^2-s\,y^2=1,\qquad(4)$$

which is Pell's equation, and $z=qsy$. Going back through Sections 4 and 5,

$$k=\frac{mx+qsy+m}{2},\qquad n=k+\frac{x-1}{2}-m.\qquad(5)$$

For $m=1$, the pivot $4$ has $x=3$, $z=4$ and $y=2$, and indeed $3^2-2\cdot2^2=1$.

7. All solutions of Pell's equation

For a positive integer $s$ that is not a square, (4) has infinitely many solutions. The solution in positive integers with the smallest $x$ is called the fundamental solution $(x_1,y_1)$, and every solution in positive integers is one of its powers: ^[6] ^[7]

$$x_j+y_j\sqrt{s}=\bigl(x_1+y_1\sqrt{s}\bigr)^{j},\qquad j=1,2,3,\dots$$

Multiplying out one more factor gives the recurrence $x_{j+1}=x_1x_j+sy_1y_j$ and $y_{j+1}=x_1y_j+y_1x_j$. For $s=2$ it produces $(3,2)$, $(17,12)$, $(99,70)$, $(577,408)$, and for $s=3$ it produces $(2,1)$, $(7,4)$, $(26,15)$, $(97,56)$.

The fundamental solution comes from the continued fraction of $\sqrt s$, the expansion $a_0+1/(a_1+1/(a_2+\cdots))$ with positive integers $a_i$. Starting from $a_0=\lfloor\sqrt s\rfloor$, $\mu_0=0$ and $d_0=1$, the integer recurrence

$$\begin{aligned}\mu_{i+1}&=d_ia_i-\mu_i,\\ d_{i+1}&=\frac{s-\mu_{i+1}^2}{d_i},\\ a_{i+1}&=\Bigl\lfloor\frac{a_0+\mu_{i+1}}{d_{i+1}}\Bigr\rfloor\end{aligned}$$

produces the $a_i$ without irrational arithmetic. The convergents $p_i/r_i$, the truncations of the expansion, follow from $p_i=a_ip_{i-1}+p_{i-2}$ and $r_i=a_ir_{i-1}+r_{i-2}$ with $p_{-1}=r_{-2}=1$ and $p_{-2}=r_{-1}=0$. ^[8] The fundamental solution is one of the convergents, ^[6] and the numerators $p_i$ increase, so the first convergent with $p_i^2-sr_i^2=1$ is $(x_1,y_1)$. For $s=14$ the convergents are $3/1$, $4/1$, $11/3$ and $15/4$, with $p_i^2-14r_i^2$ equal to $-5$, $2$, $-5$ and $1$, so $(x_1,y_1)=(15,4)$.

This search is always short here. The pair $(2m+1,2q)$ solves (4), because $(2m+1)^2-s(2q)^2$${}=4m^2+4m+1-4m(m+1)=1$; hence $x_1\le2m+1$. For $s=14$, which comes from $m=7$ with $56=14\cdot2^2$, the fundamental solution is this pair itself. The numerators satisfy $p_i\ge p_{i-1}+p_{i-2}$, so they grow at least like the Fibonacci numbers, and $x_1\le2\cdot70710+1=141421$, which is below $F_{27}=196418$, is reached within $25$ convergents. For $L=10^{10}$ no $s$ needs more than $16$.

8. Which solutions are pivots

Theorem. Fix $m\ge1$ and write $m(m+1)=sq^2$ as in Section 6. Formula (5) is a one-to-one correspondence between the pairs $(k,n)$ that satisfy the definition with this $m$ and the solutions $(x,y)$ of (4) in positive integers with $x$ odd and $x\ge2m+1$.

Proof. Let $(k,n)$ satisfy the definition. Then $x=2(n-k+m)+1$ is odd and at least $2m+1$, and by Sections 4 to 6 there is an integer $y$ with (4) and $z=qsy$. Here $y\ne0$, since $x^2\gt1$. If $y$ were negative, then $z\lt0$, and (3) would give $z^2-m^2x^2=m(x^2-m-1)\gt0$, so $|z|\gt mx$ and $2k=mx+z+m\lt m$, which contradicts Section 3. So $y\gt0$.

Conversely, let $(x,y)$ be such a solution, put $z=qsy$ and define $k$ and $n$ by (5). Then $n-k=(x-1)/2-m$ is an integer because $x$ is odd, and $2k=m(x+1)+qsy$ is even because $x+1$ and $qs$ are even. So $k$ and $n$ are integers, $k\gt0$, and $n\ge k$ because $x\ge2m+1$. Since $z^2=sq^2(x^2-1)=m(m+1)(x^2-1)$, the steps of Sections 2 to 5 run backwards, and the two sums are equal. Every step is invertible, so different solutions give different pairs. $\square$

Parity is easy to read off. If $s$ is even, $x^2=1+sy^2$ is odd, so every solution has odd $x$. If $s$ is odd, $x^2-y^2\equiv x^2-sy^2=1\pmod 2$, so $x$ is odd exactly when $y$ is even; for $s=3$ the solutions above alternate between even and odd $x$. The smallest admissible solution is always $(2m+1,2q)$ from Section 7, and (5) turns it into $k=n=2m(m+1)$, the family of Section 3.

9. Worked example: all pivots up to 1000

For $L=1000$ we need $m\le21$, since $2\cdot21\cdot22=924$ and $2\cdot22\cdot23=1012$. Every such $m$ gives $2m(m+1)$ through $(2m+1,2q)$, and the later solutions give:

$m=1$, $s=2$, $q=1$, $k=(x+2y+1)/2$: $(17,12)$, $(99,70)$ and $(577,408)$ give $21$, $120$ and $697$, with $n+1=29$, $169$ and $985$, the hypotenuses of Section 1; the next solution gives $4060$;
$m=2$, $s=6$, $q=1$: $(49,20)$ gives $110$ with $n=132$, and $(485,198)$ already gives $1080$;
$m=3$, $s=3$, $q=2$: $(26,15)$ is skipped because $x$ is even, and (5) would give $k=171/2$; $(97,56)$ gives $315$ with $n=360$;
$m=4$, $s=5$, $q=2$: $(161,72)$ gives $684$ with $n=760$;
$m=8$, $s=2$, $q=6$: $(3,2)$ is skipped because $3\lt17$. It gives $k=28$ and $n=21$, and $20^2+\cdots+28^2=22^2+\cdots+29^2$ is true, but the two blocks overlap; removing the common squares $22^2,\dots,28^2$ leaves $20^2+21^2=29^2$. Then $(17,12)$ gives $144=2\cdot8\cdot9$, and $(99,70)$ gives $820$ with $n=861$.

For every other $m\le21$, the second admissible solution already gives a pivot above $1000$; for $m=5$, $(241,44)$ gives $1265$. That makes $21+7=28$ values, and since $684=2\cdot18\cdot19$ appears twice, there are $27$ distinct pivots with the sum $9167$, exactly the brute-force list of Section 1.

10. Shared Pell equations and repeated pivots

Different values of $m$ can have the same $s$, and then they share the solutions of (4). For $s=2$, the equation $m(m+1)=2q^2$ says that the triangular number $m(m+1)/2$ is the square $q^2$; up to $70710$ this happens for $m=1$, $8$, $49$, $288$, $1681$, $9800$ and $57121$. ^[7] Such sharing is rare: the $70710$ values of $m$ for $L=10^{10}$ have $70561$ different values of $s$, and only $275$ of them share their $s$ with another $m$.

A pivot can also arise from several values of $m$, like $684$. Up to $10^{10}$ this happens exactly three times, each time for two values of $m$: $684$ ($m=4$ and $18$), $231880$ ($m=80$ and $340$) and $74872084$ ($m=1444$ and $6118$); in each case the pivot is $2m(m+1)$ for the larger $m$. So repeated pivots must be removed before the sum is formed.

How the Code Works

The C++ program ^[10] reads --limit= (default $10^{10}$) and --skip-checkpoints in parse_arguments, runs run_checkpoints unless told not to, and prints solve(limit), the sum of the pivots returned by generate_square_pivots.

generate_square_pivots first computes the bound $\lfloor(\sqrt{2L+1}-1)/2\rfloor$ of Section 3 in long double arithmetic, which gives exactly $70710$ for the target. build_primes(m_max + 2) is a linear sieve, which crosses out every composite number once, by its smallest prime factor; ^[9] it fills a smallest-prime-factor table internally but returns only the list of primes, $7004$ of them. squarefree_and_square_part then splits $m(m+1)=sq^2$ by trial division: a prime $p$ with exponent $e$ contributes $p$ to $s$ if $e$ is odd and $p^{\lfloor e/2\rfloor}$ to $q$, and the loop stops as soon as $p^2$ exceeds the remaining cofactor, which is then $1$ or a prime; a prime cofactor goes into $s$. Primes up to $m_{\max}+2$ suffice, because the number being divided is below $(m+1)^2$.

For each $m$ the program also stores x_max${}=\lfloor(2L-m)/m\rfloor$: by (5), $2k\ge mx+m$, so $k\le L$ forces $x\le(2L-m)/m$. Since $qsy$ is about $mx$, this bound is about twice as large as necessary, which is harmless. The triples $(m,q,x_{\max})$ go into an unordered_map keyed by $s$, so that the values of $m$ of Section 10 with the same $s$ share one list of Pell solutions.

For every key $s$, pell_solutions_up_to takes the largest x_max of its group. find_fundamental_pell_limited runs the recurrence of Section 7 in 128-bit integers (its variable m is the $\mu_i$ of Section 7, and p, q are our $p_i$, $r_i$) and returns the first convergent with $p_i^2-sr_i^2=1$. It would give up if a numerator exceeded the bound or after $2\,000\,000$ steps, but neither happens, because $x_1\le2m+1\le$x_max. The recurrence of Section 7 then lists the solutions up to the bound, with 128-bit products, so no step overflows before the comparison with the bound.

Each entry $m$ of the group then runs through this list and applies the theorem of Section 8: it skips $x\lt2m+1$, where $n\lt k$, and even $x$; it forms the numerator $mx+qsy+m$ of (5), skips it if it is odd, and skips $k\gt L$, which occurs for solutions computed for a smaller $m$ of the group or within the margin of x_max. By the parity argument of Section 8 the odd-numerator test never rejects a candidate, and the final test $k\gt0$ always passes; both are harmless. Finally std::sort and std::unique remove the repeated pivots of Section 10, and solve adds the rest in a 64-bit unsigned integer.

The checks use brute_pivots, which scans $m=1,2,\dots$ below $k$ for each $k\le5000$ and stops at the first $m$ that works: it requires $m$ to divide $(m+1)k(k-m)$, forms $c=n(n+m+1)$ from (1), tests with a corrected integer square root whether the discriminant $(m+1)^2+4c$ of $n^2+(m+1)n-c=0$ is a perfect square, and recovers $n$, which must be an integer with $n\ge k$. run_checkpoints requires that this list equals generate_square_pivots(5000), $62$ pivots, and that $4$, $21$, $24$ and $110$ appear among the pivots up to $120$; a failure ends the program with exit code $2$, an unknown argument with exit code $1$.

The Python and Java programs mirror these functions. Python integers are unbounded, and dict.fromkeys removes the repeated pivots from the sorted list. Java uses long, a HashMap for the groups, BigInteger for the step of the Pell recurrence and a LinkedHashSet for the repeated pivots. Its continued-fraction loop works in long, which is safe because no numerator exceeds $x_1\le141421$ there, and the largest numerator $mx+qsy+m$ for the target is about $5.1\cdot10^{14}$, far below $2^{63}$.

Complexity and Verification

Let $M\approx\sqrt{L/2}$ be the largest $m$, so $M=70710$ for the target. The sieve takes $O(M)$ time. Trial division is the most expensive step: for one $m$ it may run through all primes up to about $m/2$, for instance $3762$ loop passes for the prime $m=70621$, where $m+1=2\cdot35311$. So the worst case is $O(M^2/\log M)=O(L/\log L)$, but in the actual run the loop makes $10\,900\,717$ passes, about $154$ per $m$. Factoring $m$ and $m+1$ separately with the smallest-prime-factor table of the sieve would cost only $O(\log m)$ each, but the C++ program already finishes in hundredths of a second.

The Pell part is small. No $s$ needs more than $16$ convergents, $207\,042$ in total; no group stores more than $17$ solutions ($s=3$), $72\,320$ in total; and the entries test $73\,014$ candidates, of which $338$ fail $x\ge2m+1$, $89$ have even $x$, $517$ give $k\gt L$ and $72\,070$ are accepted. Sorting them takes $O(P\log P)$ time for $P=72\,070$, and $72\,067$ distinct pivots remain, $70\,710$ of them of the form $2m(m+1)$. Each group has $O(\log L)$ solutions, because $x$ at least doubles from one solution to the next, so this part costs $O(\sqrt L\log L)$ time. The number $P$ of accepted candidates is even $O(M)$: for fixed $m$, the value $2k-m=mx+qsy$ starts at $4m^2+3m$ and at least doubles from one solution to the next, so $m$ contributes at most $1+\log_2\bigl(L/(2m^2)\bigr)$ pivots, and since $L/2\lt(M+2)^2$ and $M!\ge(M/e)^M$, these bounds add up to $O(M)$. So the primes, the groups and the pivots take $O(\sqrt L)$ memory, and only one group's solutions are stored at a time. The sum is below $72\,067\cdot10^{10}\lt2^{63}$, so 64-bit integers suffice. On our test machine the C++ program, compiled with -O2, finishes in about $0.06$ seconds including its checks, the Java program in about $0.2$ seconds and the Python program in about $1$ second.

The programs check themselves with the brute force and the four sample pivots described above. For this article we also ran checks that are not part of the published programs:

three brute forces agree with an independent implementation of Sections 6 to 8, pivot by pivot: one compares the explicit sums for every $m\le2k$, without the bound of Section 3, for $k\le1200$ ($31$ pivots); one tests (1) through its discriminant for every $m\le2k$ and $k\le20\,000$ ($120$ pivots); one moves two pointers along the explicit sums for every $m$ with $2m(m+1)\le3\cdot10^6$ and every $k\le3\cdot10^6$ ($1337$ pivots);
an instrumented copy of the method confirmed for $L=10^{10}$ that $x_1\le2m+1$ for every $m$, that the first accepted solution of every $m$ is $(2m+1,2q)$, that no group is empty and that the odd-numerator test never rejects a candidate, and it produced the counts above;
the three published programs print the same number for $L=10^{10}$, and for $L=120$, $1000$, $5000$, $10^6$ and $10^8$ they all print $587$, $9167$, $106277$, $255978835$ and $243309539837$, in agreement with the independent implementation.

Footnotes and References

The references below supply the background facts used above; the reduction to Pell's equation, the correspondence of Section 8, the bounds and the operation counts are derived in this article.

Project Euler 261 — Pivotal Square Sums. The official statement calls a positive integer k a square-pivot if there are integers m > 0 and n ≥ k with (k − m)² + ⋯ + k² = (n + 1)² + ⋯ + (n + m)², gives the examples 4, 21, 24 and 110, and asks for the sum of all distinct square-pivots ≤ 10¹⁰.
OEIS A046090. For the Pythagorean triples (X, X + 1, Z) ordered by Z, the sequence lists the values X + 1: 1, 4, 21, 120, 697, 4060, 23661, …, with the recurrence a(n) = 6a(n − 1) − a(n − 2) − 2. Apart from the degenerate first term, these are the pivots with m = 1.
OEIS A046092. The sequence a(n) = 2n(n + 1), four times the triangular numbers, starts 0, 4, 12, 24, 40, 60; a comment notes that the squares of the n + 1 consecutive numbers from a(n) − n to a(n) have the same sum as the squares of the n numbers that follow a(n), for example 10² + 11² + 12² = 13² + 14².
Wikipedia — Square pyramidal number. The sum of the first n squares is Pₙ = n(n + 1)(2n + 1)/6, a special case of Faulhaber's formula that can be proved by induction.
Wikipedia — Square-free integer. The section on square-free factors states that every positive integer can be written in a unique way as the product of a largest possible square and a square-free integer, its square-free part.
Wikipedia — Pell's equation. For a positive non-square n, x² − ny² = 1 has infinitely many integer solutions (Lagrange). The fundamental solution, the positive solution with the smallest x, corresponds to a convergent of the continued fraction of √n, and all remaining solutions follow from xₖ + yₖ√n = (x₁ + y₁√n)ᵏ, which gives the recurrence used in Section 7.
K. Conrad — Pell's Equation, I. Theorem 5.3 proves that, with (x₁, y₁) the positive solution of x² − dy² = 1 with the smallest y (equivalently the smallest x), every solution in positive integers is given by x + y√d = (x₁ + y₁√d)ᵏ with k ≥ 1. Section 3 shows that the triangular numbers that are squares correspond to the solutions of x² − 2y² = 1, with x = 2m + 1 for the triangular number m(m + 1)/2.
Wikipedia — Periodic continued fraction. The section on the canonical form gives the integer recurrence for the continued fraction of √S used in Section 7 and notes that, for every positive non-square integer n, the continued fraction of √n repeats after its first term a₀. Wikipedia — Simple continued fraction gives the recurrences hₙ = aₙhₙ₋₁ + hₙ₋₂ and kₙ = aₙkₙ₋₁ + kₙ₋₂ for the convergents, with h₋₁ = k₋₂ = 1 and h₋₂ = k₋₁ = 0.
cp-algorithms — Linear Sieve. The linear sieve computes the smallest prime factor of every number up to n and the list of primes in O(n) time, and it reaches every composite number exactly once.
C++ — Euler261.cpp. The linked immutable C++ revision contains generate_square_pivots, find_fundamental_pell_limited, pell_solutions_up_to, brute_pivots and run_checkpoints. It is the implementation record for this article.

Problem 261 source code

C++

#include <algorithm>
#include <cstdint>
#include <cmath>
#include <iostream>
#include <string>
#include <unordered_map>
#include <utility>
#include <vector>

namespace {

using i64 = long long;
using u64 = std::uint64_t;
using i128 = __int128_t;

struct Options {
    u64 limit = 10'000'000'000ULL;
    bool run_checkpoints = true;
};

struct MEntry {
    int m = 0;
    u64 q = 0;
    u64 x_max = 0;
};

struct PellSolution {
    u64 x = 0;
    u64 y = 0;
};

bool parse_u64_after_prefix(const std::string& arg, const std::string& prefix, u64& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }
    u64 parsed = 0;
    for (char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        parsed = parsed * 10ULL + static_cast<u64>(c - '0');
    }
    value = parsed;
    return true;
}

bool parse_arguments(int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        if (parse_u64_after_prefix(arg, "--limit=", options.limit)) {
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.limit >= 1;
}

std::vector<int> build_primes(const int n) {
    std::vector<int> spf(static_cast<std::size_t>(n + 1), 0);
    std::vector<int> primes;
    primes.reserve(n / 10);
    for (int i = 2; i <= n; ++i) {
        if (spf[static_cast<std::size_t>(i)] == 0) {
            spf[static_cast<std::size_t>(i)] = i;
            primes.push_back(i);
        }
        for (int p : primes) {
            const i64 v = static_cast<i64>(i) * p;
            if (v > n || p > spf[static_cast<std::size_t>(i)]) {
                break;
            }
            spf[static_cast<std::size_t>(v)] = p;
        }
    }
    return primes;
}

std::pair<u64, u64> squarefree_and_square_part(u64 n, const std::vector<int>& primes) {
    u64 s = 1;
    u64 q = 1;
    u64 x = n;

    for (int p : primes) {
        const u64 pp = static_cast<u64>(p);
        if (pp * pp > x) {
            break;
        }
        if (x % pp != 0) {
            continue;
        }
        int e = 0;
        while (x % pp == 0) {
            x /= pp;
            ++e;
        }
        if ((e & 1) != 0) {
            s *= pp;
        }
        for (int i = 0; i < e / 2; ++i) {
            q *= pp;
        }
    }

    if (x > 1) {
        s *= x;
    }

    return {s, q};
}

u64 floor_sqrt_u64(u64 x) {
    u64 r = static_cast<u64>(std::sqrt(static_cast<long double>(x)));
    while ((r + 1) * (r + 1) <= x) {
        ++r;
    }
    while (r * r > x) {
        --r;
    }
    return r;
}

bool find_fundamental_pell_limited(u64 s, u64 x_limit, u64& x1, u64& y1) {
    const u64 a0 = floor_sqrt_u64(s);
    if (a0 * a0 == s) {
        return false;
    }

    i128 m = 0;
    i128 d = 1;
    i128 a = static_cast<i128>(a0);

    i128 p_prev2 = 0;
    i128 p_prev1 = 1;
    i128 q_prev2 = 1;
    i128 q_prev1 = 0;

    for (int iter = 0; iter < 2'000'000; ++iter) {
        const i128 p = a * p_prev1 + p_prev2;
        const i128 q = a * q_prev1 + q_prev2;

        const i128 diff = p * p - static_cast<i128>(s) * q * q;
        if (diff == 1 && p > 1 && p <= static_cast<i128>(x_limit)) {
            x1 = static_cast<u64>(p);
            y1 = static_cast<u64>(q);
            return true;
        }

        if (p > static_cast<i128>(x_limit)) {
            return false;
        }

        p_prev2 = p_prev1;
        p_prev1 = p;
        q_prev2 = q_prev1;
        q_prev1 = q;

        m = d * a - m;
        d = (static_cast<i128>(s) - m * m) / d;
        a = (static_cast<i128>(a0) + m) / d;
    }

    return false;
}

std::vector<PellSolution> pell_solutions_up_to(u64 s, u64 x_limit) {
    u64 x1 = 0;
    u64 y1 = 0;
    if (!find_fundamental_pell_limited(s, x_limit, x1, y1)) {
        return {};
    }

    std::vector<PellSolution> out;
    u64 x = x1;
    u64 y = y1;

    while (x <= x_limit) {
        out.push_back({x, y});

        const i128 nx = static_cast<i128>(x1) * x + static_cast<i128>(s) * y1 * y;
        const i128 ny = static_cast<i128>(x1) * y + static_cast<i128>(y1) * x;
        if (nx > static_cast<i128>(x_limit)) {
            break;
        }
        x = static_cast<u64>(nx);
        y = static_cast<u64>(ny);
    }

    return out;
}

std::vector<u64> generate_square_pivots(const u64 limit) {
    const int m_max = static_cast<int>((std::sqrt(2.0L * static_cast<long double>(limit) + 1.0L) - 1.0L) / 2.0L);
    const std::vector<int> primes = build_primes(m_max + 2);

    std::unordered_map<u64, std::vector<MEntry>> grouped;
    grouped.reserve(static_cast<std::size_t>(m_max * 2));

    for (int m = 1; m <= m_max; ++m) {
        const u64 n = static_cast<u64>(m) * static_cast<u64>(m + 1);
        const auto [s, q] = squarefree_and_square_part(n, primes);
        const u64 x_max = (2ULL * limit - static_cast<u64>(m)) / static_cast<u64>(m);
        grouped[s].push_back({m, q, x_max});
    }

    std::vector<u64> pivots;
    pivots.reserve(1 << 20);

    for (auto& kv : grouped) {
        const u64 s = kv.first;
        std::vector<MEntry>& entries = kv.second;

        u64 x_max = 0;
        for (const MEntry& e : entries) {
            x_max = std::max(x_max, e.x_max);
        }

        const std::vector<PellSolution> sols = pell_solutions_up_to(s, x_max);
        if (sols.empty()) {
            continue;
        }

        for (const MEntry& e : entries) {
            const u64 min_d = static_cast<u64>(2 * e.m + 1);
            const i128 m128 = static_cast<i128>(e.m);
            const i128 qs128 = static_cast<i128>(e.q) * static_cast<i128>(s);

            for (const PellSolution& sol : sols) {
                if (sol.x < min_d) {
                    continue;
                }
                if ((sol.x & 1ULL) == 0ULL) {
                    continue;
                }

                const i128 numerator = m128 * static_cast<i128>(sol.x) +
                                       qs128 * static_cast<i128>(sol.y) + m128;
                if ((numerator & 1) != 0) {
                    continue;
                }

                const i128 k = numerator / 2;
                if (k > static_cast<i128>(limit)) {
                    continue;
                }
                if (k > 0) {
                    pivots.push_back(static_cast<u64>(k));
                }
            }
        }
    }

    std::sort(pivots.begin(), pivots.end());
    pivots.erase(std::unique(pivots.begin(), pivots.end()), pivots.end());
    return pivots;
}

u64 solve(const u64 limit) {
    const std::vector<u64> pivots = generate_square_pivots(limit);
    u64 sum = 0;
    for (u64 k : pivots) {
        sum += k;
    }
    return sum;
}

std::vector<u64> brute_pivots(const int limit) {
    std::vector<u64> pivots;
    for (int k = 1; k <= limit; ++k) {
        bool ok = false;
        for (int m = 1; m < k && !ok; ++m) {
            const i128 lhs = static_cast<i128>(m + 1) * k * (k - m);
            if (lhs % m != 0) {
                continue;
            }
            const i128 c = lhs / m;
            const i128 d = static_cast<i128>(m + 1) * (m + 1) + 4 * c;
            const i64 s = static_cast<i64>(std::sqrt(static_cast<long double>(d)));
            i64 root = s;
            while (static_cast<i128>(root + 1) * (root + 1) <= d) {
                ++root;
            }
            while (static_cast<i128>(root) * root > d) {
                --root;
            }
            if (static_cast<i128>(root) * root != d) {
                continue;
            }
            const i64 numer = -static_cast<i64>(m + 1) + root;
            if ((numer & 1) != 0) {
                continue;
            }
            const i64 n = numer / 2;
            if (n >= k) {
                ok = true;
            }
        }
        if (ok) {
            pivots.push_back(static_cast<u64>(k));
        }
    }
    return pivots;
}

bool run_checkpoints() {
    const auto fast_small = generate_square_pivots(5000);
    const auto brute_small = brute_pivots(5000);
    if (fast_small != brute_small) {
        std::cerr << "Checkpoint failed for limit=5000" << '\n';
        return false;
    }

    const auto sample = generate_square_pivots(120);
    const std::vector<u64> required = {4, 21, 24, 110};
    for (u64 x : required) {
        if (!std::binary_search(sample.begin(), sample.end(), x)) {
            std::cerr << "Checkpoint failed: missing sample pivot " << x << '\n';
            return false;
        }
    }

    return true;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }
    if (options.run_checkpoints && !run_checkpoints()) {
        return 2;
    }
    std::cout << solve(options.limit) << '\n';
    return 0;
}

Python

#!/usr/bin/env python3
import sys
import math

def build_primes(n):
    spf = [0] * (n + 1)
    primes = []
    for i in range(2, n + 1):
        if spf[i] == 0:
            spf[i] = i
            primes.append(i)
        for p in primes:
            v = i * p
            if v > n or p > spf[i]:
                break
            spf[v] = p
    return primes

def squarefree_and_square_part(n, primes):
    s = 1
    q = 1
    x = n
    for p in primes:
        pp = p
        if pp * pp > x:
            break
        if x % pp != 0:
            continue
        e = 0
        while x % pp == 0:
            x //= pp
            e += 1
        if e % 2 != 0:
            s *= pp
        for _ in range(e // 2):
            q *= pp
    if x > 1:
        s *= x
    return s, q

def floor_sqrt_u64(x):
    r = int(math.sqrt(float(x)))
    while (r + 1) * (r + 1) <= x:
        r += 1
    while r * r > x:
        r -= 1
    return r

def find_fundamental_pell_limited(s, x_limit):
    a0 = floor_sqrt_u64(s)
    if a0 * a0 == s:
        return None
    m = 0
    d = 1
    a = a0
    p_prev2 = 0
    p_prev1 = 1
    q_prev2 = 1
    q_prev1 = 0
    for _ in range(2_000_000):
        p = a * p_prev1 + p_prev2
        q = a * q_prev1 + q_prev2
        diff = p * p - s * q * q
        if diff == 1 and p > 1 and p <= x_limit:
            return p, q
        if p > x_limit:
            return None
        p_prev2 = p_prev1
        p_prev1 = p
        q_prev2 = q_prev1
        q_prev1 = q
        m = d * a - m
        d = (s - m * m) // d
        a = (a0 + m) // d
    return None

def pell_solutions_up_to(s, x_limit):
    res = find_fundamental_pell_limited(s, x_limit)
    if res is None:
        return []
    x1, y1 = res
    solutions = []
    x = x1
    y = y1
    while x <= x_limit:
        solutions.append((x, y))
        nx = x1 * x + s * y1 * y
        ny = x1 * y + y1 * x
        if nx > x_limit:
            break
        x = nx
        y = ny
    return solutions

def generate_square_pivots(limit):
    m_max = int((math.sqrt(2.0 * float(limit) + 1.0) - 1.0) / 2.0)
    primes = build_primes(m_max + 2)
    grouped = {}
    for m in range(1, m_max + 1):
        n = m * (m + 1)
        s, q = squarefree_and_square_part(n, primes)
        x_max = (2 * limit - m) // m
        if s not in grouped:
            grouped[s] = []
        grouped[s].append((m, q, x_max))
    pivots = []
    for s, entries in grouped.items():
        x_max = max(e[2] for e in entries)
        sols = pell_solutions_up_to(s, x_max)
        if not sols:
            continue
        for m, q, _ in entries:
            min_d = 2 * m + 1
            m128 = m
            qs128 = q * s
            for x, y in sols:
                if x < min_d:
                    continue
                if x % 2 == 0:
                    continue
                numerator = m128 * x + qs128 * y + m128
                if numerator % 2 != 0:
                    continue
                k = numerator // 2
                if k > limit:
                    continue
                if k > 0:
                    pivots.append(k)
    pivots.sort()
    return list(dict.fromkeys(pivots))

def solve(limit):
    pivots = generate_square_pivots(limit)
    return sum(pivots)

def brute_pivots(limit):
    pivots = []
    for k in range(1, limit + 1):
        ok = False
        for m in range(1, k):
            lhs = (m + 1) * k * (k - m)
            if lhs % m != 0:
                continue
            c = lhs // m
            d = (m + 1) * (m + 1) + 4 * c
            s = int(math.sqrt(float(d)))
            root = s
            while (root + 1) * (root + 1) <= d:
                root += 1
            while root * root > d:
                root -= 1
            if root * root != d:
                continue
            numer = -(m + 1) + root
            if numer % 2 != 0:
                continue
            n = numer // 2
            if n >= k:
                ok = True
                break
        if ok:
            pivots.append(k)
    return pivots

def run_checkpoints():
    fast_small = generate_square_pivots(5000)
    brute_small = brute_pivots(5000)
    if fast_small != brute_small:
        print("Checkpoint failed for limit=5000", file=sys.stderr)
        return False
    sample = generate_square_pivots(120)
    required = [4, 21, 24, 110]
    for x in required:
        if x not in sample:
            print(f"Checkpoint failed: missing sample pivot {x}", file=sys.stderr)
            return False
    return True

def parse_arguments(args):
    options = {
        "limit": 10_000_000_000,
        "run_checkpoints": True
    }
    for arg in args:
        if arg == "--skip-checkpoints":
            options["run_checkpoints"] = False
        elif arg.startswith("--limit="):
            try:
                options["limit"] = int(arg[8:])
            except ValueError:
                print(f"Unknown argument: {arg}", file=sys.stderr)
                return None
        else:
            print(f"Unknown argument: {arg}", file=sys.stderr)
            return None
    if options["limit"] < 1:
        return None
    return options

def main():
    args = sys.argv[1:]
    options = parse_arguments(args)
    if options is None:
        return 1
    if options["run_checkpoints"] and not run_checkpoints():
        return 2
    print(solve(options["limit"]))
    return 0

if __name__ == "__main__":
    sys.exit(main())

Java

import java.util.*;
import java.math.*;

class Euler261 {
    private static class Options {
        long limit = 10_000_000_000L;
        boolean runCheckpoints = true;
    }

    private static class MEntry {
        int m;
        long q;
        long xMax;
        
        MEntry(int m, long q, long xMax) {
            this.m = m;
            this.q = q;
            this.xMax = xMax;
        }
    }

    private static class PellSolution {
        long x;
        long y;
        
        PellSolution(long x, long y) {
            this.x = x;
            this.y = y;
        }
    }

    private static List<Integer> buildPrimes(int n) {
        int[] spf = new int[n + 1];
        List<Integer> primes = new ArrayList<>();
        for (int i = 2; i <= n; i++) {
            if (spf[i] == 0) {
                spf[i] = i;
                primes.add(i);
            }
            for (int p : primes) {
                long v = (long)i * p;
                if (v > n || p > spf[i]) break;
                spf[(int)v] = p;
            }
        }
        return primes;
    }

    private static long[] squarefreeAndSquarePart(long n, List<Integer> primes) {
        long s = 1;
        long q = 1;
        long x = n;
        for (int p : primes) {
            long pp = p;
            if (pp * pp > x) break;
            if (x % pp != 0) continue;
            int e = 0;
            while (x % pp == 0) {
                x /= pp;
                e++;
            }
            if ((e & 1) != 0) s *= pp;
            for (int i = 0; i < e/2; i++) q *= pp;
        }
        if (x > 1) s *= x;
        return new long[]{s, q};
    }

    private static long floorSqrtU64(long x) {
        long r = (long)Math.sqrt((double)x);
        while ((r + 1) * (r + 1) <= x) r++;
        while (r * r > x) r--;
        return r;
    }

    private static long[] findFundamentalPellLimited(long s, long xLimit) {
        long a0 = floorSqrtU64(s);
        if (a0 * a0 == s) return null;
        
        long m = 0;
        long d = 1;
        long a = a0;
        
        long pPrev2 = 0;
        long pPrev1 = 1;
        long qPrev2 = 1;
        long qPrev1 = 0;
        
        for (int iter = 0; iter < 2_000_000; iter++) {
            long p = a * pPrev1 + pPrev2;
            long q = a * qPrev1 + qPrev2;
            
            long diff = p * p - s * q * q;
            if (diff == 1 && p > 1 && p <= xLimit) {
                return new long[]{p, q};
            }
            
            if (p > xLimit) return null;
            
            pPrev2 = pPrev1;
            pPrev1 = p;
            qPrev2 = qPrev1;
            qPrev1 = q;
            
            m = d * a - m;
            d = (s - m * m) / d;
            a = (a0 + m) / d;
        }
        return null;
    }

    private static List<PellSolution> pellSolutionsUpTo(long s, long xLimit) {
        long[] fundamental = findFundamentalPellLimited(s, xLimit);
        if (fundamental == null) return new ArrayList<>();
        
        long x1 = fundamental[0];
        long y1 = fundamental[1];
        
        List<PellSolution> solutions = new ArrayList<>();
        long x = x1;
        long y = y1;
        
        while (x <= xLimit) {
            solutions.add(new PellSolution(x, y));
            BigInteger nx = BigInteger.valueOf(x1).multiply(BigInteger.valueOf(x))
                          .add(BigInteger.valueOf(s).multiply(BigInteger.valueOf(y1)).multiply(BigInteger.valueOf(y)));
            BigInteger ny = BigInteger.valueOf(x1).multiply(BigInteger.valueOf(y))
                          .add(BigInteger.valueOf(y1).multiply(BigInteger.valueOf(x)));
            if (nx.compareTo(BigInteger.valueOf(xLimit)) > 0) break;
            x = nx.longValue();
            y = ny.longValue();
        }
        return solutions;
    }

    private static List<Long> generateSquarePivots(long limit) {
        int mMax = (int)((Math.sqrt(2.0 * limit + 1.0) - 1.0) / 2.0);
        List<Integer> primes = buildPrimes(mMax + 2);
        Map<Long, List<MEntry>> grouped = new HashMap<>();
        
        for (int m = 1; m <= mMax; m++) {
            long n = (long)m * (m + 1);
            long[] sqParts = squarefreeAndSquarePart(n, primes);
            long s = sqParts[0];
            long q = sqParts[1];
            long xMax = (2 * limit - m) / m;
            grouped.computeIfAbsent(s, k -> new ArrayList<>()).add(new MEntry(m, q, xMax));
        }
        
        List<Long> pivots = new ArrayList<>();
        for (Map.Entry<Long, List<MEntry>> entry : grouped.entrySet()) {
            long s = entry.getKey();
            List<MEntry> entries = entry.getValue();
            
            long xMax = 0;
            for (MEntry e : entries) {
                xMax = Math.max(xMax, e.xMax);
            }
            
            List<PellSolution> sols = pellSolutionsUpTo(s, xMax);
            if (sols.isEmpty()) continue;
            
            for (MEntry e : entries) {
                long minD = 2L * e.m + 1;
                long m128 = e.m;
                long qs128 = e.q * s;
                
                for (PellSolution sol : sols) {
                    if (sol.x < minD) continue;
                    if ((sol.x & 1L) == 0L) continue;
                    
                    long numerator = m128 * sol.x + qs128 * sol.y + m128;
                    if ((numerator & 1L) != 0L) continue;
                    
                    long k = numerator / 2;
                    if (k > limit) continue;
                    if (k > 0) pivots.add(k);
                }
            }
        }
        
        Collections.sort(pivots);
        return new ArrayList<>(new LinkedHashSet<>(pivots));
    }

    private static long solve(long limit) {
        List<Long> pivots = generateSquarePivots(limit);
        long sum = 0;
        for (Long k : pivots) {
            sum += k;
        }
        return sum;
    }

    private static List<Long> brutePivots(int limit) {
        List<Long> pivots = new ArrayList<>();
        for (int k = 1; k <= limit; k++) {
            boolean ok = false;
            for (int m = 1; m < k && !ok; m++) {
                long lhs = (long)(m + 1) * k * (k - m);
                if (lhs % m != 0) continue;
                long c = lhs / m;
                long d = (long)(m + 1) * (m + 1) + 4 * c;
                long s = (long)Math.sqrt((double)d);
                long root = s;
                while ((root + 1) * (root + 1) <= d) root++;
                while (root * root > d) root--;
                if (root * root != d) continue;
                long numer = -(m + 1) + root;
                if ((numer & 1L) != 0) continue;
                long n = numer / 2;
                if (n >= k) ok = true;
            }
            if (ok) pivots.add((long)k);
        }
        return pivots;
    }

    private static boolean runCheckpoints() {
        List<Long> fastSmall = generateSquarePivots(5000);
        List<Long> bruteSmall = brutePivots(5000);
        if (!fastSmall.equals(bruteSmall)) {
            System.err.println("Checkpoint failed for limit=5000");
            return false;
        }
        
        List<Long> sample = generateSquarePivots(120);
        long[] required = {4, 21, 24, 110};
        for (long x : required) {
            if (!sample.contains(x)) {
                System.err.println("Checkpoint failed: missing sample pivot " + x);
                return false;
            }
        }
        return true;
    }

    private static Options parseArguments(String[] args) {
        Options options = new Options();
        for (String arg : args) {
            if (arg.equals("--skip-checkpoints")) {
                options.runCheckpoints = false;
            } else if (arg.startsWith("--limit=")) {
                try {
                    options.limit = Long.parseLong(arg.substring(8));
                } catch (NumberFormatException e) {
                    System.err.println("Unknown argument: " + arg);
                    return null;
                }
            } else {
                System.err.println("Unknown argument: " + arg);
                return null;
            }
        }
        if (options.limit < 1) return null;
        return options;
    }

    public static void main(String[] args) {
        Options options = parseArguments(args);
        if (options == null) {
            System.exit(1);
            return;
        }
        if (options.runCheckpoints && !runCheckpoints()) {
            System.exit(2);
            return;
        }
        System.out.println(solve(options.limit));
    }
}