Problem 342: The Totient of a Square Is a Cube

Project Euler Problem 342 Solution

EulerSolve provides an optimized solution for Project Euler Problem 342, The Totient of a Square Is a Cube, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary We need the value of $$\sum_{\substack{1 \le n \lt 10^{10}\\ \varphi(n^2)\text{ is a cube}}} n.$$ A direct scan up to $10^{10}$ would be hopelessly slow, even with a fast totient routine. The key is to describe exactly which prime-exponent patterns make $\varphi(n^2)$ a perfect cube, and then construct only those $n$. The numerical Project Euler answer itself is intentionally omitted here. Mathematical Approach Write the prime factorization of $n$ as $$n=\prod_{p} p^{e_p},\qquad e_p\ge 0.$$ Only finitely many exponents are nonzero. The entire solution is built around turning the cube condition on $\varphi(n^2)$ into congruence constraints modulo $3$. Totient Formula for a Square Euler's totient function is multiplicative on coprime arguments, and for prime powers we have $$\varphi(p^k)=p^k-p^{k-1}=p^{k-1}(p-1).$$ Applying this to $n^2=\prod_{p \mid n} p^{2e_p}$ yields $$\varphi(n^2)=\prod_{p \mid n}\varphi\!\left(p^{2e_p}\right)=\prod_{p \mid n} p^{2e_p-1}(p-1).$$ Equivalently, since $\varphi(n)=n\prod_{p \mid n}\left(1-\frac1p\right)$, we also obtain the compact identity $$\varphi(n^2)=n\varphi(n).$$ This factorization is what exposes the prime exponents that must be checked modulo $3$. Cube Criterion via Prime Valuations An integer is a perfect cube if and only if every prime exponent in its factorization is divisible by $3$....

Detailed mathematical approach

Problem Summary

Euler's totient function $\varphi(m)$ counts the integers in $\{1,\dots,m\}$ coprime to $m$. The statement's example is $n=50$: here $50^2=2500=2^2\times 5^4$, so $\varphi(2500)=2\times 4\times 5^3$, which is $2^3\times 5^3=10^3$. Thus $2500$ is a square whose totient is a cube. We must find the sum of all $n$ with $1\lt n\lt 10^{10}$ for which $\varphi(n^2)$ is a perfect cube. ^[1]

The whole solution in one sentence: because every prime factor of $p-1$ is smaller than $p$, we can choose the prime factors of $n$ from the largest downwards so that every exponent in $\varphi(n^2)$ becomes a multiple of $3$, and then reach the remaining solutions by raising the chosen exponents in steps of $3$.

Mathematical Approach

1. The totient of a square

Two standard facts about $\varphi$ suffice: it is multiplicative, so $\varphi(ab)=\varphi(a)\varphi(b)$ when $\gcd(a,b)=1$, and $\varphi(p^k)=p^{k-1}(p-1)$ on a prime power. ^[2] Writing $n=\prod_p p^{e_p}$, so that $n^2=\prod_p p^{2e_p}$, they give

$$\varphi(n^2)=\prod_{p\mid n}p^{2e_p-1}(p-1).$$

The same calculation gives the compact identity $\varphi(n^2)=n\,\varphi(n)$, the case $k=2$ of $\varphi(n^k)=n^{k-1}\varphi(n)$. ^[2]

2. Cubes and prime exponents

In $72=2^3\cdot 3^2$ the prime $2$ carries the exponent $3$ and the prime $3$ carries $2$: we write $\nu_2(72)=3$, $\nu_3(72)=2$. In general $\nu_q(m)$ is the exponent of $q$ in $m$, zero when $q\nmid m$, so that $m=\prod_q q^{\nu_q(m)}$ and $\nu_q(xy)=\nu_q(x)+\nu_q(y)$. ^[3] By unique factorisation, ^[4] $m$ is a cube if and only if $3\mid\nu_q(m)$ for every $q$: if $m=c^3$ then $\nu_q(m)=3\nu_q(c)$, and conversely $c=\prod_q q^{\nu_q(m)/3}$ is then an integer with $c^3=m$. The problem is therefore a system of congruences modulo $3$, one per prime, on the exponents of $n$.

3. The exponent of a prime in the totient

Apply $\nu_q$ to the product of Section 1. The factor $p^{2e_p-1}$ contributes only when $p=q$, and each factor $p-1$ contributes $\nu_q(p-1)$. Writing $[q\mid n]$ for $1$ when $q$ divides $n$ and $0$ otherwise,

$$\nu_q\bigl(\varphi(n^2)\bigr)=[q\mid n]\,(2e_q-1)+s_q,$$

$$s_q=\sum_{p\mid n}\nu_q(p-1).$$

The first term is what the $q$-part of $n$ contributes itself; the pending sum $s_q$ collects what the chosen primes add through their factors $p-1$. By Section 2, $n$ is a solution precisely when

