Problem 501: Eight Divisors

View on Project Euler

Project Euler Problem 501 Solution

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

Problem Summary We want the number of integers \(x\le n\) that have exactly eight positive divisors. For the target scale \(n=10^{12}\), testing each integer separately would be far too slow, so the solution classifies every valid number by its prime-factor pattern and turns the problem into a small set of prime-counting formulas. Mathematical Approach Let $$f(n)=\#\{x\in \mathbb{Z}_{>0}:x\le n,\ d(x)=8\},$$ where \(d(x)\) is the divisor-counting function. If $$x=\prod_i p_i^{a_i},$$ then $$d(x)=\prod_i (a_i+1).$$ So the entire problem is to find all exponent patterns whose contribution to \(d(x)\) is exactly \(8\), and then count how many numbers of each pattern lie below \(n\). Step 1: Classify all numbers with exactly eight divisors The multiplicative partitions of \(8\) are $$8,\qquad 4\cdot 2,\qquad 2\cdot 2\cdot 2.$$ Therefore every integer with exactly eight divisors has one of the three mutually exclusive shapes $$x=p^7,\qquad x=p^3q,\qquad x=pqr,$$ where \(p,q,r\) are primes, \(p\ne q\) in the middle case, and \(p\lt q\lt r\) in the last case. No other factorization type is possible, because any other exponent list would make \(\prod (a_i+1)\neq 8\)....

Detailed mathematical approach

Problem Summary

We want the number of integers \(x\le n\) that have exactly eight positive divisors. For the target scale \(n=10^{12}\), testing each integer separately would be far too slow, so the solution classifies every valid number by its prime-factor pattern and turns the problem into a small set of prime-counting formulas.

Mathematical Approach

Let

$$f(n)=\#\{x\in \mathbb{Z}_{>0}:x\le n,\ d(x)=8\},$$

where \(d(x)\) is the divisor-counting function. If

$$x=\prod_i p_i^{a_i},$$

then

$$d(x)=\prod_i (a_i+1).$$

So the entire problem is to find all exponent patterns whose contribution to \(d(x)\) is exactly \(8\), and then count how many numbers of each pattern lie below \(n\).

Step 1: Classify all numbers with exactly eight divisors

The multiplicative partitions of \(8\) are

$$8,\qquad 4\cdot 2,\qquad 2\cdot 2\cdot 2.$$

Therefore every integer with exactly eight divisors has one of the three mutually exclusive shapes

$$x=p^7,\qquad x=p^3q,\qquad x=pqr,$$

where \(p,q,r\) are primes, \(p\ne q\) in the middle case, and \(p\lt q\lt r\) in the last case. No other factorization type is possible, because any other exponent list would make \(\prod (a_i+1)\neq 8\).

Step 2: Count the three-prime case \(pqr\)

Fix ordered primes \(p\lt q\lt r\) with

$$pqr\le n.$$

For a fixed pair \((p,q)\), the third prime must satisfy

$$q\lt r\le \left\lfloor\frac{n}{pq}\right\rfloor.$$

Hence the number of admissible \(r\) values is

$$\pi\!\left(\left\lfloor\frac{n}{pq}\right\rfloor\right)-\pi(q),$$

where \(\pi(x)\) is the prime-counting function. Also, once \(q^2>n/p\), even the smallest possible \(r\) is too large, so \(q\) only needs to run up to \(\left\lfloor\sqrt{n/p}\right\rfloor\). This gives

$$A(n)=\sum_{\substack{p\lt q\\ pq^2\le n}}\left(\pi\!\left(\left\lfloor\frac{n}{pq}\right\rfloor\right)-\pi(q)\right).$$

Step 3: Count the mixed-power case \(p^3q\)

Now fix a prime \(p\). Any prime \(q\) with

$$q\le \left\lfloor\frac{n}{p^3}\right\rfloor$$

produces a candidate \(p^3q\le n\). If we temporarily ignore the restriction \(q\ne p\), the contribution from this \(p\) is simply \(\pi(\lfloor n/p^3\rfloor)\). Summing over all possible base primes yields the raw count

$$B_{\mathrm{raw}}(n)=\sum_{p\le \lfloor \sqrt[3]{n}\rfloor}\pi\!\left(\left\lfloor\frac{n}{p^3}\right\rfloor\right).$$

This is exactly the quantity the implementation accumulates before correcting the bad diagonal case.

Step 4: Correct the overcount and add the \(p^7\) family

The raw sum \(B_{\mathrm{raw}}(n)\) incorrectly includes the choice \(q=p\). That choice gives

