Problem 302: Strong Achilles Numbers

Project Euler Problem 302 Solution

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

Problem Summary A strong Achilles number is an integer $n$ such that both $n$ and $\varphi(n)$ are Achilles numbers. An integer is Achilles when it is powerful, but not a perfect power. The problem asks for the number of such $n$ below $10^{18}$. Mathematical Approach Write $$n=\prod_{i=1}^r p_i^{a_i},\qquad a_i\ge 1.$$ 1) Achilles condition in exponent language The number $n$ is powerful exactly when every exponent is at least $2$: $$a_i\ge 2\qquad \text{for all }i.$$ It is a perfect $k$-th power exactly when all exponents are divisible by the same $k \ge 2$, equivalently when $$\gcd(a_1,a_2,\dots,a_r)\ge 2.$$ Therefore $$n\text{ is Achilles}\iff a_i\ge 2\ \forall i\quad\text{and}\quad \gcd(a_1,\dots,a_r)=1.$$ This turns the first half of the problem into a statement about exponent vectors rather than about the decimal value of $n$. 2) What the totient contributes Euler's product formula gives $$\varphi(n)=\prod_{i=1}^r p_i^{a_i-1}(p_i-1).$$ Fix a prime $q$. Its exponent inside $\varphi(n)$ is $$v_q(\varphi(n))=\sum_{i=1}^r v_q(p_i-1)+\sum_{\substack{1\le i\le r\\p_i=q}}(a_i-1).$$ So every exponent of $\varphi(n)$ comes from two sources: 1. factors already forced by the numbers $p_i-1$; 2. the self-contribution $a_i-1$ when the prime $q$ itself appears in $n$....

Detailed mathematical approach

Problem Summary

A positive integer $n$ is powerful if $p^2$ divides $n$ for every prime $p$ dividing $n$, and a perfect power if it is a power of another positive integer; an Achilles number is powerful but not a perfect power, as $864=2^5\cdot3^3$ and $1800=2^3\cdot3^2\cdot5^2$ are. It is a Strong Achilles number when $n$ and $\varphi(n)$ are both Achilles numbers, $\varphi$ being Euler's totient: $864$ qualifies, since $\varphi(864)=288=2^5\cdot3^2$, while $1800$ does not, since $\varphi(1800)=480=2^5\cdot3\cdot5$. There are $7$ such numbers below $10^4$ and $656$ below $10^8$; we must count those below $10^{18}$. ^[1]

The whole solution in one sentence: instead of testing numbers we build their prime factorisations, picking the primes of $n$ from the largest down, because each factor $p-1$ added to $\varphi(n)$ is made of smaller primes only — so a prime's exponent in $\varphi(n)$ is settled the instant the search drops below it, and both Achilles tests survive as two running gcds.

Mathematical Approach

1. Both conditions in the language of exponents

For $n\gt1$ write $n=\prod_{i=1}^{r}p_i^{a_i}$ with distinct primes $p_i$, $r\ge1$ and $a_i\ge1$. Powerful means precisely that every $a_i$ is at least $2$. ^[3] And if $n=m^k$ with $k\ge2$ then every exponent of $n$ is $k$ times the matching exponent of $m$, so $k$ divides $\gcd(a_1,\dots,a_r)$; conversely if $d=\gcd(a_1,\dots,a_r)\ge2$ then $n=\bigl(\prod_i p_i^{a_i/d}\bigr)^{d}$. So $n$ is a perfect power exactly when that gcd exceeds $1$, ^[4] and an Achilles number is described completely by

$$a_i\ge2\ \ (1\le i\le r),\qquad \gcd(a_1,\dots,a_r)=1.$$

Two things follow. One prime cannot do it, since the gcd is then $a_1\ge2$, so every Achilles number has at least two distinct prime factors; and the smallest is $72=2^3\cdot3^2$, the only one below $100$. ^[2]

2. The exponents of the totient

Euler's totient $\varphi(n)$ counts the integers from $1$ to $n$ whose gcd with $n$ is $1$; for $n=12$ these are $1,5,7,11$, so $\varphi(12)=4$. Exactly $p^{a-1}$ of the numbers up to $p^a$ are multiples of $p$, so $\varphi(p^a)=p^{a-1}(p-1)$; with multiplicativity on coprime factors that gives Euler's product formula ^[5]

$$\varphi(n)=\prod_{i=1}^{r}p_i^{a_i-1}(p_i-1).$$

Let $v_q(m)$ denote the exponent of the prime $q$ in $m$. Reading that formula factor by factor,

$$v_q(\varphi(n))=\sum_{i=1}^{r}v_q(p_i-1)+e_q,$$

where $e_q=a_j-1$ if $q=p_j$ is a prime of $n$ and $e_q=0$ otherwise. Every prime of $n$ donates the whole factorisation of $p_i-1$, plus a self-contribution $a_i-1$ at $p_i$ itself. By Section 1, $\varphi(n)$ is Achilles exactly when its nonzero exponents are all at least $2$ and their gcd is $1$. The smallest Strong Achilles number is $500=2^2\cdot5^3$, with $\varphi=200=2^3\cdot5^2$; the two ways the test on $\varphi$ can fail are both visible nearby:

$1800=2^3\cdot3^2\cdot5^2$, $\varphi=480=2^5\cdot3\cdot5$: here $v_3=v_5=1$, so $\varphi$ is not powerful;
$108=2^2\cdot3^3$, $\varphi=36=2^2\cdot3^2=6^2$: $\varphi$ is powerful but a perfect square.

3. The largest prime factor needs exponent at least three

Let $P$ be the largest prime factor of $n$. Every other prime $q$ of $n$ has $q\lt P$, and each prime divisor of $q-1$ is at most $q-1$, hence below $P$, so no $q-1$ holds a factor $P$. Section 2 then leaves $v_P(\varphi(n))=a_P-1\ge1$, and a powerful $\varphi(n)$ needs it to be at least $2$:

$$a_P\ge3.$$

With Section 1 — one more prime at least, with exponent at least $2$ — that gives $n\ge4P^3$, so for $n\le10^{18}$

$$P\le\left\lfloor\sqrt[3]{10^{18}/4}\right\rfloor=629960.$$

