Problem 1013: Sum 25 Divisors

Project Euler Problem 1013 Solution

Reduce both conditions to n = (pq)^4 with primes p < q that are both 1 modulo 5, then add the 33 such fourth powers below 10^15 found with a small sieve. Implementations are available in C++, Python and Java.

Detailed mathematical approach

Problem Summary

For a positive integer $n$ let $d(n)$ be the number of its divisors and $\sigma(n)$ their sum. The statement notes that $13521270961$ is the first number with $d(n)=25$ and $25\mid\sigma(n)$, and asks for the sum of all $n\lt10^{15}$ with these two properties. ^[1]

Mathematical Approach

Both conditions can be read off the prime factorization. The divisor count forces $n$ to be a fourth power of a special shape, the divisibility of $\sigma(n)$ by $25$ depends only on the residues of its primes modulo $5$, and what remains is a short sum over pairs of primes. Only $33$ numbers below $10^{15}$ qualify.

1. Exactly 25 divisors

Write $n=p_1^{a_1}p_2^{a_2}\cdots p_r^{a_r}$ with distinct primes $p_i$. Every divisor of $n$ is $p_1^{b_1}\cdots p_r^{b_r}$ with $0\le b_i\le a_i$, which gives ^[2]

$$d(n)=\prod_{i=1}^{r}(a_i+1),$$

$$\sigma(n)=\prod_{i=1}^{r}\left(1+p_i+\dots+p_i^{a_i}\right).$$

Every factor $a_i+1$ is at least $2$, and $25=5\cdot5$ splits into such factors only as $25$ or as $5\cdot5$. Hence $d(n)=25$ exactly when $n=p^{24}$ for a prime $p$, or $n=p^4q^4$ for primes $p\lt q$. In both cases $n$ is a fourth power: $p^{24}=(p^6)^4$ and $p^4q^4=(pq)^4$.

2. When 25 divides the divisor sum

For a prime $p$ put $s(p)=\sigma(p^4)$${}=1+p+p^2+p^3+p^4$. Its residue modulo $5$ depends only on the residue of $p$:

If $p\equiv1\pmod5$, write $p=1+5t$. The binomial theorem gives $p^i\equiv1+5ti\pmod{25}$, so $s(p)\equiv5+5t(0+1+2+3+4)$${}=5+50t\equiv5\pmod{25}$. Thus $5$ divides $s(p)$, but $25$ does not.
If $p=5$, then $s(p)=781\equiv1\pmod5$.
Otherwise $p\not\equiv0,1\pmod5$. By Fermat's little theorem $p^5\equiv p\pmod5$, ^[3] hence $(p-1)\,s(p)=p^5-1$${}\equiv p-1\pmod5$. Since $p-1$ is invertible modulo $5$, $s(p)\equiv1\pmod5$.

The divisor sum is multiplicative, so $\sigma(p^4q^4)=s(p)\,s(q)$. Each factor contains the prime $5$ at most once, and only when its prime is $\equiv1\pmod5$. Therefore

$$25\mid\sigma(p^4q^4)\Leftrightarrow p\equiv q\equiv1\ (\mathrm{mod}\ 5).$$

The shape $n=p^{24}$ contributes nothing: $5^{24}\approx5.96\cdot10^{16}$ already exceeds $10^{15}$, so only $p=2$ and $p=3$ remain, and $\sigma(2^{24})=2^{25}-1\equiv6\pmod{25}$ while $\sigma(3^{24})=\frac{3^{25}-1}{2}\equiv21\pmod{25}$.

3. The sum over prime pairs

The requested sum is therefore the sum of $(pq)^4$ over all pairs of primes $p\lt q$ with $p\equiv q\equiv1\pmod5$ and $(pq)^4\lt10^{15}$. Since $5623^4\approx9.997\cdot10^{14}$ and $5624^4\approx1.0004\cdot10^{15}$, the last condition means $pq\le5623$. The smallest prime $\equiv1\pmod5$ is $11$, so $q\le\lfloor5623/11\rfloor=511$, and only the primes $11,31,41,61,71,101,\dots,491$ occur. With $p=11$ there are $21$ admissible partners $q$, with $p=31$ there are $7$, with $p=41$ there are $4$, and with $p=61$ there is $1$, namely $q=71$; $p=71$ has none because $71\cdot101=7171\gt5623$. In total $33$ numbers qualify, the smallest being $341^4$, $451^4$, $671^4$, $781^4$ and $1111^4$. The first is the number from the statement: $341^4=(11\cdot31)^4$${}=13521270961$, with $\sigma(11^4)=16105=5\cdot3221$ and $\sigma(31^4)=954305$${}=5\cdot11\cdot17351$.