$$p^3q=p^4,$$

and \(p^4\) has only

$$d(p^4)=5$$

divisors, not \(8\). This invalid term appears exactly when

$$p^4\le n,$$

so the number of bad terms is

$$\pi\!\left(\left\lfloor n^{1/4}\right\rfloor\right).$$

Therefore the correct count for the \(p^3q\) family is

$$B(n)=B_{\mathrm{raw}}(n)-\pi\!\left(\left\lfloor n^{1/4}\right\rfloor\right).$$

The remaining family is \(p^7\le n\), whose contribution is

$$C(n)=\pi\!\left(\left\lfloor n^{1/7}\right\rfloor\right).$$

Putting the three disjoint cases together, we obtain

$$\boxed{f(n)=A(n)+B(n)+C(n).}$$

Worked Example: \(n=100\)

For the \(pqr\) family, only \(p=2\) contributes. With \((p,q)=(2,3)\), we get

$$\pi\!\left(\left\lfloor\frac{100}{6}\right\rfloor\right)-\pi(3)=\pi(16)-\pi(3)=6-2=4,$$

corresponding to \(30,42,66,78\). With \((p,q)=(2,5)\), we get

$$\pi\!\left(\left\lfloor\frac{100}{10}\right\rfloor\right)-\pi(5)=\pi(10)-\pi(5)=4-3=1,$$

which gives \(70\). Hence

$$A(100)=5.$$

For the \(p^3q\) family, the raw count is

$$B_{\mathrm{raw}}(100)=\pi(12)+\pi(3)=5+2=7.$$

The invalid choices \(q=p\) are counted by

$$\pi\!\left(\left\lfloor 100^{1/4}\right\rfloor\right)=\pi(3)=2,$$

so

$$B(100)=7-2=5.$$

These five numbers are \(24,40,54,56,88\). Finally, \(100\lt 2^7\), so

$$C(100)=0.$$

Therefore

$$f(100)=5+5+0=10,$$

which matches the checkpoint used by the implementation.

How the Code Works

The C++, Python, and Java implementations first compute exact integer square, cube, and seventh roots, with small corrective adjustments so every boundary case is handled safely. They then sieve all primes up to \(\lfloor\sqrt{n}\rfloor\). That range is enough for explicit iteration, because every prime that appears in the outer loops lies at most in that interval.

Next, the implementation builds a combinatorial prime-counting structure on two domains: direct arguments up to \(\sqrt{n}\), and quotient arguments of the form \(\lfloor n/i\rfloor\). After a prime-by-prime update pass, this structure can answer every needed value of \(\pi(x)\) in constant time.

With those prime counts available, the implementation evaluates the same three contributions as the mathematics: prime pairs for the \(pqr\) family, the raw sum for \(p^3q\), the subtraction of the forbidden \(p^4\) terms, and finally the \(p^7\) contribution. The sum of those parts is the final answer.

Complexity Analysis

Let \(v=\lfloor\sqrt{n}\rfloor\). The sieve and the prime-counting preprocessing use \(O(v)\) memory and about \(O(v\log\log v)\) arithmetic on the explicitly sieved range. The \(p^3q\) stage iterates over primes \(p\le n^{1/3}\), so it requires only \(O(\pi(n^{1/3}))\) prime-count queries.

The dominant work is the \(pqr\) stage, which visits each prime pair \((p,q)\) satisfying \(p\lt q\) and \(pq^2\le n\). Its loop count is

$$O\!\left(\sum_{p\le n^{1/3}}\pi\!\left(\sqrt{\frac{n}{p}}\right)\right),$$

and each visit performs only \(O(1)\) table lookups. In practice this is vastly smaller than scanning all integers up to \(n\), and it is easily fast enough for \(n=10^{12}\).

Footnotes and References

  1. Problem page: https://projecteuler.net/problem=501
  2. Divisor function: Wikipedia — Divisor function
  3. Prime-counting function: Wikipedia — Prime-counting function
  4. Sieve of Eratosthenes: Wikipedia — Sieve of Eratosthenes
  5. Fundamental theorem of arithmetic: Wikipedia — Fundamental theorem of arithmetic

Problem 501 source code

C++

#include <cmath>
#include <cstdint>
#include <iostream>
#include <string>
#include <vector>

namespace {

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

struct Options {
    u64 n = 1'000'000'000'000ULL;
    bool run_checkpoints = true;
};

bool parse_u64_after_prefix(const std::string& arg, const std::string& prefix, u64& value) {
    if (arg.rfind(prefix, 0U) != 0U) return false;
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) return false;
    u64 parsed = 0ULL;
    for (char ch : tail) {
        if (ch < '0' || ch > '9') return false;
        parsed = parsed * 10ULL + static_cast<u64>(ch - '0');
    }
    value = parsed;
    return true;
}

