Problem 439: Sum of Sum of Divisors

Project Euler Problem 439 Solution

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

Problem Summary Let $\sigma(n)$ denote the sum of all positive divisors of $n$. The problem asks for $$S(N)=\sum_{x=1}^{N}\sum_{y=1}^{N}\sigma(xy)$$ with $N=10^{11}$, reduced modulo $10^9$. A direct double loop is hopeless, so the solution rewrites $\sigma(xy)$ into a form that separates the common part of $x$ and $y$, then evaluates the resulting sums by quotient blocks. Mathematical Approach Step 1: A Möbius Identity for $\sigma(xy)$ The key identity used by the implementation is $$\sigma(xy)=\sum_{d\mid \gcd(x,y)} \mu(d)\,d\,\sigma\!\left(\frac{x}{d}\right)\sigma\!\left(\frac{y}{d}\right),$$ where $\mu$ is the Möbius function. Both sides are multiplicative in the pair $(x,y)$, so it is enough to check prime powers. Let $x=p^a$ and $y=p^b$. If $a=0$ or $b=0$, then $\gcd(x,y)=1$ and the identity is immediate. For $a,b\ge 1$, only $d=1$ and $d=p$ contribute, because $\mu(p^k)=0$ for $k\ge 2$. Thus $$\begin{aligned} \sum_{d\mid \gcd(p^a,p^b)} \mu(d)\,d\,\sigma\!\left(\frac{p^a}{d}\right)\sigma\!\left(\frac{p^b}{d}\right) &=\sigma(p^a)\sigma(p^b)-p\,\sigma(p^{a-1})\sigma(p^{b-1})\\ &=\frac{(p^{a+1}-1)(p^{b+1}-1)-p(p^a-1)(p^b-1)}{(p-1)^2}\\ &=\frac{p^{a+b+1}-1}{p-1} =\sigma(p^{a+b}). \end{aligned}$$ Since $p^{a+b}=xy$, the identity holds for prime powers and therefore for all $x,y$....

Detailed mathematical approach

Problem Summary

Let $\sigma(n)$ denote the sum of all positive divisors of $n$. The problem asks for

$$S(N)=\sum_{x=1}^{N}\sum_{y=1}^{N}\sigma(xy)$$

with $N=10^{11}$, reduced modulo $10^9$. A direct double loop is hopeless, so the solution rewrites $\sigma(xy)$ into a form that separates the common part of $x$ and $y$, then evaluates the resulting sums by quotient blocks.

Mathematical Approach

Step 1: A Möbius Identity for $\sigma(xy)$

The key identity used by the implementation is

$$\sigma(xy)=\sum_{d\mid \gcd(x,y)} \mu(d)\,d\,\sigma\!\left(\frac{x}{d}\right)\sigma\!\left(\frac{y}{d}\right),$$

where $\mu$ is the Möbius function. Both sides are multiplicative in the pair $(x,y)$, so it is enough to check prime powers. Let $x=p^a$ and $y=p^b$.

If $a=0$ or $b=0$, then $\gcd(x,y)=1$ and the identity is immediate. For $a,b\ge 1$, only $d=1$ and $d=p$ contribute, because $\mu(p^k)=0$ for $k\ge 2$. Thus

$$\begin{aligned} \sum_{d\mid \gcd(p^a,p^b)} \mu(d)\,d\,\sigma\!\left(\frac{p^a}{d}\right)\sigma\!\left(\frac{p^b}{d}\right) &=\sigma(p^a)\sigma(p^b)-p\,\sigma(p^{a-1})\sigma(p^{b-1})\\ &=\frac{(p^{a+1}-1)(p^{b+1}-1)-p(p^a-1)(p^b-1)}{(p-1)^2}\\ &=\frac{p^{a+b+1}-1}{p-1} =\sigma(p^{a+b}). \end{aligned}$$

Since $p^{a+b}=xy$, the identity holds for prime powers and therefore for all $x,y$.

Step 2: Collapse the Double Sum

Insert the identity into the definition of $S(N)$:

$$S(N)=\sum_{x=1}^{N}\sum_{y=1}^{N}\sum_{d\mid \gcd(x,y)} \mu(d)\,d\,\sigma\!\left(\frac{x}{d}\right)\sigma\!\left(\frac{y}{d}\right).$$

