Problem 678: Fermat-like Equations

Project Euler Problem 678 Solution

EulerSolve provides an optimized solution for Project Euler Problem 678, Fermat-like Equations, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Let $F(n)$ be the number of tuples $(a,b,c,e,f)$ of positive integers satisfying $$a^e+b^e=c^f\le n,\qquad a<b,\qquad e\ge 2,\qquad f\ge 3.$$ The goal is to compute $F(10^{18})$. The challenge is that the same right-hand side value can be a perfect power in more than one way, while the left-hand side behaves very differently for $e=2$, $e=3$, and $e\ge 4$. Mathematical Approach The implementation separates the problem into three parts. First it catalogs all perfect powers $c^f\le n$ with $f\ge 3$. Then it counts representations $a^e+b^e=N$ using a special method for squares, another for cubes, and a direct scan for higher exponents. Step 1: Catalog Every Perfect Power on the Right-Hand Side For every base $c\ge 2$ and exponent $f\ge 3$, generate $$N=c^f\le n.$$ This produces a catalogue of perfect-power descriptions, not merely a set of distinct values. That distinction matters because one numerical value may support several exponents. If a number is both a cube and a fifth power, then a single identity $a^e+b^e=N$ contributes once for each admissible exponent $f$. For each numerical value $N$, the implementation stores the full set of exponents $f$ with $N=c^f$. The largest possible exponent is $$f_{\max}=\left\lfloor \log_2 n \right\rfloor,$$ because $2^f\le n$ is the smallest perfect power having exponent $f$....

Detailed mathematical approach

Problem Summary

Let $F(n)$ be the number of tuples $(a,b,c,e,f)$ of positive integers satisfying

$$a^e+b^e=c^f\le n,\qquad a<b,\qquad e\ge 2,\qquad f\ge 3.$$

The goal is to compute $F(10^{18})$. The challenge is that the same right-hand side value can be a perfect power in more than one way, while the left-hand side behaves very differently for $e=2$, $e=3$, and $e\ge 4$.

Mathematical Approach

The implementation separates the problem into three parts. First it catalogs all perfect powers $c^f\le n$ with $f\ge 3$. Then it counts representations $a^e+b^e=N$ using a special method for squares, another for cubes, and a direct scan for higher exponents.

Step 1: Catalog Every Perfect Power on the Right-Hand Side

For every base $c\ge 2$ and exponent $f\ge 3$, generate

$$N=c^f\le n.$$

This produces a catalogue of perfect-power descriptions, not merely a set of distinct values. That distinction matters because one numerical value may support several exponents. If a number is both a cube and a fifth power, then a single identity $a^e+b^e=N$ contributes once for each admissible exponent $f$.

For each numerical value $N$, the implementation stores the full set of exponents $f$ with $N=c^f$. The largest possible exponent is

$$f_{\max}=\left\lfloor \log_2 n \right\rfloor,$$

because $2^f\le n$ is the smallest perfect power having exponent $f$.

Step 2: Treat $e=2$ with the Sum-of-Two-Squares Theorem

Fix a perfect power $N=c^f$, and write

$$N=\prod p^{\alpha_p}.$$

The classical formula for the number of integer solutions of $x^2+y^2=N$ is

$$r_2(N)=4\prod_{p\equiv 1 \pmod 4}(\alpha_p+1),$$

provided every prime $p\equiv 3\pmod 4$ appears with even exponent. If one such prime has odd exponent, then $r_2(N)=0$.

This count includes signs and order, so it is not yet the desired number of positive pairs with $a<b$. Two special families must be removed first:

$$\text{axis solutions }(t,0)\text{ or }(0,t),\qquad \text{equal solutions }(t,t).$$

If $N$ is a square, there are four axis solutions. If $N=2t^2$, there are four equal solutions. Every remaining nonzero unequal representation occurs in exactly $8$ signed-and-ordered forms, so the required count is