How the Code Works

The C++ program ^[6] follows Section 3. fourth_root(n) finds by binary search the largest $b$ with $b^4\le n$, and solve(limit) calls it with limit - 1, so that $b$ is the largest base whose fourth power lies below the limit; for the limit $10^{15}$ this is $b=5623$. A sieve of Eratosthenes up to $b$ keeps the primes $\equiv1\pmod5$. ^[4] Two nested loops over these primes add $(pq)^4$ for $p\lt q$ with $pq\le b$; the loops stop as soon as $p\gt b/p$ or $q\gt b/p$, comparisons written with integer division so that no product can overflow. Every term is below $10^{15}$ and there are $33$ of them, so the total fits easily in unsigned 64-bit arithmetic. The program rejects limits above $10^{15}$, returns $0$ for limits up to $1$, and without arguments prints the sum for the limit $10^{15}$.

qualifies(n) is an independent brute-force test: it factors $n$ by trial division, multiplies the factors $a_i+1$ and $1+p_i+\dots+p_i^{a_i}$, and checks $d(n)=25$ and $25\mid\sigma(n)$. --self-test runs run_tests. It confirms that $13521270961$ qualifies, that the bound is strict (the sum below $13521270961$ is $0$ and the sum below $13521270962$ is $13521270961$) and that empty ranges give $0$. Then it applies qualifies to every fourth power $b^4\lt10^{15}$. Since every number with $25$ divisors is a fourth power, this enumerates all candidates; at each qualifying $n$ the test requires solve(n) and solve(n + 1) to equal the running sums, and at the end it compares the full sum. The Python and Java versions mirror the C++ program function by function; Java uses Math.addExact and Math.multiplyExact, so an overflow would raise an exception instead of wrapping around.

Complexity and Verification

With $B=5623$, the sieve takes $O(B\log\log B)$ steps and the pair loop visits $33$ pairs, so the sum appears instantly. The self-test is the expensive part: trial division of $b^4$ needs at most about $b$ steps, $O(B^2)$ in total, which still takes only a fraction of a second in each language. All three programs pass the self-test with $33$ qualifying numbers and print the same sum. As an independent check, SymPy's divisor_count and divisor_sigma, applied to every fourth power below $10^{15}$, find the same $33$ numbers and the same sum, ^[5] and these numbers are exactly the products $(pq)^4$ of Section 3.

Footnotes and References

The references below supply the background facts used above. The reduction to the products $(pq)^4$ and the residues of $\sigma(p^4)$ modulo $5$ and $25$ are derived explicitly in this article.

Project Euler 1013 — Sum 25 Divisors. The official statement gives 13521270961 as the first number with exactly 25 divisors whose divisor sum is divisible by 25, and asks for the sum of all such numbers below 10¹⁵. The example tests the implementation without disclosing the requested value.
Wikipedia — Divisor function. For n = p₁^a₁ ⋯ pᵣ^aᵣ the article gives σ₀(n) = ∏(aᵢ + 1) for the number of divisors and σ₁(n) as the product of 1 + pᵢ + ⋯ + pᵢ^aᵢ for their sum, and it notes that σₓ is multiplicative. Sections 1 and 2 use these formulas.
Wikipedia — Fermat's little theorem. For a prime p and any integer a, aᵖ ≡ a (mod p). Section 2 applies it with p = 5.
Wikipedia — Sieve of Eratosthenes. The sieve lists the primes up to a bound by crossing out the multiples of each prime, starting from its square; the program uses it for the primes up to 5623.
SymPy — Number theory. divisor_count(n) returns the number of divisors of n, and divisor_sigma(n), which this page lists as moved to sympy.functions.combinatorial.numbers, returns their sum. Both were used for the independent check, which is not part of the published solutions.
C++ — Euler1013.cpp. The linked immutable C++ revision contains solve, the brute-force check qualifies and the checks run by --self-test. It is the implementation record for this article.

Problem 1013 source code

C++

#include <cstdint>
#include <cstdlib>
#include <iostream>
#include <stdexcept>
#include <string>
#include <vector>