$$[q\mid n]\,(2e_q-1)+s_q\equiv 0 \pmod 3$$

for every prime $q$ — the congruence saying that $3$ divides the left-hand side.

4. The largest prime settles itself

One observation makes the problem finite and small. If $q$ divides $p-1$ then $q\le p-1\lt p$, so a chosen prime $p$ can only add to the pending sums of primes strictly below $p$; in particular $\nu_q(q-1)=0$, so no prime feeds its own pending sum.

Let $P$ be the largest prime factor of a solution $n\gt 1$. Every $p\mid n$ has $p\le P$, hence $p-1\lt P$ and $\nu_P(p-1)=0$. So $s_P=0$, and the congruence at $q=P$ reads $2e_P-1\equiv 0\pmod 3$, that is

$$e_P\equiv 2 \pmod 3 .$$

In particular $e_P\ge 2$, so $P^2\le P^{e_P}\le n\lt 10^{10}$ and $P\lt 10^5$. The other prime factors are smaller still, so every prime factor of every solution is at most $99999$, the integer part of $\sqrt{10^{10}-1}$. There are $9592$ such primes, and that finite list is the whole search space. The same congruence appears as a comment in the OEIS entry for this sequence. ^[5]

5. Three exponent classes

Fix $q$ and suppose its pending sum is known, $s_q\equiv c\pmod 3$ with $c\in\{0,1,2\}$. Section 3 then leaves two kinds of choice.

If $q\nmid n$, the condition is $c\equiv 0$: a prime may be left out only when its pending residue already vanishes.
If $q\mid n$ with exponent $e_q\ge 1$, the condition is $c+2e_q-1\equiv 0\pmod 3$. Multiplying $2e_q\equiv 1-c$ by $2$ and using $4\equiv 1$ gives $e_q\equiv 2+c\pmod 3$.

The admissible exponents thus form one arithmetic progression of step $3$, whose smallest member depends only on $c$: for $c=0$ they are $2,5,8,\dots$; for $c=1$, $3,6,9,\dots$; for $c=2$, $1,4,7,\dots$. Residue $0$ leaves a genuine branch — leave the prime out, or take it with exponent $2$ — while residues $1$ and $2$ force the prime into $n$, with minimal exponent $3$ and $1$ respectively.

6. Walking down the primes

Now visit the primes $99991,\dots,5,3,2$ in decreasing order. On arrival at $q$ its pending sum is already final: by Section 4 only larger primes can add to it, and all of those have been decided. Every decision is therefore local, and the search is ordinary backtracking. ^[6] Two tests cut the tree down: a prime with residue $0$ whose square no longer fits can only be left out, and is skipped without branching; a prime with nonzero residue whose minimal forced power does not fit kills the branch.

7. Lifting a solution by cubes

The walk gives every chosen prime its smallest admissible exponent, so it finds only part of the solutions; the rest come from a lifting rule. If $n$ satisfies the cube condition and $p\mid n$, so does $p^3n$, because raising $e_p$ by $3$ changes $2e_p-1$ by $6\equiv 0$ while the set of primes dividing $n$, hence every pending sum, is unchanged; the lifted number is a solution as long as it stays below $10^{10}$. ^[5]

Conversely every solution arises this way exactly once. For each prime $p\mid n$ let $m_p$ be the unique element of $\{1,2,3\}$ with $m_p\equiv e_p\pmod 3$; then $m_p\le e_p$, since $e_p\equiv 0\pmod 3$ forces $e_p\ge 3$. The base $n_0=\prod_{p\mid n}p^{m_p}$ has the same prime divisors as $n$, so its pending sums are the same while each $2e_p-1$ drops by a multiple of $6$; by Section 3 it is a solution too, and $n_0\le n$. Hence

$$n=n_0\prod_{p\mid n_0}p^{3t_p},\qquad t_p\ge 0,$$

with $t_p=(e_p-m_p)/3$ uniquely determined. The solutions split into disjoint families, one per base, and the walk of Section 6 produces precisely the bases. Summing a family is a nested sum: with $R=\lfloor(10^{10}-1)/n_0\rfloor$ we need the sum of all products $\prod p^{3t_p}$ at most $R$, which a short recursion over the primes of $n_0$ evaluates exactly.

8. A complete example by hand

Take the primes $5$, $3$ and $2$ and follow Section 6; all larger primes have residue $0$ and are left out. Since $5-1=2^2$ and $3-1=2$, taking $5$ adds $2$ to the residue at the prime $2$ and taking $3$ adds $1$ there, while $2-1=1$ adds nothing anywhere.

