Problem 211: Divisor Square Sum

Project Euler Problem 211 Solution

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

Problem Summary For this problem we set $L=64{,}000{,}000$ and study the divisor-square-sum function $$\sigma_2(n)=\sum_{d\mid n} d^2.$$ The goal is to compute $$\sum_{\substack{1 \le n \lt L \\ \exists r\in \mathbb{Z}_{\ge 0}:\ \sigma_2(n)=r^2}} n.$$ The condition is arithmetical rather than combinatorial: we are not counting divisors, but adding their squares and then testing whether the resulting integer is itself a square. Since $\sigma_2(1)=1$, the sum already includes $n=1$; the real work is to find all larger $n$ below the limit that survive the same test. Mathematical Approach The common structure behind all three implementations is the prime-factor description of $\sigma_2$. Once that formula is available, the problem becomes: generate or evaluate integers $n\lt L$, compute $\sigma_2(n)$ exactly, and keep only those for which the value is a square. Multiplicativity of $\sigma_2$ If $a$ and $b$ are coprime, then every divisor of $ab$ can be written uniquely as $d_1d_2$ with $d_1\mid a$ and $d_2\mid b$. Therefore $$\sigma_2(ab)=\sum_{d\mid ab} d^2=\sum_{d_1\mid a}\sum_{d_2\mid b}(d_1d_2)^2=\left(\sum_{d_1\mid a} d_1^2\right)\left(\sum_{d_2\mid b} d_2^2\right)=\sigma_2(a)\sigma_2(b).$$ So $\sigma_2$ is multiplicative....

Detailed mathematical approach

Problem Summary

For a positive integer $n$, let $\sigma_2(n)$ be the sum of the squares of its positive divisors; the statement gives the example $\sigma_2(10)=1+4+25+100=130$. We must find the sum of all $n$ with $0\lt n\lt 64\,000\,000$ for which $\sigma_2(n)$ is a perfect square, the square of an integer. ^[1] The number $n=1$ counts, since $\sigma_2(1)=1=1^2$.

The whole solution in one sentence: we build every number below the limit exactly once from its prime factorisation, smallest prime first, so that $\sigma_2$ of each new number is $\sigma_2$ of the number it was built from times $\sigma_2$ of one prime power, and we test each value exactly for being a square.

Mathematical Approach

No divisors need to be listed: Section 2 splits $\sigma_2$ into a product over prime powers, and Section 6 obtains each new value from an earlier one with two multiplications.

1. Small cases

For $n=1,2,\dots,12$ the values of $\sigma_2(n)$ are

$$\begin{aligned}&1,\ 5,\ 10,\ 21,\ 26,\ 50,\\&50,\ 85,\ 91,\ 130,\ 122,\ 210.\end{aligned}$$

Apart from the first, none of them is a square. The next square appears at $n=42$, whose divisors $1,2,3,6,7,14,21,42$ give

$$\begin{aligned}\sigma_2(42)&=1+4+9+36+49\\&\quad+196+441+1764\\&=2500=50^2.\end{aligned}$$

Below $1000$ exactly five numbers qualify, $1$, $42$, $246$, $287$ and $728$, with the sum $1304$. Below $10^4$ there are $12$ of them and below $10^6$ there are $41$; they form the sequence A046655 of the OEIS. ^[2]

2. Multiplicativity

Suppose $a$ and $b$ share no prime factor. By the fundamental theorem of arithmetic, ^[3] every divisor $d$ of $ab$ then splits in exactly one way as $d=d_1d_2$ with $d_1\mid a$ and $d_2\mid b$: the prime powers of $d$ whose primes divide $a$ form $d_1$, the others form $d_2$. Conversely every such product divides $ab$. Hence

$$\begin{aligned}\sigma_2(ab)&=\sum_{d_1\mid a}\ \sum_{d_2\mid b}d_1^2\,d_2^2\\&=\sigma_2(a)\,\sigma_2(b).\end{aligned}$$

