Project Euler

Sum of Squares

Consider equations of the form a^2 + b^2 = n with 0 <= a <= b and a, b, n integers. Let P = {p prime: p equiv 1 (mod 4), p < 150}. For each squarefree N that is a product of a non-empty subset of P...

Source sync Apr 19, 2026

Problem #0273

Level Level 12

Solved By 1,630

Languages C++, Python

Answer 2032447591196869022

Length 399 words

number_theorysearchmodular_arithmetic

Consider equations of the form: $a^2 + b^2 = N$, $0 \le a \le b$, $a$, $b$ and $N$ integer.

For $N=65$ there are two solutions: $a=1$, $b=8$ and $a=4$, $b=7$.

We call $S(N)$ the sum of the values of $a$ of all solutions of $a^2 + b^2 = N$, $0 \le a \le b$, $a$, $b$ and $N$ integer.

Thus $S(65) = 1 + 4 = 5$.

Find $\sum S(N)$, for all squarefree $N$ only divisible by primes of the form $4k+1$ with $4k+1 < 150$.

Problem 273: Sum of Squares

Mathematical Foundation

Theorem 1 (Fermat’s Two-Square Theorem). A prime $p$ is representable as a sum of two squares if and only if $p = 2$ or $p \equiv 1 \pmod{4}$ .

Proof. ( $\Longleftarrow$ ) If $p \equiv 1 \pmod{4}$ , then $-1$ is a quadratic residue mod $p$ , so $p \mid (m^2 + 1)$ for some $m$ . In the Gaussian integers $\mathbb{Z}[i]$ , $p \mid (m+i)(m-i)$ but $p \nmid (m \pm i)$ , so $p$ is not a Gaussian prime. Hence $p = \pi \bar{\pi}$ with $|\pi|^2 = a^2 + b^2 = p$ . ( $\Longrightarrow$ ) If $a^2 + b^2 = p$ with $p$ odd, then $a^2 \equiv -b^2 \pmod{p}$ , so $(ab^{-1})^2 \equiv -1 \pmod{p}$ , forcing $p \equiv 1 \pmod{4}$ . $\square$

Theorem 2 (Gaussian Integer Representation). Let $N = p_1 p_2 \cdots p_k$ be a squarefree product of primes $p_j \equiv 1 \pmod{4}$ , with Gaussian factorisation $p_j = \pi_j \bar{\pi}_j$ . The representations $a^2 + b^2 = N$ with $a, b \ge 0$ are in bijection with

\alpha = u \prod_{j=1}^{k} \pi_j^{\epsilon_j} \bar{\pi}_j^{1 - \epsilon_j}, \quad \epsilon_j \in \{0, 1\},\; u \in \{1, i, -1, -i\}

where $\alpha = a + bi$ and we identify $(a, b) = (|\operatorname{Re}(\alpha)|, |\operatorname{Im}(\alpha)|)$ .

Proof. $\mathbb{Z}[i]$ is a PID, so $N = \prod_j \pi_j \bar{\pi}_j$ is the unique factorisation (up to units). Any Gaussian integer $\alpha$ with $|\alpha|^2 = N$ must be a product of exactly one element from each pair $\{\pi_j, \bar{\pi}_j\}$ , times a unit. There are $2^k$ choices for the $\epsilon_j$ and 4 units, giving $4 \cdot 2^k$ Gaussian integers. Identifying $(a, b)$ with $(|a|, |b|)$ and $a \le b$ reduces the count by a factor of 8 (generically), yielding $2^{k-1}$ distinct representations when $k \ge 1$ (with possible boundary cases when $a = 0$ or $a = b$ ). $\square$

Lemma 1 (Brahmagupta—Fibonacci Identity). $(a^2 + b^2)(c^2 + d^2) = (ac - bd)^2 + (ad + bc)^2 = (ac + bd)^2 + (ad - bc)^2$ .

Theorem 3 (Incremental Construction). If $N_1$ has representations $\{(a_i, b_i)\}$ and $p = c^2 + d^2$ is a new prime, then $N_1 \cdot p$ has representations obtained by applying the Brahmagupta—Fibonacci identity to each $(a_i, b_i)$ with $(c, d)$ :

(a_i, b_i) \otimes (c, d) \;\longmapsto\; (|a_i c - b_i d|,\; a_i d + b_i c) \text{ and } (|a_i c + b_i d|,\; |a_i d - b_i c|)

Each existing representation produces exactly two new ones (possibly after swapping coordinates to ensure $a \le b$ ).

Proof. This follows from Lemma 1 and the Gaussian integer product: $(a_i + b_i \, i)(c + d\, i)$ and $(a_i + b_i \, i)(c - d\, i)$ yield the two new representations. $\square$