At $5$: the residue is $0$, so $5$ is optional, with exponent $2$ if taken.
At $3$: the residue is still $0$ either way, so $3$ is optional too, again with exponent $2$.
At $2$: the residue $c$ now follows from the earlier decisions. Neither $5$ nor $3$ taken: $c=0$, so $2$ is optional, giving $n=1$, excluded by the range, or $n=2^2=4$. Only $3$: $c=1$ forces $e_2=3$ and $n=2^3 3^2=72$. Only $5$: $c=2$ forces $e_2=1$ and $n=2\cdot 5^2=50$. Both: $c=2+1\equiv 0$, so $2$ is optional again, giving $n=3^2 5^2=225$ or $n=2^2 3^2 5^2=900$.

The five leaves give the bases $4$, $50$, $72$, $225$ and $900$, and the totients of their squares are cubes: $\varphi(4^2)=8=2^3$, $\varphi(50^2)=1000=10^3$, $\varphi(72^2)=1728=12^3$, $\varphi(225^2)=27000=30^3$ and $\varphi(900^2)=216000=60^3$. Lifting each base by cubes of its own primes and staying below $1000$ gives $4,32,256$ from $4$; $50,400$ from $50$; $72,576$ from $72$; and the bases $225$ and $900$ alone, as $225\cdot 3^3=6075$ and $900\cdot 2^3=7200$ are already too large. Together:

$$4,\,32,\,50,\,72,\,225,\,256,\,400,\,576,\,900,$$

exactly the list a brute-force scan of $n\lt 1000$ returns, with sum $2515$; the two lifts above, and $225\cdot 5^3=28125$, are solutions further out.

9. An example with a larger prime

Start at $19$ with residue $0$ and take it with exponent $2$. As $19-1=18=2\cdot 3^2$, this adds $2$ at the prime $3$ and $1$ at the prime $2$. At $3$ the residue $2$ forces $e_3=1$, and $3-1=2$ lifts the residue at $2$ to $1+1=2$, which forces $e_2=1$. The base is $2\cdot 3\cdot 19^2=2166$, with $\varphi(2166^2)=1481544=114^3$. The largest prime occurring in any base here is $62501$: as $62501-1=2^2\cdot 5^6$ and $3\mid 6$, the prime $5$ is untouched while $2$ is forced with $e_2=1$, giving the base $2\cdot 62501^2=7812750002$; the bare square $62501^2$ is not a solution, its residue at $2$ being $2\ne 0$.

How the Code Works

All three programs implement the same algorithm in a Solver class built with the exclusive limit $10^{10}$; the identifiers below are the C++ and Python ones, which Java spells buildPrimesAndSpf, applyFactors, addSolution and multiplierSumRec. The field max_prime is the integer square root of $10^{10}-1$, namely $99999$ — the bound of Section 4. build_primes_and_spf sets spf[j] = j, then for each prime $i$ up to $\sqrt{99999}$ walks the multiples of $i$ from $i^2$ upwards, overwriting spf[j] only where it still equals $j$; ^[7] afterwards spf[j] is the smallest prime factor of $j$, and the fixed points from $2$ on are the $9592$ primes, which primes_desc holds reversed as in Section 6. The same array factorises every $p-1$, and factors_mod3_of_p_minus_1[p] keeps the pairs $(q,\nu_q(p-1)\bmod 3)$ with nonzero residue — $29786$ out of $31638$ prime powers, the rest being invisible modulo $3$. The residues of Section 5 live in residue_mod3, and apply_factors(p, sign) adds those pairs when the walk takes $p$ and subtracts them on backtracking.

dfs(idx, current_n) is the walk, and its first loop is the skip of Section 6. If that loop runs off the end of primes_desc, the partial product is a complete base and goes to add_solution; otherwise the prime $p$ at the current index is decided by its residue $c$. For $c=0$ the C++ and Python versions first recurse with $p$ left out, then, if $p^2$ fits, apply the factors of $p$, record $p$ in included_primes, recurse with $current\_n\cdot p^2$ and undo both. For $c\ne0$ the minimal exponent is m = (2 + c) % 3, read as $3$ when it comes out $0$, exactly as in Section 5, and the single forced branch is taken, or dropped when that power does not fit.

add_solution(base_n) discards $base\_n\le1$, the empty choice $n=1$, and otherwise sets the budget $R=\lfloor(limit-1)/base\_n\rfloor$ and calls multiplier_sum_rec(0, R), the family sum of Section 7. Once every recorded prime has been passed it returns $1$, the multiplier of the empty product; at each earlier position it runs that prime's multiplier through $1,p^3,p^6,\dots$ while this stays within the budget, adds it times the sum returned for the later primes with budget remaining / mul, and stops when another $p^3$ would overshoot. For the base $72$ under a limit of $1000$ the budget is $R=13$, the multipliers are $1$ and $8$, and $72\cdot 9=648=72+576$; in general $base\_n$ times the returned sum goes into answer.