Every prime the search can need lies below $630000$, the size of the prime table in the code. The largest that occurs is $629857$: here $n=2^2\cdot629857^3$ equals $999\,507\,074\,182\,743\,172$ and $629857-1=2^5\cdot3^9$, so $\varphi(n)=2^6\cdot3^9\cdot629857^2$, exponents $6,9,2$ — all at least $2$, gcd $1$. Nothing larger works: ten primes lie between $629857$ and $629960$, each forced to exponent $3$, leaving a budget of $4$, so the one extra factor can only be $2^2$ — and each of the ten has an odd prime of multiplicity $1$ in $P-1$, which $2^2$ cannot repair.

4. Choosing the primes from the top down

The observation that organises the search is the one just used: $p-1$ is built from primes smaller than $p$. So pick the primes of $n$ in strictly decreasing order. Once the search is below a prime $q$, nothing can change $v_q(\varphi(n))$: a later prime $p'\lt q$ donates the factors of $p'-1\lt q$, none of which is $q$, and the self-contribution at $q$ would need $q$ itself, no longer allowed. That exponent is final. So the state stays small:

a budget $L$: if the chosen prime powers multiply to $M$, then $L=\lfloor10^{18}/M\rfloor$;
a multiset $A$ — repeated entries kept — of the primes below the last choice contributed by the factors $p_i-1$, so $12=2^2\cdot3$ adds $2,2,3$;
$g_n$, the gcd of the exponents chosen so far in $n$;
$g_\varphi$, the gcd of the already final exponents of $\varphi(n)$.

5. The exponent floor of the next prime

Let $p$ be the next prime, smaller than everything chosen, with exponent $a$, and let $t\ge0$ be its multiplicity in $A$. By Sections 2 and 4 the final value is $v_p(\varphi(n))=t+(a-1)$, so a powerful $\varphi(n)$ demands

$$t+(a-1)\ge2.$$

That splits the search in two. If $p$ is new, $t=0$ and $a\ge3$: a fresh prime costs at least $p^3$. If $p$ already lies in $A$, then $t\ge1$ and $a\ge2$ suffices — which Section 1 demands anyway, so it costs nothing extra.

Two examples make the difference concrete. In $1125=3^2\cdot5^3$ the prime $3$ is new, since $5-1=4$ donates no $3$, and $v_3(\varphi)=1$ spoils it; exponent $3$ gives the cube $3375=15^3$; the first exponent that works is $4$, and $10125=3^4\cdot5^3$, with $\varphi=2^3\cdot3^3\cdot5^2$, is the smallest odd Strong Achilles number. By contrast $12348=2^2\cdot3^2\cdot7^3$ manages with exponent $2$ twice, because $7-1=6$ and $3-1=2$ have put copies of $2$ and $3$ into $A$: indeed $\varphi(12348)=2^3\cdot3^2\cdot7^2$.

6. A prime seen once sets a lower bound

Suppose a prime $q$ occurs in $A$ with multiplicity exactly $1$. Its exponent in $\varphi(n)$ must still grow, in one of two ways: choose $q$ itself, or choose a prime $p\gt q$ with $q\mid p-1$. Both occur. After $13^3$ the factor $13-1=12$ leaves a single $3$ in $A$, and $7^4$ adds a second copy through $7-1=6$: $13^3\cdot7^4=5\,274\,997$ is a Strong Achilles number with $\varphi=2^3\cdot3^2\cdot7^3\cdot13^2$, although $3$ does not divide it. But if the next prime were below $q$, every later one would be too, and neither repair would remain. Hence:

If $A$ contains a prime of multiplicity $1$, the next prime must be at least the largest such prime.

That prunes hard, and it does more than save time. When the search steps down, every prime of $A$ above the new choice leaves the state for good — and by the bound each has multiplicity at least $2$, so its final exponent goes straight into $g_\varphi$, with no separate test anywhere in the descent.

7. When a partial choice is already an answer

Every node already describes a complete candidate: $n$ is the product of the prime powers chosen so far, and the recursion only decides whether to extend it. It is a Strong Achilles number exactly when $g_n=1$, so $n$ is no perfect power; every multiplicity in $A$ is at least $2$, those primes getting no self-contribution, so their exponent in $\varphi(n)$ is their multiplicity; and the gcd of $g_\varphi$ with those multiplicities is $1$, so $\varphi(n)$ is no perfect power either. One scan of $A$ and a few gcds decide it. The two gcd conditions prune nothing, though: one more exponent can pull a running gcd down to $1$, so the descent continues whatever the test says — the prefix $5^3$ has $g_n=3$, $g_\varphi=2$ and is not counted, yet extending it by $2^2$ gives $500$.

8. The whole search below ten thousand, by hand

With $L=10^4$, Section 3 allows $P^3\le L/4$, so $P\le13$. Take $P=5$ with $a=3$: the prefix is $125$, $5-1=4$ puts two copies of $2$ into $A$, $g_n=3$, and $v_5(\varphi)=2$ is already final, so $g_\varphi=2$; the budget is $80$. One next prime is the fresh prime $3$, where only exponent $3$ fits: that gives the cube $3375=15^3$, and a further $2^2$ would reach $13500$. The other is $2$, with $t=2$, so $a\ge2$ and $v_2(\varphi)=2+(a-1)$:

$a=2$: $n=500$, $\gcd(3,2)=1$, $v_2(\varphi)=3$, $\gcd(2,3)=1$ — counted;
$a=3$: $n=1000=10^3$, $\gcd(3,3)=3$ — a perfect cube;
$a=4$: $n=2000$, $v_2(\varphi)=5$, $\gcd(2,5)=1$ — counted;
$a=5$: $n=4000$, $v_2(\varphi)=6$, $\gcd(2,6)=2$; indeed $\varphi=2^6\cdot5^2=40^2$;
$a=6$: $n=8000=20^3$; $a=7$ already exceeds $10^4$.

