Problem 530: GCD of Divisors

View on Project Euler

Project Euler Problem 530 Solution

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

Problem Summary The divisor-based quantity in this problem can be reindexed as an ordered factor-pair sum: $$S(N)=\sum_{n=1}^{N}\sum_{d\mid n}\gcd\!\left(d,\frac{n}{d}\right)=\sum_{ab\le N}\gcd(a,b).$$ So the task is to add \(\gcd(a,b)\) over all ordered pairs \((a,b)\) whose product is at most \(N\). At \(N=10^{15}\), direct enumeration is impossible, so the solution rewrites the sum using Euler's totient function and the divisor summatory function. Mathematical Approach We work with the equivalent form $$S(N)=\sum_{ab\le N}\gcd(a,b).$$ The efficient formula comes from separating the common-divisor contribution, then grouping repeated floor values so the large range can be handled in blocks. Step 1: Reindex the Divisor Pairs For each \(n\) and each divisor \(d\mid n\), the pair \((d,n/d)\) is an ordered factor pair whose product is exactly \(n\). Summing over all \(n\le N\) and all divisors therefore gives exactly the same set of terms as summing over all ordered pairs \((a,b)\) with \(ab\le N\). This is the first important simplification: the original divisor language becomes a two-variable summation problem with a clean product bound....

Detailed mathematical approach

Problem Summary

The divisor-based quantity in this problem can be reindexed as an ordered factor-pair sum:

$$S(N)=\sum_{n=1}^{N}\sum_{d\mid n}\gcd\!\left(d,\frac{n}{d}\right)=\sum_{ab\le N}\gcd(a,b).$$

So the task is to add \(\gcd(a,b)\) over all ordered pairs \((a,b)\) whose product is at most \(N\). At \(N=10^{15}\), direct enumeration is impossible, so the solution rewrites the sum using Euler's totient function and the divisor summatory function.

Mathematical Approach

We work with the equivalent form

$$S(N)=\sum_{ab\le N}\gcd(a,b).$$

The efficient formula comes from separating the common-divisor contribution, then grouping repeated floor values so the large range can be handled in blocks.

Step 1: Reindex the Divisor Pairs

For each \(n\) and each divisor \(d\mid n\), the pair \((d,n/d)\) is an ordered factor pair whose product is exactly \(n\). Summing over all \(n\le N\) and all divisors therefore gives exactly the same set of terms as summing over all ordered pairs \((a,b)\) with \(ab\le N\).

This is the first important simplification: the original divisor language becomes a two-variable summation problem with a clean product bound.

Step 2: Expand \(\gcd\) with Euler's Totient Function

Use the classical identity

$$\sum_{t\mid m}\varphi(t)=m.$$

Applying it with \(m=\gcd(a,b)\) gives

$$\gcd(a,b)=\sum_{t\mid \gcd(a,b)}\varphi(t).$$

Substitute this into the sum and swap the order of summation:

$$S(N)=\sum_{ab\le N}\sum_{t\mid a,\ t\mid b}\varphi(t).$$

Now write \(a=tx\) and \(b=ty\). The condition becomes

$$t^2xy\le N.$$

Hence

$$S(N)=\sum_{t\le \sqrt N}\varphi(t)\,D\!\left(\left\lfloor\frac{N}{t^2}\right\rfloor\right),$$

where \(D(m)\) counts ordered pairs \((x,y)\) with \(xy\le m\).

Step 3: Introduce the Divisor Summatory Function

The pair-counting function is

$$D(m)=\#\{(x,y)\in\mathbb{Z}_{>0}^2:xy\le m\}.$$

For a fixed \(x\), there are exactly \(\left\lfloor m/x\right\rfloor\) valid choices for \(y\). Therefore

$$D(m)=\sum_{x=1}^{m}\left\lfloor\frac{m}{x}\right\rfloor.$$

Equivalently, every integer \(n\le m\) contributes once for each divisor pair \((x,y)\) with \(xy=n\), so

$$D(m)=\sum_{n=1}^{m}\tau(n),$$