Exchange the order of summation and write $x=da$, $y=db$. Then $a,b\le \lfloor N/d\rfloor$, so

$$S(N)=\sum_{d=1}^{N}\mu(d)\,d\left(\sum_{a\le N/d}\sigma(a)\right)\left(\sum_{b\le N/d}\sigma(b)\right).$$

Define the summatory divisor function

$$A(M)=\sum_{n=1}^{M}\sigma(n).$$

This gives the single outer sum

$$\boxed{S(N)=\sum_{d=1}^{N}\mu(d)\,d\,A\!\left(\left\lfloor\frac{N}{d}\right\rfloor\right)^2.}$$

Step 3: Fast Evaluation of $A(M)$

We start from the divisor definition:

$$A(M)=\sum_{n=1}^{M}\sum_{e\mid n} e=\sum_{e=1}^{M} e\left\lfloor\frac{M}{e}\right\rfloor.$$

The implementations use an equivalent form that is better suited to quotient blocking. Let

$$T(q)=\frac{q(q+1)}{2}.$$

Then

$$A(M)=\sum_{t=1}^{M} T\!\left(\left\lfloor\frac{M}{t}\right\rfloor\right),$$

because for each $t$, the inner sum $\sum_{e\le M/t} e$ is exactly $T(\lfloor M/t\rfloor)$. Now $\left\lfloor M/t\right\rfloor$ is constant on intervals $l\le t\le r$ with

$$q=\left\lfloor\frac{M}{l}\right\rfloor,\qquad r=\left\lfloor\frac{M}{q}\right\rfloor.$$

So one whole block contributes

$$\sum_{t=l}^{r} T\!\left(\left\lfloor\frac{M}{t}\right\rfloor\right)=(r-l+1)\,T(q).$$

There are only $O(\sqrt{M})$ such blocks, which is the standard harmonic-sum speedup.

Step 4: Prefix Sums of $d\,\mu(d)$

The outer formula needs interval sums of $\mu(d)\,d$. Define

$$W(n)=\sum_{k=1}^{n} k\,\mu(k).$$

The implementation computes small values by sieve and large values by memoized recursion. The recurrence comes from the identity

$$\sum_{k=1}^{n} k\,W\!\left(\left\lfloor\frac{n}{k}\right\rfloor\right)=1.$$

Indeed, after expanding $W$,

$$\sum_{k=1}^{n} k\sum_{m\le n/k} m\,\mu(m)=\sum_{km\le n} km\,\mu(m)=\sum_{t\le n} t\sum_{m\mid t}\mu(m)=1,$$

because $\sum_{m\mid t}\mu(m)=0$ unless $t=1$. Therefore

$$W(n)=1-\sum_{k=2}^{n} k\,W\!\left(\left\lfloor\frac{n}{k}\right\rfloor\right).$$

Again, equal quotients are grouped into intervals. For any block $l\le k\le r$ with $\lfloor n/k\rfloor=v$, we replace the whole block by

$$\left(\sum_{k=l}^{r} k\right)W(v).$$

Step 5: Final Harmonic Decomposition

The same block idea is applied to the final sum for $S(N)$. If $\lfloor N/d\rfloor=q$ on $l\le d\le r$, then the entire interval contributes

$$A(q)^2\sum_{d=l}^{r} d\,\mu(d)=A(q)^2\bigl(W(r)-W(l-1)\bigr).$$

So the computation walks over the $O(\sqrt N)$ distinct quotient values of $\lfloor N/d\rfloor$, evaluates $A(q)$ once per block, and weights it by the corresponding difference of weighted Möbius prefixes.

Worked Example: $N=3$

For $N=3$, the divisor sums are $\sigma(1)=1$, $\sigma(2)=3$, $\sigma(3)=4$, so

$$A(3)=1+3+4=8,\qquad A(1)=1.$$

The outer formula gives

$$S(3)=\mu(1)\cdot 1\cdot A(3)^2+\mu(2)\cdot 2\cdot A(1)^2+\mu(3)\cdot 3\cdot A(1)^2,$$

hence

$$S(3)=1\cdot 64-2\cdot 1-3\cdot 1=59,$$

which matches the small validation target used by the implementations.