The sweep $P=5,\ a=4$ (prefix $625$, $g_\varphi=3$) adds $5000$ and rejects $2500=50^2$ and $10000=10^4$; $a=5$ leaves no room. For $P=3$ write $n=3^a2^b$, so $\varphi(n)=3^{a-1}2^b$, and the conditions reduce to $a\ge3$, $b\ge2$, $\gcd(a,b)=1$, $\gcd(a-1,b)=1$. The budget permits $2\le b\le8$ at $a=3$, where $b=5,7$ pass and give $864$ and $3456$; $2\le b\le6$ at $a=4$, where only $b=5$ passes, giving $2592$; $2\le b\le5$ at $a=5$, where only $b=3$ passes, giving $1944$; $b\in\{2,3\}$ at $a=6$, both failing $\gcd(a,b)=1$ at $2916=54^2$ and $5832=18^3$; and $b=2$ at $a=7$, where $8748$ has totient $2916$. Larger $a$ leaves no room. For $P=7$ the set $A=\{2,3\}$ forces the next prime to be at least $3$, and the budget $29$ excludes the fresh prime $5$: $3^2\cdot7^3=3087$ has totient $42^2$, $3^3\cdot7^3=9261=21^3$, higher exponents of $3$ exceed the budget, and the smallest extension $2^2\cdot3^2\cdot7^3=12348$ is too large, as is $7^4\cdot3^2$. For $P=11$, $A=\{2,5\}$ forces $11^3\cdot5^2=33275$; for $P=13$, $A=\{2,2,3\}$ forces $13^3\cdot3^2=19773$; both are too large.

The complete list below $10^4$ is $500$, $864$, $1944$, $2000$, $2592$, $3456$, $5000$ — the $7$ numbers of the statement, and the first seven terms of the sequence of Strong Achilles numbers. ^[6]

9. Each number is produced exactly once

Because the primes decrease strictly, a set of prime powers can be assembled in one order only, so no candidate is reached twice. Nothing is missed either. Follow the primes of any Strong Achilles number $n\le10^{18}$ downwards: the largest exponent obeys Section 3, every later one Section 5, every prefix fits the budget, and whenever a singleton $q$ sits in $A$ the next prime of $n$ cannot be below $q$, since $v_q(\varphi(n))$ would then stay $1$. So Section 6 never prunes this path, and Section 7 counts $n$ once its last prime is taken.

10. Why the tree stays small, and brute force does not

Four bounds cut the recursion down: $p^3\le L$ for a new prime, $p^2\le L$ for an active one, the lower bound of Section 6, and the global $630000$ of Section 3. Against that, the count $k(x)$ of powerful numbers up to $x$ satisfies ^[3]

$$k(x)=\frac{\zeta(3/2)}{\zeta(3)}\sqrt{x}+O\!\left(x^{1/3}\right),$$

a result of Golomb, in which $\zeta(s)=\sum_{j\ge1}j^{-s}$ and $\zeta(3/2)/\zeta(3)\approx2.1733$. Listing all powerful numbers is a fine check at moderate bounds — there are $2\,158\,391$ up to $10^{12}$ — but below $10^{18}$ it would visit some $2.17\cdot10^9$ candidates, most still needing the exponents of their totient.

How the Code Works

C++ builds two tables once in prepare_data, called from main; Java does the same in prepareData from solve, Python in Solver.__init__, with solve_main starting the recursion. A sieve on the 2–3 wheel keeps only the residues $1$ and $5$ modulo $6$, which covers every prime above $3$; ^[7] for each sieving prime $p$ one chain starts at $p^2$ and the other at $p(p+4)$ if $p\equiv1\pmod6$, at $p(p+2)$ if $p\equiv5\pmod6$, both stepping by $6p$ so that each stays in one residue class. With kPrimeSieveN equal to $6\cdot105000=630000$, the bound of Section 3, that leaves $51341$ primes in primes. Then the index of each prime is appended to divisors[j] for every multiple $j$ of each of its powers: $2\,262\,483$ writes, after which divisors[m] is the ascending multiset of prime indices of $m$ with multiplicity and divisors_count[m] its length. Writing $\Omega(m)$ for the number of prime factors of $m$ with multiplicity, that product is at least $2^{\Omega(m)}$, so $\Omega(m)\le\lfloor\log_2m\rfloor\le19$ for $m\le630000$, with equality at $2^{19}$; the $24$ slots per row suffice.

count_recursive(max_prime_idx, active_factors, active_count, gcd_exp_n, gcd_exp_phi, lim) — countRecursive in Java — carries the state of Section 4: lim is $L$, active_factors is $A$ as a sorted index array, gcd_exp_n and gcd_exp_phi are $g_n$ and $g_\varphi$, and max_prime_idx caps the next index. The root passes an empty $A$, both gcds $0$ — since $\gcd(0,x)=x$, the first exponent initialises the running value — and $L=10^{18}$. The first block is the test of Section 7: with gcd_exp_n equal to $1$ it walks the equal runs of active_factors, stops at any run shorter than $2$, folds each run length into a local gcd and increments answer if that gcd ends at $1$.

Then the two index bounds. min_idx is the lower bound of Section 6, the largest index occurring exactly once in $A$; the scan inspects the top run and the interior positions, never the first, which it need not: the prime $2$ joins $A$ as soon as any odd prime is chosen, so a non-empty $A$ starts at the index of $2$, where the fallback $0$ is correct anyway. For max_idx, cube_root_floor takes std::pow as a guess and corrects it both ways with exact integer multiplications, so no floating-point error survives; its result bounds the next fresh prime but need not be prime, so a binary search over primes turns it into an index, and a clamp to max_prime_idx - 1 keeps the primes strictly decreasing.

The first loop extends the number by an active prime, the branch $t\ge1$ of Section 5, walking the runs of $A$ downwards. For a run of length overdeg at the prime $p$ it merges $A$ with divisors[p-1] into merged in one pass, pops every run above the chosen index into g2 — legitimate by Section 6 — then drops the overdeg copies of that index, which the self-contribution replaces. Exponents run deg = 2, 3, ... while pow <= lim / p, and the child receives gcd(gcd_exp_n, deg) and gcd(g2, deg - 1 + overdeg), which is $t+(a-1)$ exactly. The loop ends with if (overdeg == 1) break, the other half of Section 6. The second loop takes a brand-new prime instead: from min_idx to max_idx, skipping indices in $A$, merging and popping the same way, with exponents from deg = 3 and gcd(g2, deg - 1) for the child.

