Project Euler

Pivotal Square Sums

A positive integer k is called a square-pivot if there exist integers m > 0 and n >= k such that the sum of (m+1) consecutive squares ending at k equals the sum of m consecutive squares starting at...

Source sync Apr 19, 2026

Problem #0261

Level Level 17

Solved By 830

Languages C++, Python

Answer 238890850232021

Length 405 words

modular_arithmeticsequencenumber_theory

Let us call a positive integer $k$ a square-pivot, if there is a pair of integers $m > 0$ and $n \ge k$, such that the sum of the $(m+1)$ consecutive squares up to $k$ equals the sum of the $m$ consecutive squares from $(n+1)$ on: $$(k - m)^2 + \cdots + k^2 = (n + 1)^2 + \cdots + (n + m)^2.$$ Some small square-pivots are

$\textcolor{red}{4}$: $3^2 + \textcolor{red}{4}^2 = 5^2$
$\textcolor{red}{21}$: $20^2 + \textcolor{red}{21}^2 = 29^2$
$\textcolor{red}{24}$: $21^2 + 22^2 + 23^2 + \textcolor{red}{24}^2 = 25^2 + 26^2 + 27^2$
$\textcolor{red}{110}$: $108^2 + 109^2 + \textcolor{red}{110}^2 = 133^2 + 134^2$

Find the sum of all distinct square-pivots $\le 10^{10}$.

Problem 261: Pivotal Square Sums

Mathematical Development

Notation

Denote the sum-of-squares function $S(t) = \sum_{i=1}^{t} i^2 = \tfrac{t(t+1)(2t+1)}{6}$ .

Theorem 1 (Reduction to generalized Pell equation). If $k$ is a square-pivot with parameters $(m, n)$ , then writing $m = s p^2$ with $s$ squarefree and $m + 1 = g r^2$ with $g$ squarefree, there exist non-negative integers $q, u$ satisfying

g u^2 - s q^2 = 1.

Equivalently, setting $D = sg$ , $X = gu$ , $Y = q$ , one obtains the generalized Pell equation $X^2 - D Y^2 = g$ .

Proof. The pivot condition is $S(k) - S(k - m - 1) = S(n + m) - S(n)$ . Expanding both sides using the closed form of $S$ and introducing the substitutions $x = 2k - m$ and $t = 2n + m + 1$ , algebraic simplification yields:

x^2(m+1) - t^2 m = m(m+1).

Dividing by $m(m+1)$ , we obtain $\frac{x^2}{m} - \frac{t^2}{m+1} = 1$ , which factors as

x^2 = m \cdot a, \qquad t^2 = (m+1)(a+1)

for some non-negative integer $a$ (where $a = \frac{x^2}{m}$ ). With the squarefree decompositions $m = s p^2$ and $m + 1 = g r^2$ , we deduce $x = s p q$ (so that $a = s q^2$ ) and $t = g r u$ (so that $a + 1 = g u^2$ ). Subtracting:

g u^2 - s q^2 = (a + 1) - a = 1.

Setting $X = gu$ and $Y = q$ gives $X^2 - D Y^2 = g^2 u^2 - sg q^2 = g(gu^2 - sq^2) = g$ . $\square$

Lemma 1 (Base solution). The pair $(X_0, Y_0) = (gr, p)$ satisfies $X_0^2 - D Y_0^2 = g$ .

Proof. Direct computation: $X_0^2 - D Y_0^2 = g^2 r^2 - sg \cdot p^2 = g(gr^2 - sp^2) = g\bigl((m+1) - m\bigr) = g$ . $\square$

Lemma 2 (Pell recurrence). Let $(x_1, y_1)$ be the fundamental solution to $X^2 - DY^2 = 1$ (obtained via the continued-fraction expansion of $\sqrt{D}$ ). Then all solutions to $X^2 - DY^2 = g$ in the orbit of $(X_0, Y_0)$ are generated by the recurrence

(X_{j+1}, Y_{j+1}) = (X_j \, x_1 + D \, Y_j \, y_1, \; X_j \, y_1 + Y_j \, x_1).

Proof. By the Brahmagupta—Fibonacci identity, if $(X_j, Y_j)$ satisfies $X_j^2 - D Y_j^2 = g$ and $(x_1, y_1)$ satisfies $x_1^2 - D y_1^2 = 1$ , then

(X_j x_1 + D Y_j y_1)^2 - D(X_j y_1 + Y_j x_1)^2 = (X_j^2 - D Y_j^2)(x_1^2 - D y_1^2) = g \cdot 1 = g.

Moreover, the new solution is strictly larger than $(X_j, Y_j)$ since $x_1 \ge 1$ and $y_1 \ge 1$ . The orbit exhausts all positive solutions in the given class (a standard result in the theory of Pell equations; see, e.g., Jacobson and Williams, Solving the Pell Equation, CMS, 2009). $\square$

Lemma 3 (Recovering $k$ and $n$ ). Given a Pell solution $(X, Y)$ with $Y = q$ , the pivot is