How the Code Works

The C++, Python, and Java implementations all follow the same structure. First they precompute Möbius values up to a practical cutoff and store the prefix sums of $k\mu(k)$. For larger arguments they reuse the recurrence for $W(n)$ with memoization, so repeated quotient values are solved only once.

Next they enumerate the quotient blocks of $\lfloor N/d\rfloor$. For each block they evaluate $A(q)$ by the triangular-number block formula, square it modulo $10^9$, multiply by the weighted Möbius difference for that block, and add the result to the running total. The large independent blocks can also be accumulated in parallel.

Complexity Analysis

The direct definition needs $N^2$ evaluations of $\sigma(xy)$, which is completely infeasible. After the Möbius decomposition, the outer sum has only $O(\sqrt N)$ distinct quotient blocks. One evaluation of $A(M)$ costs $O(\sqrt M)$, and summing this over the distinct values $M=\lfloor N/d\rfloor$ gives about $O(N^{3/4})$ total block work. The weighted Möbius prefix is also handled sublinearly by combining a sieve for small values with memoized quotient recursion for large ones. Memory usage is linear in the chosen sieve cutoff, plus $O(\sqrt N)$ cached block data.

Footnotes and References

Problem page: https://projecteuler.net/problem=439
Divisor function $\sigma(n)$: Wikipedia — Divisor function
Möbius function and inversion: Wikipedia — Möbius function
Dirichlet convolution: Wikipedia — Dirichlet convolution
Harmonic-lemma style floor blocking: cp-algorithms — Number of divisors / sum of divisors

Problem 439 source code

C++

#include <algorithm>
#include <atomic>
#include <chrono>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <string>
#include <thread>
#include <unordered_map>
#include <vector>

using int64 = long long;
using i128 = __int128_t;
using u128 = __uint128_t;

static constexpr int64 kMod = 1000000000LL;
static constexpr int64 kDefaultN = 100000000000LL;
static constexpr int kSieveLimit = 20000000;

static inline int64 mod_norm(int64 x) {
  x %= kMod;
  if (x < 0)
    x += kMod;
  return x;
}

static inline int64 mod_mul(int64 a, int64 b) {
  return static_cast<int64>((static_cast<i128>(a) * b) % kMod);
}

static inline int64 sum_arith_mod(int64 l, int64 r) {
  u128 t = static_cast<u128>(l + r) * static_cast<u128>(r - l + 1) / 2;
  return static_cast<int64>(t % kMod);
}

// A(M) = sum_{n<=M} sigma(n) = sum_{k<=M} k * floor(M/k).
static int64 sum_sigma_prefix(int64 M) {
  int64 total = 0;
  for (int64 l = 1, r; l <= M; l = r + 1) {
    int64 q = M / l;
    r = M / q;
    u128 t = static_cast<u128>(q) * static_cast<u128>(q + 1) / 2;
    t *= static_cast<u128>(r - l + 1);
    total += static_cast<int64>(t % kMod);
    if (total >= kMod)
      total %= kMod;
  }
  return total % kMod;
}

// Summatory of k * mu(k) with memoization.
class MobiusPrefix {
public:
  explicit MobiusPrefix(int64 max_n, size_t cache_hint) {
    limit_ = static_cast<int>(max_n);
    mu_.assign(limit_ + 1, 0);
    prefix_.assign(limit_ + 1, 0);
    cache_.reserve(cache_hint);
    build();
  }

  int64 sum_prefix(int64 n) {
    if (n <= 0)
      return 0;
    if (n <= limit_)
      return prefix_[static_cast<size_t>(n)];
    auto it = cache_.find(n);
    if (it != cache_.end())
      return it->second;

    int64 ans = 1;
    for (int64 l = 2, r; l <= n; l = r + 1) {
      int64 v = n / l;
      r = n / v;
      int64 sum_lr = sum_arith_mod(l, r);
      int64 sub = mod_mul(sum_lr, sum_prefix(v));
      ans -= sub;
      if (ans < 0)
        ans += kMod;
    }
    cache_[n] = static_cast<int32_t>(ans);
    return ans;
  }

private:
  int limit_ = 0;
  std::vector<int8_t> mu_;
  std::vector<int32_t> prefix_;
  std::unordered_map<int64, int32_t> cache_;