The implementations differ in integer arithmetic, recursion and self-checks. C++ accumulates in unsigned __int128, a compiler extension, ^[8] prints through to_string_u128, and carries the only self-checks: a totient-sieve brute_sum, an exact is_perfect_cube whose integer_cuberoot_floor corrects a cbrt estimate up and down, and run_checkpoints, disabled by --skip-checkpoints. Java works in signed 64-bit long, with quotient guards before the products in dfs but plain arithmetic in the family sum; that is safe because the $2395$ solutions are each below $10^{10}$, so the total, and every quantity accumulated on the way to it, stays under $2395\cdot 10^{10}$. Its dfs also takes the inclusion branch first and handles exclusion by continuing its enclosing loop with idx = i + 1, which keeps its recursion depth at $7$ rather than $9593$. Python needs no overflow care but mirrors the recursive C++ shape, so the file lifts the interpreter's stack limit from its default $1000$ with sys.setrecursionlimit(20000), without which it would stop with a RecursionError. ^[9]

Complexity and Verification

Preparation is cheap and bounded: the sieve costs $O(M\log\log M)$ for $M=99999$, ^[7] and factorising every $p-1$ through spf takes $O(\log p)$ divisions. Storage is $O(M)$ machine words with a large constant: the C++ containers come to about $3.1$ MB, of which $2.4$ MB is the $10^5$ factor-list objects and only $0.5$ MB the sieve and residue tables, while Python's lists of references cost several times that. C++ and Python also hold $9593$ recursion frames at the deepest point.

The search has no clean closed form: it is a backtracking tree, small in practice only because residues $1$ and $2$ leave no choice and the limit kills branches early. On the target input it makes $34683$ calls to dfs, of which $15669$ die at once because a forced power does not fit, and $36192$ calls to apply_factors; only $106$ of the $9592$ primes ever enter a base, at most $6$ at a time. The real work is the skip loop, whose body runs $42479502$ times, $42445737$ of those iterations advancing past a prime that can only be left out. Jumping straight to the first prime whose square still fits would be wrong, because the scan must also stop at the next prime with nonzero residue, which forces a decision; all three versions keep the linear scan. Measured once on one machine, the run times are $0.04$ s for C++ with --skip-checkpoints, $0.23$ s with the checkpoint, $0.13$ s for Java and $2.4$ s for Python — against the $10^{10}$ totient values a plain scan would need.

The C++ file packages two checks. The first asserts that $50\cdot\varphi(50)=50\cdot20$ is a cube, the statement's example; since it uses the identity $\varphi(n^2)=n\varphi(n)$ rather than $\varphi(2500)$ itself, it tests the cube test and not the identity. The second is a genuine cross-check: brute_sum(200000) sieves the totients below $200000$, keeps the $n$ with $n\varphi(n)$ a cube, and agrees with the fast solver at $3673238$.

This article adds the following checks. The identity $\varphi(n^2)=n\varphi(n)$ and the product form of Section 1 hold for all $n\le3000$, and the rule of Section 3 reproduces the true factorisation of $\varphi(n^2)$ for all $n\le2000$. A brute-force totient sieve with an exact integer cube test agrees with the enumeration element for element at the limits $10^3$, $10^4$, $2\cdot10^5$, $10^6$ and $10^7$, whose sums are $2515$, $62061$, $3673238$, $23294055$ and $563098362$. At the full limit the $917$ bases expand to $2395$ distinct values with no repetition, each re-tested by factorising it from scratch against Section 3, and all satisfy the lemma of Section 4. The first terms match the OEIS entry, ^[5] and the three programs print the same total.

Footnotes and References

The sources below supply the background facts used above; the valuation congruence, the bound on the prime factors, the exponent classes and the cube-lifted families are derived in this article.

Project Euler 342 — The Totient of a Square Is a Cube. The official statement works out n = 50, where 50² = 2500 = 2² × 5⁴ and φ(2500) = 2³ × 5³, and asks for the sum of all n with 1 < n < 10¹⁰ such that φ(n²) is a cube.
Wikipedia — Euler's totient function. It states that φ(mn) = φ(m)φ(n) for relatively prime m and n, gives φ(p^k) = p^(k−1)(p − 1) for a prime p and k ≥ 1, and lists the identity φ(n^k) = n^(k−1)φ(n); Section 1 uses all three.
Wikipedia — p-adic valuation. It defines the valuation of n as the exponent of the highest power of the prime p dividing n and states the product rule ν_p(rs) = ν_p(r) + ν_p(s), both used in Sections 2 and 3.
Wikipedia — Fundamental theorem of arithmetic. Every integer greater than 1 is a prime or is uniquely a product of primes, up to the order of the factors. The article does not discuss perfect powers; the cube criterion of Section 2 is proved from that uniqueness.
OEIS A114076 — numbers k such that k·φ(k) is a cube. Since φ(k²) = kφ(k) this is the sequence of solutions, and its terms 1, 4, 32, 50, 72, 225, 256, 400, 576, 900, 1944, 2048, 2166, … agree with the enumeration above. Two of its comments are the facts proved in Sections 4 and 7: for every term above 1 the exponent of the largest prime factor, A071178, is 2 modulo 3; and p³ times a term is again a term when p divides it.
Wikipedia — Backtracking. It describes building candidates incrementally and abandoning one as soon as it cannot possibly be completed, pruning the whole subtree below it, which is the scheme of Section 6.
Wikipedia — Sieve of Eratosthenes. It describes enumerating the multiples of each prime i from i² upwards and gives the running time O(n log log n); the code stores the smallest prime factor instead of a plain mark, which turns the table into a factorising device for the numbers p − 1.
GCC manual — 128-bit integers. The manual states that __int128 is supported as an extension on targets with a wide enough integer mode, unsigned __int128 being its unsigned form.
Python documentation — sys.setrecursionlimit. The function sets the maximum depth of the interpreter stack, a limit that stops runaway recursion from overflowing the C stack, and a value too low for the current depth raises RecursionError.