Functions with this property are called multiplicative. ^[4] For example $\sigma_2(10)=\sigma_2(2)\,\sigma_2(5)$${}=5\cdot26=130$, the value of the statement. The condition matters: $\sigma_2(4)=1+4+16=21$, not $\sigma_2(2)^2=25$.

3. Prime powers

The divisors of a prime power $p^e$ are $1,p,p^2,\dots,p^e$, so $\sigma_2(p^e)$ is a geometric series: ^[4]

$$\begin{aligned}G(p,e)&=1+p^2+p^4+\dots+p^{2e}\\&=\frac{p^{2e+2}-1}{p^2-1}.\end{aligned}$$

Multiplying the sum by $p^2$ and adding $1$ gives the sum with one more term, which is the update rule of the C++ and Python programs:

$$\begin{aligned}G(p,0)&=1,\\G(p,e+1)&=p^2\,G(p,e)+1.\end{aligned}$$

For $p=2$ it produces $1,5,21,85,341$; the Java program instead keeps the newest term $p^{2j}$ and adds it to the sum. Together with Section 2, for $n=p_1^{e_1}p_2^{e_2}\cdots p_k^{e_k}$ with distinct primes $p_i$,

$$\begin{aligned}\sigma_2(n)&=G(p_1,e_1)\,G(p_2,e_2)\\&\quad\cdots G(p_k,e_k).\end{aligned}$$

For instance $728=2^3\cdot7\cdot13$ gives $\sigma_2(728)=85\cdot50\cdot170$${}=722500=850^2$.

4. When a product is a square

In this example no factor is a square, yet the product is. Every positive integer can be written in exactly one way as $x=q\,k^2$ with $q$ squarefree, that is, divisible by no square other than $1$; this $q$, the squarefree part of $x$, is the product of the primes with an odd exponent in $x$. ^[5] Exponents add when numbers are multiplied, so a product is a square exactly when the product of the squarefree parts of its factors is a square. For $\sigma_2(728)$ the factors $85=5\cdot17$, $50=2\cdot5^2$ and $170=2\cdot5\cdot17$ have the squarefree parts $85$, $2$ and $170$, and each of the primes $2$, $5$ and $17$ divides exactly two of them.

The same bookkeeping explains several small solutions. The values $\sigma_2(6)=\sigma_2(7)=50=2\cdot5^2$, $\sigma_2(41)=1682=2\cdot29^2$ and $\sigma_2(239)=57122=2\cdot169^2$ all have the squarefree part $2$, and no two of $6$, $7$, $41$, $239$ have a common prime factor, so for any two of them $\sigma_2(ab)=2s^2\cdot2t^2=(2st)^2$. This gives $42=6\cdot7$, with $\sigma_2(42)=50\cdot50$, and likewise $246$, $287$, $1434$, $1673$ and $9799$; the product of all four, $411558$, qualifies too.

So the search cannot keep a yes-or-no mark per factor: it must carry the exact value of $\sigma_2$, or at least its squarefree part, and test only the complete product.

5. How large the values become

The map $d\mapsto n/d$ runs through the divisors of $n$ in reverse order, so

$$\begin{aligned}\sigma_2(n)&=\sum_{d\mid n}\Bigl(\frac nd\Bigr)^2\\&=n^2\sum_{d\mid n}\frac1{d^2}.\end{aligned}$$

The last sum consists of some of the terms of $1+\frac14+\frac19+\dots$, whose total is $\pi^2/6\approx1.6449$ by Euler's solution of the Basel problem, ^[6] so $\sigma_2(n)\lt\frac{\pi^2}6\,n^2$. Below the limit $L=64\cdot10^6$ this gives $\sigma_2(n)\lt6.74\cdot10^{15}$, and the largest value that actually occurs is $\sigma_2(63977760)\approx6.56\cdot10^{15}$. Both are below $2^{53}\approx9.01\cdot10^{15}$, far below the bound $2^{63}\approx9.22\cdot10^{18}$ of signed 64-bit integers, so every value fits into a 64-bit integer and is even held exactly by a double-precision number, which represents every integer up to $2^{53}$ without rounding. ^[7] The final total is at most $1+2+\dots+(L-1)$${}=L(L-1)/2\lt2.05\cdot10^{15}$, so it fits as well.