$$\frac{r_2(N)-\text{axis}-\text{equal}}{8}.$$

Because $N=c^f$, it is enough to factor $c$ and multiply every prime exponent by $f$.

Step 3: Treat $e=3$ by Factoring $a^3+b^3$

For cubes we use the identity

$$a^3+b^3=(a+b)(a^2-ab+b^2).$$

Set

$$u=a+b,\qquad v=a^2-ab+b^2,\qquad N=uv.$$

Now enumerate divisors $u\mid N$, set $v=N/u$, and recover $a$ and $b$ from symmetric data. Since

$$u^2-v=3ab,$$

we obtain

$$ab=\frac{u^2-v}{3}.$$

Then $a$ and $b$ are roots of

$$t^2-ut+ab=0,$$

so the discriminant

$$\Delta=u^2-4ab$$

must be a perfect square. If $\Delta=s^2$ and $u\pm s$ are even, then

$$a=\frac{u-s}{2},\qquad b=\frac{u+s}{2}.$$

Only solutions with $a>0$ and $a<b$ are kept. The implementation also skips right-hand-side exponents divisible by $3$, because then $c^f$ is itself a cube and $a^3+b^3=z^3$ would contradict Fermat's Last Theorem for positive integers.

Step 4: Scan Directly for $e\ge 4$

For a fixed exponent $e\ge 4$, precompute

$$1^e,2^e,\dots,\left\lfloor n^{1/e}\right\rfloor^e.$$

For each $b$, admissible values of $a$ satisfy

$$1\le a\le \min\!\left(b-1,\left\lfloor (n-b^e)^{1/e}\right\rfloor\right).$$

Each candidate sum

$$s=a^e+b^e$$

is looked up in the perfect-power catalogue. If $s=c^f$ for several exponents $f$, each such exponent contributes except the multiples of $e$. Indeed, if $e\mid f$, then

$$s=c^f=\left(c^{f/e}\right)^e,$$

which would produce a nontrivial solution of $a^e+b^e=z^e$, impossible for $e\ge 3$.

There is also a parity pruning for even $e$. If both $a$ and $b$ are odd, then

$$a^e+b^e\equiv 1+1\equiv 2 \pmod 4,$$

but a perfect power $c^f$ with $f\ge 3$ is either odd or divisible by $8$, never $2\bmod 4$. So when $e$ is even and $b$ is odd, only even values of $a$ need to be tested.

Step 5: Why Duplicate Perfect-Power Values Still Matter

The count is over exponent choices as well as numerical values. A single value $N$ may correspond to several admissible right-hand-side exponents, so a solution $a^e+b^e=N$ may contribute more than once. This is why the implementation keeps both a list of all perfect-power descriptions and a compact record of which exponents belong to each value.

Worked Example: $n=1000$

The checkpoint built into the implementation is

$$F(1000)=7.$$

The seven counted tuples come from exactly four right-hand-side values:

$$125=2^2+11^2=5^2+10^2=5^3,$$

$$243=3^3+6^3=3^5,$$

$$625=7^2+24^2=15^2+20^2=5^4,$$

$$1000=10^2+30^2=18^2+26^2=10^3.$$

That is $2+1+2+2=7$. The square cases are detected by the two-squares theorem, the cubic case is reconstructed from divisors and a discriminant test, and no exponent $e\ge 4$ contributes below $1000$.

How the Code Works

The C++, Python, and Java implementations all follow the same counting strategy. They first enumerate all perfect powers $c^f\le 10^{18}$ with $f\ge 3$, preserving every individual $(c,f)$ description while also recording, for each value $N$, the set of exponents that produce it.

Next, the implementation builds a smallest-prime-factor sieve up to the largest base appearing in the catalogue. The $e=2$ branch uses that sieve to factor bases quickly and then applies the sum-of-two-squares theorem. The $e=3$ branch enumerates divisors of the perfect power, reconstructs possible pairs $(a,b)$, and caches the cubic count for each numerical value so repeated values are not recomputed.