namespace {

using u64 = std::uint64_t;
constexpr u64 LIMIT = 1'000'000'000'000'000;

constexpr u64 fourth_power(const u64 n) {
    const u64 square = n * n;
    return square * square;
}

unsigned fourth_root(const u64 n) {
    unsigned low = 0, high = 65'536;
    while (high - low > 1) {
        const unsigned middle = low + (high - low) / 2;
        if (fourth_power(middle) <= n) low = middle;
        else high = middle;
    }
    return low;
}

u64 solve(const u64 limit) {
    if (limit > LIMIT) throw std::invalid_argument("Supported limit is at most 10^15");
    if (limit <= 1) return 0;
    const unsigned bound = fourth_root(limit - 1);
    std::vector<bool> composite(bound + 1);
    std::vector<u64> primes;
    for (unsigned p = 2; p <= bound; ++p) {
        if (composite[p]) continue;
        if (p % 5 == 1) primes.push_back(p);
        for (unsigned multiple = p * p; multiple <= bound; multiple += p) composite[multiple] = true;
    }
    u64 sum = 0;
    for (std::size_t i = 0; i < primes.size(); ++i) {
        const u64 p = primes[i];
        if (p > bound / p) break;
        for (std::size_t j = i + 1; j < primes.size(); ++j) {
            const u64 q = primes[j];
            if (q > bound / p) break;
            sum += fourth_power(p * q);
        }
    }
    return sum;
}

bool qualifies(u64 n) {
    u64 count = 1, sum = 1;
    for (u64 p = 2; p <= n / p; ++p) {
        if (n % p != 0) continue;
        unsigned exponent = 0;
        u64 power = 1, factor_sum = 1;
        do {
            n /= p;
            ++exponent;
            power *= p;
            factor_sum += power;
        } while (n % p == 0);
        count *= exponent + 1;
        sum *= factor_sum;
    }
    if (n > 1) {
        count *= 2;
        sum *= n + 1;
    }
    return count == 25 && sum % 25 == 0;
}

void require(const bool condition, const std::string& description) {
    if (!condition) throw std::runtime_error("Check failed: " + description);
}

void run_tests() {
    constexpr u64 FIRST = 13'521'270'961;
    require(qualifies(FIRST), "given first number");
    require(solve(FIRST) == 0 && solve(FIRST + 1) == FIRST, "first number and strict cutoff");
    require(solve(0) == 0 && solve(1) == 0, "empty ranges");
    u64 expected = 0;
    unsigned count = 0;
    for (u64 base = 1; fourth_power(base) < LIMIT; ++base) {
        const u64 n = fourth_power(base);
        if (!qualifies(n)) continue;
        require(solve(n) == expected, "exclusive cutoff at each qualifying number");
        expected += n;
        ++count;
        require(solve(n + 1) == expected, "inclusive successor at each qualifying number");
    }
    require(solve(LIMIT) == expected, "full-range factorization check");
    std::cout << "All checks passed; " << count << " qualifying numbers.\n";
}

}

int main(int argc, char* argv[]) {
    try {
        if (argc == 2 && std::string(argv[1]) == "--self-test") {
            run_tests();
            return EXIT_SUCCESS;
        }
        if (argc != 1) throw std::invalid_argument("Usage: Euler1013 [--self-test]");
        std::cout << solve(LIMIT) << '\n';
    } catch (const std::exception& error) {
        std::cerr << error.what() << '\n';
        return EXIT_FAILURE;
    }
    return EXIT_SUCCESS;
}

Python

#!/usr/bin/env python3
"""Project Euler Problem 1013 - Sum 25 Divisors."""

import sys

LIMIT = 10**15


def fourth_power(n):
    square = n * n
    return square * square


def fourth_root(n):
    low, high = 0, 65_536
    while high - low > 1:
        middle = low + (high - low) // 2
        if fourth_power(middle) <= n:
            low = middle
        else:
            high = middle
    return low


def solve(limit):
    if limit > LIMIT:
        raise ValueError("Supported limit is at most 10^15")
    if limit <= 1:
        return 0
    bound = fourth_root(limit - 1)
    composite = bytearray(bound + 1)
    primes = []
    for p in range(2, bound + 1):
        if composite[p]:
            continue
        if p % 5 == 1:
            primes.append(p)
        for multiple in range(p * p, bound + 1, p):
            composite[multiple] = 1
    total = 0
    for i, p in enumerate(primes):
        if p > bound // p:
            break
        for q in primes[i + 1:]:
            if q > bound // p:
                break
            total += fourth_power(p * q)
    return total


def qualifies(n):
    count = total = 1
    p = 2
    while p <= n // p:
        if n % p == 0:
            exponent = 0
            power = factor_sum = 1
            while True:
                n //= p
                exponent += 1
                power *= p
                factor_sum += power
                if n % p != 0:
                    break
            count *= exponent + 1
            total *= factor_sum
        p += 1
    if n > 1:
        count *= 2
        total *= n + 1
    return count == 25 and total % 25 == 0