6. Building every number from its factorisation

Arrange the numbers below the limit $L$ in a tree with the root $1$. The children of a node $m=p_1^{e_1}\cdots p_k^{e_k}$ with $p_1\lt\dots\lt p_k$ are the numbers $m\,q^f\lt L$ with a prime $q\gt p_k$ and an exponent $f\ge1$; the children of the root are the prime powers below $L$. The prime $q$ does not divide $m$, so Sections 2 and 3 give

$$\sigma_2(m\,q^f)=\sigma_2(m)\,G(q,f),$$

so, starting from $\sigma_2(1)=1$, each node costs at most two multiplications: one extends $G(q,f-1)$ to $G(q,f)$ by the rule of Section 3, and one multiplies by $\sigma_2(m)$.

Claim. Every $n$ with $2\le n\lt L$ is a node of the tree exactly once.

Proof. Write $n=p_1^{e_1}\cdots p_k^{e_k}$ with $p_1\lt\dots\lt p_k$. The numbers $p_1^{e_1}$, $p_1^{e_1}p_2^{e_2}$, …, $n$ divide $n$, so they are below $L$, and each of them is a child of the one before; hence $n$ is a node. Conversely, the path from the root to a node multiplies in prime powers with increasing primes, so it spells out a factorisation of the node with increasing primes. By the uniqueness of the factorisation, ^[3] only one path ends at $n$. $\square$

A depth-first search walks through the tree: at a node $m$ it tries the primes $q\gt p_k$ in increasing order and, for each, the exponents $f=1,2,\dots$. Once $m\,q^f\ge L$, larger exponents give larger products, so the exponent loop can stop; once even $m\,q\ge L$, every later prime $q'$ gives $m\,q'^f\gt m\,q\ge L$, so the prime loop can stop. Neither stop loses a node. The tree is shallow: the product of the first nine primes, $2\cdot3\cdots23=223092870$, exceeds the limit, so no number below it has more than $8$ distinct prime factors, and no path has more than $8$ steps.

7. Worked example: the limit 20

For $L=20$ the search visits the numbers

$$\begin{aligned}&2,\ 6,\ 18,\ 10,\ 14,\ 4,\\&12,\ 8,\ 16,\ 3,\ 15,\ 9,\\&5,\ 7,\ 11,\ 13,\ 17,\ 19\end{aligned}$$

in this order. Below the node $2$, with $\sigma_2(2)=5$, the prime $3$ gives $6=2\cdot3$ with $\sigma_2(6)=5\cdot10=50$ and $18=2\cdot3^2$ with $5\cdot91=455$, and its exponent loop stops at $2\cdot3^3=54$. Neither $6$ nor $18$ has a child, as $6\cdot5=30$ and $18\cdot5=90$ exceed $19$. The primes $5$ and $7$ give $10$ with $5\cdot26=130$ and $14$ with $5\cdot50=250$, and the prime $11$ ends the loop, since $2\cdot11=22\gt19$. The node $4$ has the single child $12$, with $21\cdot10=210$; the nodes $8$ and $16$ have none, since $8\cdot3=24\gt19$; and below $3$ only $15$ appears, with $10\cdot26=260$. All $18$ numbers from $2$ to $19$ occur once, and none of their values

$$\begin{aligned}&5,\ 50,\ 455,\ 130,\ 250,\ 21,\\&210,\ 85,\ 341,\ 10,\ 260,\ 91,\\&26,\ 50,\ 122,\ 170,\ 290,\ 362\end{aligned}$$

is a square, so the sum below $20$ is just $1$. Below $1000$ the same search finds $42$, $246$, $287$ and $728$ and returns $1+42+246+287+728=1304$.

8. An exact square test