bool parse_arguments(int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        if (parse_u64_after_prefix(arg, "--n=", options.n)) continue;
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.n >= 1ULL;
}

u64 isqrt_u64(u64 n) {
    u64 x = static_cast<u64>(std::sqrt(static_cast<long double>(n)));
    while ((x + 1ULL) <= n / (x + 1ULL)) ++x;
    while (x > 0ULL && x > n / x) --x;
    return x;
}

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

bool pow_leq(u64 base, int exp, u64 limit) {
    u128 v = 1;
    for (int i = 0; i < exp; ++i) {
        v *= base;
        if (v > static_cast<u128>(limit)) return false;
    }
    return true;
}

u64 iroot7_u64(u64 n) {
    long double x = static_cast<long double>(n);
    u64 r = static_cast<u64>(std::pow(x, 1.0L / 7.0L));
    while (pow_leq(r + 1ULL, 7, n)) ++r;
    while (r > 0ULL && !pow_leq(r, 7, n)) --r;
    return r;
}

std::vector<int> prime_sieve(u64 n) {
    if (n < 2ULL) return {};
    const std::size_t m = static_cast<std::size_t>(n + 1ULL);
    std::vector<unsigned char> is_prime(m, 1);
    is_prime[0] = 0;
    is_prime[1] = 0;
    const u64 r = isqrt_u64(n);
    for (u64 p = 2; p <= r; ++p) {
        if (!is_prime[static_cast<std::size_t>(p)]) continue;
        for (u64 q = p * p; q <= n; q += p) is_prime[static_cast<std::size_t>(q)] = 0;
    }
    std::vector<int> primes;
    primes.reserve(m / 10);
    for (u64 i = 2; i <= n; ++i) {
        if (is_prime[static_cast<std::size_t>(i)]) primes.push_back(static_cast<int>(i));
    }
    return primes;
}

struct PiTables {
    u64 n = 0;
    u64 v = 0;
    std::vector<u64> smalls;
    std::vector<u64> larges;
};

PiTables tabulate_pis(u64 n, const std::vector<int>& primes) {
    PiTables t;
    t.n = n;
    t.v = isqrt_u64(n);
    t.smalls.assign(static_cast<std::size_t>(t.v + 1ULL), 0ULL);
    t.larges.assign(static_cast<std::size_t>(t.v + 1ULL), 0ULL);

    for (u64 i = 1; i <= t.v; ++i) {
        t.smalls[static_cast<std::size_t>(i)] = i - 1ULL;
        t.larges[static_cast<std::size_t>(i)] = n / i - 1ULL;
    }

    for (int p_int : primes) {
        const u64 p = static_cast<u64>(p_int);
        if (p > t.v) break;
        const u64 p_cnt = t.smalls[static_cast<std::size_t>(p - 1ULL)];
        const u64 q = p * p;
        if (q > n) break;
        const u64 end = std::min<u64>(t.v, n / q);

        for (u64 i = 1; i <= end; ++i) {
            const u64 d = i * p;
            if (d <= t.v) {
                t.larges[static_cast<std::size_t>(i)] -=
                    t.larges[static_cast<std::size_t>(d)] - p_cnt;
            } else {
                t.larges[static_cast<std::size_t>(i)] -=
                    t.smalls[static_cast<std::size_t>(n / d)] - p_cnt;
            }
        }

        for (u64 i = t.v; i >= q; --i) {
            t.smalls[static_cast<std::size_t>(i)] -=
                t.smalls[static_cast<std::size_t>(i / p)] - p_cnt;
            if (i == q) break;
        }
    }

    return t;
}