k = \frac{s p (p + q)}{2}, \qquad n = \frac{g r u - m - 1}{2}

where $u = X / g$ . The pivot is valid if and only if $k$ is a positive integer and $n \ge k$ .

Proof. From $x = spq$ and $x = 2k - m$ , we get $k = \frac{spq + m}{2} = \frac{sp(q + p)}{2}$ (using $m = sp^2$ ). Similarly, $t = gru$ and $t = 2n + m + 1$ gives $n = \frac{gru - m - 1}{2}$ . $\square$

Lemma 4 (Bound on $m$ ). For a valid pivot $k \le 10^{10}$ , we have $m \le 70710$ .

Proof. The constraint $n \ge k$ combined with the expressions for $k$ and $n$ forces $k \ge 2m(m+1)$ . Indeed, the smallest valid $q$ satisfies $q \ge p + 1$ , which after substitution yields $k \ge sp^2(p + p + 1)/2 > m$ , and the refined inequality $k \ge 2m(m+1)$ follows from the requirement $n \ge k$ . Then $2m(m+1) \le 10^{10}$ implies $m^2 < 5 \times 10^9$ , so $m \le \lfloor\sqrt{5 \times 10^9}\rfloor = 70710$ . $\square$

Editorial

Reduction: the pivot equation reduces to a generalized Pell equation X^2 - DY^2 = g where D = sg, with m = sp^2 (s squarefree) and m+1 = gr^2 (g squarefree). Iterating Pell solutions recovers k values. We loop.

Pseudocode

pivots = empty set
For m from 1 to 70710:
    compute s, p such that m = s * p^2 (s squarefree)
    compute g, r such that m + 1 = g * r^2 (g squarefree)
    D = s * g
    if D is a perfect square: skip // Pell equation X^2 - D*Y^2 = 1 has no solution
    (x1, y1) = fundamental_pell_solution(D) // via continued fractions
    (X, Y) = (g * r, p) // base solution from Lemma 1
    loop:
        q = Y
        k = s * p * (p + q) / 2
        If k > 10^10 then stop this loop
        If k is a positive integer and k >= 2 * m * (m + 1) then
            u = X / g
            n_val = (g * r * u - m - 1) / 2
            If n_val >= k then
                pivots.add(k)
        (X, Y) = (X * x1 + D * Y * y1, X * y1 + Y * x1) // Lemma 2
Return sum(pivots)

Complexity Analysis

Time. The outer loop runs $M = 70710$ iterations. For each $m$ :

The squarefree decomposition uses a smallest-prime-factor sieve, computed once in $O(M \log \log M)$ .
Finding the fundamental Pell solution via continued fractions takes $O(\sqrt{D})$ arithmetic operations, where $D = sg \le m(m+1)$ . With caching, each distinct $D$ is solved only once.
The inner Pell iteration produces at most $O(\log(10^{10}))$ solutions per orbit, since $X_j$ grows exponentially.

Total: $O(M \log \log M + M \cdot (\sqrt{D_{\max}} + \log N))$ where $N = 10^{10}$ and $D_{\max} = O(M^2)$ .

Space. $O(M)$ for the sieve and Pell cache; $O(|\text{pivots}|)$ for the result set.

Answer

\boxed{238890850232021}

C++ project_euler/problem_261/solution.cpp

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef __int128 lll;

/*
 * Problem 261: Pivotal Square Sums
 *
 * Reduction to generalized Pell equation X^2 - D*Y^2 = g,
 * where m = s*p^2, m+1 = g*r^2, D = s*g.
 * Iterate Pell solutions to recover pivots k = s*p*(p+q)/2.
 */

const ll LIMIT = 10000000000LL;
int spf_arr[100002];

void build_spf(int n) {
    for (int i = 0; i <= n; i++) spf_arr[i] = i;
    for (int i = 2; (ll)i * i <= n; i++)
        if (spf_arr[i] == i)
            for (int j = i * i; j <= n; j += i)
                if (spf_arr[j] == j) spf_arr[j] = i;
}

pair<int, int> sqfree_sqrt(int n) {
    int sf = 1, sq = 1;
    while (n > 1) {
        int p = spf_arr[n], e = 0;
        while (n % p == 0) { n /= p; e++; }
        if (e & 1) sf *= p;
        for (int i = 0; i < e / 2; i++) sq *= p;
    }
    return {sf, sq};
}

unordered_map<ll, pair<ll, ll>> pcache;

pair<ll, ll> pell_fund(ll D) {
    auto it = pcache.find(D);
    if (it != pcache.end()) return it->second;
    ll a0 = (ll)sqrtl((long double)D);
    while (a0 * a0 > D) a0--;
    while ((a0 + 1) * (a0 + 1) <= D) a0++;
    ll m = 0, d = 1, a = a0;
    ll n1 = 1, n0 = a0, d1 = 0, d0 = 1;
    for (int i = 0; i < 100000000; i++) {
        m = d * a - m;
        d = (D - m * m) / d;
        if (d == 0) break;
        a = (a0 + m) / d;
        ll n2 = n1, d2 = d1;
        n1 = n0; d1 = d0;
        n0 = a * n1 + n2; d0 = a * d1 + d2;
        if ((lll)n0 * n0 - (lll)D * d0 * d0 == 1) {
            pcache[D] = {n0, d0};
            return {n0, d0};
        }
        if (n0 > 4000000000000000000LL) break;
    }
    pcache[D] = {0, 0};
    return {0, 0};
}