A value $x$ is a square exactly when $r^2=x$ for $r=\lfloor\sqrt x\rfloor$, the largest integer with $r^2\le x$. The C++ and Java programs truncate a floating-point square root to obtain a first guess for $r$: Java applies Math.sqrt to a double, which holds $x$ exactly (Section 5), and gets the correctly rounded root, ^[8] while C++ applies std::sqrt to a long double. Even a correctly rounded root can make the guess one too large just below a square. For $k=8\cdot10^7$ the root $\sqrt{k^2-1}$ is about $k-6.25\cdot10^{-9}$, while neighbouring doubles near $k$ are $2^{-26}\approx1.49\cdot10^{-8}$ apart, ^[7] so the nearest double is $k$ itself, although $k^2\gt k^2-1$. The programs therefore correct the guess: they increase $r$ while $(r+1)^2\le x$ and then decrease it while $r^2\gt x$. Each decrease keeps $(r+1)^2\gt x$, so in the end $r^2\le x\lt(r+1)^2$, whatever the rounding was. Python's math.isqrt computes $\lfloor\sqrt x\rfloor$ with integer arithmetic alone. ^[8]

9. Two other ways to compute all values

The Java program computes $\sigma_2$ for all numbers in increasing order. A sieve first stores the smallest prime factor $p$ of every $n$. ^[9] If $p=n$, the number is prime and $\sigma_2(n)=1+n^2$. Otherwise dividing out the full power $p^e$ leaves $m=n/p^e$, which is smaller than $n$ and not divisible by $p$, so

$$\sigma_2(n)=\sigma_2(m)\,G(p,e)$$

with $\sigma_2(m)$ already known; for example $12=2^2\cdot3$ gives $m=3$ and $\sigma_2(12)=10\cdot21=210$. A third method needs no factorisation at all: the divisor sieve adds $d^2$ to the entries of all multiples of $d$, for every $d\lt L$. It makes $\sum_{d\lt L}\lfloor(L-1)/d\rfloor$${}\approx L\ln L$ additions, $1\,160\,244\,999$ for the problem, and serves below as an independent check.

How the Code Works

The C++ program ^[10] reads the options --limit=, by default $64\,000\,000$, and --skip-checkpoints. Unless the checkpoints are skipped, run_checkpoints requires the sums below $10$ and $1000$ to be $1$ and $1304$ (Sections 1 and 7) before main prints the result of solve for the requested limit.

solve builds the prime list with prime_list_upto(limit - 1), the sieve of Eratosthenes on a byte table: ^[9] for every $p\le\sqrt{L-1}$ still marked as prime it crosses off $p^2,p^2+p,\dots$, and the $3\,785\,086$ survivors for the default limit go into a vector of 32-bit integers. contribution_root(primes, idx, limit) handles the children of the root that are powers of $p$ = primes[idx]: for pwr $=p,p^2,\dots$ up to $L-1$ it applies s = s * p2 + 1, the rule of Section 3, adds pwr if s is a square, and passes the node on to contribution.

contribution(primes, from, s0, n0, limit) searches below the node n0, with s0 $=\sigma_2(n_0)$, as in Section 6. Its prime loop starts at index from and stops when n0 > max_n / p, the test $n_0\,p\gt L-1$ written without the product. Its exponent loop runs while n0 <= max_n / pwr: it forms n = n0 * pwr, updates s and forms s2 = s0 * s $=\sigma_2(n)$, adds n if s2 is a square and recurses with the next prime index, until pwr > max_n / p. By Section 5 all these numbers fit into the unsigned 64-bit type u64. is_perfect_square truncates a long double square root and runs the correction loops of Section 8 with products in unsigned __int128. Finally, solve adds the contributions of all possible smallest primes to the initial total $1$, which stands for $n=1$; by Section 6 their subtrees are disjoint and together contain every integer from $2$ to $L-1$.

