Problem 561: Divisor Pairs

Project Euler Problem 561 Solution

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

Problem Summary Let $p_m\#=\prod_{i=1}^{m} p_i$ be the primorial of the first $m$ primes, let $\tau(n)$ be the divisor-counting function, and define $$S(N)=\sum_{d \mid N}\tau(d)-\tau(N),\qquad E(m,n)=v_2\!\left(S\!\left((p_m\#)^n\right)\right).$$ The required quantity is the cumulative sum $$Q(m,N)=\sum_{n=1}^{N}E(m,n).$$ For the actual problem instance, $m=904961$ and $N=10^{12}$. Because that $m$ is odd and $N$ is enormous, the solution must turn the divisor expression into a closed formula rather than evaluating every term separately. Mathematical Approach The key observation is that a primorial power has exactly $m$ distinct prime factors and all of them appear with the same exponent $n$. That symmetry makes both the divisor sum and the $2$-adic valuation collapse into elementary arithmetic. Step 1: Expand the divisor expression on a primorial power For $N=(p_m\#)^n$, each of the $m$ primes has exponent $n$, so $$\tau\!\left((p_m\#)^n\right)=(n+1)^m.$$ Every divisor is obtained by choosing one exponent from $0$ to $n$ independently for each prime. Therefore $$\sum_{d \mid (p_m\#)^n}\tau(d)=\left(\sum_{k=0}^{n}(k+1)\right)^m=\left(\frac{(n+1)(n+2)}{2}\right)^m.$$ Substituting this into the definition of $S$ gives $$S\!\left((p_m\#)^n\right)=\left(\frac{(n+1)(n+2)}{2}\right)^m-(n+1)^m.$$ Step 2: Evaluate the odd case Let $n=2k-1$....

Detailed mathematical approach

Problem Summary

Let $p_m\#=\prod_{i=1}^{m} p_i$ be the primorial of the first $m$ primes, let $\tau(n)$ be the divisor-counting function, and define