where \(\tau(n)\) is the divisor-count function. The whole problem is now a single weighted sum of totients times divisor-summatory values.

Step 4: Evaluate \(D(m)\) by Quotient Blocks

A naive loop up to \(m\) is still too slow when \(m\) is large. The key observation is that the quotient

$$q=\left\lfloor\frac{m}{\ell}\right\rfloor$$

stays constant on an interval \([\ell,r]\), where

$$r=\left\lfloor\frac{m}{q}\right\rfloor.$$

So one whole block contributes

$$q\,(r-\ell+1).$$

Advancing from one block to the next yields

$$D(m)=\sum_{\text{blocks }[\ell,r]}\left\lfloor\frac{m}{\ell}\right\rfloor(r-\ell+1),$$

and the number of blocks is only \(O(\sqrt m)\). This is the hyperbola-style speedup used whenever \(D(m)\) has to be evaluated for a large argument.

Step 5: Split the Outer Sum at \(N^{1/3}\)

Let

$$K=\left\lfloor N^{1/3}\right\rfloor,\qquad M=\left\lfloor\frac{N}{K^2}\right\rfloor.$$

Then the outer sum naturally splits into two regimes.

For \(t>K\), we have \(\left\lfloor N/t^2\right\rfloor\le M\), so all needed \(D(m)\) values are small. The implementation precomputes \(\tau(1),\dots,\tau(M)\), takes prefix sums, and obtains every small \(D(m)\) by table lookup.

For \(t\le K\), the argument \(\left\lfloor N/t^2\right\rfloor\) is large, but many consecutive \(t\) share the same value. If

$$v=\left\lfloor\frac{N}{t^2}\right\rfloor,$$

then the largest index with the same quotient is

$$r=\left\lfloor\sqrt{\frac{N}{v}}\right\rfloor.$$

That whole interval contributes

$$D(v)\sum_{u=t}^{r}\varphi(u).$$

Prefix sums of \(\varphi\) make the interval sum an \(O(1)\) lookup, so the expensive work is concentrated only in the distinct large \(D(v)\) calls.

Worked Example: \(N=10\)

Here \(\lfloor\sqrt{10}\rfloor=3\), so the transformed formula becomes

$$S(10)=\sum_{t=1}^{3}\varphi(t)\,D\!\left(\left\lfloor\frac{10}{t^2}\right\rfloor\right).$$

The required totients are

$$\varphi(1)=1,\qquad \varphi(2)=1,\qquad \varphi(3)=2.$$

The required divisor-summatory values are

$$D(10)=\sum_{x=1}^{10}\left\lfloor\frac{10}{x}\right\rfloor=27,\qquad D(2)=3,\qquad D(1)=1.$$

Therefore

$$S(10)=1\cdot 27+1\cdot 3+2\cdot 1=32,$$

which matches the checkpoint used by the implementation.

How the Code Works

The C++, Python, and Java implementations all follow the same number-theoretic formula. They first compute the integer square-root limit \(\lfloor\sqrt N\rfloor\) and the cube-root split point \(\lfloor N^{1/3}\rfloor\). Next they build a totient sieve up to \(\lfloor\sqrt N\rfloor\) and store prefix sums so that \(\sum_{u=l}^{r}\varphi(u)\) can be recovered instantly.

For the small-argument region, the implementation precomputes divisor counts up to \(M=\lfloor N/K^2\rfloor\) by visiting multiples of each divisor, then turns those counts into prefix sums to obtain \(D(1),D(2),\dots,D(M)\). For the large-argument region, it groups equal values of \(\left\lfloor N/t^2\right\rfloor\), evaluates each large \(D(v)\) with quotient blocks, and multiplies by the corresponding totient-interval sum.

The C++ and Java implementations parallelize those independent large-group contributions before adding the small tail. The Python implementation is a thin wrapper: it builds the compiled solver when needed, runs it, and returns the printed result. The compiled version also checks the two checkpoints \(S(10)=32\) and \(S(1000)=12776\) before evaluating the full target.

Complexity Analysis