Both loops shrink the budget as lim / pow / p, not lim / (pow * p): two successive integer divisions give the same floor and cannot overflow on the way. C++ holds the budget and the prime powers in unsigned std::uint64_t and Java in signed long, which also reaches $10^{18}$; both leave an exponent loop when one more multiplication of pow would overflow that type, while Python needs no guard and keeps divisors as a list of lists. One further difference is a genuine defect: after the prime $2$ is chosen the child is entered with max_prime_idx equal to $0$, and C++ evaluates primes[max_prime_idx - 1] before the clamp that follows — an out-of-bounds read. Java guards it with maxPrimeIdx > 0, Python's negative index merely yields the last prime, and in all three the candidate range is clamped to empty just afterwards, so the counts are unaffected. Python's recursion limit of $20000$ is far more than needed: Sections 3 and 5 force $n\ge P^3\prod q^2$, allowing at most $8$ distinct primes below $10^{18}$, and the deepest path here selects $7$. Only the C++ program checks the statement's checkpoints, with assert. ^[8]

Complexity and Verification

Let $N=630000$. The wheel sieve costs $O(N\log\log N)$, and the divisor table makes one write per multiple of each prime power,

$$D_N=\sum_{p^k\le N}\left\lfloor\frac{N}{p^k}\right\rfloor=2\,262\,483,$$

whose $k=1$ terms are again $O(N\log\log N)$ while the higher powers add $O(N)$. The table reserves $630001\cdot24$ integer slots, $60\,480\,096$ bytes at four bytes each; Java adds a header per row and Python stores lists of references, so resident size is language-dependent.

The search has no closed-form cost, since the tree's shape depends on how the numbers $p-1$ factorise. Let $V$ count the calls, $C$ the candidate primes merged and $F$ bound the list lengths. A call scans $A$ and may run the binary search, $O(F+\log N)$; each candidate prime costs one merge, $O(F)$, reused for all its exponents; each child costs two gcds, on operands below $60$. So the search takes

$$O\bigl(V(F+\log N)+CF\bigr)$$

operations, with $O(DF)$ recursion storage at depth $D$. Measured, $V$ is $1468$ at $L=10^8$, $26\,671$ at $10^{12}$, $214\,384$ at $10^{15}$ and $1\,608\,314$ at $10^{18}$, the last traversal merging $C=756\,529$ candidate primes. On the machine used here the run takes $8$ ms of table building and $52$ ms of search in C++, about half a second in Java and $0.4+4.9$ s in Python.

The state arrays hold $70$ prime indices, and that is provably enough: every chosen prime has $a_i\ge2$, so $\Omega(p_i-1)\le\log_2 p_i\le\tfrac12 a_i\log_2 p_i$ and

$$\sum_i\Omega(p_i-1)\le\tfrac12\log_2 n\lt30,$$

while a merge adds at most $\Omega(p-1)\le16$ entries, the maximum over the primes below $630000$ and reached at the Fermat prime $65537$. So the merged length never passes $46$; an instrumented run reaches $19$ for the length of $A$ and $20$ after a merge.

Those assertions, $7$ below $10^4$ and $656$ below $10^8$, are active unless the program is compiled with NDEBUG; the three programs printed the same value at every bound tried here. One detail is worth stating: the recursion treats the bound as inclusive, as the counts $1$ for $L=500$ and $0$ for $L=499$ show, but $10^4$, $10^8$ and $10^{18}$ are perfect squares and hence never Achilles numbers, so "below" and "up to" agree.

The check written for this article shares the exponent tests of Section 1 but none of the state or the pruning: it enumerates every powerful number up to $x$ by a separate recursion over primes with exponents at least $2$, computes $\varphi$ from that factorisation and tests it directly, using neither $A$ nor the bound of Section 6. It agrees with the published solution at $x=10^4$, $10^5$, $10^6$, $10^7$, $10^8$, $10^{10}$ and $10^{12}$, with counts $7$, $34$, $99$, $260$, $656$, $3424$ and $16\,176$; its first twelve terms match the published sequence. ^[6] The largest exponent anywhere is $55$, in $2^{55}\cdot3^3=972\,777\,519\,512\,027\,136$ with totient $2^{55}\cdot3^2$: no prime above $2$ can carry exponent $56$ inside the bound, and beside $2^a$ a second prime is needed whose largest member gives at least $3^3$, so $n\ge27\cdot2^a$ while $27\cdot2^{56}\gt10^{18}$.

Footnotes and References

The references below supply the background facts used above. The exponent form of both conditions, the exponent floors, the bound on the next prime and the counts are derived here.

Project Euler 302 — Strong Achilles Numbers. The official statement defines powerful numbers, perfect powers and Achilles numbers, calls S a Strong Achilles number when S and φ(S) are both Achilles numbers, illustrates that with 864 and the failure of 1800, and gives 7 below 10⁴ and 656 below 10⁸ before asking for the count below 10¹⁸.
Wikipedia — Achilles number. An Achilles number is one that is powerful but not a perfect power; the article lists 72, 108, 200, 288, 392, 432 as the smallest, notes that 72 is the only one from 1 to 100, and introduces strong Achilles numbers, the smallest being 500 = 2²·5³ with totient 200 = 2³·5², and 864.
Wikipedia — Powerful number. The article defines a powerful number as one in which every prime factor has exponent at least 2, gives the representation m = a²b³ with b squarefree, and states Golomb's 1970 result: the count from 1 to x lies between c·√x − 3·∛x and c·√x, with c = ζ(3/2)/ζ(3) ≈ 2.173.
Wikipedia — Perfect power. A perfect power is a number of the form mᵏ with m > 1 and k > 1, and the article states the criterion of Section 1: n with exponents α₁, …, α_r is a perfect power if and only if their greatest common divisor exceeds 1.
Wikipedia — Euler's totient function. The article gives the product formula φ(n) = p₁^(k₁−1)(p₁ − 1) ⋯ p_r^(k_r−1)(p_r − 1), states φ(pᵏ) = p^(k−1)(p − 1) and records that φ is multiplicative on coprime arguments, which is what Section 2 reads the exponents off.
OEIS A194085 — Strong Achilles numbers. The sequence collects the Achilles numbers m for which φ(m) is also an Achilles number; its data begin 500, 864, 1944, 2000, 2592, 3456, 5000, 10125, 10368, 12348, 12500, 16875, with which the independent enumeration was compared term by term.
Wikipedia — Wheel factorization. The article explains the wheel with basis 2 and 3: all numbers coprime to 2 and 3 come from the set containing 1 and 5 by repeated additions of 6, leaving a third of all integers as candidates — the wheel the sieve runs on.
C++ — Euler302.cpp. The linked immutable revision contains prepare_data, cube_root_floor, count_recursive and the two assertions on the checkpoints. It is the implementation record for this article.