The Python program defines is_perfect_square, prime_list_upto and contribution inside solve, with the limit fixed and without options or checkpoints. Its sieve crosses off the multiples of each prime with one slice assignment to a bytearray, the root loop is written out in solve, which returns the total as a string, and is_perfect_square uses math.isqrt. Python integers have no size limit, and the recursion is at most eight calls deep (Section 6).

The Java program follows Section 9 in main. Starting from spf[i] = i, for every i with $i^2\lt L$ and spf[i] == i it sets spf[j] = i on the multiples $j=i^2,i^2+i,\dots$ not yet marked, so that the int array spf holds the smallest prime factors. The long array sigma2 then gets $1+n^2$ for a prime $n$ and sigma2[m] * s otherwise; the inner loop divides m by $p$ while it can and adds pk $=p^2,p^4,\dots$ to s, which ends as $G(p,e)$. A last loop tests every value with Math.sqrt and the two correction loops and adds the $n$ that pass, $n=1$ included through sigma2[1] = 1. The casts in (long) n * n and (long) p * p make Java multiply in long, which by Section 5 suffices for every value and for the total.

Complexity and Verification

The sieve uses the $1007$ primes up to $7999$ and makes $153\,592\,611$ crossings, in line with its $O(L\log\log L)$ cost. ^[9] The search visits each of the $L-2$ numbers from $2$ to $L-1$ once and does constant work per visit: the exponent step that creates the node, one square test, and one call whose prime loop creates a child in every step but the last. An instrumented copy counted $63\,999\,998$ visits and $127\,833\,814$ prime-loop steps, root included, fewer than two per visit. So the search costs $O(L)$ and the whole method $O(L\log\log L)$ time. The Java program runs the same sieve and $O(L)$ further steps: its exponent loop runs once for each composite $n$ and once more for every further factor $p$ of $n$, and at most $L/p^e$ numbers below $L$ are divisible by $p^e$, so the extra passes number at most

$$L\sum_p\sum_{e\ge2}\frac1{p^e}=L\sum_p\frac1{p(p-1)}\lt L,$$

since the differences $\frac1{p(p-1)}=\frac1{p-1}-\frac1p$ add up to less than $1$. At the default limit the loop makes $99\,415\,816$ passes. All three programs use $O(L)$ memory: the C++ program peaked at about $95$ MB, mostly the $64$ MB byte table and $15$ MB of primes, the Python program at about $330$ MB, and Java's two arrays take $12$ bytes per number, $768$ MB, which the Java heap must provide. On an Apple M4 Max laptop a run takes about $0.35$ s in C++ (clang++ -O2), $0.65$ s in Java and $21$ s in CPython.

The C++ program checks itself with the two checkpoints, and all three programs print the same total. For this article the C++ program agreed with a direct summation over the divisors for every limit from $1$ to $3001$, set with --limit=, and copies of the Python and Java programs with the limit replaced by $10^3,10^4,\dots,10^7$ agreed with it, for example on $211234335$ below $10^7$. The instrumented copy confirmed at full size that every number from $2$ to $L-1$ is visited exactly once, at depth at most $8$. An independent C program computed all values with the divisor sieve of Section 9 and exact integer square roots. It found $5$, $12$, $22$, $41$, $82$ and $138$ qualifying numbers below $10^3$, $10^4$, $10^5$, $10^6$, $10^7$ and $64\cdot10^6$, the first $37$ equal to the terms listed in A046655, ^[2] confirmed the maximum of Section 5, never needed to correct the floating-point guess of Section 8, and gave the same final sum as the three published programs. These checks are not part of the published solutions.

Footnotes and References

The references below supply the background facts used above; the factorisation tree, its correctness and cost, the bound of Section 5 and the worked examples are derived in this article.