  void build() {
    std::vector<int> primes;
    primes.reserve(limit_ / 10);
    std::vector<uint8_t> is_comp(limit_ + 1, 0);
    mu_[1] = 1;
    for (int i = 2; i <= limit_; ++i) {
      if (!is_comp[i]) {
        primes.push_back(i);
        mu_[i] = -1;
      }
      for (int p : primes) {
        int64 v = static_cast<int64>(i) * p;
        if (v > limit_)
          break;
        is_comp[static_cast<size_t>(v)] = 1;
        if (i % p == 0) {
          mu_[static_cast<size_t>(v)] = 0;
          break;
        }
        mu_[static_cast<size_t>(v)] = static_cast<int8_t>(-mu_[i]);
      }
    }
    int64 running = 0;
    for (int i = 1; i <= limit_; ++i) {
      running += static_cast<int64>(i) * mu_[i];
      running = mod_norm(running);
      prefix_[static_cast<size_t>(i)] = static_cast<int32_t>(running);
    }
  }
};

struct RangeData {
  int64 l;
  int64 r;
  int64 m;
  int64 mu_sum;
};

// Uses identity:
//   sigma(xy) = sum_{d|gcd(x,y)} mu(d) * d * sigma(x/d) * sigma(y/d).
// Summing over i,j<=N gives S(N)=sum_{d<=N} mu(d)*d*A(N/d)^2.
static int64 solve(int64 N, MobiusPrefix &pref, int threads) {
  std::vector<RangeData> ranges;
  ranges.reserve(
      static_cast<size_t>(2.0 * std::sqrt(static_cast<long double>(N)) + 10));
  for (int64 l = 1, r; l <= N; l = r + 1) {
    int64 q = N / l;
    r = N / q;
    ranges.push_back({l, r, q, 0});
  }

  int64 prev_prefix = 0;
  for (auto &range : ranges) {
    int64 pref_r = pref.sum_prefix(range.r);
    int64 pref_lm1 = (range.l == 1) ? 0 : prev_prefix;
    int64 mu_sum = pref_r - pref_lm1;
    if (mu_sum < 0)
      mu_sum += kMod;
    range.mu_sum = mu_sum;
    prev_prefix = pref_r;
  }

  if (threads < 1)
    threads = 1;
  threads = std::min<int>(threads, static_cast<int>(ranges.size()));
  std::atomic<size_t> next_idx(0);
  std::vector<int64> partial(static_cast<size_t>(threads), 0);

  auto worker = [&](int tid) {
    int64 local = 0;
    while (true) {
      size_t idx = next_idx.fetch_add(1, std::memory_order_relaxed);
      if (idx >= ranges.size())
        break;
      const auto &range = ranges[idx];
      if (range.mu_sum == 0)
        continue;
      int64 a = sum_sigma_prefix(range.m);
      int64 term = mod_mul(range.mu_sum, mod_mul(a, a));
      local += term;
      if (local >= kMod)
        local %= kMod;
    }
    partial[static_cast<size_t>(tid)] = local % kMod;
  };

  std::vector<std::thread> pool;
  pool.reserve(static_cast<size_t>(threads));
  for (int t = 0; t < threads; ++t)
    pool.emplace_back(worker, t);
  for (auto &th : pool)
    th.join();

  int64 total = 0;
  for (int64 v : partial) {
    total += v;
    if (total >= kMod)
      total %= kMod;
  }
  return total % kMod;
}

static bool validate(MobiusPrefix &pref, int threads) {
  struct Test {
    int64 N;
    int64 expected_mod;
  };
  const std::vector<Test> tests = {
      {3, 59},
      {1000, 563576517282LL % kMod},
      {100000, 215766508},
  };
  bool ok = true;
  for (const auto &test : tests) {
    int64 got = solve(test.N, pref, threads);
    std::cout << "Validation S(" << test.N << ") mod 1e9: " << got;
    if (got == test.expected_mod) {
      std::cout << " [PASS]";
    } else {
      std::cout << " [FAIL] Expected " << test.expected_mod;
      ok = false;
    }
    std::cout << "\n";
  }
  return ok;
}