Problem 302 source code

C++

#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <iostream>

namespace {

using u64 = std::uint64_t;

constexpr int kPrimeSieveN = 6 * 105000;
int prime_count = 0;
int primes[kPrimeSieveN / 4];
unsigned char erat[kPrimeSieveN + 1];

int divisors[kPrimeSieveN + 1][24];
unsigned char divisors_count[kPrimeSieveN + 1];

u64 answer = 0;

int gcd_small(int a, int b) {
    while (b != 0) {
        const int c = a % b;
        a = b;
        b = c;
    }
    return a;
}

int cube_root_floor(const u64 n) {
    if (n <= 1ULL) {
        return static_cast<int>(n);
    }
    int r = static_cast<int>(std::pow(static_cast<long double>(n), 1.0L / 3.0L));
    while (static_cast<u64>(r + 1) * static_cast<u64>(r + 1) * static_cast<u64>(r + 1) <= n) ++r;
    while (static_cast<u64>(r) * static_cast<u64>(r) * static_cast<u64>(r) > n) --r;
    return r;
}

void count_recursive(const int max_prime_idx,
                     const int* active_factors,
                     const int active_count,
                     int gcd_exp_n,
                     int gcd_exp_phi,
                     const u64 lim) {
    if (lim == 0ULL) {
        return;
    }

    if (gcd_exp_n == 1) {
        int g = gcd_exp_phi;
        int i = 0;
        while (i < active_count) {
            int j = i + 1;
            while (j < active_count && active_factors[j] == active_factors[i]) ++j;
            const int cnt = j - i;
            if (cnt < 2) {
                g = 0;
                break;
            }
            g = gcd_small(g, cnt);
            i = j;
        }
        if (g == 1) {
            ++answer;
        }
    }

    int min_idx = 0;
    if (active_count >= 2) {
        if (active_factors[active_count - 1] != active_factors[active_count - 2]) {
            min_idx = active_factors[active_count - 1];
        } else {
            for (int k = active_count - 2; k >= 1; --k) {
                if (active_factors[k] != active_factors[k + 1] && active_factors[k] != active_factors[k - 1]) {
                    min_idx = active_factors[k];
                    break;
                }
            }
        }
    } else if (active_count == 1) {
        min_idx = active_factors[0];
    }

    int max_idx = cube_root_floor(lim);
    if (primes[max_prime_idx - 1] > max_idx) {
        int a = 0;
        int b = prime_count - 1;
        while (a < b - 1) {
            const int c = (a + b) / 2;
            if (primes[c] > max_idx) {
                b = c;
            } else {
                a = c;
            }
        }
        max_idx = a;
    }
    if (max_idx > max_prime_idx - 1) {
        max_idx = max_prime_idx - 1;
    }

    int merged[70];
    for (int j = active_count - 1; j >= 0;) {
        int overdeg = 0;
        for (int y = j; y >= 0 && active_factors[y] == active_factors[j]; --y) ++overdeg;

        if (active_factors[j] < max_prime_idx) {
            int k1 = 0;
            int k2 = 0;
            int merged_count = 0;
            const int idx = active_factors[j];
            const int p_minus_one = primes[idx] - 1;
            const int div_cnt = divisors_count[p_minus_one];

            while (k1 < active_count && k2 < div_cnt) {
                if (active_factors[k1] < divisors[p_minus_one][k2]) {
                    merged[merged_count++] = active_factors[k1++];
                } else {
                    merged[merged_count++] = divisors[p_minus_one][k2++];
                }
            }
            while (k2 < div_cnt) merged[merged_count++] = divisors[p_minus_one][k2++];
            while (k1 < active_count) merged[merged_count++] = active_factors[k1++];

            int g2 = gcd_exp_phi;
            while (merged_count > 0 && merged[merged_count - 1] > active_factors[j]) {
                int y = merged_count - 1;
                while (y >= 0 && merged[y] == merged[merged_count - 1]) --y;
                g2 = gcd_small(g2, merged_count - y - 1);
                merged_count = y + 1;
            }
            if (merged_count > 0 && merged[merged_count - 1] == active_factors[j]) {
                merged_count -= overdeg;
            }

            const u64 p = static_cast<u64>(primes[idx]);
            u64 pow = p;
            int deg = 2;
            while (pow <= lim / p) {
                const int new_g1 = gcd_small(gcd_exp_n, deg);
                const int new_g2 = gcd_small(g2, deg - 1 + overdeg);
                count_recursive(idx, merged, merged_count, new_g1, new_g2, lim / pow / p);

                if (pow > (std::numeric_limits<u64>::max() / p)) break;
                pow *= p;
                ++deg;
            }
        }

        j -= overdeg;
        if (overdeg == 1) break;
    }

    int now = 0;
    while (now < active_count && active_factors[now] < min_idx) ++now;

    for (int idx = min_idx; idx <= max_idx; ++idx) {
        if (now < active_count && idx == active_factors[now]) {
            while (now < active_count && active_factors[now] == idx) ++now;
            continue;
        }

        int k1 = 0;
        int k2 = 0;
        int merged_count = 0;
        const int p_minus_one = primes[idx] - 1;
        const int div_cnt = divisors_count[p_minus_one];
        while (k1 < active_count && k2 < div_cnt) {
            if (active_factors[k1] < divisors[p_minus_one][k2]) {
                merged[merged_count++] = active_factors[k1++];
            } else {
                merged[merged_count++] = divisors[p_minus_one][k2++];
            }
        }
        while (k2 < div_cnt) merged[merged_count++] = divisors[p_minus_one][k2++];
        while (k1 < active_count) merged[merged_count++] = active_factors[k1++];

        int g2 = gcd_exp_phi;
        while (merged_count > 0 && merged[merged_count - 1] > idx) {
            int y = merged_count - 1;
            while (y >= 0 && merged[y] == merged[merged_count - 1]) --y;
            g2 = gcd_small(g2, merged_count - y - 1);
            merged_count = y + 1;
        }

        const u64 p = static_cast<u64>(primes[idx]);
        u64 pow = p * p;
        int deg = 3;
        while (pow <= lim / p) {
            const int new_g1 = gcd_small(gcd_exp_n, deg);
            const int new_g2 = gcd_small(g2, deg - 1);
            count_recursive(idx, merged, merged_count, new_g1, new_g2, lim / pow / p);

            if (pow > (std::numeric_limits<u64>::max() / p)) break;
            pow *= p;
            ++deg;
        }
    }
}

void prepare_data() {
    primes[0] = 2;
    primes[1] = 3;
    prime_count = 2;
    erat[1] = 1;

    const int sq = static_cast<int>(std::sqrt(static_cast<long double>(kPrimeSieveN)));
    for (int i = 0; i <= sq; i += 6) {
        if (erat[i + 1] == 0) {
            const int p = i + 1;
            const int step = p * 6;
            for (int j = p * p; j <= kPrimeSieveN; j += step) erat[j] = 1;
            for (int j = p * (p + 4); j <= kPrimeSieveN; j += step) erat[j] = 1;
            primes[prime_count++] = p;
        }
        if (erat[i + 5] == 0) {
            const int p = i + 5;
            const int step = p * 6;
            for (int j = p * p; j <= kPrimeSieveN; j += step) erat[j] = 1;
            for (int j = p * (p + 2); j <= kPrimeSieveN; j += step) erat[j] = 1;
            primes[prime_count++] = p;
        }
    }

    const int start = sq - (sq % 6) + 6;
    for (int i = start; i < kPrimeSieveN; i += 6) {
        if (erat[i + 1] == 0) primes[prime_count++] = i + 1;
        if (erat[i + 5] == 0) primes[prime_count++] = i + 5;
    }

    for (int i = 0; i < prime_count; ++i) {
        for (u64 deg = static_cast<u64>(primes[i]); deg <= static_cast<u64>(kPrimeSieveN);) {
            for (int j = static_cast<int>(deg); j <= kPrimeSieveN; j += static_cast<int>(deg)) {
                divisors[j][divisors_count[j]++] = i;
            }
            if (deg <= static_cast<u64>(kPrimeSieveN / primes[i])) {
                deg *= static_cast<u64>(primes[i]);
            } else {
                break;
            }
        }
    }
}

u64 solve(const u64 limit) {
    answer = 0;
    int root_factors[70] = {0};
    count_recursive(prime_count, root_factors, 0, 0, 0, limit);
    return answer;
}

}  // namespace