For $e\ge 4$, the implementation precomputes the table of $e$-th powers, scans the admissible range $a<b$, applies the parity pruning for even exponents, and consults the perfect-power catalogue for every sum. The C++ implementation parallelizes large outer loops; the Python implementation keeps the same logic in a direct serial form; the Java implementation delegates to the same computational core.

Complexity Analysis

Let

$$M=\sum_{f=3}^{\lfloor \log_2 n\rfloor}\left\lfloor n^{1/f}\right\rfloor.$$

Generating the perfect-power catalogue costs $O(M)$, dominated by cubes and therefore behaving like $O(n^{1/3})$. Building the smallest-prime-factor sieve up to the largest base costs $O(n^{1/3}\log\log n)$ time and $O(n^{1/3})$ memory.

The $e=2$ branch is essentially linear in the number of catalogue entries after sieve preprocessing. The $e=3$ branch depends on divisor enumeration for the perfect powers that survive the $3\nmid f$ filter, but caching makes repeated values cheap. The dominant direct scan is

$$\sum_{e=4}^{\lfloor \log_2 n\rfloor} O\!\left(n^{2/e}\right),$$

with the largest contribution coming from $e=4$. In practice, the special handling of $e=2$ and $e=3$, parity pruning, perfect-power lookup, and parallel work splitting make the full $n=10^{18}$ computation feasible.

Footnotes and References

Project Euler problem page: https://projecteuler.net/problem=678
Perfect power: Wikipedia - Perfect power
Fermat's theorem on sums of two squares: Wikipedia - Fermat's theorem on sums of two squares
Sum of two cubes: Wikipedia - Sum of two cubes
Fermat's Last Theorem: Wikipedia - Fermat's Last Theorem

Problem 678 source code

C++

#include <pthread.h>
#include <bit>
#include <atomic>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <unordered_map>
#include <unistd.h>
#include <utility>
#include <vector>