Let \(L=\lfloor\sqrt N\rfloor\), \(K=\lfloor N^{1/3}\rfloor\), and \(M=\lfloor N/K^2\rfloor\). The totient sieve costs \(O(L\log\log L)\) time and \(O(L)\) memory. The small-table preparation for divisor counts costs

$$\sum_{d=1}^{M}\left\lfloor\frac{M}{d}\right\rfloor=O(M\log M)$$

time and \(O(M)\) memory.

In the large-value region there are only \(O(K)\) distinct quotient groups, and one call to the block method for \(D(v)\) costs \(O(\sqrt v)\). Summing those costs over \(t\le K\) gives

$$O\!\left(\sum_{t=1}^{K}\sqrt{\frac{N}{t^2}}\right)=O(\sqrt N\log K).$$

So the overall running time is \(O(\sqrt N\log N)\) in the worst case, while memory is dominated by the totient arrays and remains \(O(\sqrt N)\). Parallel execution improves wall-clock time for the grouped large-value part but does not change the asymptotic bound.

Footnotes and References

  1. Problem page: https://projecteuler.net/problem=530
  2. Euler's totient function and the identity \(\sum_{d\mid n}\varphi(d)=n\): Wikipedia — Euler's totient function
  3. Divisor summatory function: Wikipedia — Divisor summatory function
  4. Dirichlet hyperbola method: Wikipedia — Dirichlet hyperbola method
  5. Divisor function \(\tau(n)\): Wikipedia — Divisor function

Problem 530 source code

C++

#include <cmath>
#include <cstdint>
#include <iostream>
#include <numeric>
#include <pthread.h>
#include <string>
#include <unistd.h>
#include <vector>
#include <algorithm>
#include <functional>