int main() {
    prepare_data();
    assert(solve(10'000ULL) == 7ULL);
    assert(solve(100'000'000ULL) == 656ULL);
    std::cout << solve(1'000'000'000'000'000'000ULL) << '\n';
    return 0;
}

Python

import sys

sys.setrecursionlimit(20000)

class Solver:
    def __init__(self):
        self.kPrimeSieveN = 6 * 105000
        self.primes = []
        
        erat = [0] * (self.kPrimeSieveN + 1)
        erat[1] = 1
        
        self.primes.append(2)
        self.primes.append(3)
        
        sq = int(self.kPrimeSieveN ** 0.5)
        for i in range(0, sq + 1, 6):
            if not erat[i + 1]:
                p = i + 1
                step = p * 6
                for j in range(p * p, self.kPrimeSieveN + 1, step): erat[j] = 1
                for j in range(p * (p + 4), self.kPrimeSieveN + 1, step): erat[j] = 1
                self.primes.append(p)
            if not erat[i + 5]:
                p = i + 5
                step = p * 6
                for j in range(p * p, self.kPrimeSieveN + 1, step): erat[j] = 1
                for j in range(p * (p + 2), self.kPrimeSieveN + 1, step): erat[j] = 1
                self.primes.append(p)
                
        start = sq - (sq % 6) + 6
        for i in range(start, self.kPrimeSieveN, 6):
            if not erat[i + 1]: self.primes.append(i + 1)
            if not erat[i + 5]: self.primes.append(i + 5)
            
        self.prime_count = len(self.primes)
        
        self.divisors = [[] for _ in range(self.kPrimeSieveN + 1)]
        for i, p in enumerate(self.primes):
            deg = p
            while deg <= self.kPrimeSieveN:
                for j in range(deg, self.kPrimeSieveN + 1, deg):
                    self.divisors[j].append(i)
                if deg <= self.kPrimeSieveN // p:
                    deg *= p
                else:
                    break
                    
        self.answer = 0

    def gcd_small(self, a, b):
        while b != 0:
            a, b = b, a % b
        return a

    def cube_root_floor(self, n):
        if n <= 1: return int(n)
        r = int(n ** (1.0 / 3.0))
        while (r + 1) * (r + 1) * (r + 1) <= n: r += 1
        while r * r * r > n: r -= 1
        return r

    def count_recursive(self, max_prime_idx, active_factors, active_count, gcd_exp_n, gcd_exp_phi, lim):
        if lim == 0: return

        if gcd_exp_n == 1:
            g = gcd_exp_phi
            i = 0
            while i < active_count:
                j = i + 1
                while j < active_count and active_factors[j] == active_factors[i]:
                    j += 1
                cnt = j - i
                if cnt < 2:
                    g = 0
                    break
                g = self.gcd_small(g, cnt)
                i = j
            if g == 1:
                self.answer += 1

        min_idx = 0
        if active_count >= 2:
            if active_factors[active_count - 1] != active_factors[active_count - 2]:
                min_idx = active_factors[active_count - 1]
            else:
                for k in range(active_count - 2, 0, -1):
                    if active_factors[k] != active_factors[k + 1] and active_factors[k] != active_factors[k - 1]:
                        min_idx = active_factors[k]
                        break
        elif active_count == 1:
            min_idx = active_factors[0]

        max_idx = self.cube_root_floor(lim)
        if self.primes[max_prime_idx - 1] > max_idx:
            a = 0
            b = self.prime_count - 1
            while a < b - 1:
                c = (a + b) // 2
                if self.primes[c] > max_idx:
                    b = c
                else:
                    a = c
            max_idx = a
        if max_idx > max_prime_idx - 1:
            max_idx = max_prime_idx - 1

        merged = [0] * 70
        j = active_count - 1
        while j >= 0:
            overdeg = 0
            y = j
            while y >= 0 and active_factors[y] == active_factors[j]:
                overdeg += 1
                y -= 1

            if active_factors[j] < max_prime_idx:
                k1 = 0
                k2 = 0
                merged_count = 0
                idx = active_factors[j]
                p_minus_one = self.primes[idx] - 1
                divs = self.divisors[p_minus_one]
                div_cnt = len(divs)

                while k1 < active_count and k2 < div_cnt:
                    if active_factors[k1] < divs[k2]:
                        merged[merged_count] = active_factors[k1]
                        k1 += 1
                    else:
                        merged[merged_count] = divs[k2]
                        k2 += 1
                    merged_count += 1
                while k2 < div_cnt:
                    merged[merged_count] = divs[k2]
                    k2 += 1
                    merged_count += 1
                while k1 < active_count:
                    merged[merged_count] = active_factors[k1]
                    k1 += 1
                    merged_count += 1

                g2 = gcd_exp_phi
                while merged_count > 0 and merged[merged_count - 1] > active_factors[j]:
                    y = merged_count - 1
                    while y >= 0 and merged[y] == merged[merged_count - 1]:
                        y -= 1
                    g2 = self.gcd_small(g2, merged_count - y - 1)
                    merged_count = y + 1
                if merged_count > 0 and merged[merged_count - 1] == active_factors[j]:
                    merged_count -= overdeg

                p = self.primes[idx]
                pow_val = p
                deg = 2
                while pow_val <= lim // p:
                    new_g1 = self.gcd_small(gcd_exp_n, deg)
                    new_g2 = self.gcd_small(g2, deg - 1 + overdeg)
                    self.count_recursive(idx, merged, merged_count, new_g1, new_g2, lim // pow_val // p)

                    pow_val *= p
                    deg += 1

            j -= overdeg
            if overdeg == 1: break

        now = 0
        while now < active_count and active_factors[now] < min_idx:
            now += 1

        for idx in range(min_idx, max_idx + 1):
            if now < active_count and idx == active_factors[now]:
                while now < active_count and active_factors[now] == idx:
                    now += 1
                continue

            k1 = 0
            k2 = 0
            merged_count = 0
            p_minus_one = self.primes[idx] - 1
            divs = self.divisors[p_minus_one]
            div_cnt = len(divs)

            while k1 < active_count and k2 < div_cnt:
                if active_factors[k1] < divs[k2]:
                    merged[merged_count] = active_factors[k1]
                    k1 += 1
                else:
                    merged[merged_count] = divs[k2]
                    k2 += 1
                merged_count += 1
            while k2 < div_cnt:
                merged[merged_count] = divs[k2]
                k2 += 1
                merged_count += 1
            while k1 < active_count:
                merged[merged_count] = active_factors[k1]
                k1 += 1
                merged_count += 1

            g2 = gcd_exp_phi
            while merged_count > 0 and merged[merged_count - 1] > idx:
                y = merged_count - 1
                while y >= 0 and merged[y] == merged[merged_count - 1]:
                    y -= 1
                g2 = self.gcd_small(g2, merged_count - y - 1)
                merged_count = y + 1

            p = self.primes[idx]
            pow_val = p * p
            deg = 3
            while pow_val <= lim // p:
                new_g1 = self.gcd_small(gcd_exp_n, deg)
                new_g2 = self.gcd_small(g2, deg - 1)
                self.count_recursive(idx, merged, merged_count, new_g1, new_g2, lim // pow_val // p)

                pow_val *= p
                deg += 1

    def solve_main(self, limit):
        self.answer = 0
        root_factors = [0] * 70
        self.count_recursive(self.prime_count, root_factors, 0, 0, 0, limit)
        return self.answer

def solve(limit=10**18):
    solver = Solver()
    return str(solver.solve_main(limit))

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

Java

public class Euler302 {
    static final int kPrimeSieveN = 6 * 105000;
    static int primeCount = 0;
    static int[] primes = new int[kPrimeSieveN / 4];
    static byte[] erat = new byte[kPrimeSieveN + 1];

    static int[][] divisors = new int[kPrimeSieveN + 1][24];
    static byte[] divisorsCount = new byte[kPrimeSieveN + 1];

    static long answer = 0;

    static int gcdSmall(int a, int b) {
        while (b != 0) {
            int c = a % b;
            a = b;
            b = c;
        }
        return a;
    }

    static int cubeRootFloor(long n) {
        if (n <= 1L) {
            return (int) n;
        }
        int r = (int) Math.pow(n, 1.0 / 3.0);
        while ((long) (r + 1) * (r + 1) * (r + 1) <= n)
            ++r;
        while ((long) r * r * r > n)
            --r;
        return r;
    }

    static void countRecursive(int maxPrimeIdx, int[] activeFactors, int activeCount, int gcdExpN, int gcdExpPhi,
            long lim) {
        if (lim == 0L)
            return;

        if (gcdExpN == 1) {
            int g = gcdExpPhi;
            int i = 0;
            while (i < activeCount) {
                int j = i + 1;
                while (j < activeCount && activeFactors[j] == activeFactors[i])
                    ++j;
                int cnt = j - i;
                if (cnt < 2) {
                    g = 0;
                    break;
                }
                g = gcdSmall(g, cnt);
                i = j;
            }
            if (g == 1) {
                ++answer;
            }
        }

        int minIdx = 0;
        if (activeCount >= 2) {
            if (activeFactors[activeCount - 1] != activeFactors[activeCount - 2]) {
                minIdx = activeFactors[activeCount - 1];
            } else {
                for (int k = activeCount - 2; k >= 1; --k) {
                    if (activeFactors[k] != activeFactors[k + 1] && activeFactors[k] != activeFactors[k - 1]) {
                        minIdx = activeFactors[k];
                        break;
                    }
                }
            }
        } else if (activeCount == 1) {
            minIdx = activeFactors[0];
        }

        int maxIdx = cubeRootFloor(lim);
        if (maxPrimeIdx > 0 && primes[maxPrimeIdx - 1] > maxIdx) {
            int a = 0;
            int b = maxPrimeIdx - 1;
            while (a < b - 1) {
                int c = (a + b) / 2;
                if (primes[c] > maxIdx) {
                    b = c;
                } else {
                    a = c;
                }
            }
            maxIdx = a;
        }
        if (maxIdx > maxPrimeIdx - 1) {
            maxIdx = maxPrimeIdx - 1;
        }

        int[] merged = new int[70];
        for (int j = activeCount - 1; j >= 0;) {
            int overdeg = 0;
            for (int y = j; y >= 0 && activeFactors[y] == activeFactors[j]; --y)
                ++overdeg;

            if (activeFactors[j] < maxPrimeIdx) {
                int k1 = 0, k2 = 0, mergedCount = 0;
                int idx = activeFactors[j];
                int pMinusOne = primes[idx] - 1;
                int divCnt = divisorsCount[pMinusOne];

                while (k1 < activeCount && k2 < divCnt) {
                    if (activeFactors[k1] < divisors[pMinusOne][k2]) {
                        merged[mergedCount++] = activeFactors[k1++];
                    } else {
                        merged[mergedCount++] = divisors[pMinusOne][k2++];
                    }
                }
                while (k2 < divCnt)
                    merged[mergedCount++] = divisors[pMinusOne][k2++];
                while (k1 < activeCount)
                    merged[mergedCount++] = activeFactors[k1++];

                int g2 = gcdExpPhi;
                while (mergedCount > 0 && merged[mergedCount - 1] > activeFactors[j]) {
                    int y = mergedCount - 1;
                    while (y >= 0 && merged[y] == merged[mergedCount - 1])
                        --y;
                    g2 = gcdSmall(g2, mergedCount - y - 1);
                    mergedCount = y + 1;
                }
                if (mergedCount > 0 && merged[mergedCount - 1] == activeFactors[j]) {
                    mergedCount -= overdeg;
                }

                long p = primes[idx];
                long pow = p;
                int deg = 2;
                while (pow <= lim / p) {
                    int newG1 = gcdSmall(gcdExpN, deg);
                    int newG2 = gcdSmall(g2, deg - 1 + overdeg);
                    countRecursive(idx, merged, mergedCount, newG1, newG2, lim / pow / p);

                    if (pow > (Long.MAX_VALUE / p))
                        break;
                    pow *= p;
                    ++deg;
                }
            }

            j -= overdeg;
            if (overdeg == 1)
                break;
        }

        int now = 0;
        while (now < activeCount && activeFactors[now] < minIdx)
            ++now;

        for (int idx = minIdx; idx <= maxIdx; ++idx) {
            if (now < activeCount && idx == activeFactors[now]) {
                while (now < activeCount && activeFactors[now] == idx)
                    ++now;
                continue;
            }

            int k1 = 0, k2 = 0, mergedCount = 0;
            int pMinusOne = primes[idx] - 1;
            int divCnt = divisorsCount[pMinusOne];

            while (k1 < activeCount && k2 < divCnt) {
                if (activeFactors[k1] < divisors[pMinusOne][k2]) {
                    merged[mergedCount++] = activeFactors[k1++];
                } else {
                    merged[mergedCount++] = divisors[pMinusOne][k2++];
                }
            }
            while (k2 < divCnt)
                merged[mergedCount++] = divisors[pMinusOne][k2++];
            while (k1 < activeCount)
                merged[mergedCount++] = activeFactors[k1++];

            int g2 = gcdExpPhi;
            while (mergedCount > 0 && merged[mergedCount - 1] > idx) {
                int y = mergedCount - 1;
                while (y >= 0 && merged[y] == merged[mergedCount - 1])
                    --y;
                g2 = gcdSmall(g2, mergedCount - y - 1);
                mergedCount = y + 1;
            }

            long p = primes[idx];
            long pow = p * p;
            int deg = 3;
            while (pow <= lim / p) {
                int newG1 = gcdSmall(gcdExpN, deg);
                int newG2 = gcdSmall(g2, deg - 1);
                countRecursive(idx, merged, mergedCount, newG1, newG2, lim / pow / p);

                if (pow > (Long.MAX_VALUE / p))
                    break;
                pow *= p;
                ++deg;
            }
        }
    }

    static void prepareData() {
        primes[0] = 2;
        primes[1] = 3;
        primeCount = 2;
        erat[1] = 1;

        int sq = (int) Math.sqrt(kPrimeSieveN);
        for (int i = 0; i <= sq; i += 6) {
            if (erat[i + 1] == 0) {
                int p = i + 1;
                int step = p * 6;
                for (int j = p * p; j <= kPrimeSieveN; j += step)
                    erat[j] = 1;
                for (int j = p * (p + 4); j <= kPrimeSieveN; j += step)
                    erat[j] = 1;
                primes[primeCount++] = p;
            }
            if (erat[i + 5] == 0) {
                int p = i + 5;
                int step = p * 6;
                for (int j = p * p; j <= kPrimeSieveN; j += step)
                    erat[j] = 1;
                for (int j = p * (p + 2); j <= kPrimeSieveN; j += step)
                    erat[j] = 1;
                primes[primeCount++] = p;
            }
        }

        int start = sq - (sq % 6) + 6;
        for (int i = start; i < kPrimeSieveN; i += 6) {
            if (erat[i + 1] == 0)
                primes[primeCount++] = i + 1;
            if (erat[i + 5] == 0)
                primes[primeCount++] = i + 5;
        }

        for (int i = 0; i < primeCount; ++i) {
            for (long deg = primes[i]; deg <= kPrimeSieveN;) {
                for (int j = (int) deg; j <= kPrimeSieveN; j += (int) deg) {
                    divisors[j][divisorsCount[j]++] = i;
                }
                if (deg <= kPrimeSieveN / primes[i]) {
                    deg *= primes[i];
                } else {
                    break;
                }
            }
        }
    }

    public static String solve() {
        prepareData();
        answer = 0;
        int[] rootFactors = new int[70];
        countRecursive(primeCount, rootFactors, 0, 0, 0, 1000000000000000000L);
        return String.valueOf(answer);
    }

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