def require(condition, description):
    if not condition:
        raise RuntimeError("Check failed: " + description)


def run_tests():
    first = 13_521_270_961
    require(qualifies(first), "given first number")
    require(solve(first) == 0 and solve(first + 1) == first, "first number and strict cutoff")
    require(solve(0) == 0 and solve(1) == 0, "empty ranges")
    expected = 0
    count = 0
    base = 1
    while fourth_power(base) < LIMIT:
        n = fourth_power(base)
        if qualifies(n):
            require(solve(n) == expected, "exclusive cutoff at each qualifying number")
            expected += n
            count += 1
            require(solve(n + 1) == expected, "inclusive successor at each qualifying number")
        base += 1
    require(solve(LIMIT) == expected, "full-range factorization check")
    print(f"All checks passed; {count} qualifying numbers.")


def main(argv):
    try:
        if len(argv) == 2 and argv[1] == "--self-test":
            run_tests()
            return 0
        if len(argv) != 1:
            raise ValueError("Usage: Euler1013.py [--self-test]")
        print(solve(LIMIT))
    except (ValueError, RuntimeError) as error:
        print(error, file=sys.stderr)
        return 1
    return 0


if __name__ == "__main__":
    sys.exit(main(sys.argv))

Java

import java.util.ArrayList;
import java.util.List;

public class Euler1013 {
    private static final long LIMIT = 1_000_000_000_000_000L;

    private static long fourthPower(long n) {
        long square = n * n;
        return square * square;
    }

    private static int fourthRoot(long n) {
        int low = 0, high = 65_536;
        while (high - low > 1) {
            int middle = low + (high - low) / 2;
            if (fourthPower(middle) <= n) low = middle;
            else high = middle;
        }
        return low;
    }

    static long solve(long limit) {
        if (limit > LIMIT) throw new IllegalArgumentException("Supported limit is at most 10^15");
        if (limit <= 1) return 0;
        int bound = fourthRoot(limit - 1);
        boolean[] composite = new boolean[bound + 1];
        List<Long> primes = new ArrayList<>();
        for (int p = 2; p <= bound; ++p) {
            if (composite[p]) continue;
            if (p % 5 == 1) primes.add((long) p);
            for (long multiple = (long) p * p; multiple <= bound; multiple += p) composite[(int) multiple] = true;
        }
        long sum = 0;
        for (int i = 0; i < primes.size(); ++i) {
            long p = primes.get(i);
            if (p > bound / p) break;
            for (int j = i + 1; j < primes.size(); ++j) {
                long q = primes.get(j);
                if (q > bound / p) break;
                sum = Math.addExact(sum, fourthPower(p * q));
            }
        }
        return sum;
    }

    static boolean qualifies(long n) {
        long count = 1, sum = 1;
        for (long p = 2; p <= n / p; ++p) {
            if (n % p != 0) continue;
            int exponent = 0;
            long power = 1, factorSum = 1;
            do {
                n /= p;
                ++exponent;
                power *= p;
                factorSum += power;
            } while (n % p == 0);
            count *= exponent + 1;
            sum = Math.multiplyExact(sum, factorSum);
        }
        if (n > 1) {
            count *= 2;
            sum = Math.multiplyExact(sum, n + 1);
        }
        return count == 25 && sum % 25 == 0;
    }

    private static void require(boolean condition, String description) {
        if (!condition) throw new IllegalStateException("Check failed: " + description);
    }

    private static void runTests() {
        final long first = 13_521_270_961L;
        require(qualifies(first), "given first number");
        require(solve(first) == 0 && solve(first + 1) == first, "first number and strict cutoff");
        require(solve(0) == 0 && solve(1) == 0, "empty ranges");
        long expected = 0;
        int count = 0;
        for (long base = 1; fourthPower(base) < LIMIT; ++base) {
            long n = fourthPower(base);
            if (!qualifies(n)) continue;
            require(solve(n) == expected, "exclusive cutoff at each qualifying number");
            expected += n;
            ++count;
            require(solve(n + 1) == expected, "inclusive successor at each qualifying number");
        }
        require(solve(LIMIT) == expected, "full-range factorization check");
        System.out.println("All checks passed; " + count + " qualifying numbers.");
    }

    public static void main(String[] args) {
        try {
            if (args.length == 1 && args[0].equals("--self-test")) {
                runTests();
                return;
            }
            if (args.length != 0) throw new IllegalArgumentException("Usage: java Euler1013 [--self-test]");
            System.out.println(solve(LIMIT));
        } catch (RuntimeException error) {
            System.err.println(error.getMessage());
            System.exit(1);
        }
    }
}