$$S(N)=\sum_{d \mid N}\tau(d)-\tau(N),\qquad E(m,n)=v_2\!\left(S\!\left((p_m\#)^n\right)\right).$$

The required quantity is the cumulative sum

$$Q(m,N)=\sum_{n=1}^{N}E(m,n).$$

For the actual problem instance, $m=904961$ and $N=10^{12}$. Because that $m$ is odd and $N$ is enormous, the solution must turn the divisor expression into a closed formula rather than evaluating every term separately.

Mathematical Approach

The key observation is that a primorial power has exactly $m$ distinct prime factors and all of them appear with the same exponent $n$. That symmetry makes both the divisor sum and the $2$-adic valuation collapse into elementary arithmetic.

Step 1: Expand the divisor expression on a primorial power

For $N=(p_m\#)^n$, each of the $m$ primes has exponent $n$, so

$$\tau\!\left((p_m\#)^n\right)=(n+1)^m.$$

Every divisor is obtained by choosing one exponent from $0$ to $n$ independently for each prime. Therefore

$$\sum_{d \mid (p_m\#)^n}\tau(d)=\left(\sum_{k=0}^{n}(k+1)\right)^m=\left(\frac{(n+1)(n+2)}{2}\right)^m.$$

Substituting this into the definition of $S$ gives

$$S\!\left((p_m\#)^n\right)=\left(\frac{(n+1)(n+2)}{2}\right)^m-(n+1)^m.$$

Step 2: Evaluate the odd case

Let $n=2k-1$. Then $n+1=2k$ and $n+2=2k+1$, so

$$S\!\left((p_m\#)^{2k-1}\right)=\bigl(k(2k+1)\bigr)^m-(2k)^m=k^m\left((2k+1)^m-2^m\right).$$

The bracketed factor is odd, because an odd number minus an even number is odd. Hence the entire $2$-adic valuation comes from the factor $k^m$:

$$E(m,2k-1)=m\,v_2(k)=m\,v_2\!\left(\frac{n+1}{2}\right).$$

Step 3: Evaluate the even case

Let $n=2r$. Then $n+1=2r+1$ is odd, and the same identity becomes

$$S\!\left((p_m\#)^{2r}\right)=\bigl((2r+1)(r+1)\bigr)^m-(2r+1)^m=(2r+1)^m\left((r+1)^m-1\right).$$

Now split according to the residue of $n$ modulo $4$.

If $r$ is odd, then $n\equiv 2 \pmod 4$ and $r+1$ is even, so $(r+1)^m-1$ is odd. Therefore

$$E(m,n)=0 \qquad \text{when } n\equiv 2 \pmod 4.$$

If $r=2t$, then $n=4t$ and $r+1=2t+1$ is odd. Since $m=904961$ is odd,

$$\begin{aligned} (2t+1)^m-1&=(2t)\left((2t+1)^{m-1}+(2t+1)^{m-2}+\cdots+1\right). \end{aligned}$$

The parenthesized sum is odd because it contains an odd number of odd terms. Hence

$$E(m,4t)=v_2(2t)=1+v_2(t)=v_2(n)-1.$$

For this odd $m$, the valuation rule is therefore

$$E(m,n)= \begin{cases} m\,v_2\!\left(\frac{n+1}{2}\right), & n \text{ odd},\\ 0, & n\equiv 2 \pmod 4,\\ v_2(n)-1, & 4\mid n. \end{cases}$$

Step 4: Sum the residue classes

We now sum $E(m,n)$ from $1$ to $N$. The class $n\equiv 2 \pmod 4$ contributes nothing, so only odd numbers and multiples of $4$ remain.

For odd $n$, write $n=2k-1$ with

$$K=\left\lfloor\frac{N+1}{2}\right\rfloor.$$

Then

$$\sum_{\substack{1\le n\le N\\ n\text{ odd}}}E(m,n)=m\sum_{k=1}^{K}v_2(k)=m\,v_2(K!).$$

For multiples of $4$, write $n=4t$ with

$$T=\left\lfloor\frac{N}{4}\right\rfloor.$$

Then

$$\sum_{\substack{1\le n\le N\\ 4\mid n}}E(m,n)=\sum_{t=1}^{T}\bigl(1+v_2(t)\bigr)=T+v_2(T!).$$

So the whole cumulative quantity reduces to

$$Q(m,N)=m\,v_2(K!)+T+v_2(T!).$$

Step 5: Replace factorial valuations with Legendre's formula

Legendre's formula for the prime $2$ says

$$v_2(x!)=\sum_{j\ge 1}\left\lfloor\frac{x}{2^j}\right\rfloor=x-\operatorname{popcount}(x).$$

This turns the final sum into direct arithmetic:

$$Q(m,N)=m\left(K-\operatorname{popcount}(K)\right)+T+\left(T-\operatorname{popcount}(T)\right).$$

No iteration up to $N$ survives the derivation.

Worked Example: $N=8$

This checkpoint is small enough to verify by hand and already shows the full mechanism. Here

$$K=\left\lfloor\frac{8+1}{2}\right\rfloor=4,\qquad T=\left\lfloor\frac{8}{4}\right\rfloor=2.$$

Using Legendre,

$$v_2(4!)=4-\operatorname{popcount}(4)=4-1=3,\qquad v_2(2!)=2-\operatorname{popcount}(2)=2-1=1.$$

Therefore

$$Q(m,8)=3m+3.$$

For the actual problem parameter $m=904961$, this becomes

$$Q(904961,8)=3\cdot 904961+3=2714886,$$

which matches the checkpoint used by the implementations.

How the Code Works

The C++, Python, and Java implementations never loop from $1$ to $N$. They apply the closed form derived above directly to the odd input $m=904961$.

Each implementation computes

$$K=\left\lfloor\frac{N+1}{2}\right\rfloor,\qquad T=\left\lfloor\frac{N}{4}\right\rfloor,$$

then evaluates $v_2(K!)$ and $v_2(T!)$ via the identity $x-\operatorname{popcount}(x)$, and finally forms

$$m\,v_2(K!)+T+v_2(T!).$$

The C++ implementation also includes two small sanity checks: one confirms that $S(6)=5$, and another confirms the checkpoint $Q(904961,8)=2714886$. After that, the exact decimal answer for $N=10^{12}$ is printed.

Complexity Analysis

For the fixed-size integer inputs used here, the algorithm runs in $O(1)$ time and $O(1)$ memory. The derivation removes every loop over $n$; only a handful of integer divisions, bit counts, and additions remain.

Footnotes and References

Problem page: https://projecteuler.net/problem=561
Primorial: Wikipedia — Primorial
Divisor function: Wikipedia — Divisor function
$p$-adic valuation: Wikipedia — $p$-adic valuation
Legendre's formula: Wikipedia — Legendre's formula

Problem 561 source code

C++

#include <algorithm>
#include <cassert>
#include <cstdint>
#include <iostream>
#include <string>

namespace {

using u64 = std::uint64_t;
using u128 = unsigned __int128;

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

static u64 v2_u64(u64 x) {
    if (x == 0) return 64;
    return static_cast<u64>(__builtin_ctzll(x));
}

// For p=2, Legendre gives v2(n!) = sum_{k>=1} floor(n/2^k) = n - popcount(n).
static u64 v2_factorial(u64 n) {
    return n - static_cast<u64>(__builtin_popcountll(n));
}

// Directly compute S((p_m#)^n) for small m,n using the divisor-structure formula:
// S(N) = sum_{d|N} tau(d) - tau(N), and for N = (p_m#)^n:
//   tau(N) = (n+1)^m
//   sum_{d|N} tau(d) = (sum_{k=0..n} (k+1))^m = ((n+1)(n+2)/2)^m
static u64 S_primorial_power_small(u64 m, u64 n) {
    const u64 t = (n + 1) * (n + 2) / 2;
    u128 A = 1, B = 1;
    for (u64 i = 0; i < m; ++i) {
        A *= static_cast<u128>(t);
        B *= static_cast<u128>(n + 1);
    }
    const u128 S = A - B;
    return static_cast<u64>(S);  // used only for tiny cases where it fits
}

// For odd m, the 2-adic valuation of S((p_m#)^n) collapses to a simple piecewise rule:
// Let m be odd and n>=1.
//   If n is odd:   E(m,n) = m * v2((n+1)/2)
//   If n≡2 (mod4): E(m,n) = 0
//   If 4|n:        E(m,n) = v2(n) - 1
static u64 E_odd_m(u64 m, u64 n) {
    assert((m & 1) == 1);
    if (n & 1) {
        return m * v2_u64((n + 1) / 2);
    }
    if ((n & 3) == 2) return 0;
    return v2_u64(n) - 1;
}

static u128 Q(u64 m, u64 N) {
    assert((m & 1) == 1);

    // Sum E(m,n) for n=1..N using the piecewise rule:
    // - odd n: contribute m*v2((n+1)/2); letting k=(n+1)/2, k=1..floor((N+1)/2) -> m*v2(K!)
    // - n≡2 (mod4): contribute 0
    // - n=4t: contribute v2(4t)-1 = 1+v2(t); sum_{t<=floor(N/4)} (1+v2(t)) = T + v2(T!)
    const u64 K = (N + 1) / 2;
    const u64 T = N / 4;

    const u128 part_odd = static_cast<u128>(m) * static_cast<u128>(v2_factorial(K));
    const u128 part_mult4 = static_cast<u128>(T) + static_cast<u128>(v2_factorial(T));
    return part_odd + part_mult4;
}

}  // namespace

int main() {
    std::ios::sync_with_stdio(false);
    std::cin.tie(nullptr);

    // Statement check: E(2,1)=0 since S(6)=5.
    {
        const u64 s = S_primorial_power_small(2, 1);
        assert(s == 5);
        assert(v2_u64(s) == 0);
    }

    // Statement check: Q(8)=2714886 for m=904961.
    {
        const u64 m = 904961;
        u128 sum = 0;
        for (u64 i = 1; i <= 8; ++i) sum += E_odd_m(m, i);
        assert(to_string_u128(sum) == "2714886");
        assert(to_string_u128(Q(m, 8)) == "2714886");
    }

    const u64 m = 904961;
    const u64 N = 1'000'000'000'000ULL;
    std::cout << to_string_u128(Q(m, N)) << '\n';
    return 0;
}

Python

def v2_factorial(n):
    return n - n.bit_count()

def Q(m, N):
    K = (N + 1) // 2
    T = N // 4
    
    part_odd = m * v2_factorial(K)
    part_mult4 = T + v2_factorial(T)
    return part_odd + part_mult4

def solve():
    m = 904961
    N = 1000000000000
    return str(Q(m, N))

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

Java

public class Euler561 {
    static long v2Factorial(long n) {
        return n - Long.bitCount(n);
    }

    static java.math.BigInteger qFunc(long m, long N) {
        long k = (N + 1) / 2;
        long t = N / 4;

        java.math.BigInteger partOdd = java.math.BigInteger.valueOf(m)
                .multiply(java.math.BigInteger.valueOf(v2Factorial(k)));

        java.math.BigInteger partMult4 = java.math.BigInteger.valueOf(t)
                .add(java.math.BigInteger.valueOf(v2Factorial(t)));

        return partOdd.add(partMult4);
    }

    public static String solve() {
        long m = 904961;
        long n = 1000000000000L;
        return qFunc(m, n).toString();
    }

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