namespace {

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

u64 isqrt_u64(const u64 x) {
    u64 r = static_cast<u64>(std::sqrt(static_cast<long double>(x)));
    while ((r + 1) != 0 && (r + 1) * (r + 1) <= x) {
        ++r;
    }
    while (r * r > x) {
        --r;
    }
    return r;
}

u64 icbrt_u64(const u64 x) {
    u64 r = static_cast<u64>(std::cbrt(static_cast<long double>(x)));
    while ((r + 1) != 0 && (r + 1) * (r + 1) * (r + 1) <= x) {
        ++r;
    }
    while (r * r * r > x) {
        --r;
    }
    return r;
}

std::string to_string_u128(u128 x) {
    if (x == 0) {
        return "0";
    }
    std::string s;
    while (x > 0) {
        const u64 digit = static_cast<u64>(x % 10);
        s.push_back(static_cast<char>('0' + digit));
        x /= 10;
    }
    std::reverse(s.begin(), s.end());
    return s;
}

// D(n) = sum_{m<=n} tau(m) = sum_{i=1..n} floor(n/i).
u64 divisor_summatory(const u64 n) {
    u128 res = 0;
    for (u64 l = 1; l <= n;) {
        const u64 q = n / l;
        const u64 r = n / q;
        res += static_cast<u128>(q) * static_cast<u128>(r - l + 1);
        l = r + 1;
    }
    return static_cast<u64>(res);
}

struct Part1Group {
    u64 v;
    u64 sum_phi;
};

struct Part1WorkerCtx {
    const std::vector<Part1Group>* groups;
    unsigned worker_id;
    unsigned worker_count;
    u128 partial;
};

void* part1_worker_main(void* ptr) {
    auto* ctx = static_cast<Part1WorkerCtx*>(ptr);
    u128 local = 0;
    for (std::size_t i = ctx->worker_id; i < ctx->groups->size(); i += ctx->worker_count) {
        const Part1Group& g = (*ctx->groups)[i];
        local += static_cast<u128>(g.sum_phi) * static_cast<u128>(divisor_summatory(g.v));
    }
    ctx->partial = local;
    return nullptr;
}

unsigned choose_thread_count(std::size_t tasks) {
    if (tasks <= 1U) {
        return 1U;
    }
    long cpu_count = sysconf(_SC_NPROCESSORS_ONLN);
    unsigned threads = (cpu_count > 0) ? static_cast<unsigned>(cpu_count) : 1U;
    if (threads > 8U) {
        threads = 8U;
    }
    if (threads > tasks) {
        threads = static_cast<unsigned>(tasks);
    }
    return threads == 0U ? 1U : threads;
}

std::vector<u32> totient_sieve(const u64 n) {
    std::vector<u32> phi(static_cast<std::size_t>(n + 1), 0);
    for (u64 i = 0; i <= n; ++i) {
        phi[static_cast<std::size_t>(i)] = static_cast<u32>(i);
    }
    if (n >= 1) {
        phi[1] = 1;
    }
    for (u64 p = 2; p <= n; ++p) {
        if (phi[static_cast<std::size_t>(p)] != p) {
            continue;
        }
        for (u64 m = p; m <= n; m += p) {
            const u32 cur = phi[static_cast<std::size_t>(m)];
            phi[static_cast<std::size_t>(m)] = static_cast<u32>(cur - cur / static_cast<u32>(p));
        }
    }
    return phi;
}

u128 compute_F(const u64 N) {
    const u64 limit = isqrt_u64(N);
    const u64 K = icbrt_u64(N);

    const auto phi = totient_sieve(limit);
    std::vector<u64> phi_prefix(static_cast<std::size_t>(limit + 1), 0ULL);
    for (u64 i = 1; i <= limit; ++i) {
        phi_prefix[static_cast<std::size_t>(i)] =
            phi_prefix[static_cast<std::size_t>(i - 1)] + static_cast<u64>(phi[static_cast<std::size_t>(i)]);
    }

    const u64 max_small = N / (K * K);  // ~= K
    std::vector<u32> tau(static_cast<std::size_t>(max_small + 1), 0);
    for (u64 d = 1; d <= max_small; ++d) {
        for (u64 m = d; m <= max_small; m += d) {
            ++tau[static_cast<std::size_t>(m)];
        }
    }
    std::vector<u64> D_small(static_cast<std::size_t>(max_small + 1), 0);
    for (u64 i = 1; i <= max_small; ++i) {
        D_small[static_cast<std::size_t>(i)] =
            D_small[static_cast<std::size_t>(i - 1)] + static_cast<u64>(tau[static_cast<std::size_t>(i)]);
    }

    // F(N) = sum_{t<=sqrt(N)} phi(t) * D(floor(N/t^2)).
    u128 ans = 0;

    // For t <= K, floor(N/t^2) is large; compute D via divisor_summatory.
    // Group equal quotients.
    std::vector<Part1Group> groups;
    groups.reserve(static_cast<std::size_t>(K));
    u64 t = 1;
    while (t <= K) {
        const u64 v = N / (t * t);
        u64 r = isqrt_u64(N / v);
        if (r > K) {
            r = K;
        }
        const u64 sum_phi =
            phi_prefix[static_cast<std::size_t>(r)] - phi_prefix[static_cast<std::size_t>(t - 1)];
        groups.push_back(Part1Group{v, sum_phi});
        t = r + 1;
    }

    const unsigned thread_count = choose_thread_count(groups.size());
    if (thread_count == 1U) {
        for (const Part1Group& g : groups) {
            ans += static_cast<u128>(g.sum_phi) * static_cast<u128>(divisor_summatory(g.v));
        }
    } else {
        std::vector<pthread_t> threads(thread_count);
        std::vector<Part1WorkerCtx> ctx(thread_count);
        unsigned created = 0U;
        bool failed = false;
        for (unsigned w = 0U; w < thread_count; ++w) {
            ctx[w] = Part1WorkerCtx{&groups, w, thread_count, 0};
            if (pthread_create(&threads[w], nullptr, part1_worker_main, &ctx[w]) != 0) {
                failed = true;
                break;
            }
            ++created;
        }
        for (unsigned w = 0U; w < created; ++w) {
            pthread_join(threads[w], nullptr);
        }
        if (failed) {
            for (const Part1Group& g : groups) {
                ans += static_cast<u128>(g.sum_phi) * static_cast<u128>(divisor_summatory(g.v));
            }
        } else {
            for (const Part1WorkerCtx& w : ctx) {
                ans += w.partial;
            }
        }
    }

    // For t > K, v <= max_small, so D can be tabled.
    for (u64 i = K + 1; i <= limit; ++i) {
        const u64 v = N / (i * i);
        ans += static_cast<u128>(phi[static_cast<std::size_t>(i)]) *
               static_cast<u128>(D_small[static_cast<std::size_t>(v)]);
    }

    return ans;
}

bool run_checkpoints() {
    if (compute_F(10) != static_cast<u128>(32)) {
        std::cerr << "Checkpoint failed: F(10)\n";
        return false;
    }
    if (compute_F(1000) != static_cast<u128>(12776)) {
        std::cerr << "Checkpoint failed: F(1000)\n";
        return false;
    }
    return true;
}

}  // namespace