Project Euler 211 — Divisor Square Sum. The official statement defines σ₂(n) as the sum of the squares of the divisors of n, gives σ₂(10) = 130 and asks for the sum of all n with 0 < n < 64 000 000 for which σ₂(n) is a perfect square.
OEIS A046655 — Numbers whose sum of the squares of divisors is also a square number. The data section lists the terms 1, 42, 246, 287, 728, … up to 741096; the example is 246 = 2·3·41 with the sum 84100 = 290².
Wikipedia — Fundamental theorem of arithmetic. Every integer greater than 1 is a product of primes in exactly one way up to the order of the factors, written as p₁^n₁ ⋯ p_k^n_k with primes p₁ < ⋯ < p_k.
Wikipedia — Divisor function. σ_x(n) is the sum of the x-th powers of the positive divisors of n; it is multiplicative, and for n = ∏ p_i^a_i it is the product of the geometric sums (p_i^((a_i+1)x) − 1)/(p_i^x − 1).
Wikipedia — Square-free integer. Every positive integer is in a unique way the product of a largest possible square and a square-free integer, its square-free part: the product of the primes whose exponent is odd.
Wikipedia — Basel problem. The sum of the reciprocals of the squares of the positive integers is π²/6 ≈ 1.644934, as Leonhard Euler showed in 1734.
Wikipedia — Double-precision floating-point format. The significand has 53 bits, integers from −2⁵³ to 2⁵³ are exact, and between 2ⁿ and 2ⁿ⁺¹ the representable numbers are the multiples of 2ⁿ⁻⁵².
Java — Math.sqrt returns the double closest to the true square root, the squareRoot operation of IEEE 754. Python — math.isqrt returns the greatest integer a with a² ≤ n.
Wikipedia — Sieve of Eratosthenes. The multiples of each prime p are crossed off from p² on, because the smaller ones are already marked; only primes up to √n are needed, and the sieve takes O(n log log n) operations.
C++ — Euler211.cpp. The linked immutable revision contains prime_list_upto, contribution_root, contribution, is_perfect_square and run_checkpoints; it is the implementation record for this article.

Problem 211 source code

C++

#include <algorithm>
#include <atomic>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <limits>
#include <pthread.h>
#include <string>
#include <unistd.h>
#include <vector>

namespace {

using u32 = std::uint32_t;
using u64 = std::uint64_t;

struct Options {
    int limit = 64000000;
    bool run_checkpoints = true;
    int threads = 0;
};

bool parse_int_after_prefix(const std::string& arg, const std::string& prefix, int& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }

    int parsed = 0;
    for (char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        parsed = parsed * 10 + static_cast<int>(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_int_after_prefix(arg, "--limit=", options.limit)) {
            continue;
        }
        if (parse_int_after_prefix(arg, "--threads=", options.threads)) {
            continue;
        }

        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }

    return options.limit >= 1 && options.threads >= 0;
}

inline bool is_perfect_square(const u64 x) {
    u64 r = static_cast<u64>(std::sqrt(static_cast<long double>(x)));
    while (static_cast<unsigned __int128>(r + 1ULL) * static_cast<unsigned __int128>(r + 1ULL) <=
           x) {
        ++r;
    }
    while (static_cast<unsigned __int128>(r) * static_cast<unsigned __int128>(r) > x) {
        --r;
    }
    return static_cast<unsigned __int128>(r) * static_cast<unsigned __int128>(r) == x;
}

std::vector<u32> prime_list_upto(const int limit) {
    if (limit < 2) {
        return {};
    }
    std::vector<std::uint8_t> is_prime(static_cast<std::size_t>(limit + 1), 1U);
    is_prime[0] = 0U;
    is_prime[1] = 0U;
    const int r = static_cast<int>(std::sqrt(static_cast<long double>(limit)));
    for (int p = 2; p <= r; ++p) {
        if (is_prime[static_cast<std::size_t>(p)] == 0U) {
            continue;
        }
        const int step = p;
        int m = p * p;
        while (m <= limit) {
            is_prime[static_cast<std::size_t>(m)] = 0U;
            m += step;
        }
    }
    std::vector<u32> primes;
    primes.reserve(static_cast<std::size_t>(limit / std::max(1.0L, std::log((long double)limit))));
    for (int p = 2; p <= limit; ++p) {
        if (is_prime[static_cast<std::size_t>(p)] != 0U) {
            primes.push_back(static_cast<u32>(p));
        }
    }
    return primes;
}