Problem 342 source code

C++

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

namespace {

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

std::string to_string_u128(u128 value) {
    if (value == 0) {
        return "0";
    }
    std::string s;
    while (value > 0) {
        const int digit = static_cast<int>(value % 10);
        s.push_back(static_cast<char>('0' + digit));
        value /= 10;
    }
    std::reverse(s.begin(), s.end());
    return s;
}

u64 integer_cuberoot_floor(const u128 x) {
    long double approx = std::cbrt(static_cast<long double>(x));
    u64 r = static_cast<u64>(approx);

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

bool is_perfect_cube(const u128 x) {
    const u64 r = integer_cuberoot_floor(x);
    return static_cast<u128>(r) * r * r == x;
}

u128 brute_sum(const int limit) {
    std::vector<int> phi(static_cast<std::size_t>(limit));
    for (int i = 0; i < limit; ++i) {
        phi[static_cast<std::size_t>(i)] = i;
    }

    for (int p = 2; p < limit; ++p) {
        if (phi[static_cast<std::size_t>(p)] != p) {
            continue;
        }
        for (int x = p; x < limit; x += p) {
            phi[static_cast<std::size_t>(x)] -= phi[static_cast<std::size_t>(x)] / p;
        }
    }

    u128 sum = 0;
    for (int n = 2; n < limit; ++n) {
        const u128 value = static_cast<u128>(n) * static_cast<u128>(phi[static_cast<std::size_t>(n)]);
        if (is_perfect_cube(value)) {
            sum += static_cast<u128>(n);
        }
    }
    return sum;
}

struct Solver {
    u64 limit_exclusive;
    int max_prime;

    std::vector<int> primes_desc;
    std::vector<int> spf;
    std::vector<std::vector<std::pair<int, int>>> factors_mod3_of_p_minus_1;

    std::vector<std::uint8_t> residue_mod3;
    std::vector<int> included_primes;

    u128 answer = 0;

    explicit Solver(const u64 limit)
        : limit_exclusive(limit), max_prime(static_cast<int>(std::sqrt(static_cast<long double>(limit - 1)))) {}

    void build_primes_and_spf() {
        spf.assign(static_cast<std::size_t>(max_prime + 1), 0);
        for (int i = 0; i <= max_prime; ++i) {
            spf[static_cast<std::size_t>(i)] = i;
        }

        for (int i = 2; static_cast<i64>(i) * i <= max_prime; ++i) {
            if (spf[static_cast<std::size_t>(i)] != i) {
                continue;
            }
            for (int j = i * i; j <= max_prime; j += i) {
                if (spf[static_cast<std::size_t>(j)] == j) {
                    spf[static_cast<std::size_t>(j)] = i;
                }
            }
        }

        std::vector<int> primes_asc;
        for (int i = 2; i <= max_prime; ++i) {
            if (spf[static_cast<std::size_t>(i)] == i) {
                primes_asc.push_back(i);
            }
        }

        primes_desc = primes_asc;
        std::reverse(primes_desc.begin(), primes_desc.end());

        factors_mod3_of_p_minus_1.assign(static_cast<std::size_t>(max_prime + 1), {});
        for (int p : primes_asc) {
            int x = p - 1;
            auto& vec = factors_mod3_of_p_minus_1[static_cast<std::size_t>(p)];
            while (x > 1) {
                const int q = spf[static_cast<std::size_t>(x)];
                int cnt = 0;
                while (x % q == 0) {
                    x /= q;
                    ++cnt;
                }
                cnt %= 3;
                if (cnt != 0) {
                    vec.emplace_back(q, cnt);
                }
            }
        }

        residue_mod3.assign(static_cast<std::size_t>(max_prime + 1), 0);
    }

    void apply_factors(const int p, const int sign) {
        for (const auto& [q, e] : factors_mod3_of_p_minus_1[static_cast<std::size_t>(p)]) {
            int v = static_cast<int>(residue_mod3[static_cast<std::size_t>(q)]);
            if (sign > 0) {
                v += e;
            } else {
                v -= e;
            }
            v %= 3;
            if (v < 0) {
                v += 3;
            }
            residue_mod3[static_cast<std::size_t>(q)] = static_cast<std::uint8_t>(v);
        }
    }

    u64 pow_u64(const u64 base, const int exp) const {
        u64 value = 1;
        for (int i = 0; i < exp; ++i) {
            value *= base;
        }
        return value;
    }

    u128 multiplier_sum_rec(const int pos, const u64 remaining) const {
        if (pos == static_cast<int>(included_primes.size())) {
            return 1;
        }

        const u64 p = static_cast<u64>(included_primes[static_cast<std::size_t>(pos)]);
        const u128 p3 = static_cast<u128>(p) * p * p;

        u128 total = 0;
        u128 mul = 1;

        while (mul <= static_cast<u128>(remaining)) {
            const u64 next_remaining = remaining / static_cast<u64>(mul);
            total += mul * multiplier_sum_rec(pos + 1, next_remaining);

            if (mul > static_cast<u128>(remaining) / p3) {
                break;
            }
            mul *= p3;
        }

        return total;
    }

    void add_solution(const u64 base_n) {
        if (base_n <= 1 || base_n >= limit_exclusive) {
            return;
        }

        const u64 remaining = (limit_exclusive - 1) / base_n;
        const u128 mul_sum = multiplier_sum_rec(0, remaining);
        answer += static_cast<u128>(base_n) * mul_sum;
    }

    void dfs(int idx, const u64 current_n) {
        int i = idx;
        const int npr = static_cast<int>(primes_desc.size());

        // Skip deterministic exclusions: residue is zero and p^2 already too large.
        while (i < npr) {
            const int p = primes_desc[static_cast<std::size_t>(i)];
            const int c = residue_mod3[static_cast<std::size_t>(p)];
            if (c != 0) {
                break;
            }

            const u128 min_include = static_cast<u128>(p) * p;
            if (static_cast<u128>(current_n) * min_include < static_cast<u128>(limit_exclusive)) {
                break;
            }
            ++i;
        }

        if (i >= npr) {
            add_solution(current_n);
            return;
        }

        const int p = primes_desc[static_cast<std::size_t>(i)];
        const int c = residue_mod3[static_cast<std::size_t>(p)];

        if (c == 0) {
            // Exclude p.
            dfs(i + 1, current_n);

            // Include p with exponent congruent to 2 (mod 3), minimal exponent 2.
            const int e0 = 2;
            const u64 power = pow_u64(static_cast<u64>(p), e0);
            if (static_cast<u128>(current_n) * static_cast<u128>(power) < static_cast<u128>(limit_exclusive)) {
                apply_factors(p, +1);
                included_primes.push_back(p);
                dfs(i + 1, current_n * power);
                included_primes.pop_back();
                apply_factors(p, -1);
            }
            return;
        }

        // Forced inclusion when residue is nonzero.
        const int m = (2 + c) % 3;
        const int e0 = (m == 0) ? 3 : m;
        const u64 power = pow_u64(static_cast<u64>(p), e0);
        if (static_cast<u128>(current_n) * static_cast<u128>(power) >= static_cast<u128>(limit_exclusive)) {
            return;
        }

        apply_factors(p, +1);
        included_primes.push_back(p);
        dfs(i + 1, current_n * power);
        included_primes.pop_back();
        apply_factors(p, -1);
    }

    u128 solve() {
        build_primes_and_spf();
        answer = 0;
        included_primes.clear();
        dfs(0, 1ULL);
        return answer;
    }
};

bool run_checkpoints() {
    // Problem statement example: n=50 qualifies.
    const u128 phi_50_sq = static_cast<u128>(50ULL) * 20ULL;  // phi(50)=20, and phi(n^2)=n*phi(n)
    if (!is_perfect_cube(phi_50_sq)) {
        std::cerr << "Checkpoint failed: n=50 should satisfy cube condition\n";
        return false;
    }

    // Cross-check on a smaller limit with brute force.
    const int small_limit = 200000;
    Solver solver_small(static_cast<u64>(small_limit));
    const u128 fast = solver_small.solve();
    const u128 slow = brute_sum(small_limit);
    if (fast != slow) {
        std::cerr << "Checkpoint failed: fast/brute mismatch for limit=" << small_limit << '\n';
        return false;
    }

    return true;
}

}  // namespace

int main(int argc, char** argv) {
    bool skip_checkpoints = false;
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            skip_checkpoints = true;
        } else {
            std::cerr << "Unknown argument: " << arg << '\n';
            return 1;
        }
    }

    if (!skip_checkpoints && !run_checkpoints()) {
        return 2;
    }

    Solver solver(10000000000ULL);
    const u128 answer = solver.solve();
    std::cout << to_string_u128(answer) << '\n';
    return 0;
}

