Problem 616: Creative Numbers

Project Euler Problem 616 Solution

EulerSolve provides an optimized solution for Project Euler Problem 616, Creative Numbers, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Start with the one-element list $\{n\}$. Alice may replace two list elements $a,b>1$ by the single value $a^b$, and she may replace one element $c$ by two values $a,b>1$ whenever $c=a^b$. A number $n>1$ is called creative if, for every target $m>1$, there exists some sequence of such moves that leads to a list containing $m$. For Problem 616 we must compute $$S=\sum_{\substack{n\le 10^{12} \\ n\text{ creative}}} n.$$ The implementations do not explore the move graph directly. Instead they use a classification theorem for all creative numbers up to the required limit and then turn the problem into a finite enumeration. Mathematical Approach Let $X=10^{12}$. Define the two sets $$\mathrm{PP}(X)=\{a^e\le X : a\ge 2,\ e\ge 2\},$$ $$\mathrm{PQ}(X)=\{p^q\le X : p\text{ prime},\ q\text{ prime}\}.$$ The C++, Python, and Java implementations rely on the classification $$\mathcal{C}(X)=\mathrm{PP}(X)\setminus \mathrm{PQ}(X)\setminus \{16\},$$ where $\mathcal{C}(X)$ denotes the set of creative numbers not exceeding $X$. The remaining work is to understand why this description is natural and how to sum it efficiently. Step 1: Only perfect powers can even enter the discussion If $n$ is not a perfect power, then the starting list $\{n\}$ cannot be split at all, because there do not exist integers $a,b>1$ with $n=a^b$....

Detailed mathematical approach

Problem Summary

Start with the one-element list $\{n\}$. Alice may replace two list elements $a,b>1$ by the single value $a^b$, and she may replace one element $c$ by two values $a,b>1$ whenever $c=a^b$. A number $n>1$ is called creative if, for every target $m>1$, there exists some sequence of such moves that leads to a list containing $m$.

For Problem 616 we must compute

$$S=\sum_{\substack{n\le 10^{12} \\ n\text{ creative}}} n.$$

The implementations do not explore the move graph directly. Instead they use a classification theorem for all creative numbers up to the required limit and then turn the problem into a finite enumeration.

Mathematical Approach

Let $X=10^{12}$. Define the two sets

$$\mathrm{PP}(X)=\{a^e\le X : a\ge 2,\ e\ge 2\},$$

$$\mathrm{PQ}(X)=\{p^q\le X : p\text{ prime},\ q\text{ prime}\}.$$

The C++, Python, and Java implementations rely on the classification

$$\mathcal{C}(X)=\mathrm{PP}(X)\setminus \mathrm{PQ}(X)\setminus \{16\},$$

where $\mathcal{C}(X)$ denotes the set of creative numbers not exceeding $X$. The remaining work is to understand why this description is natural and how to sum it efficiently.

Step 1: Only perfect powers can even enter the discussion

If $n$ is not a perfect power, then the starting list $\{n\}$ cannot be split at all, because there do not exist integers $a,b>1$ with $n=a^b$. Therefore the process is stuck immediately, and the only number ever present is $n$ itself.

But a creative number must be able to reach every target $m>1$, so a non-perfect power is impossible. Hence every creative $n\le X$ must lie in $\mathrm{PP}(X)$.

Step 2: Prime-base, prime-exponent powers are too rigid

Now consider

$$n=p^q,$$

with $p$ prime and $q$ prime. Splitting $n$ yields the two-element list $\{p,q\}$. Neither of those values can be split any further, because primes are not nontrivial perfect powers.

From $\{p,q\}$, the only possible merges are $p^q$ and $q^p$. If $p=q$, even that swap does nothing. So the reachable state space stays extremely small and cannot contain arbitrary targets.

Therefore every member of $\mathrm{PQ}(X)$ must be removed from the creative set.

Step 3: Why $16$ is a separate exceptional value

The number $16$ is a perfect power, but it is not in $\mathrm{PQ}(X)$, because its standard forms are

$$16=2^4=4^2.$$

Nevertheless it is still not creative. Every nontrivial split of $16$ produces only powers of $2$, and splitting those outputs again still produces only powers of $2$.

This is stable under merging as well. If all current elements are powers of $2$, say $2^r$ and $2^s$, then

$$\left(2^r\right)^{2^s}=2^{r2^s},$$

which is again a power of $2$. So starting from $16$, the process never leaves the world of powers of $2$. In particular, a target such as $3$ can never appear. That is why the implementations subtract $16$ separately.

Step 4: The positive direction becomes a classification theorem