u64 contribution(const std::vector<u32>& primes, const std::size_t from, const u64 s0,
                 const u64 n0, const u64 limit) {
    u64 result = 0ULL;
    const u64 max_n = limit - 1ULL;
    for (std::size_t i = from; i < primes.size(); ++i) {
        const u64 p = primes[i];
        if (n0 > max_n / p) {
            break;
        }

        const u64 p2 = p * p;
        u64 s = 1ULL;
        u64 pwr = p;

        while (n0 <= max_n / pwr) {
            const u64 n = n0 * pwr;
            s = s * p2 + 1ULL;
            const u64 s2 = s0 * s;

            if (is_perfect_square(s2)) {
                result += n;
            }
            result += contribution(primes, i + 1, s2, n, limit);

            if (pwr > max_n / p) {
                break;
            }
            pwr *= p;
        }
    }
    return result;
}

u64 contribution_root(const std::vector<u32>& primes, const std::size_t idx, const u64 limit) {
    const u64 max_n = limit - 1ULL;
    const u64 p = primes[idx];
    if (p > max_n) {
        return 0ULL;
    }

    const u64 p2 = p * p;
    u64 s = 1ULL;
    u64 pwr = p;
    u64 result = 0ULL;

    while (pwr <= max_n) {
        s = s * p2 + 1ULL;
        if (is_perfect_square(s)) {
            result += pwr;
        }
        result += contribution(primes, idx + 1, s, pwr, limit);
        if (pwr > max_n / p) {
            break;
        }
        pwr *= p;
    }

    return result;
}

struct WorkerArgs {
    const std::vector<u32>* primes = nullptr;
    u64 limit = 0ULL;
    std::atomic<std::size_t>* next_idx = nullptr;
    std::size_t max_idx = 0;
    std::size_t chunk = 1;
    u64 result = 0ULL;
};

void* worker_entry(void* raw) {
    auto* args = static_cast<WorkerArgs*>(raw);
    const auto& primes = *args->primes;
    u64 local = 0ULL;
    while (true) {
        const std::size_t start = args->next_idx->fetch_add(args->chunk, std::memory_order_relaxed);
        if (start >= args->max_idx) {
            break;
        }
        const std::size_t end = std::min(args->max_idx, start + args->chunk);
        for (std::size_t i = start; i < end; ++i) {
            local += contribution_root(primes, i, args->limit);
        }
    }
    args->result = local;
    return nullptr;
}

u64 solve(const int limit, int requested_threads) {
    if (limit <= 1) {
        return 0ULL;
    }

    const std::vector<u32> primes = prime_list_upto(limit - 1);
    u64 total = 1ULL;

    if (requested_threads <= 0) {
        long hw = sysconf(_SC_NPROCESSORS_ONLN);
        if (hw < 1) {
            hw = 1;
        }
        requested_threads = static_cast<int>(hw);
    }
    if (requested_threads < 1) {
        requested_threads = 1;
    }
    const std::size_t max_idx = primes.size();
    if (max_idx == 0) {
        return total;
    }

    std::size_t thread_count = static_cast<std::size_t>(requested_threads);
    if (thread_count > max_idx) {
        thread_count = max_idx;
    }
    if (thread_count <= 1) {
        for (std::size_t i = 0; i < max_idx; ++i) {
            total += contribution_root(primes, i, static_cast<u64>(limit));
        }
        return total;
    }

    std::atomic<std::size_t> next_idx(0);
    constexpr std::size_t chunk = 8;
    std::vector<pthread_t> tids(thread_count);
    std::vector<WorkerArgs> args(thread_count);

    bool create_failed = false;
    std::size_t started = 0;
    for (std::size_t t = 0; t < thread_count; ++t) {
        args[t].primes = &primes;
        args[t].limit = static_cast<u64>(limit);
        args[t].next_idx = &next_idx;
        args[t].max_idx = max_idx;
        args[t].chunk = chunk;
        args[t].result = 0ULL;

        if (pthread_create(&tids[t], nullptr, worker_entry, &args[t]) != 0) {
            create_failed = true;
            break;
        }
        ++started;
    }

    for (std::size_t t = 0; t < started; ++t) {
        pthread_join(tids[t], nullptr);
        total += args[t].result;
    }

    if (create_failed) {
        total = 1ULL;
        for (std::size_t i = 0; i < max_idx; ++i) {
            total += contribution_root(primes, i, static_cast<u64>(limit));
        }
    }

    return total;
}