Python

import math
import sys

# Increase recursion depth if needed
sys.setrecursionlimit(20000)

class Solver:
    def __init__(self, limit):
        self.limit_exclusive = limit
        self.max_prime = int(math.isqrt(limit - 1))
        
        self.primes_desc = []
        self.spf = []
        self.factors_mod3_of_p_minus_1 = []
        self.residue_mod3 = []
        self.included_primes = []
        self.answer = 0

    def build_primes_and_spf(self):
        self.spf = list(range(self.max_prime + 1))
        for i in range(2, math.isqrt(self.max_prime) + 1):
            if self.spf[i] != i:
                continue
            for j in range(i * i, self.max_prime + 1, i):
                if self.spf[j] == j:
                    self.spf[j] = i
                    
        primes_asc = [i for i in range(2, self.max_prime + 1) if self.spf[i] == i]
        self.primes_desc = primes_asc[::-1]
        
        self.factors_mod3_of_p_minus_1 = [[] for _ in range(self.max_prime + 1)]
        for p in primes_asc:
            x = p - 1
            vec = self.factors_mod3_of_p_minus_1[p]
            while x > 1:
                q = self.spf[x]
                cnt = 0
                while x % q == 0:
                    x //= q
                    cnt += 1
                cnt %= 3
                if cnt != 0:
                    vec.append((q, cnt))
                    
        self.residue_mod3 = [0] * (self.max_prime + 1)

    def apply_factors(self, p, sign):
        for q, e in self.factors_mod3_of_p_minus_1[p]:
            v = self.residue_mod3[q]
            if sign > 0:
                v += e
            else:
                v -= e
            v %= 3
            if v < 0:
                v += 3
            self.residue_mod3[q] = v

    def multiplier_sum_rec(self, pos, remaining):
        if pos == len(self.included_primes):
            return 1
            
        p = self.included_primes[pos]
        p3 = p * p * p
        
        total = 0
        mul = 1
        
        while mul <= remaining:
            next_remaining = remaining // mul
            total += mul * self.multiplier_sum_rec(pos + 1, next_remaining)
            
            if mul > remaining // p3:
                break
            mul *= p3
            
        return total

    def add_solution(self, base_n):
        if base_n <= 1 or base_n >= self.limit_exclusive:
            return
            
        remaining = (self.limit_exclusive - 1) // base_n
        mul_sum = self.multiplier_sum_rec(0, remaining)
        self.answer += base_n * mul_sum

    def dfs(self, idx, current_n):
        i = idx
        npr = len(self.primes_desc)
        
        while i < npr:
            p = self.primes_desc[i]
            c = self.residue_mod3[p]
            if c != 0:
                break
                
            min_include = p * p
            if current_n * min_include < self.limit_exclusive:
                break
            i += 1
            
        if i >= npr:
            self.add_solution(current_n)
            return
            
        p = self.primes_desc[i]
        c = self.residue_mod3[p]
        
        if c == 0:
            self.dfs(i + 1, current_n)
            
            e0 = 2
            power = p ** e0
            if current_n * power < self.limit_exclusive:
                self.apply_factors(p, +1)
                self.included_primes.append(p)
                self.dfs(i + 1, current_n * power)
                self.included_primes.pop()
                self.apply_factors(p, -1)
            return
            
        m = (2 + c) % 3
        e0 = 3 if m == 0 else m
        power = p ** e0
        if current_n * power >= self.limit_exclusive:
            return
            
        self.apply_factors(p, +1)
        self.included_primes.append(p)
        self.dfs(i + 1, current_n * power)
        self.included_primes.pop()
        self.apply_factors(p, -1)

    def run(self):
        self.build_primes_and_spf()
        self.answer = 0
        self.included_primes = []
        self.dfs(0, 1)
        return self.answer

