Project Euler

Cube-full Divisors

Let s(n) be the number of cube-full divisors of n, where a positive integer is cube-full if every prime factor appears with exponent at least 3. Define S(N) = sum_(n=1)^N s(n). Find S(10^18).

Source sync Apr 19, 2026

Problem #0694

Level Level 14

Solved By 1,227

Languages C++, Python

Answer 1339784153569958487

Length 242 words

number_theoryrecursionlinear_algebra

A positive integer $n$ is considered cube-full, if for every prime $p$ that divides $n$, so does $p^3$. Note that $1$ is considered cube-full.

Let $s(n)$ be the function that counts the number of cube-full divisors of $n$. For example, $1$, $8$ and $16$ are the three cube-full divisors of $16$. Therefore, $s(16)=3$.

Let $S(n)$ represent the summatory function of $s(n)$, that is $S(n)=\displaystyle \sum _{i=1}^n s(i)$.

You are given $S(16) = 19$, $S(100) = 126$ and $S(10000) = 13344$.

Find $S(10^{18})$.

Problem 694: Cube-full Divisors

Mathematical Foundation

Lemma 1. For every positive integer $N$ ,

S(N) = \sum_{\substack{c \le N \\ c \text{ cube-full}}} \left\lfloor \frac{N}{c} \right\rfloor.

Proof. Expand $s(n)$ as

s(n) = \sum_{\substack{c \mid n \\ c \text{ cube-full}}} 1.

Summing over $n \le N$ and exchanging the order of summation gives

S(N) = \sum_{n \le N}\sum_{\substack{c \mid n \\ c \text{ cube-full}}} 1 = \sum_{\substack{c \le N \\ c \text{ cube-full}}}\#\{n \le N : c \mid n\} = \sum_{\substack{c \le N \\ c \text{ cube-full}}} \left\lfloor \frac{N}{c} \right\rfloor.

$\square$

Lemma 2. Every cube-full number has a unique factorisation

c = \prod_{i=1}^{r} p_i^{e_i}, \qquad e_i \ge 3,

with distinct primes $p_1 < \cdots < p_r$ .

Proof. This is just the uniqueness of prime factorisation, together with the defining condition $e_i \ge 3$ for every prime divisor. $\square$

Theorem. Recursively enumerating primes in increasing order and, for each chosen prime, allowing exponents $3,4,5,\dots$ , visits every cube-full number at most $N$ exactly once.

Proof. By Lemma 2 every cube-full number corresponds to a unique strictly increasing list of primes and a unique exponent attached to each prime. The recursion chooses the next prime only from larger primes, so the same set of prime factors cannot be generated in two different orders. Hence each cube-full number appears exactly once. $\square$

Editorial

This is exactly the method used by the existing Python and C++ implementations. We recursively enumerate cube-full numbers. We then iterate over the next prime $p$ , try $p^3, p^4, p^5, \dots$ as long as the product stays $\le N$ ,. Finally, recurse only on larger primes.

Pseudocode

Sieve all primes up to $N^{1/3} = 10^6$
Recursively enumerate cube-full numbers:
start with `cur = 1`,
for the next prime $p$, try $p^3, p^4, p^5, \dots$ as long as the product stays $\le N$,
recurse only on larger primes
For each generated cube-full number $c$, add $\lfloor N/c \rfloor$

Complexity Analysis

The running time is proportional to the number of cube-full integers up to $N$ , which is sparse compared with all integers up to $N$ . The sieve phase costs $O(N^{1/3}\log\log N)$ with $N=10^{18}$ .

Time: dominated by the recursive enumeration.
Space: $O(\pi(N^{1/3}))$ for the prime list plus recursion stack.

Answer

\boxed{1339784153569958487}

C++ project_euler/problem_694/solution.cpp

#include <bits/stdc++.h>
using namespace std;

/*
 * Project Euler Problem 694: Cube-full Divisors
 *
 * s(n) = number of cube-full divisors of n (where p|d => p^3|d for all primes p).
 * S(n) = sum_{i=1}^{n} s(i).
 * Find S(10^18).
 *
 * Key insight: S(n) = sum over all cube-full numbers d of floor(n/d).
 * A cube-full number has the form: product of p_i^{a_i} where all a_i >= 3.
 * We can write any cube-full number as a^3 * b where b is square-free and
 * gcd(a, b) ... no, that's not right.
 *
 * Actually a cube-full number c can be written as c = a^3 * b^4 * d^5 ...
 * More precisely, enumerate cube-full numbers by their prime factorization.
 *
 * Approach: enumerate all cube-full numbers up to N = 10^18.
 * A cube-full number c = p1^a1 * p2^a2 * ... where all ai >= 3.
 * We can write c = m^3 * k where k is also cube-full, but this overcounts.
 *
 * Better: recursively enumerate cube-full numbers.
 * For each prime p, try exponents 3, 4, 5, 6, ..., and recurse.
 * Since p^3 <= N, we need primes up to N^(1/3) = 10^6.
 *
 * S(N) = sum over all cube-full c <= N of floor(N / c).
 */

typedef long long ll;
typedef unsigned long long ull;

const ll N = 1000000000000000000LL; // 10^18

// Sieve primes up to 10^6
vector<int> primes;
vector<bool> is_prime;

void sieve(int limit) {
    is_prime.assign(limit + 1, true);
    is_prime[0] = is_prime[1] = false;
    for (int i = 2; i <= limit; i++) {
        if (is_prime[i]) {
            primes.push_back(i);
            for (long long j = (long long)i * i; j <= limit; j += i)
                is_prime[j] = false;
        }
    }
}

ll answer = 0;

// Recursively enumerate cube-full numbers.
// cur = current cube-full number built so far
// pidx = index of next prime to consider (to avoid permutations)
// We try each prime p = primes[pidx] with exponent e >= 3,
// as long as cur * p^3 <= N.
void enumerate(ll cur, int pidx) {
    // cur is a cube-full number; add floor(N / cur) to answer.
    answer += N / cur;

    // Try next primes
    for (int i = pidx; i < (int)primes.size(); i++) {
        ll p = primes[i];
        ll p3 = p * p * p;
        if (cur > N / p3) break; // cur * p^3 > N

        ll val = cur * p3;
        for (int e = 3; val <= N; e++) {
            enumerate(val, i + 1);
            if (val > N / p) break;
            val *= p;
        }
    }
}

int main() {
    // N^(1/3) = 10^6
    sieve(1000001);

    enumerate(1, 0);

    cout << answer << endl;
    return 0;
}

Python project_euler/problem_694/solution.py

"""
Project Euler Problem 694: Cube-full Divisors

S(n) = sum_{i=1}^{n} s(i) where s(i) counts cube-full divisors of i.
S(n) = sum over all cube-full c <= n of floor(n/c).

Enumerate cube-full numbers recursively using primes up to n^(1/3).
"""

def solve():
    N = 10**18
    limit = 10**6 + 1

    # Sieve of Eratosthenes
    is_prime = [True] * (limit + 1)
    is_prime[0] = is_prime[1] = False
    primes = []
    for i in range(2, limit + 1):
        if is_prime[i]:
            primes.append(i)
            for j in range(i * i, limit + 1, i):
                is_prime[j] = False

    answer = [0]

    def enumerate(cur, pidx):
        answer[0] += N // cur
        for i in range(pidx, len(primes)):
            p = primes[i]
            p3 = p * p * p
            if cur * p3 > N:
                break
            val = cur * p3
            while val <= N:
                enumerate(val, i + 1)
                if val > N // p:
                    break
                val *= p

    enumerate(1, 0)
    print(answer[0])

if __name__ == "__main__":
    solve()