int main() {
    if (!run_checkpoints()) {
        return 1;
    }

    constexpr u64 N = 1'000'000'000'000'000ULL;
    const u128 ans = compute_F(N);
    std::cout << to_string_u128(ans) << '\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.List;

public class Euler530 {

    static long isqrt(long x) {
        if (x == 0)
            return 0;
        long r = (long) Math.sqrt(x);
        while ((r + 1) * (r + 1) <= x)
            r++;
        while (r * r > x)
            r--;
        return r;
    }

    static long icbrt(long x) {
        if (x == 0)
            return 0;
        long r = (long) Math.cbrt(x);
        while ((r + 1) * (r + 1) * (r + 1) <= x)
            r++;
        while (r * r * r > x)
            r--;
        return r;
    }

    static long divisorSummatory(long n) {
        long res = 0;
        for (long l = 1; l <= n;) {
            long q = n / l;
            long r = n / q;
            res += q * (r - l + 1);
            l = r + 1;
        }
        return res;
    }

    static class Part1Group {
        long v;
        long sumPhi;

        Part1Group(long v, long sumPhi) {
            this.v = v;
            this.sumPhi = sumPhi;
        }
    }

    public static void main(String[] args) {
        long N = 1000000000000000L;
        long limit = isqrt(N);
        long K = icbrt(N);

        int[] phi = new int[(int) limit + 1];
        for (int i = 0; i <= limit; i++)
            phi[i] = i;
        if (limit >= 1)
            phi[1] = 1;

        for (int p = 2; p <= limit; p++) {
            if (phi[p] == p) {
                for (int m = p; m <= limit; m += p) {
                    phi[m] -= phi[m] / p;
                }
            }
        }

        long[] phiPrefix = new long[(int) limit + 1];
        long sum = 0;
        for (int i = 1; i <= limit; i++) {
            sum += phi[i];
            phiPrefix[i] = sum;
        }

        long maxSmall = N / (K * K);
        int[] tau = new int[(int) maxSmall + 1];
        for (int d = 1; d <= maxSmall; d++) {
            for (int m = d; m <= maxSmall; m += d) {
                tau[m]++;
            }
        }

        long[] DSmall = new long[(int) maxSmall + 1];
        sum = 0;
        for (int i = 1; i <= maxSmall; i++) {
            sum += tau[i];
            DSmall[i] = sum;
        }

        List<Part1Group> groups = new ArrayList<>();
        long t = 1;
        while (t <= K) {
            long v = N / (t * t);
            if (v == 0)
                break;
            long r = isqrt(N / v);
            if (r > K)
                r = K;
            long sumPhi = phiPrefix[(int) r] - phiPrefix[(int) (t - 1)];
            groups.add(new Part1Group(v, sumPhi));
            t = r + 1;
        }

        long ansPart1 = groups.parallelStream().mapToLong(g -> {
            return g.sumPhi * divisorSummatory(g.v);
        }).sum();

        long ansPart2 = 0;
        for (int i = (int) (K + 1); i <= limit; i++) {
            long v = N / ((long) i * i);
            ansPart2 += (long) phi[i] * DSmall[(int) v];
        }

        System.out.println(ansPart1 + ansPart2);
    }
}