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

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$. Using the $p$-adic valuation $\nu_q(m)$, this becomes

$$\nu_q\bigl(\varphi(n^2)\bigr)\equiv 0\pmod 3\qquad\text{for every prime }q.$$

From the product formula above, each prime $q$ receives contributions from two sources:

$$\nu_q\bigl(\varphi(n^2)\bigr)=\mathbf{1}_{q\mid n}(2e_q-1)+\sum_{p \mid n}\nu_q(p-1).$$

The first term is the exponent contributed by the $q$-part of $n$ itself. The second term records how often $q$ appears inside the factors $p-1$ contributed by every chosen prime $p$.

Why a Descending Prime DFS Works

The code processes primes in descending order. This is crucial because every prime divisor of $p-1$ is strictly smaller than $p$. Therefore, once the search reaches a prime $q$, no later unprocessed prime can create a new $q$-contribution through a factor of the form $p-1$.

So when $q$ is examined, the residue

$$c_q\equiv\sum_{\substack{p \mid n\\ p>q}}\nu_q(p-1)\pmod 3$$

is already final. The array residue_mod3 stores exactly these pending residues. This makes the decision at each prime local rather than global, which is why the DFS can prune so aggressively.

Forced and Optional Exponent Classes

Suppose the current residue at prime $q$ is $c\in\{0,1,2\}$.

If $q$ is excluded from $n$, then its own contribution $2e_q-1$ does not exist. Since no smaller future prime can add more $q$-valuation, exclusion is legal only when $c=0$.

If $q$ is included with exponent $e_q$, then we must satisfy

$$c+(2e_q-1)\equiv 0\pmod 3.$$

Solving modulo $3$ gives

$$e_q\equiv c+2\pmod 3.$$

Hence the minimal admissible exponents are

$$c=0 \Rightarrow e_q=2,\qquad c=1 \Rightarrow e_q=3,\qquad c=2 \Rightarrow e_q=1.$$

This is exactly the branch structure used in the implementation. Residue $0$ means the prime is optional; residues $1$ and $2$ force inclusion with the smallest matching exponent class.

Worked Example: $n = 72$

A small nontrivial valid example is

$$72=2^3\cdot 3^2.$$

Using the prime-power formula,

$$\varphi(72^2)=2^{2\cdot 3-1}(2-1)\cdot 3^{2\cdot 2-1}(3-1)=2^5\cdot 1\cdot 3^3\cdot 2=2^6\cdot 3^3=12^3.$$

The residue interpretation is illuminating. Prime $3$ may start with exponent $2$ because residue $0$ requires the class $e_3\equiv 2\pmod 3$. But $3-1=2$ contributes one pending factor of $2$, so when the DFS later reaches prime $2$ it sees residue $c=1$. Then

$$1+(2e_2-1)\equiv 0\pmod 3\Rightarrow e_2\equiv 0\pmod 3,$$

and the minimal choice is $e_2=3$. This is how the search reconstructs $72$ from congruence data alone.

Generating Whole Families from a Minimal Base

Once a minimal valid base solution $n_0$ is fixed, every included prime exponent may be increased by multiples of $3$:

$$n=n_0\prod_{p\in\mathcal{P}} p^{3t_p},\qquad t_p\ge 0.$$

This preserves the cube condition because adding $3$ to an exponent changes $2e_p-1$ by $6$, which is $0\pmod 3$, while the factors $p-1$ remain unchanged. The program therefore enumerates each minimal base only once, and then sums the entire cube-power family attached to it.

How the Code Works

The solver first builds a smallest-prime-factor sieve up to $\left\lfloor\sqrt{10^{10}-1}\right\rfloor$. That range is enough for minimal base factors: whenever a new prime is introduced from residue $0$, its smallest admissible exponent is $2$, so it must satisfy $p^2 \lt 10^{10}$.

Next, the code pre-factors every $p-1$ and stores only the exponents modulo $3$ in factors_mod3_of_p_minus_1. During DFS, choosing a prime $p$ simply adds those stored contributions into residue_mod3.

The DFS proceeds from large primes to small primes. At each step it reads the current residue, decides whether the prime is excluded, optional, or forced, multiplies the current base by the smallest valid power, and updates the residue table. A useful prune appears when the residue is $0$ but even $p^2$ would exceed the limit: then that prime cannot start a new branch and is skipped deterministically.

When the search reaches the end, it has produced one minimal base solution. The auxiliary recursion multiplier_sum_rec then sums all admissible multipliers $1,p^3,p^6,\dots$ for the included primes, so the C++, Java, and Python versions all count whole families without rediscovering them one by one.

Complexity Analysis

The preprocessing stage is essentially a sieve plus factorizations of $p-1$ for primes up to $\sqrt{L}$, where $L=10^{10}$. The DFS itself does not have a neat closed-form worst-case bound; in the abstract it is exponential in the number of candidate primes, but in practice the tree is much smaller because residues frequently force a unique decision and the limit cuts off large branches early.

Compared with brute-force totient evaluation for every $n \lt 10^{10}$, the savings are enormous. Memory usage stays modest: mainly the SPF table, the precomputed factor lists, the residue array, and the recursion stack.

Footnotes and References

Problem page: https://projecteuler.net/problem=342
Euler's totient function: https://en.wikipedia.org/wiki/Euler%27s_totient_function
Multiplicative function: https://en.wikipedia.org/wiki/Multiplicative_function
$p$-adic valuation: https://en.wikipedia.org/wiki/P-adic_valuation
Hardy, G. H.; Wright, E. M. An Introduction to the Theory of Numbers, 6th ed., Oxford University Press. Chapters on multiplicative functions and prime factorization.

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());
    }
}