def solve():
    limit = 10000000000
    solver = Solver(limit)
    ans = solver.run()
    return str(ans)

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

Java

import java.util.*;

public class Euler342 {

    static class Solver {
        long limit_exclusive;
        int max_prime;

        List<Integer> primes_desc = new ArrayList<>();
        int[] spf;
        List<List<int[]>> factorsMod3OfPMinus1;
        byte[] residueMod3;
        List<Integer> includedPrimes = new ArrayList<>();
        long answer = 0;

        Solver(long limit) {
            this.limit_exclusive = limit;
            this.max_prime = (int) Math.sqrt(limit - 1);
        }

        void buildPrimesAndSpf() {
            spf = new int[max_prime + 1];
            for (int i = 0; i <= max_prime; i++)
                spf[i] = i;

            for (int i = 2; i <= Math.sqrt(max_prime); i++) {
                if (spf[i] != i)
                    continue;
                for (int j = i * i; j <= max_prime; j += i) {
                    if (spf[j] == j)
                        spf[j] = i;
                }
            }

            List<Integer> primesAsc = new ArrayList<>();
            for (int i = 2; i <= max_prime; i++) {
                if (spf[i] == i)
                    primesAsc.add(i);
            }

            for (int i = primesAsc.size() - 1; i >= 0; i--) {
                primes_desc.add(primesAsc.get(i));
            }

            factorsMod3OfPMinus1 = new ArrayList<>(max_prime + 1);
            for (int i = 0; i <= max_prime; i++) {
                factorsMod3OfPMinus1.add(new ArrayList<>());
            }

            for (int p : primesAsc) {
                int x = p - 1;
                List<int[]> vec = factorsMod3OfPMinus1.get(p);
                while (x > 1) {
                    int q = spf[x];
                    int cnt = 0;
                    while (x % q == 0) {
                        x /= q;
                        cnt++;
                    }
                    cnt %= 3;
                    if (cnt != 0) {
                        vec.add(new int[] { q, cnt });
                    }
                }
            }

            residueMod3 = new byte[max_prime + 1];
        }