Editorial

For squarefree N that is a product of primes p = 1 mod 4 with p < 150, let S(N) = sum of all a values where a^2 + b^2 = N, 0 <= a <= b. Find sum of S(N) over all such N. Primes = 1 mod 4 less than 150: 5,13,17,29,37,41,53,61,73,89,97,101,109,113,137,149 (16 primes) Uses Gaussian integers and Brahmagupta-Fibonacci identity. We 16 primes. We then iterate over each prime p, find (c, d) with c^2 + d^2 = p, c ≤ d. Finally, we run a depth-first search over the subsets while maintaining the current list of $(a,b)$ representations.

Pseudocode

16 primes
For each prime p, find (c, d) with c^2 + d^2 = p, c ≤ d
DFS over subsets, maintaining list of (a, b) representations
reps = list of (a, b) with a ≤ b for current product
if reps is not empty
for i from index to 15
Two new representations
Bootstrap: each single prime p = c^2 + d^2 gives one representation

Complexity Analysis

Time: Each subset of size $k$ produces $2^{k-1}$ representations. The total work over all $2^{16} - 1$ non-empty subsets is $\sum_{k=1}^{16} \binom{16}{k} \cdot 2^{k-1} = \frac{1}{2}(3^{16} - 1) \approx 2.15 \times 10^7$ . Each representation requires $O(1)$ arithmetic. Total: $O(3^{16})$ .
Space: The DFS stack depth is at most 16. At each level, the representation list has at most $2^{15}$ entries. Space: $O(2^{16})$ .

Answer

\boxed{2032447591196869022}

C++ project_euler/problem_273/solution.cpp

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

// Problem 273: Sum of Squares
//
// For a squarefree N that is a product of primes p = 1 mod 4 with p < 150,
// let S(N) = sum of all values of a where a^2 + b^2 = N, 0 <= a <= b.
// Find the sum of S(N) over all such N.
//
// Primes p = 1 mod 4 less than 150: 5, 13, 17, 29, 37, 41, 53, 61, 73, 89, 97, 101, 109, 113, 137, 149
// There are 16 such primes, giving 2^16 - 1 = 65535 squarefree products.

typedef pair<ll, ll> gi;

gi multiply(gi a, gi b) {
    return {a.first * b.first - a.second * b.second,
            a.first * b.second + a.second * b.first};
}

gi conj(gi a) {
    return {a.first, -a.second};
}

void print128(lll v) {
    if (v == 0) { cout << "0" << endl; return; }
    bool neg = false;
    if (v < 0) { neg = true; v = -v; }
    string s;
    while (v > 0) { s += ('0' + (int)(v % 10)); v /= 10; }
    if (neg) s += '-';
    reverse(s.begin(), s.end());
    cout << s << endl;
}

