Solution

Answer: 238890850232021

Let

$$S(a,b)=a^2+(a+1)^2+\cdots+b^2$$

and suppose $k$ is a square-pivot. Then there exist integers $m>0$ and $n\ge k$ such that

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

We must find the sum of all distinct such $k\le 10^{10}$.

1. Mathematical analysis

The sum of squares formula is

$$\sum_{i=1}^{t} i^2=\frac{t(t+1)(2t+1)}6.$$

Therefore

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

Writing this with prefix sums and simplifying gives

$$k^2(m+1)-km(m+1)-mn(m+n+1)=0.$$

Dividing by $m(m+1)$ gives the more symmetric form

$$\frac{k^2}{m}-k-n-\frac{n^2}{m+1}=0. \tag{1}$$

Bounding $m$

Since $n\ge k$,

$$\frac{k^2}{m}-2k-\frac{k^2}{m+1}\ge0.$$

After simplification:

$$k\ge 2m(m+1).$$

Hence for $k\le N=10^{10}$,

$$2m(m+1)\le N.$$

$$m \le \frac{-1+\sqrt{1+2N}}2 \approx 70710.$$

Thus only about $7\times 10^4$ values of $m$ need to be checked.

2. Transforming to a Pell equation

Introduce integers $x,t$ by

$$\begin{cases} 2k = m(x+1)+t,\ 2n = (m+1)(x-1)+t. \end{cases}$$

Substituting into (1) and simplifying yields

$$m(m+1)(x^2-1)=t^2. \tag{2}$$

Now write

$$m=r p^2,\qquad m+1=s q^2,$$

where $r,s$ are squarefree.

Because consecutive integers are coprime, $r,s$ are coprime too.

Equation (2) becomes

$$r s p^2 q^2 (x^2-1)=t^2.$$

Hence $rpsq\mid t$. Let

$$t=rpsq,y.$$

Then

$$x^2-rs,y^2=1. \tag{3}$$

This is a Pell equation.

So for each $m$, we:

Compute the squarefree parts $r,s$ of $m,m+1$.
Solve the Pell equation

$$x^2-Dy^2=1, \qquad D=rs.$$ 3. Generate all Pell solutions. 4. Recover

$$k=\frac{m(x+1)+t}{2}, \qquad t=rpsq,y.$$ 5. Keep all distinct $k\le 10^{10}$.

3. Python implementation

from math import isqrt

N = 10**10

# ---------------------------------------------------------
# Compute squarefree part and square part of every integer
# n = sf[n] * sq[n]^2
# ---------------------------------------------------------

M = (-1 + isqrt(1 + 2 * N)) // 2

sf = [0] * (M + 3)   # squarefree part
sq = [0] * (M + 3)   # square root of square part

for n in range(1, M + 3):

    x = 1
    y = 1
    t = n

    d = 2
    while d * d <= t:

        if t % d == 0:

            e = 0
            while t % d == 0:
                t //= d
                e += 1

            # odd exponent contributes to squarefree part
            if e & 1:
                x *= d

            # floor(e/2) contributes to square part
            y *= d ** (e >> 1)

        d += 1 if d == 2 else 2

    if t > 1:
        x *= t

    sf[n] = x
    sq[n] = y

# ---------------------------------------------------------
# Fundamental Pell solution x^2 - D y^2 = 1
# via continued fractions
# ---------------------------------------------------------

def pell_fundamental(D):

    a0 = isqrt(D)

    if a0 * a0 == D:
        return None

    m = 0
    d = 1
    a = a0

    num_prev, num = 1, a
    den_prev, den = 0, 1

    while num * num - D * den * den != 1:

        m = d * a - m
        d = (D - m * m) // d
        a = (a0 + m) // d

        num_prev, num = num, num_prev + a * num
        den_prev, den = den, den_prev + a * den

    return num, den

# ---------------------------------------------------------
# Enumerate all square-pivots
# ---------------------------------------------------------

ans_set = set()

for m in range(1, M + 1):

    D = sf[m] * sf[m + 1]

    x1, y1 = pell_fundamental(D)

    # Generate all Pell solutions
    x, y = x1, y1

    # t = sf[m]*sf[m+1]*sq[m]*sq[m+1]*y
    mult = sf[m] * sf[m + 1] * sq[m] * sq[m + 1]

    while True:

        t = mult * y

        kk = m * (x + 1) + t

        if kk > 2 * N:
            break

        nn = (m + 1) * (x - 1) + t

        # parity checks + n >= k
        if kk % 2 == 0 and nn % 2 == 0 and nn >= kk:
            ans_set.add(kk // 2)

        # Next Pell solution
        x, y = (
            x1 * x + D * y1 * y,
            x1 * y + y1 * x
        )

answer = sum(ans_set)

print(answer)

4. Code walkthrough

Step 1 — Bound the search

M = (-1 + isqrt(1 + 2 * N)) // 2

From

$$2m(m+1)\le N,$$

this gives the maximum possible $m$.

Step 2 — Decompose integers

For every integer $n$,

$$n = \text{squarefree}(n)\times (\text{square part})^2.$$

We store:

sf[n] = squarefree part
sq[n] = square-root of square part

These are needed to form the Pell discriminant.

Step 3 — Solve Pell equations

For each $m$,

$$D = \operatorname{sf}(m)\operatorname{sf}(m+1).$$

We solve

$$x^2 - D y^2 = 1.$$

The function pell_fundamental computes the fundamental solution using continued fractions.

Step 4 — Generate all solutions

All Pell solutions come from powers of the fundamental unit:

$$x+y\sqrt D=(x_1+y_1\sqrt D)^n.$$

The recurrence

x, y = (
    x1 * x + D * y1 * y,
    x1 * y + y1 * x
)

generates them successively.

Step 5 — Recover $k$

Using

$$k=\frac{m(x+1)+t}{2},$$

with

$$t=\operatorname{sf}(m)\operatorname{sf}(m+1) \operatorname{sq}(m)\operatorname{sq}(m+1),y.$$

We verify parity and $n\ge k$, then store distinct $k$.

5. Verification on the examples

The program generates:

$4$
$21$
$24$
$110$

matching the problem statement.

The total number of distinct pivots found below $10^{10}$ is $72067$, and their sum is:

$$238890850232021.$$

This magnitude is reasonable given the density of Pell-generated solutions up to $10^{10}$.

Answer: 238890850232021

Project Euler Problem 261