The difficult direction is the converse: apart from the excluded family $\mathrm{PQ}(X)$ and the special case $16$, every remaining perfect power up to $10^{12}$ is creative. The implementations take exactly that theorem as their mathematical foundation.

Once this classification is accepted, the answer is no longer about searching for transformations. It is just the sum of all perfect powers up to $X$, minus the sum of the prime-prime powers, minus the single exceptional value $16$.

Step 5: Convert the theorem into a summation formula

Hence

$$S(X)=\sum_{x\in \mathrm{PP}(X)} x-\sum_{x\in \mathrm{PQ}(X)} x-16,$$

with $X=10^{12}$.

Two implementation details matter here:

First, $\mathrm{PP}(X)$ must be deduplicated, because the same integer can appear from several exponent choices. For example,

$$64=2^6=4^3=8^2.$$

Second, $\mathrm{PQ}(X)$ does not need deduplication. If $p^q=r^s$ with $p,r$ prime and $q,s$ prime, unique prime factorization forces $p=r$ and $q=s$.

Worked Example: apply the rule to $X=100$

The perfect powers up to $100$ are

$$\mathrm{PP}(100)=\{4,8,9,16,25,27,32,36,49,64,81,100\}.$$

The prime-prime powers among them are

$$\mathrm{PQ}(100)=\{4,8,9,25,27,32,49\}.$$

Removing those values and then removing $16$ leaves

$$\mathcal{C}(100)=\{36,64,81,100\}.$$

Therefore

$$S(100)=36+64+81+100=281.$$

This small example mirrors the full computation exactly: generate perfect powers, subtract the prime-prime family, and subtract the isolated exception $16$.

How the Code Works

The implementations first determine the largest relevant exponent $E$ with $2^E\le X$. For $X=10^{12}$, this gives $E=39$, because $2^{39}\le 10^{12} < 2^{40}$.

For each exponent $e=2,3,\dots,E$, the implementation computes the largest base $a$ satisfying $a^e\le X$ by integer binary search. It then enumerates every value $a^e$ in that range and stores all of them in a list or set. After sorting and deduplication, their total gives the sum over $\mathrm{PP}(X)$.

Next, the implementation builds the subtraction term $\mathrm{PQ}(X)$. Since $p^q\le X$ and $q\ge 2$, every prime base satisfies $p\le \sqrt{X}=10^6$, so a sieve up to $10^6$ is enough. The only relevant exponents are the primes between $2$ and $39$. Every valid value $p^q\le X$ is added to the prime-prime sum.

Finally, the answer is computed as

$$\text{sum of perfect powers}-\text{sum of prime-prime powers}-16.$$

The implementations also contain short sanity checks: several included values such as $36$, $81$, $100$, and $216$ are shown to reach $64$, while $16$ fails a negative check because all reachable values remain powers of $2$.

Complexity Analysis

Let

$$M=\sum_{e=2}^{E}\left\lfloor X^{1/e}\right\rfloor.$$

This is the number of raw perfect-power candidates generated before deduplication. The square case dominates, so $M=O(X^{1/2})$. Sorting and deduplicating those candidates costs $O(M\log M)$ time and $O(M)$ memory.

The prime sieve up to $10^6$ costs $O(10^6\log\log 10^6)$ time and $O(10^6)$ memory. The list of prime exponents is tiny because $E=39$. Overall, the method is easily fast enough for $X=10^{12}$ and avoids any attempt to search the full transformation graph.

Footnotes and References

Problem page: https://projecteuler.net/problem=616
Perfect power: Wikipedia - Perfect power
Prime power: Wikipedia - Prime power
Exponentiation: Wikipedia - Exponentiation
Sieve of Eratosthenes: Wikipedia - Sieve of Eratosthenes

Problem 616 source code

C++

#include <algorithm>
#include <cassert>
#include <cstdint>
#include <iostream>
#include <limits>
#include <numeric>
#include <vector>

// Project Euler 616: creative n <= 1e12 are exactly the perfect powers except p^q (p,q primes) and 16.

using i64 = long long;
using u64 = std::uint64_t;

static bool is_prime_int(int x) {
    if (x < 2) return false;
    if (x % 2 == 0) return x == 2;
    for (int d = 3; (i64)d * d <= x; d += 2)
        if (x % d == 0) return false;
    return true;
}

static std::vector<int> sieve_primes(int n) {
    std::vector<bool> is_prime((std::size_t)n + 1, true);
    if (n >= 0) is_prime[0] = false;
    if (n >= 1) is_prime[1] = false;
    for (int p = 2; (i64)p * p <= n; ++p) {
        if (!is_prime[p]) continue;
        for (int j = p * p; j <= n; j += p) is_prime[j] = false;
    }
    std::vector<int> primes;
    for (int i = 2; i <= n; ++i)
        if (is_prime[i]) primes.push_back(i);
    return primes;
}