int main(int argc, char **argv) {
  int64 N = kDefaultN;
  int threads = static_cast<int>(std::thread::hardware_concurrency());
  if (threads == 0)
    threads = 4;
  bool run_validation = true;

  for (int i = 1; i < argc; ++i) {
    std::string arg = argv[i];
    if (arg.rfind("--threads=", 0) == 0) {
      threads = std::max(1, std::stoi(arg.substr(10)));
    } else if (arg == "--no-validate") {
      run_validation = false;
    } else {
      N = std::stoll(arg);
    }
  }

  int64 sieve_limit = std::min<int64>(kSieveLimit, N);
  size_t cache_hint =
      static_cast<size_t>(2.0 * std::sqrt(static_cast<long double>(N)) + 1000);
  MobiusPrefix pref(sieve_limit, cache_hint);

  if (run_validation) {
    std::cout << "Running validations...\n";
    if (!validate(pref, std::min(threads, 4))) {
      std::cout << "Validation failed. Aborting.\n";
      return 1;
    }
    std::cout << "--------------------------------\n";
  }

  std::cout << "Solving for N = " << N << "...\n";
  auto start = std::chrono::high_resolution_clock::now();
  int64 result = solve(N, pref, threads);
  auto end = std::chrono::high_resolution_clock::now();
  std::chrono::duration<double> elapsed = end - start;

  std::cout << "Result S(N) mod 1e9: " << result << "\n";
  std::cout << "Answer: " << result << "\n";
  // std::cout << "Time elapsed: " << elapsed.count() << "s\n";
  return 0;
}

Python

from __future__ import annotations

import re
import shutil
import subprocess
from pathlib import Path

ANSWER_RE = re.compile(r"answer\s*:\s*(.+)$", re.IGNORECASE)
EQUAL_RE = re.compile(r"=\s*(.+)$")


def parse_output(stdout: str) -> str:
    lines = [line.strip() for line in stdout.splitlines() if line.strip()]
    if not lines:
        return ""
    answers = []
    equals = []
    for line in lines:
        m1 = ANSWER_RE.search(line)
        if m1:
            answers.append(m1.group(1).strip())
        m2 = EQUAL_RE.search(line)
        if m2:
            equals.append(m2.group(1).strip())
    if answers:
        return answers[-1]
    if equals:
        return equals[-1]
    return lines[-1]


def should_skip_cpp_checkpoints(src: Path) -> bool:
    try:
        text = src.read_text(encoding="utf-8", errors="ignore")
    except OSError:
        return False
    return "--skip-checkpoints" in text


def run_cpp(binary: Path, src: Path, root: Path) -> str:
    cmd = [str(binary)]
    if should_skip_cpp_checkpoints(src):
        cmd.append("--skip-checkpoints")

    try:
        return subprocess.check_output(cmd, text=True, cwd=root)
    except subprocess.CalledProcessError:
        return subprocess.check_output(cmd, text=True, cwd=src.parent)


def solve() -> str:
    problem_id = __file__.split("Euler")[-1].split(".")[0]
    root = Path(__file__).resolve().parent.parent
    src = root / "solutionsCpp" / f"Euler{problem_id}.cpp"
    binary = root / "solutionsCpp" / f".euler{problem_id}_py_bridge"

    if not binary.exists() or src.stat().st_mtime > binary.stat().st_mtime:
        compiler = shutil.which("clang++") or shutil.which("g++")
        if not compiler:
            raise RuntimeError("No C++ compiler found (clang++/g++).")
        subprocess.check_call([compiler, "-std=c++17", "-O2", str(src), "-o", str(binary)])

    output = run_cpp(binary=binary, src=src, root=root)
    parsed = parse_output(output)
    if not parsed:
        raise RuntimeError(f"Euler{problem_id} bridge produced empty output.")
    return parsed


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

Java

import java.util.ArrayList;
import java.util.HashMap;
import java.util.List;
import java.util.Map;
import java.util.stream.IntStream;

public class Euler439 {
    static final long MOD = 1000000000L;
    static final long N = 100000000000L;
    static final int SIEVE_LIMIT = 5000000;

    static long sumArithMod(long l, long r) {
        long count = r - l + 1;
        long sum = l + r;
        long val;
        if (count % 2 == 0) {
            val = ((count / 2) % MOD) * (sum % MOD) % MOD;
        } else {
            val = (count % MOD) * ((sum / 2) % MOD) % MOD;
        }
        return val;
    }