int main() {
    // Find primes p = 1 mod 4 less than 150
    vector<int> primes;
    for (int p = 2; p < 150; p++) {
        bool is_p = true;
        for (int d = 2; d * d <= p; d++)
            if (p % d == 0) { is_p = false; break; }
        if (is_p && p % 4 == 1) primes.push_back(p);
    }
    int np = primes.size(); // 16

    // For each prime, find Gaussian integer factor pi = a + bi with a^2 + b^2 = p, a <= b
    vector<gi> prime_gi(np);
    for (int i = 0; i < np; i++) {
        int p = primes[i];
        for (int a = 1; a * a < p; a++) {
            int b2 = p - a * a;
            int b = (int)round(sqrt(b2));
            if (b * b == b2 && a <= b) {
                prime_gi[i] = {a, b};
                break;
            }
        }
    }

    // For each non-empty subset of the 16 primes (2^16 - 1 subsets),
    // compute all representations a^2 + b^2 = N and sum the a values (with a <= b).
    //
    // For a product of k primes, there are 2^(k-1) essentially distinct representations
    // (the factor of 2 reduction is because conjugating ALL factors gives (a,-b) -> same {a,b}).
    //
    // We use the Gaussian integer product: for each subset, multiply the Gaussian factors
    // choosing conjugate or not for each. alpha = product -> a = Re(alpha), b = Im(alpha).
    // We want 0 <= a <= b.
    //
    // Optimization: enumerate subsets and for each, iterate over 2^(k-1) sign masks.
    // With k up to 16, 2^15 = 32768 per subset, and 65535 subsets, total ~ 2 * 10^9.
    // That's too slow for brute force.
    //
    // Better: build incrementally. Use DP where state is the current set of representations.
    // For each new prime, combine existing representations with the new prime's representation.
    //
    // For a single prime p = a^2 + b^2: representations are {(a,b)} (just one with a<=b).
    // When multiplying two representations:
    // (a^2+b^2)(c^2+d^2) = (ac-bd)^2 + (ad+bc)^2 = (ac+bd)^2 + (|ad-bc|)^2
    //
    // So from one existing representation (a,b) and one new prime representation (c,d),
    // we get two new representations: (|ac-bd|, ad+bc) and (|ad-bc|, ac+bd)
    // (taking absolute values and ensuring a <= b).
    //
    // We process primes one by one. We maintain for the "current product" a list of
    // all (a,b) representations. Adding a new prime doubles the number of representations.
    //
    // But we need to process ALL subsets. The key: we want sum of S(N) over all N.
    // S(N) = sum of a values over all representations (a,b) of N with a <= b.
    //
    // We can compute the total by DFS over subsets of primes:
    // At each step, we either include or exclude the next prime.
    // If included, combine each existing representation with the prime's representation.
    // Track the sum of a values across all representations.
    //
    // To avoid storing all representations (which grow exponentially), we process
    // the inclusion of each prime and aggregate the sum of a values.
    //
    // Actually, for 16 primes, max representations per subset = 2^15 = 32768.
    // Total representations across all subsets = sum over k=1..16 of C(16,k) * 2^(k-1)
    // = (1/2) * sum of C(16,k) * 2^k = (1/2) * (3^16 - 1) ~ 21.5 million.
    // That's manageable!

    lll total = 0;

    // DFS: idx = current prime index, reps = list of current representations
    function<void(int, vector<gi>&)> dfs = [&](int idx, vector<gi>& reps) {
        // Try including each remaining prime
        for (int i = idx; i < np; i++) {
            gi pg = prime_gi[i];
            vector<gi> new_reps;

            if (reps.empty()) {
                // First prime: just one representation
                new_reps.push_back(pg);
            } else {
                // Combine each existing rep with the new prime
                for (auto& r : reps) {
                    // Two new representations from Brahmagupta-Fibonacci
                    ll a1 = abs(r.first * pg.first - r.second * pg.second);
                    ll b1 = abs(r.first * pg.second + r.second * pg.first);
                    ll a2 = abs(r.first * pg.first + r.second * pg.second);
                    ll b2 = abs(r.first * pg.second - r.second * pg.first);
                    // Ensure a <= b
                    if (a1 > b1) swap(a1, b1);
                    if (a2 > b2) swap(a2, b2);
                    new_reps.push_back({a1, b1});
                    new_reps.push_back({a2, b2});
                }
            }

            // Compute S(N) for this subset: sum of a values
            for (auto& r : new_reps) {
                total += r.first;
            }

            // Continue with more primes
            dfs(i + 1, new_reps);
        }
    };

    vector<gi> empty;
    dfs(0, empty);

    print128(total);
    return 0;
}

Python project_euler/problem_273/solution.py

"""
Problem 273: Sum of Squares

For squarefree N that is a product of primes p = 1 mod 4 with p < 150,
let S(N) = sum of all a values where a^2 + b^2 = N, 0 <= a <= b.
Find sum of S(N) over all such N.

Primes = 1 mod 4 less than 150: 5,13,17,29,37,41,53,61,73,89,97,101,109,113,137,149
(16 primes)

Uses Gaussian integers and Brahmagupta-Fibonacci identity.
"""

def solve():
    # Find primes p = 1 mod 4, p < 150
    primes = []
    for p in range(2, 150):
        if all(p % d != 0 for d in range(2, int(p**0.5) + 1)) and p > 1 and p % 4 == 1:
            primes.append(p)

    # Gaussian integer representation for each prime
    prime_gi = []
    for p in primes:
        for a in range(1, p):
            b2 = p - a * a
            if b2 <= 0:
                break
            b = int(b2**0.5 + 0.5)
            if b * b == b2 and a <= b:
                prime_gi.append((a, b))
                break

    total = 0

    def dfs(idx, reps):
        nonlocal total
        for i in range(idx, len(primes)):
            pg = prime_gi[i]
            if not reps:
                new_reps = [pg]
            else:
                new_reps = []
                for (a, b) in reps:
                    c, d = pg
                    # Brahmagupta-Fibonacci: two new representations
                    a1 = abs(a * c - b * d)
                    b1 = abs(a * d + b * c)
                    a2 = abs(a * c + b * d)
                    b2 = abs(a * d - b * c)
                    if a1 > b1: a1, b1 = b1, a1
                    if a2 > b2: a2, b2 = b2, a2
                    new_reps.append((a1, b1))
                    new_reps.append((a2, b2))

            # Add S(N) for this subset
            for (a, b) in new_reps:
                total += a

            dfs(i + 1, new_reps)

    dfs(0, [])
    print(total)

if __name__ == "__main__":
    solve()