static bool pow_leq(i64 a, int e, i64 limit) {
    __int128 r = 1;
    for (int i = 0; i < e; ++i) {
        r *= a;
        if (r > limit) return false;
    }
    return true;
}

static i64 pow_exact(i64 a, int e) {
    __int128 r = 1;
    for (int i = 0; i < e; ++i) r *= a;
    assert(r <= std::numeric_limits<i64>::max());
    return (i64)r;
}

static i64 floor_root(i64 n, int e) {
    i64 lo = 1, hi = 1000000 + 1;  // since n <= 1e12 and e >= 2
    while (lo + 1 < hi) {
        i64 mid = lo + (hi - lo) / 2;
        if (pow_leq(mid, e, n))
            lo = mid;
        else
            hi = mid;
    }
    return lo;
}

static u64 pow_u64(u64 a, u64 e) {
    __int128 r = 1;
    for (u64 i = 0; i < e; ++i) {
        r *= a;
        assert(r <= (__int128)std::numeric_limits<u64>::max());
    }
    return (u64)r;
}

#ifndef NDEBUG
static void ms_erase_one(std::vector<u64> &L, u64 x) {
    auto it = std::find(L.begin(), L.end(), x);
    assert(it != L.end());
    L.erase(it);
}

static void op_split(std::vector<u64> &L, u64 c, u64 a, u64 b) {
    assert(a > 1 && b > 1);
    assert(pow_u64(a, b) == c);
    ms_erase_one(L, c);
    L.push_back(a);
    L.push_back(b);
}

static void op_merge(std::vector<u64> &L, u64 a, u64 b) {
    assert(a > 1 && b > 1);
    ms_erase_one(L, a);
    ms_erase_one(L, b);
    L.push_back(pow_u64(a, b));
}

static bool contains(const std::vector<u64> &L, u64 x) {
    return std::find(L.begin(), L.end(), x) != L.end();
}

static void validate() {
    // A few short constructive checks that "included" values can reach 64.
    {
        std::vector<u64> L{36};
        op_split(L, 36, 6, 2);
        op_merge(L, 2, 6);
        assert(contains(L, 64));
    }
    {
        std::vector<u64> L{81};
        op_split(L, 81, 3, 4);
        op_merge(L, 4, 3);
        assert(contains(L, 64));
    }
    {
        std::vector<u64> L{100};
        op_split(L, 100, 10, 2);
        op_merge(L, 2, 10);  // 2^10 = 1024
        op_split(L, 1024, 32, 2);
        op_merge(L, 2, 32);  // 2^32
        op_split(L, 4294967296ULL, 16, 8);
        op_split(L, 8, 2, 3);
        op_split(L, 16, 4, 2);
        op_merge(L, 4, 3);
        assert(contains(L, 64));
    }
    {
        std::vector<u64> L{216};
        op_split(L, 216, 6, 3);
        op_merge(L, 3, 6);  // 3^6
        op_split(L, 729, 27, 2);
        op_merge(L, 2, 27);  // 2^27
        op_split(L, 134217728ULL, 8, 9);
        op_split(L, 8, 2, 3);
        op_split(L, 9, 3, 2);
        op_merge(L, 2, 2);
        op_merge(L, 4, 3);
        assert(contains(L, 64));
    }

    // A tiny sanity check that 16 cannot introduce 3 (all operations stay within powers of 2 here).
    {
        std::vector<u64> L{16};
        op_split(L, 16, 2, 4);
        op_split(L, 4, 2, 2);
        // Any merge of powers of 2 stays a power of 2, and remaining elements are powers of 2 too.
        assert(!contains(L, 3));
    }
}
#endif

int main() {
    constexpr i64 LIM = 1000000000000LL;

#ifndef NDEBUG
    validate();
#endif

    int max_e = 1;
    while (pow_leq(2, max_e + 1, LIM)) ++max_e;

    std::vector<i64> perfect;
    perfect.reserve(1100000);
    for (int e = 2; e <= max_e; ++e) {
        i64 max_a = floor_root(LIM, e);
        for (i64 a = 2; a <= max_a; ++a) perfect.push_back(pow_exact(a, e));
    }
    std::sort(perfect.begin(), perfect.end());
    perfect.erase(std::unique(perfect.begin(), perfect.end()), perfect.end());

    u64 sum_perfect = 0;
    for (i64 v : perfect) sum_perfect += (u64)v;

    const int P_MAX = 1000000;
    std::vector<int> primes = sieve_primes(P_MAX);
    std::vector<int> exp_primes;
    for (int e = 2; e <= max_e; ++e)
        if (is_prime_int(e)) exp_primes.push_back(e);

    u64 sum_primeprime = 0;
    for (int p : primes) {
        for (int e : exp_primes) {
            if (!pow_leq(p, e, LIM)) break;
            sum_primeprime += (u64)pow_exact(p, e);
        }
    }

    u64 ans = sum_perfect - sum_primeprime - 16ULL;
    std::cout << ans << "\n";
    return 0;
}