namespace {

using u64 = std::uint64_t;
using u128 = __uint128_t;

struct Entry {
    u64 n;
    int c;
    int f;
};

int detect_thread_count(std::size_t work_items) {
    long cores = ::sysconf(_SC_NPROCESSORS_ONLN);
    int threads = (cores > 0) ? static_cast<int>(cores) : 4;
    if (threads < 1) threads = 1;
    if (threads > 4) threads = 4;
    if (work_items > 0 && static_cast<std::size_t>(threads) > work_items) {
        threads = static_cast<int>(work_items);
    }
    return threads;
}

struct EWorkerTask {
    const std::vector<u64>* pw = nullptr;
    const std::vector<int>* amax_by_b = nullptr;
    const std::unordered_map<u64, u64>* perfect_mask = nullptr;
    u64 mask_filter = 0ULL;
    bool even_e = false;
    std::atomic<int>* next_b = nullptr;
    int max_b = 0;
    u64 partial = 0ULL;
};

void* e_worker_entry(void* raw) {
    auto* task = static_cast<EWorkerTask*>(raw);
    u64 local = 0ULL;

    while (true) {
        const int b = task->next_b->fetch_add(1, std::memory_order_relaxed);
        if (b > task->max_b) break;

        const int amax = (*task->amax_by_b)[static_cast<std::size_t>(b)];
        if (amax <= 0) continue;

        const u64 pb = (*task->pw)[static_cast<std::size_t>(b)];
        if (task->even_e && ((b & 1) != 0)) {
            for (int a = 2; a <= amax; a += 2) {
                const u64 s = pb + (*task->pw)[static_cast<std::size_t>(a)];
                const auto it = task->perfect_mask->find(s);
                if (it != task->perfect_mask->end()) {
                    local += static_cast<u64>(__builtin_popcountll(static_cast<unsigned long long>(it->second & task->mask_filter)));
                }
            }
        } else {
            for (int a = 1; a <= amax; ++a) {
                const u64 s = pb + (*task->pw)[static_cast<std::size_t>(a)];
                const auto it = task->perfect_mask->find(s);
                if (it != task->perfect_mask->end()) {
                    local += static_cast<u64>(__builtin_popcountll(static_cast<unsigned long long>(it->second & task->mask_filter)));
                }
            }
        }
    }

    task->partial = local;
    return nullptr;
}

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

std::vector<Entry> generate_entries(u64 n, int& max_f, int& max_c) {
    max_f = 0;
    max_c = 0;
    while ((1ULL << max_f) <= n && max_f < 62) {
        ++max_f;
    }

    std::vector<Entry> out;
    for (int f = 3; f <= max_f; ++f) {
        for (int c = 2;; ++c) {
            const u64 v = pow_limit(static_cast<u64>(c), f, n);
            if (v > n) {
                break;
            }
            out.push_back({v, c, f});
            if (c > max_c) {
                max_c = c;
            }
        }
    }
    return out;
}

std::vector<int> build_spf(int n) {
    std::vector<int> spf(static_cast<std::size_t>(n + 1), 0);
    std::vector<int> primes;
    primes.reserve(static_cast<std::size_t>(n / 10));

    for (int i = 2; i <= n; ++i) {
        if (spf[static_cast<std::size_t>(i)] == 0) {
            spf[static_cast<std::size_t>(i)] = i;
            primes.push_back(i);
        }
        for (int p : primes) {
            const u64 v = static_cast<u64>(p) * static_cast<u64>(i);
            if (v > static_cast<u64>(n) || p > spf[static_cast<std::size_t>(i)]) {
                break;
            }
            spf[static_cast<std::size_t>(v)] = p;
        }
    }

    return spf;
}

std::vector<std::pair<int, int>> factorize_int(int x, const std::vector<int>& spf) {
    std::vector<std::pair<int, int>> fac;
    int v = x;
    while (v > 1) {
        const int p = spf[static_cast<std::size_t>(v)];
        int e = 0;
        do {
            v /= p;
            ++e;
        } while (v % p == 0);
        fac.push_back({p, e});
    }
    return fac;
}

u64 count_sum_two_squares_from_c_f(int c, int f, const std::vector<int>& spf) {
    const auto fac = factorize_int(c, spf);

    u64 prod = 1ULL;
    bool all_even_other = true;
    int e2 = 0;

    for (const auto& pe : fac) {
        const int p = pe.first;
        const int e = pe.second * f;

        if (p == 2) {
            e2 = e;
        } else if ((e & 1) != 0) {
            all_even_other = false;
        }

        if ((p & 3) == 3 && (e & 1) != 0) {
            return 0ULL;
        }
        if ((p & 3) == 1) {
            prod *= static_cast<u64>(e + 1);
        }
    }

    const u64 r2 = 4ULL * prod;
    const u64 zero_repr = (all_even_other && ((e2 & 1) == 0)) ? 4ULL : 0ULL;
    const u64 equal_repr = (all_even_other && e2 >= 1 && (((e2 - 1) & 1) == 0)) ? 4ULL : 0ULL;

    return (r2 - zero_repr - equal_repr) / 8ULL;
}

u64 isqrt_u128(u128 x) {
    long double ld = std::sqrt(static_cast<long double>(x));
    u64 r = static_cast<u64>(ld);

    while (static_cast<u128>(r + 1ULL) * static_cast<u128>(r + 1ULL) <= x) {
        ++r;
    }
    while (static_cast<u128>(r) * static_cast<u128>(r) > x) {
        --r;
    }
    return r;
}

std::vector<u64> divisors_from_factorization(const std::vector<std::pair<u64, int>>& fac) {
    std::vector<u64> divs{1ULL};
    for (const auto& pe : fac) {
        const u64 p = pe.first;
        const int e = pe.second;
        const std::size_t old = divs.size();
        u64 mul = 1ULL;
        for (int i = 1; i <= e; ++i) {
            mul *= p;
            for (std::size_t j = 0; j < old; ++j) {
                divs.push_back(divs[j] * mul);
            }
        }
    }
    return divs;
}

u64 count_sum_two_cubes_from_factorization(u64 n, const std::vector<std::pair<u64, int>>& fac) {
    const std::vector<u64> divs = divisors_from_factorization(fac);
    u64 count = 0ULL;

    for (u64 u : divs) {
        const u64 v = n / u;
        const u128 u2 = static_cast<u128>(u) * static_cast<u128>(u);

        if (u2 <= static_cast<u128>(v)) {
            continue;
        }

        const u128 t = u2 - static_cast<u128>(v);
        if (t % 3U != 0U) {
            continue;
        }
        const u128 p = t / 3U;

        if (u2 < 4U * p) {
            continue;
        }
        const u128 delta = u2 - 4U * p;

        const u64 s = isqrt_u128(delta);
        if (static_cast<u128>(s) * static_cast<u128>(s) != delta) {
            continue;
        }

        if (s >= u || ((u - s) & 1ULL) != 0ULL) {
            continue;
        }

        const u64 a = (u - s) / 2ULL;
        const u64 b = (u + s) / 2ULL;
        if (a == 0ULL || a >= b) {
            continue;
        }

        const u128 lhs = static_cast<u128>(a) * a * a + static_cast<u128>(b) * b * b;
        if (lhs == static_cast<u128>(n)) {
            ++count;
        }
    }

    return count;
}

u64 solve(u64 n) {
    int max_f = 0;
    int max_c = 0;
    const std::vector<Entry> entries = generate_entries(n, max_f, max_c);

    std::unordered_map<u64, u64> perfect_mask;
    perfect_mask.reserve(entries.size() * 2U);
    for (const auto& e : entries) {
        perfect_mask[e.n] |= (1ULL << e.f);
    }

    std::vector<u64> filter(static_cast<std::size_t>(max_f + 1), 0ULL);
    for (int e = 0; e <= max_f; ++e) {
        u64 m = 0ULL;
        for (int f = 3; f <= max_f; ++f) {
            if (e == 0 || (f % e) != 0) {
                m |= (1ULL << f);
            }
        }
        filter[static_cast<std::size_t>(e)] = m;
    }

    const std::vector<int> spf = build_spf(max_c);

    u64 ans = 0ULL;

    for (const auto& e : entries) {
        ans += count_sum_two_squares_from_c_f(e.c, e.f, spf);
    }

    std::unordered_map<u64, u64> cache3;
    cache3.reserve(1 << 16);

    for (const auto& e : entries) {
        if ((e.f % 3) == 0) {
            continue;
        }

        auto it = cache3.find(e.n);
        if (it != cache3.end()) {
            ans += it->second;
            continue;
        }

        const auto fac_c = factorize_int(e.c, spf);
        std::vector<std::pair<u64, int>> fac_n;
        fac_n.reserve(fac_c.size());
        for (const auto& pe : fac_c) {
            fac_n.push_back({static_cast<u64>(pe.first), pe.second * e.f});
        }

        const u64 val = count_sum_two_cubes_from_factorization(e.n, fac_n);
        cache3.emplace(e.n, val);
        ans += val;
    }

    for (int e = 4; e <= max_f; ++e) {
        std::vector<u64> pw;
        pw.push_back(0ULL);
        for (u64 b = 1;; ++b) {
            const u64 v = pow_limit(b, e, n);
            if (v > n) {
                break;
            }
            pw.push_back(v);
        }

        const int max_b = static_cast<int>(pw.size()) - 1;
        if (max_b < 2) {
            continue;
        }

        const bool even_e = ((e & 1) == 0);
        const u64 mask_filter = filter[static_cast<std::size_t>(e)];

        std::vector<int> amax_by_b(static_cast<std::size_t>(max_b + 1), 0);
        int amax_limit = max_b;
        for (int b = 2; b <= max_b; ++b) {
            const u64 pb = pw[static_cast<std::size_t>(b)];
            while (amax_limit > 0 && pw[static_cast<std::size_t>(amax_limit)] > n - pb) {
                --amax_limit;
            }
            int amax = amax_limit;
            if (amax >= b) amax = b - 1;
            amax_by_b[static_cast<std::size_t>(b)] = amax;
        }

        const int threads = detect_thread_count(static_cast<std::size_t>(max_b - 1));
        if (threads <= 1 || max_b < 3000) {
            for (int b = 2; b <= max_b; ++b) {
                const int amax = amax_by_b[static_cast<std::size_t>(b)];
                if (amax <= 0) continue;
                const u64 pb = pw[static_cast<std::size_t>(b)];
                if (even_e && ((b & 1) != 0)) {
                    for (int a = 2; a <= amax; a += 2) {
                        const u64 s = pb + pw[static_cast<std::size_t>(a)];
                        const auto it = perfect_mask.find(s);
                        if (it != perfect_mask.end()) {
                            ans += static_cast<u64>(__builtin_popcountll(static_cast<unsigned long long>(it->second & mask_filter)));
                        }
                    }
                } else {
                    for (int a = 1; a <= amax; ++a) {
                        const u64 s = pb + pw[static_cast<std::size_t>(a)];
                        const auto it = perfect_mask.find(s);
                        if (it != perfect_mask.end()) {
                            ans += static_cast<u64>(__builtin_popcountll(static_cast<unsigned long long>(it->second & mask_filter)));
                        }
                    }
                }
            }
        } else {
            std::vector<pthread_t> handles(static_cast<std::size_t>(threads));
            std::vector<EWorkerTask> tasks(static_cast<std::size_t>(threads));
            std::atomic<int> next_b{2};

            for (int t = 0; t < threads; ++t) {
                auto& task = tasks[static_cast<std::size_t>(t)];
                task.pw = &pw;
                task.amax_by_b = &amax_by_b;
                task.perfect_mask = &perfect_mask;
                task.mask_filter = mask_filter;
                task.even_e = even_e;
                task.next_b = &next_b;
                task.max_b = max_b;
                task.partial = 0ULL;
                pthread_create(&handles[static_cast<std::size_t>(t)], nullptr, e_worker_entry, &task);
            }
            for (int t = 0; t < threads; ++t) {
                pthread_join(handles[static_cast<std::size_t>(t)], nullptr);
                ans += tasks[static_cast<std::size_t>(t)].partial;
            }
        }
    }

    return ans;
}

}  // namespace