int main() {
    int mmax = 70710;
    while (2LL * (mmax + 1) * ((ll)(mmax + 1) + 1) <= LIMIT) mmax++;

    build_spf(mmax + 2);
    unordered_set<ll> pivots;
    pivots.reserve(1 << 17);

    for (int m = 1; m <= mmax; m++) {
        auto [s, p] = sqfree_sqrt(m);
        auto [g, r] = sqfree_sqrt(m + 1);
        ll D = (ll)s * g;
        auto [x1, y1] = pell_fund(D);
        if (x1 == 0) continue;

        ll x = (ll)g * r, y = p;
        for (int iter = 0; iter < 100; iter++) {
            ll q = y;
            lll num = (lll)s * p * ((ll)p + q);
            if (num / 2 > LIMIT && q > p) break;
            if (num % 2 == 0) {
                ll k = (ll)(num / 2);
                if (k <= LIMIT && k >= 2LL * m * (m + 1)) {
                    ll u = x / g;
                    lll tv = (lll)g * r * u;
                    lll rem = tv - m - 1;
                    if (rem >= 0 && rem % 2 == 0) {
                        ll n_val = (ll)(rem / 2);
                        if (n_val >= k)
                            pivots.insert(k);
                    }
                }
            }
            lll xn = (lll)x * x1 + (lll)y * y1 * D;
            lll yn = (lll)x * y1 + (lll)y * x1;
            if (xn > (lll)4e18 || yn > (lll)4e18) break;
            x = (ll)xn;
            y = (ll)yn;
        }
    }

    ll ans = 0;
    for (ll k : pivots) ans += k;
    cout << ans << endl;
    return 0;
}

Python project_euler/problem_261/solution.py

"""
Problem 261: Pivotal Square Sums

Find the sum of all distinct square-pivots k <= 10^10.

Reduction: the pivot equation reduces to a generalized Pell equation
X^2 - D*Y^2 = g where D = s*g, with m = s*p^2 (s squarefree) and
m+1 = g*r^2 (g squarefree). Iterating Pell solutions recovers k values.
"""

import math

LIMIT = 10**10


def build_spf(n):
    """Smallest prime factor sieve up to n."""
    spf = list(range(n + 1))
    for i in range(2, int(math.isqrt(n)) + 1):
        if spf[i] == i:
            for j in range(i * i, n + 1, i):
                if spf[j] == j:
                    spf[j] = i
    return spf


def squarefree_and_sqrt(n, spf):
    """Decompose n = sf * sq^2 with sf squarefree."""
    sf, sq = 1, 1
    while n > 1:
        p = spf[n]
        e = 0
        while n % p == 0:
            n //= p
            e += 1
        if e & 1:
            sf *= p
        sq *= p ** (e // 2)
    return sf, sq


def pell_fundamental(D, cache):
    """Fundamental solution (x, y) to x^2 - D*y^2 = 1 via continued fractions."""
    if D in cache:
        return cache[D]
    a0 = int(math.isqrt(D))
    m, d, a = 0, 1, a0
    num1, num = 1, a
    den1, den = 0, 1
    while num * num - D * den * den != 1:
        m = d * a - m
        d = (D - m * m) // d
        a = (a0 + m) // d
        num2, num1 = num1, num
        den2, den1 = den1, den
        num = a * num1 + num2
        den = a * den1 + den2
    cache[D] = (num, den)
    return num, den


def solve(limit=LIMIT):
    mmax = (math.isqrt(1 + 2 * limit) - 1) // 2
    spf = build_spf(mmax + 2)
    pell_cache = {}
    pivots = set()

    for m in range(1, mmax + 1):
        s, p = squarefree_and_sqrt(m, spf)
        g, r = squarefree_and_sqrt(m + 1, spf)
        D = s * g
        x1, y1 = pell_fundamental(D, pell_cache)
        X, Y = g * r, p

        while True:
            q = Y
            numerator = s * p * (p + q)
            if numerator // 2 > limit and q > p:
                break
            if numerator % 2 == 0:
                k = numerator // 2
                if k <= limit and k >= 2 * m * (m + 1):
                    u = X // g
                    t = g * r * u
                    if (t - m - 1) % 2 == 0:
                        n = (t - m - 1) // 2
                        if n >= k:
                            pivots.add(k)
            X, Y = X * x1 + Y * y1 * D, X * y1 + Y * x1

    return sum(pivots)


if __name__ == "__main__":
    answer = solve()
assert answer == 238890850232021
print(answer)