bool run_checkpoints(const int threads) {
    if (solve(10, threads) != 1ULL) {
        std::cerr << "Checkpoint failed for limit=10" << '\n';
        return false;
    }
    if (solve(1000, threads) != 1304ULL) {
        std::cerr << "Checkpoint failed for limit=1000" << '\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(options.threads)) {
        return 2;
    }

    std::cout << solve(options.limit, options.threads) << '\n';
    return 0;
}

Python

import math

def solve():
    limit = 64000000

    def is_perfect_square(x):
        r = math.isqrt(x)
        return r * r == x

    def prime_list_upto(lim):
        if lim < 2:
            return []
        is_p = bytearray(b'\x01' * (lim + 1))
        is_p[0] = 0
        is_p[1] = 0
        r = math.isqrt(lim)
        for p in range(2, r + 1):
            if is_p[p]:
                is_p[p*p::p] = bytearray(len(is_p[p*p::p]))
        return [p for p in range(2, lim + 1) if is_p[p]]

    primes = prime_list_upto(limit - 1)

    def contribution(from_idx, s0, n0, lim):
        result = 0
        max_n = lim - 1
        for i in range(from_idx, len(primes)):
            p = primes[i]
            if n0 > max_n // p:
                break
            p2 = p * p
            s = 1
            pwr = p
            while n0 <= max_n // pwr:
                n = n0 * pwr
                s = s * p2 + 1
                s2 = s0 * s
                if is_perfect_square(s2):
                    result += n
                result += contribution(i + 1, s2, n, lim)
                if pwr > max_n // p:
                    break
                pwr *= p
        return result

    total = 1
    for i in range(len(primes)):
        p = primes[i]
        max_n = limit - 1
        if p > max_n:
            break
        p2 = p * p
        s = 1
        pwr = p
        while pwr <= max_n:
            s = s * p2 + 1
            if is_perfect_square(s):
                total += pwr
            total += contribution(i + 1, s, pwr, limit)
            if pwr > max_n // p:
                break
            pwr *= p
    return str(total)

if __name__ == '__main__':
    print(solve())

Java

import java.util.*;

public class Euler211 {
    public static void main(String[] args) {
        int limit = 64000000;
        long[] sigma2 = new long[limit];
        sigma2[1] = 1;
        int[] spf = new int[limit];
        for (int i = 0; i < limit; i++)
            spf[i] = i;
        for (int i = 2; (long) i * i < limit; i++)
            if (spf[i] == i)
                for (int j = i * i; j < limit; j += i)
                    if (spf[j] == j)
                        spf[j] = i;
        for (int n = 2; n < limit; n++) {
            int p = spf[n];
            if (p == n) {
                sigma2[n] = 1 + (long) n * n;
            } else {
                int m = n;
                long pk = 1, s = 1;
                while (m % p == 0) {
                    m /= p;
                    pk *= (long) p * p;
                    s += pk;
                }
                sigma2[n] = sigma2[m] * s;
            }
        }
        long total = 0;
        for (int n = 1; n < limit; n++) {
            long s2 = sigma2[n];
            long r = (long) Math.sqrt(s2);
            while ((r + 1) * (r + 1) <= s2)
                r++;
            while (r * r > s2)
                r--;
            if (r * r == s2)
                total += n;
        }
        System.out.println(total);
    }
}