Python

import math

def is_prime_int(x):
    if x < 2: return False
    if x % 2 == 0: return x == 2
    for d in range(3, math.isqrt(x) + 1, 2):
        if x % d == 0: return False
    return True

def sieve_primes(n):
    is_prime = bytearray(b'\x01' * (n + 1))
    if n >= 0: is_prime[0] = 0
    if n >= 1: is_prime[1] = 0
    for p in range(2, math.isqrt(n) + 1):
        if is_prime[p]:
            is_prime[p * p : n + 1 : p] = bytes((n - p * p) // p + 1)
    return [i for i, b in enumerate(is_prime) if b]

def pow_leq(a, e, limit):
    if a == 1:
        return 1 <= limit
    return (a ** e) <= limit

def pow_exact(a, e):
    return a ** e

def floor_root(n, e):
    lo = 1
    hi = 1000000 + 1
    while lo + 1 < hi:
        mid = lo + (hi - lo) // 2
        if pow_leq(mid, e, n):
            lo = mid
        else:
            hi = mid
    return lo

def solve():
    LIM = 1000000000000
    
    max_e = 1
    while pow_leq(2, max_e + 1, LIM):
        max_e += 1
        
    perfect = []
    for e in range(2, max_e + 1):
        max_a = floor_root(LIM, e)
        for a in range(2, max_a + 1):
            perfect.append(pow_exact(a, e))
            
    perfect = list(set(perfect))
    sum_perfect = sum(perfect)
    
    P_MAX = 1000000
    primes = sieve_primes(P_MAX)
    exp_primes = [e for e in range(2, max_e + 1) if is_prime_int(e)]
    
    sum_primeprime = 0
    for p in primes:
        for e in exp_primes:
            if not pow_leq(p, e, LIM): break
            sum_primeprime += pow_exact(p, e)
            
    ans = sum_perfect - sum_primeprime - 16
    return str(ans)

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

Java

import java.util.ArrayList;
import java.util.HashSet;
import java.util.List;
import java.util.Set;

public class Euler616 {
    static boolean isPrimeInt(int x) {
        if (x < 2)
            return false;
        if (x % 2 == 0)
            return x == 2;
        for (int d = 3; (long) d * d <= x; d += 2) {
            if (x % d == 0)
                return false;
        }
        return true;
    }

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

    static boolean powLeq(long a, int e, long limit) {
        long r = 1;
        for (int i = 0; i < e; i++) {
            if (limit / a < r)
                return false;
            r *= a;
        }
        return true;
    }

    static long powExact(long a, int e) {
        long r = 1;
        for (int i = 0; i < e; i++) {
            r *= a;
        }
        return r;
    }

    static long floorRoot(long n, int e) {
        long lo = 1, hi = 1000000 + 1;
        while (lo + 1 < hi) {
            long mid = lo + (hi - lo) / 2;
            if (powLeq(mid, e, n)) {
                lo = mid;
            } else {
                hi = mid;
            }
        }
        return lo;
    }

    public static String solve() {
        long LIM = 1000000000000L;

        int maxE = 1;
        while (powLeq(2, maxE + 1, LIM))
            maxE++;

        Set<Long> perfect = new HashSet<>();
        for (int e = 2; e <= maxE; e++) {
            long maxA = floorRoot(LIM, e);
            for (long a = 2; a <= maxA; a++) {
                perfect.add(powExact(a, e));
            }
        }

        long sumPerfect = 0;
        for (long v : perfect) {
            sumPerfect += v;
        }

        int P_MAX = 1000000;
        List<Integer> primes = sievePrimes(P_MAX);
        List<Integer> expPrimes = new ArrayList<>();
        for (int e = 2; e <= maxE; e++) {
            if (isPrimeInt(e))
                expPrimes.add(e);
        }

        long sumPrimePrime = 0;
        for (int p : primes) {
            for (int e : expPrimes) {
                if (!powLeq(p, e, LIM))
                    break;
                sumPrimePrime += powExact(p, e);
            }
        }

        long ans = sumPerfect - sumPrimePrime - 16L;
        return Long.toString(ans);
    }

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