int main() {
    assert(solve(1'000ULL) == 7ULL);
    assert(solve(100'000ULL) == 53ULL);
    assert(solve(10'000'000ULL) == 287ULL);

    std::cout << solve(1'000'000'000'000'000'000ULL) << "\n";
    return 0;
}

Python

import math

def solve():
    N = 10**18

    def isqrt(n):
        r=int(math.isqrt(n))
        while (r+1)*(r+1)<=n: r+=1
        while r*r>n: r-=1
        return r

    def pow_lim(b,e):
        v=1
        for _ in range(e):
            if v>N//b: return N+1
            v*=b
        return v

    def spf_sieve(n):
        s=list(range(n+1)); s[0]=s[1]=0
        for i in range(2,int(n**0.5)+1):
            if s[i]==i:
                for j in range(i*i,n+1,i):
                    if s[j]==j: s[j]=i
        return s

    def factorize(x,spf):
        f=[]; v=x
        while v>1:
            p=spf[v]; e=0
            while v%p==0: v//=p; e+=1
            f.append((p,e))
        return f

    def cnt_sq(c,f_exp,spf):
        fac=factorize(c,spf); prod=1; all_even=True; e2=0
        for p,e0 in fac:
            e=e0*f_exp
            if p==2: e2=e
            elif e&1: all_even=False
            if p%4==3 and e&1: return 0
            if p%4==1: prod*=e+1
        r2=4*prod
        zero=4 if (all_even and e2%2==0) else 0
        eq=4 if (all_even and e2>=1 and (e2-1)%2==0) else 0
        return (r2-zero-eq)//8

    def divisors(fac):
        d=[1]
        for p,e in fac:
            old=len(d); mul=1
            for i in range(1,e+1):
                mul*=p
                for j in range(old): d.append(d[j]*mul)
        return d

    def cnt_cubes(n2,fac_n):
        divs=divisors(fac_n); count=0
        for u in divs:
            v=n2//u; u2=u*u
            if u2<=v: continue
            t=u2-v
            if t%3!=0: continue
            p=t//3
            if u2<4*p: continue
            delta=u2-4*p; s=isqrt(delta)
            if s*s!=delta: continue
            if s>=u or (u-s)&1: continue
            a=(u-s)//2; b=(u+s)//2
            if a==0 or a>=b: continue
            if a**3+b**3==n2: count+=1
        return count

    # Generate entries
    max_f=0
    while (1<<max_f)<=N and max_f<62: max_f+=1
    entries=[]; max_c=0
    for f in range(3,max_f+1):
        for c in range(2,10**7):
            v=pow_lim(c,f)
            if v>N: break
            entries.append((v,c,f))
            if c>max_c: max_c=c

    perfect_mask={}
    for v,c,f in entries:
        perfect_mask[v]=perfect_mask.get(v,0)|(1<<f)

    filt=[0]*(max_f+1)
    for e in range(max_f+1):
        m=0
        for f in range(3,max_f+1):
            if e==0 or f%e!=0: m|=1<<f
        filt[e]=m

    spf=spf_sieve(max_c)
    ans=0

    # Sum-of-two-squares
    for v,c,f in entries:
        ans+=cnt_sq(c,f,spf)

    # Sum-of-two-cubes
    cache3={}
    for v,c,f in entries:
        if f%3==0: continue
        if v in cache3: ans+=cache3[v]; continue
        fac_c=factorize(c,spf)
        fac_n=[(p,e*f) for p,e in fac_c]
        val=cnt_cubes(v,fac_n); cache3[v]=val; ans+=val

    # Higher powers (e>=4)
    for e in range(4,max_f+1):
        pw=[0]
        for b in range(1,10**7):
            v=pow_lim(b,e)
            if v>N: break
            pw.append(v)
        mb=len(pw)-1
        if mb<2: continue
        mf=filt[e]; even_e=(e&1)==0
        for b in range(2,mb+1):
            pb=pw[b]; amax=b-1
            while amax>0 and pw[amax]>N-pb: amax-=1
            if amax<=0: continue
            start=2 if (even_e and b&1) else 1
            step=2 if (even_e and b&1) else 1
            for a in range(start,amax+1,step):
                s=pb+pw[a]
                if s in perfect_mask:
                    ans+=bin(perfect_mask[s]&mf).count('1')
    return str(ans)

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

Java

import java.nio.file.*;
import java.util.*;
import java.util.regex.*;

public class Euler678 {
    private static final Pattern ANSWER_RE = Pattern.compile("answer\\s*:\\s*(.+)$", Pattern.CASE_INSENSITIVE);
    private static final Pattern EQUAL_RE = Pattern.compile("=\\s*(.+)$");

    private static String parseOutput(String stdout) {
        String[] lines = stdout.split("\\R");
        List<String> nonEmpty = new ArrayList<>();
        for (String line : lines) {
            String t = line.trim();
            if (!t.isEmpty()) {
                nonEmpty.add(t);
            }
        }
        if (nonEmpty.isEmpty()) {
            return "";
        }

        List<String> answers = new ArrayList<>();
        List<String> equals = new ArrayList<>();
        for (String line : nonEmpty) {
            Matcher m1 = ANSWER_RE.matcher(line);
            if (m1.find()) {
                answers.add(m1.group(1).trim());
            }
            Matcher m2 = EQUAL_RE.matcher(line);
            if (m2.find()) {
                equals.add(m2.group(1).trim());
            }
        }

        if (!answers.isEmpty()) {
            return answers.get(answers.size() - 1);
        }
        if (!equals.isEmpty()) {
            return equals.get(equals.size() - 1);
        }
        return nonEmpty.get(nonEmpty.size() - 1);
    }

    private static String pickCompiler() throws Exception {
        for (String compiler : List.of("clang++", "g++")) {
            Process probe = new ProcessBuilder("bash", "-lc", "command -v " + compiler)
                    .redirectErrorStream(true)
                    .start();
            String out = new String(probe.getInputStream().readAllBytes());
            int rc = probe.waitFor();
            if (rc == 0 && !out.trim().isEmpty()) {
                return compiler;
            }
        }
        throw new RuntimeException("No C++ compiler found (clang++/g++).");
    }

    private static Path cppSource(Path root) {
        return root.resolve("solutionsCpp").resolve("Euler678.cpp");
    }

    private static boolean shouldSkipCheckpoints(Path root) {
        Path src = cppSource(root);
        try {
            String text = Files.readString(src);
            return text.contains("--skip-checkpoints");
        } catch (Exception ex) {
            return false;
        }
    }

    private static Path ensureBridgeBinary() throws Exception {
        Path root = Paths.get(System.getProperty("user.dir"));
        Path src = cppSource(root);
        Path bin = root.resolve("solutionsCpp").resolve(".euler678_java_bridge");

        boolean rebuild = Files.notExists(bin)
                || Files.getLastModifiedTime(src).compareTo(Files.getLastModifiedTime(bin)) > 0;

        if (rebuild) {
            String compiler = pickCompiler();
            Process compile = new ProcessBuilder(
                    compiler,
                    "-std=c++17",
                    "-O2",
                    src.toString(),
                    "-o",
                    bin.toString())
                    .inheritIO()
                    .start();
            if (compile.waitFor() != 0) {
                throw new RuntimeException("Failed to compile Euler678 C++ bridge.");
            }
        }

        return bin;
    }

    private static String runBridge(Path bin, Path root, Path srcDir) throws Exception {
        List<String> cmd = new ArrayList<>();
        cmd.add(bin.toString());
        if (shouldSkipCheckpoints(root)) {
            cmd.add("--skip-checkpoints");
        }

        Process first = new ProcessBuilder(cmd)
                .directory(root.toFile())
                .redirectErrorStream(true)
                .start();
        String out = new String(first.getInputStream().readAllBytes());
        int rc = first.waitFor();
        if (rc == 0) {
            return out;
        }

        Process second = new ProcessBuilder(cmd)
                .directory(srcDir.toFile())
                .redirectErrorStream(true)
                .start();
        String out2 = new String(second.getInputStream().readAllBytes());
        int rc2 = second.waitFor();
        if (rc2 == 0) {
            return out2;
        }

        throw new RuntimeException("Euler678 C++ bridge failed.\n" + out + "\n" + out2);
    }

    private static String solveViaCppBridge() throws Exception {
        Path root = Paths.get(System.getProperty("user.dir"));
        Path src = cppSource(root);
        Path bin = ensureBridgeBinary();
        String out = runBridge(bin, root, src.getParent());
        String parsed = parseOutput(out);
        if (parsed.isEmpty()) {
            throw new RuntimeException("Euler678 C++ bridge produced empty output.");
        }
        return parsed;
    }

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