u64 solve(u64 n) {
    const u64 v = isqrt_u64(n);
    const std::vector<int> primes = prime_sieve(v);
    const PiTables pi = tabulate_pis(n, primes);

    u64 ans = 0ULL;
    const u64 cbrt_n = icbrt_u64(n);
    const u64 pi_cbrt_n = pi.smalls[static_cast<std::size_t>(cbrt_n)];

    for (u64 pi_idx = 0; pi_idx < pi_cbrt_n; ++pi_idx) {
        const u64 p = static_cast<u64>(primes[static_cast<std::size_t>(pi_idx)]);
        if (p * p * p >= n) break;
        const u64 m = n / p;
        const u64 q_lim = pi.smalls[static_cast<std::size_t>(isqrt_u64(m))];

        for (u64 pj_idx = pi_idx + 1ULL; pj_idx < q_lim; ++pj_idx) {
            const u64 q = static_cast<u64>(primes[static_cast<std::size_t>(pj_idx)]);
            const u64 r = m / q;
            const u64 pi_r = (r <= pi.v)
                                 ? pi.smalls[static_cast<std::size_t>(r)]
                                 : pi.larges[static_cast<std::size_t>(p * q)];
            ans += pi_r - pi.smalls[static_cast<std::size_t>(q)];
        }
    }

    for (u64 pi_idx = 0; pi_idx < pi_cbrt_n; ++pi_idx) {
        const u64 p = static_cast<u64>(primes[static_cast<std::size_t>(pi_idx)]);
        const u64 p3 = p * p * p;
        const u64 r = n / p3;
        if (r <= 1ULL) break;
        ans += (r <= pi.v)
                   ? pi.smalls[static_cast<std::size_t>(r)]
                   : pi.larges[static_cast<std::size_t>(p3)];
    }

    const u64 root4 = isqrt_u64(isqrt_u64(n));
    ans -= pi.smalls[static_cast<std::size_t>(root4)];
    ans += pi.smalls[static_cast<std::size_t>(iroot7_u64(n))];

    return ans;
}

u64 brute_count(u64 n) {
    u64 count = 0ULL;
    for (u64 x = 1ULL; x <= n; ++x) {
        u64 y = x;
        int d = 1;
        for (u64 p = 2ULL; p * p <= y; ++p) {
            if (y % p != 0ULL) continue;
            int e = 0;
            while (y % p == 0ULL) {
                y /= p;
                ++e;
            }
            d *= (e + 1);
        }
        if (y > 1ULL) d *= 2;
        if (d == 8) ++count;
    }
    return count;
}

bool run_checkpoints() {
    if (solve(100ULL) != 10ULL) {
        std::cerr << "Checkpoint failed: f(100)=10" << '\n';
        return false;
    }
    if (solve(1000ULL) != 180ULL) {
        std::cerr << "Checkpoint failed: f(1000)=180" << '\n';
        return false;
    }
    if (solve(1'000'000ULL) != 224'427ULL) {
        std::cerr << "Checkpoint failed: f(1e6)=224427" << '\n';
        return false;
    }
    if (solve(20'000ULL) != brute_count(20'000ULL)) {
        std::cerr << "Checkpoint failed: brute-force cross-check for n=20000" << '\n';
        return false;
    }
    return true;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) return 1;
    if (options.run_checkpoints && !run_checkpoints()) return 2;
    std::cout << solve(options.n) << '\n';
    return 0;
}

Python

import math

