Problem 211: Divisor Square Sum

Mathematical Development

Definition 1. For $k \in \mathbb{Z}_{\geq 0}$, the divisor power sum function is $\sigma_k(n) = \sum_{d \mid n} d^k$. The case $k = 2$ gives the sum-of-squares-of-divisors function.

Theorem 1 (Multiplicativity of $\sigma_2$). The function $\sigma_2$ is multiplicative: if $\gcd(a, b) = 1$, then $\sigma_2(ab) = \sigma_2(a) \cdot \sigma_2(b)$.

Proof. Let $a, b \in \mathbb{Z}_{>0}$ with $\gcd(a, b) = 1$. We claim the map

is a bijection.

Surjectivity. Let $d \mid ab$. Write $d = \prod_p p^{e_p}$. Since $\gcd(a, b) = 1$, each prime $p$ divides at most one of $a$ or $b$. Set

Then $d_1 \mid a$, $d_2 \mid b$, and $d = d_1 d_2$.

Injectivity. If $d_1 d_2 = d_1' d_2'$ with $d_1, d_1' \mid a$ and $d_2, d_2' \mid b$, then $\gcd(d_1, d_2') = 1$, so $d_1 \mid d_1'$. By symmetry $d_1' \mid d_1$, hence $d_1 = d_1'$ and then $d_2 = d_2'$.

Therefore

Lemma 1 (Geometric sum for prime powers). For a prime $p$ and integer $a \geq 0$,

Proof. The divisors of $p^a$ are precisely $1, p, p^2, \ldots, p^a$. Hence

which is a finite geometric series with ratio $p^2$. $\square$

Proposition 1 (Overflow bound). For $n < N = 6.4 \times 10^7$, unsigned 64-bit integers suffice to store $\sigma_2(n)$ in the sieve implementation.

Theorem 2 (Segmented factorization formula). Let

be the prime factorization of $n$. If we process each prime factor $p_i$ once, multiply an accumulator by

and remove $p_i^{a_i}$ from a working copy of $n$, then the final accumulator equals $\sigma_2(n)$. If a residual factor remains after all primes up to $\sqrt{n}$ have been removed, that residual factor is prime and contributes a final factor $1 + q^2$.

Proof. The product formula

follows immediately from Theorem 1 and Lemma 1.

Now remove the prime powers from a working copy of $n$. When all primes up to $\sqrt{n}$ have been processed, any remaining factor $q > 1$ cannot be composite, because a composite number has a prime divisor at most its square root. Therefore the remaining factor is either $1$ or a single prime, and multiplying by $1 + q^2$ completes the product for $\sigma_2(n)$. $\square$

Lemma 2 (Perfect square detection). A nonnegative integer $s$ is a perfect square if and only if $\lfloor \sqrt{s} \rfloor^2 = s$. Integer square root can be computed exactly via Newton’s method or a built-in isqrt.

Editorial

The useful structural fact is that $\sigma_2$ is multiplicative, so once the prime factorization of $n$ is known we can recover $\sigma_2(n)$ from the geometric-series factors

Factoring each integer from scratch would still be too slow, so the implementation works in blocks. Inside one block we keep the current values of the integers and an accumulator for their $\sigma_2$ product. Then, for each prime $p \le \sqrt{N}$, we walk through the multiples of $p$ inside that block, strip off the full power of $p$, and multiply the accumulator by the corresponding prime-power contribution.

After all small primes are processed, each number has at most one prime factor left. Multiplying by $1 + q^2$ finishes $\sigma_2(n)$, and then a single integer-square-root check tells us whether that value is a square. This keeps the same mathematics as the direct formula, but avoids storing a length-$64{,}000{,}000$ divisor table in Python.

Pseudocode

Set N = 64,000,000 and choose a block size B.
Generate all primes up to sqrt(N - 1).

answer = 0
For each block [L, R):
    values[i] = L + i for every index in the block
    sigma2[i] = 1 for every index in the block

    For each prime p in the prime list:
        Find the first multiple of p inside [L, R)

        For each multiple m of p inside the block:
            contribution = 1
            power = 1

            While values[m - L] is divisible by p:
                Divide values[m - L] by p
                power = power * p^2
                contribution = contribution + power

            Multiply sigma2[m - L] by contribution

    For each position i in the block:
        If values[i] > 1:
            Multiply sigma2[i] by 1 + values[i]^2

        r = integer square root of sigma2[i]
        If r * r = sigma2[i]:
            answer += L + i

Return answer

Complexity Analysis

• Time: Roughly $O(N \log \log N)$ arithmetic over all blocks. Each prime visits its multiples, and each prime factor is removed only as many times as its exponent occurs.
• Space: $O(B + \pi(\sqrt{N}))$, where $B$ is the chosen block size.

Answer

#include <bits/stdc++.h>
using namespace std;

int main() {
    const int N = 64000000;
    vector<unsigned long long> sigma2(N, 0);

    for (long long d = 1; d < N; d++) {
        long long d2 = d * d;
        for (long long k = d; k < N; k += d)
            sigma2[k] += d2;
    }

    long long answer = 0;
    for (int n = 1; n < N; n++) {
        unsigned long long s = sigma2[n];
        unsigned long long r = (unsigned long long)sqrt((double)s);
        while (r * r > s) r--;
        while ((r + 1) * (r + 1) <= s) r++;
        if (r * r == s)
            answer += n;
    }

    cout << answer << endl;
    return 0;
}

from math import isqrt


def primes_up_to(limit):
    sieve = bytearray(b"\x01") * (limit + 1)
    sieve[0:2] = b"\x00\x00"
    for p in range(2, isqrt(limit) + 1):
        if sieve[p]:
            sieve[p * p : limit + 1 : p] = b"\x00" * (((limit - p * p) // p) + 1)
    return [p for p in range(2, limit + 1) if sieve[p]]


def solve():
    limit = 64_000_000
    block_size = 1_000_000
    primes = primes_up_to(isqrt(limit - 1))

    answer = 0

    for block_start in range(1, limit, block_size):
        block_end = min(limit, block_start + block_size)
        size = block_end - block_start

        values = [block_start + i for i in range(size)]
        sigma2 = [1] * size

        for p in primes:
            p2 = p * p
            first = ((block_start + p - 1) // p) * p

            for multiple in range(first, block_end, p):
                idx = multiple - block_start
                term = 1
                power = 1

                while values[idx] % p == 0:
                    values[idx] //= p
                    power *= p2
                    term += power

                sigma2[idx] *= term

        for offset, remaining in enumerate(values):
            if remaining > 1:
                sigma2[offset] *= 1 + remaining * remaining

            total = sigma2[offset]
            root = isqrt(total)
            if root * root == total:
                answer += block_start + offset

    print(answer)


if __name__ == "__main__":
    solve()