        void applyFactors(int p, int sign) {
            for (int[] factor : factorsMod3OfPMinus1.get(p)) {
                int q = factor[0];
                int e = factor[1];
                int v = residueMod3[q];
                if (sign > 0) {
                    v += e;
                } else {
                    v -= e;
                }
                v %= 3;
                if (v < 0) {
                    v += 3;
                }
                residueMod3[q] = (byte) v;
            }
        }

        long multiplierSumRec(int pos, long remaining) {
            if (pos == includedPrimes.size())
                return 1;

            long p = includedPrimes.get(pos);
            long p3 = p * p * p;

            long total = 0;
            long mul = 1;

            while (mul <= remaining) {
                long next_remaining = remaining / mul;
                total += mul * multiplierSumRec(pos + 1, next_remaining);

                if (mul > remaining / p3)
                    break;
                mul *= p3;
            }
            return total;
        }

        void addSolution(long baseN) {
            if (baseN <= 1 || baseN >= limit_exclusive)
                return;
            long remaining = (limit_exclusive - 1) / baseN;
            long mulSum = multiplierSumRec(0, remaining);
            answer += baseN * mulSum;
        }

        void dfs(int idx, long currentN) {
            while (true) {
                int i = idx;
                int npr = primes_desc.size();

                while (i < npr) {
                    int p = primes_desc.get(i);
                    int c = residueMod3[p];
                    if (c != 0)
                        break;

                    long minInclude = (long) p * p;
                    // If it overflows limit_exclusive or equals it, we just continue loop to
                    // exclude it
                    // Using BigInteger or just division to check overflow
                    if (currentN < limit_exclusive / minInclude + 1) { // currentN * minInclude < limit
                        // Need to check exact currentN * minInclude < limit
                        long product = currentN * minInclude;
                        if (currentN != 0 && product / currentN == minInclude && product < limit_exclusive) {
                            break;
                        }
                    }
                    i++;
                }

                if (i >= npr) {
                    addSolution(currentN);
                    return;
                }

                int p = primes_desc.get(i);
                int c = residueMod3[p];

                if (c == 0) {
                    long power = (long) p * p;
                    if (currentN < limit_exclusive / power + 1) {
                        long product = currentN * power;
                        if (currentN != 0 && product / currentN == power && product < limit_exclusive) {
                            applyFactors(p, +1);
                            includedPrimes.add(p);
                            dfs(i + 1, product);
                            includedPrimes.remove(includedPrimes.size() - 1);
                            applyFactors(p, -1);
                        }
                    }
                    idx = i + 1;
                    continue;
                }

                int m = (2 + c) % 3;
                int e0 = (m == 0) ? 3 : m;
                long power = 1;
                for (int j = 0; j < e0; j++)
                    power *= p;

                if (currentN >= limit_exclusive / power + 1)
                    return;
                long product = currentN * power;
                if (currentN != 0 && product / currentN != power)
                    return;
                if (product >= limit_exclusive)
                    return;

                applyFactors(p, +1);
                includedPrimes.add(p);
                dfs(i + 1, product);
                includedPrimes.remove(includedPrimes.size() - 1);
                applyFactors(p, -1);
                return;
            }
        }

        long solve() {
            buildPrimesAndSpf();
            answer = 0;
            includedPrimes.clear();
            dfs(0, 1L);
            return answer;
        }
    }

    public static String solve() {
        Solver solver = new Solver(10000000000L);
        long ans = solver.solve();
        return String.valueOf(ans);
    }

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