    static long sumSigmaPrefix(long M) {
        long total = 0;
        long l = 1;
        while (l <= M) {
            long q = M / l;
            long r = M / q;
            long q_sum;
            if (q % 2 == 0) {
                q_sum = ((q / 2) % MOD) * ((q + 1) % MOD) % MOD;
            } else {
                q_sum = (q % MOD) * (((q + 1) / 2) % MOD) % MOD;
            }

            long len = (r - l + 1) % MOD;
            total = (total + q_sum * len % MOD) % MOD;
            l = r + 1;
        }
        return total;
    }

    static class MobiusPrefix {
        int limit;
        byte[] mu;
        int[] prefix;
        Map<Long, Long> cache = new HashMap<>();

        MobiusPrefix(int limit) {
            this.limit = limit;
            mu = new byte[limit + 1];
            prefix = new int[limit + 1];
            build();
        }

        void build() {
            byte[] isComp = new byte[limit + 1];
            List<Integer> primes = new ArrayList<>();
            mu[1] = 1;
            for (int i = 2; i <= limit; i++) {
                if (isComp[i] == 0) {
                    primes.add(i);
                    mu[i] = -1;
                }
                for (int p : primes) {
                    long v = (long) i * p;
                    if (v > limit)
                        break;
                    isComp[(int) v] = 1;
                    if (i % p == 0) {
                        mu[(int) v] = 0;
                        break;
                    }
                    mu[(int) v] = (byte) -mu[i];
                }
            }
            long running = 0;
            for (int i = 1; i <= limit; i++) {
                long term = ((long) i * mu[i]) % MOD;
                if (term < 0)
                    term += MOD;
                running = (running + term) % MOD;
                prefix[i] = (int) running;
            }
        }

        long sumPrefix(long n) {
            if (n <= 0)
                return 0;
            if (n <= limit)
                return prefix[(int) n];
            if (cache.containsKey(n))
                return cache.get(n);

            long ans = 1;
            long l = 2;
            while (l <= n) {
                long v = n / l;
                long r = n / v;
                long sumLr = sumArithMod(l, r);
                long sub = (sumLr * sumPrefix(v)) % MOD;
                ans = (ans - sub + MOD) % MOD;
                l = r + 1;
            }

            cache.put(n, ans);
            return ans;
        }
    }

    static class RangeData {
        long l, r, m, muSum;

        RangeData(long l, long r, long m) {
            this.l = l;
            this.r = r;
            this.m = m;
        }
    }

    public static String solve() {
        MobiusPrefix pref = new MobiusPrefix(SIEVE_LIMIT);

        List<RangeData> ranges = new ArrayList<>();
        long l = 1;
        while (l <= N) {
            long q = N / l;
            long r = N / q;
            ranges.add(new RangeData(l, r, q));
            l = r + 1;
        }

        long prevPrefix = 0;
        List<RangeData> validRanges = new ArrayList<>();
        for (RangeData rng : ranges) {
            long prefR = pref.sumPrefix(rng.r);
            long prefLm1 = (rng.l == 1) ? 0 : prevPrefix;
            long muSum = (prefR - prefLm1 + MOD) % MOD;
            prevPrefix = prefR;

            if (muSum != 0) {
                rng.muSum = muSum;
                validRanges.add(rng);
            }
        }

        int numThreads = Math.min(Runtime.getRuntime().availableProcessors(), validRanges.size());
        if (numThreads < 1)
            numThreads = 1;

        final int threads = numThreads;
        long total = IntStream.range(0, threads).parallel().mapToLong(i -> {
            int begin = i * validRanges.size() / threads;
            int end = (i + 1) * validRanges.size() / threads;
            long localSum = 0;
            for (int k = begin; k < end; k++) {
                RangeData rng = validRanges.get(k);
                long a = sumSigmaPrefix(rng.m);
                long term = rng.muSum * a % MOD * a % MOD;
                localSum = (localSum + term) % MOD;
            }
            return localSum;
        }).sum();

        return Long.toString(total % MOD);
    }

    public static void main(String[] args) {
        System.out.println(solve());
    }
}