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 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. If

$$n=\prod_{i=1}^r p_i^{e_i},$$

then

$$\sigma_2(n)=\prod_{i=1}^r \sigma_2\!\left(p_i^{e_i}\right)=\prod_{i=1}^r \left(1+p_i^2+p_i^4+\cdots+p_i^{2e_i}\right).$$

This is the key formula used everywhere in the page: the square test for $\sigma_2(n)$ is controlled completely by the prime-power factors of $n$.

Prime-power contribution and the useful recurrence

For a fixed prime $p$, define

$$G_p(e)=1+p^2+p^4+\cdots+p^{2e}=\frac{p^{2(e+1)}-1}{p^2-1}.$$

Then the factorization above becomes

$$\sigma_2(n)=\prod_i G_{p_i}(e_i).$$

The closed form is mathematically convenient, but the implementations use an even simpler iterative identity:

$$G_p(0)=1,\qquad G_p(e+1)=p^2\,G_p(e)+1.$$

Indeed, multiplying $1+p^2+\cdots+p^{2e}$ by $p^2$ and then adding 1 produces $1+p^2+\cdots+p^{2e}+p^{2e+2}$. This recurrence is exactly what lets the recursive search extend an exponent from $e$ to $e+1$ without recomputing a geometric series from scratch.

Worked example: $n=42$

A concrete example shows why the multiplicative viewpoint is so effective. Since

$$42=2\cdot 3\cdot 7,$$

we get

$$\sigma_2(42)=\left(1+2^2\right)\left(1+3^2\right)\left(1+7^2\right)=5\cdot 10\cdot 50=2500=50^2.$$

So $42$ is one of the integers that must be included in the final sum. The important point is that no divisor list had to be generated explicitly: prime factorization plus multiplicativity already gives the answer.

Why the recursive factorization search is unique

The C++ and Python implementations maintain a state of the form

$$n_0=\prod_{i=1}^k p_i^{e_i},\qquad s_0=\sigma_2(n_0),$$

with strictly increasing primes $p_1<p_2<\cdots<p_k$. To extend the state, they choose a larger prime $q$ and an exponent $f\ge 1$, then form

$$n=n_0q^f,\qquad \sigma_2(n)=s_0\,G_q(f).$$

Because each new prime is introduced only after all smaller chosen primes are fixed, every integer below the limit is generated exactly once, namely in the order dictated by its unique prime factorization. The pruning condition $n<L$ is monotone, so as soon as a prime power pushes $n$ past the limit, all larger exponents for that prime can be discarded immediately.

The sieve recurrence used by the full-array approach

The Java implementation uses the same mathematics in a different way. Write

$$n=mp^e,\qquad p\nmid m,$$

where $p$ is the smallest prime dividing $n$. Then

$$\sigma_2(n)=\sigma_2(m)\,G_p(e).$$

This gives a recurrence from a smaller integer $m$ to the current integer $n$. Once the smallest prime factor of every number is known, the algorithm can strip off the full $p$-power, build the geometric factor $G_p(e)$, and recover $\sigma_2(n)$ from a previously computed value. In other words, one family of implementations explores valid factorizations directly, while the other fills the entire table of $\sigma_2(n)$ values from $1$ up to $L-1$.

How the Code Works

C++ and Python: enumerate factorizations

These implementations first build the prime list up to $L-1$. They then recurse over increasing primes, carrying two exact integers: the current candidate $n_0$ and the already-known value $\sigma_2(n_0)$. For each newly chosen prime $p$, they try exponents $1,2,3,\dots$ as long as the product stays below the limit, update the prime-power factor with the recurrence $G_p(e+1)=p^2G_p(e)+1$, and multiply it into the current $\sigma_2$ value.

Whenever the new value of $\sigma_2(n)$ is a perfect square, the corresponding $n$ is added to the running total. Since the recursion visits each admissible factorization once, no deduplication structure is needed. The C++ version also splits the outermost prime choices across several worker threads so that separate branches of the factorization tree can be processed independently.

Java: build $\sigma_2(n)$ for every $n$

The Java implementation takes an array-based route. It first computes the smallest prime factor of every integer below the limit. Then it processes $n=2,3,\dots,L-1$ in order. If $n$ is prime, the formula is immediate:

$$\sigma_2(n)=1+n^2.$$

If $n$ is composite, the implementation extracts the full power $p^e$ of its smallest prime factor, forms $G_p(e)=1+p^2+\cdots+p^{2e}$, and multiplies that by the previously known value of $\sigma_2(m)$ for the cofactor $m$. After the array is filled, one final pass tests each $\sigma_2(n)$ for being a square and accumulates the valid $n$.

Exact perfect-square testing

The square test must be exact. Python uses the integer square-root routine directly. The C++ and Java implementations begin with a floating-point square-root estimate, but then adjust the candidate root upward or downward until $r^2\le \sigma_2(n) < (r+1)^2$. This correction step avoids false positives from rounding and guarantees that the acceptance test really means

$$\sigma_2(n)=r^2\quad\text{for some integer }r.$$

Complexity Analysis

The three implementations share the same number theory but not the same cost profile. The recursive C++ and Python approach is output-sensitive: its running time is proportional to the number of prime-exponent states visited under the bound $n<L$. The recursion depth is the number of distinct prime factors currently chosen, so it is small in practice and at worst $O(\log L)$. Memory is dominated by the prime list plus the recursion stack.

The Java approach uses more memory but has a very regular control flow. It stores arrays indexed by all integers below $L$, so its space usage is $O(L)$. Time is dominated by the smallest-prime-factor sieve, the pass that reconstructs each $\sigma_2(n)$, and the final square test sweep. For this fixed Project Euler limit, both strategies are practical: one saves memory by exploring only factorization states, while the other spends memory to make the arithmetic on each $n$ very direct.

Footnotes and References

Problem page: Project Euler 211
Divisor function: Wikipedia - Divisor function
Multiplicative function: Wikipedia - Multiplicative function
Prime factorization: Wikipedia - Prime factor
Perfect square: Wikipedia - Square number
Sieve of Eratosthenes: Wikipedia - Sieve of Eratosthenes

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);
    }
}