def solve():
    N = 10**12

    def isqrt(x): return math.isqrt(x)
    def icbrt(x):
        r = round(x ** (1/3))
        while (r+1)**3 <= x: r += 1
        while r**3 > x: r -= 1
        return r

    def iroot7(n):
        r = round(n ** (1/7))
        while (r+1)**7 <= n: r += 1
        while r**7 > n: r -= 1
        return r

    v = isqrt(N)
    # Sieve primes up to v
    sieve = bytearray([1]) * (v + 1)
    sieve[0] = sieve[1] = 0
    for p in range(2, isqrt(v) + 1):
        if sieve[p]:
            for q in range(p*p, v+1, p):
                sieve[q] = 0
    primes = [i for i in range(2, v+1) if sieve[i]]

    # Lucy pi tables
    smalls = list(range(v + 1))  # smalls[i] = i - 1 for i >= 1
    for i in range(v + 1):
        smalls[i] = i - 1
    larges = [0] * (v + 1)
    for i in range(1, v + 1):
        larges[i] = N // i - 1

    for p in primes:
        if p * p > N: break
        p_cnt = smalls[p - 1]
        q = p * p
        end = min(v, N // q)
        for i in range(1, end + 1):
            d = i * p
            if d <= v:
                larges[i] -= larges[d] - p_cnt
            else:
                larges[i] -= smalls[N // d] - p_cnt
        i = v
        while i >= q:
            smalls[i] -= smalls[i // p] - p_cnt
            if i == q: break
            i -= 1

    # Count numbers with exactly 8 divisors
    ans = 0
    cbrt_n = icbrt(N)
    pi_cbrt = smalls[cbrt_n]

    # Case 1: p*q*r (three distinct primes)
    for pi_idx in range(pi_cbrt):
        p = primes[pi_idx]
        if p**3 >= N: break
        m = N // p
        q_lim = smalls[isqrt(m)]
        for pj_idx in range(pi_idx + 1, q_lim):
            q = primes[pj_idx]
            r = m // q
            pi_r = smalls[r] if r <= v else larges[p * q]
            ans += pi_r - smalls[q]

    # Case 2: p^3 * q
    for pi_idx in range(pi_cbrt):
        p = primes[pi_idx]
        p3 = p**3
        r = N // p3
        if r <= 1: break
        ans += smalls[r] if r <= v else larges[p3]

    # Case 3: p^7 (subtract overcounting)
    root4 = isqrt(isqrt(N))
    ans -= smalls[root4]

    # Case 4: p^7 (add back)
    ans += smalls[iroot7(N)]

    return str(ans)

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

Java

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

public class Euler501 {

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

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

    private static long iroot7(long n) {
        if (n == 0)
            return 0;
        long x = (long) Math.pow(n, 1.0 / 7.0);
        while (Math.pow(x + 1, 7) <= n)
            x++;
        while (x > 0 && Math.pow(x, 7) > n)
            x--;
        return x;
    }

    private static List<Integer> primeSieve(long n) {
        if (n < 2)
            return new ArrayList<>();
        int m = (int) (n + 1);
        boolean[] isPrime = new boolean[m];
        for (int i = 2; i < m; i++)
            isPrime[i] = true;
        long r = isqrt(n);
        for (int p = 2; p <= r; p++) {
            if (isPrime[p]) {
                for (int q = p * p; q < m; q += p) {
                    isPrime[q] = false;
                }
            }
        }
        List<Integer> primes = new ArrayList<>();
        for (int i = 2; i < m; i++) {
            if (isPrime[i])
                primes.add(i);
        }
        return primes;
    }

    static class PiTables {
        long n, v;
        long[] smalls, larges;

        PiTables(long n, long v, long[] smalls, long[] larges) {
            this.n = n;
            this.v = v;
            this.smalls = smalls;
            this.larges = larges;
        }
    }

    private static PiTables tabulatePis(long n, List<Integer> primes) {
        long v = isqrt(n);
        long[] smalls = new long[(int) (v + 1)];
        long[] larges = new long[(int) (v + 1)];

        for (int i = 1; i <= v; i++) {
            smalls[i] = i - 1;
            larges[i] = n / i - 1;
        }

        for (int pInt : primes) {
            long p = pInt;
            if (p > v)
                break;
            long pCnt = smalls[(int) p - 1];
            long q = p * p;
            if (q > n)
                break;
            long end = Math.min(v, n / q);

            for (int i = 1; i <= end; i++) {
                long d = i * p;
                if (d <= v) {
                    larges[i] -= larges[(int) d] - pCnt;
                } else {
                    larges[i] -= smalls[(int) (n / d)] - pCnt;
                }
            }

            for (int i = (int) v; i >= q; i--) {
                smalls[i] -= smalls[(int) (i / p)] - pCnt;
            }
        }
        return new PiTables(n, v, smalls, larges);
    }

    public static void main(String[] args) {
        long n = 1000000000000L;
        long v = isqrt(n);
        List<Integer> primes = primeSieve(v);
        PiTables pi = tabulatePis(n, primes);

        long ans = 0;
        long cbrtN = icbrt(n);
        long piCbrtN = pi.smalls[(int) cbrtN];

        for (int piIdx = 0; piIdx < piCbrtN; piIdx++) {
            long p = primes.get(piIdx);
            if (p * p * p >= n)
                break;
            long m = n / p;
            long qLim = pi.smalls[(int) isqrt(m)];

            for (int pjIdx = piIdx + 1; pjIdx < qLim; pjIdx++) {
                long q = primes.get(pjIdx);
                long r = m / q;
                long piR;
                if (r <= pi.v) {
                    piR = pi.smalls[(int) r];
                } else {
                    piR = pi.larges[(int) (p * q)];
                }
                ans += piR - pi.smalls[(int) q];
            }
        }

        for (int piIdx = 0; piIdx < piCbrtN; piIdx++) {
            long p = primes.get(piIdx);
            long p3 = p * p * p;
            long r = n / p3;
            if (r <= 1)
                break;
            if (r <= pi.v) {
                ans += pi.smalls[(int) r];
            } else {
                ans += pi.larges[(int) p3];
            }
        }

        long root4 = isqrt(isqrt(n));
        ans -= pi.smalls[(int) root4];
        ans += pi.smalls[(int) iroot7(n)];

        System.out.println(ans);
    }
}