Problem 362: Squarefree Factors

Project Euler Problem 362 Solution

EulerSolve provides an optimized solution for Project Euler Problem 362, Squarefree Factors, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For each integer $n \ge 2$, let $\operatorname{fsf}(n)$ denote the number of unordered factorizations of $n$ into squarefree factors greater than 1. Repeated factors are allowed as long as each factor itself is squarefree. The Project Euler quantity is $$S(N)=\sum_{n=2}^{N}\operatorname{fsf}(n),\qquad N=10^{10}.$$ As a small example, $12=2^2\cdot 3$ has exactly two valid factorizations: $$12=2\cdot 2\cdot 3=2\cdot 6,$$ so $\operatorname{fsf}(12)=2$. A direct scan up to $10^{10}$ is impossible, so the implementation groups integers by their exponent pattern and counts entire families at once. Mathematical Approach Write the prime factorization of $n$ as $$n=\prod_{i=1}^{r} p_i^{e_i},\qquad e_i\ge 1.$$ The multiset of exponents $(e_1,\dots,e_r)$ is the only part that matters for the squarefree-factor count. Step 1: Replace Unordered Factorizations by a System on Subsets Every squarefree factor of $n$ is obtained by choosing a nonempty subset $J\subseteq\{1,\dots,r\}$ and taking $$f_J=\prod_{i\in J} p_i.$$ Because the factorization is unordered, all we need is the multiplicity $x_J\in\mathbb{Z}_{\ge 0}$ of each subset type....

Detailed mathematical approach

Problem Summary

The number $54$ is a product of one or more factors larger than $1$ in seven ways: $54$, $2\times27$, $3\times18$, $6\times9$, $3\times3\times6$, $2\times3\times9$ and $2\times3\times3\times3$. ^[4] A number is squarefree when no prime divides it twice, so $27=3^3$, $18=2\cdot3^2$ and $9=3^2$ are not, and only $3\times3\times6$ and $2\times3\times3\times3$ are squarefree throughout. Writing $\operatorname{Fsf}(n)$ for the number of factorisations of $n$ into one or more squarefree factors larger than $1$, this says $\operatorname{Fsf}(54)=2$. The order of the factors never matters and a factor may repeat. With

$$S(n)=\sum_{k=2}^{n}\operatorname{Fsf}(k),$$

the statement gives $S(100)=193$ and asks for $S(N)$ at $N=10^{10}$. ^[1]

The whole solution in one sentence: $\operatorname{Fsf}(n)$ depends only on the multiset of exponents in the prime factorisation of $n$, so instead of visiting $10^{10}$ integers we list the $1961$ exponent patterns occurring below the limit, count the factorisations once per pattern, and multiply by how many integers carry that pattern.

Mathematical Approach

1. What a squarefree factor can be

A squarefree number is divisible by no square other than $1$, equivalently no prime occurs in it with exponent $2$ or more, so it is a product of distinct primes. ^[2] Every factor of $n$ divides $n$, and if

$$n=\prod_{i=1}^{r}p_i^{e_i},\qquad e_i\ge1,$$

then every squarefree divisor of $n$ larger than $1$ has the form

$$f_J=\prod_{i\in J}p_i$$

for a non-empty $J\subseteq\{1,\dots,r\}$, and distinct subsets give distinct divisors. Both $12=2^2\cdot3$ and $54=2\cdot3^3$ therefore offer the three factors $2$, $3$, $6$.

2. A factorisation is a multiset of subsets

Since the order is irrelevant, a factorisation is determined by how often each $f_J$ is used; let $x_J\ge0$ be that multiplicity. The prime $p_i$ occurs once in every chosen factor whose subset contains $i$, and $e_i$ times in all, so

$$\sum_{J\ni i}x_J=e_i\qquad(1\le i\le r).$$

Conversely every non-negative solution gives a factorisation: all factors exceed $1$, and at least one is used because $e_1\ge1$. Hence:

Claim 1. $\operatorname{Fsf}(n)$ is the number of non-negative integer solutions $(x_J)$ of that system, $J$ running over the $2^r-1$ non-empty subsets of $\{1,\dots,r\}$.

For $12=2^2\cdot3$ the system reads $x_1+x_{12}=2$, $x_2+x_{12}=1$, so only $x_{12}=0$ and $x_{12}=1$ are possible, giving $(2,1,0)$, that is $2\cdot2\cdot3$, and $(1,0,1)$, that is $2\cdot6$: $\operatorname{Fsf}(12)=2$. For $54=2\cdot3^3$ the system $x_1+x_{12}=1$, $x_2+x_{12}=3$ likewise has the two solutions $2\cdot3\cdot3\cdot3$ and $3\cdot3\cdot6$, the products the statement keeps.

3. Only the exponents matter

The system above contains the exponents but not the primes themselves. Renaming the primes permutes the coordinates $1,\dots,r$ and carries solutions bijectively to solutions, so their number is unchanged. Hence $\operatorname{Fsf}(n)$ depends only on the multiset $\{e_1,\dots,e_r\}$, the prime signature of $n$. We write it as a non-increasing tuple $\sigma=(e_1,\dots,e_r)$, $e_1\ge\dots\ge e_r\ge1$, and call the common value $F(\sigma)$. ^[3] This is the whole saving: $10^{10}$ integers, but only $1961$ signatures.

4. Three signatures we can settle by hand

A prime power. Here $r=1$ and the only equation is $x_1=e_1$, so $F(e_1)=1$: the one factorisation of $p^{e}$ is $p\cdot p\cdots p$. ^[3]
Two primes. The subsets are $\{1\}$, $\{2\}$, $\{1,2\}$, and $x_{12}$ already fixes $x_1=e_1-x_{12}$ and $x_2=e_2-x_{12}$. Both must stay non-negative, so $x_{12}$ runs from $0$ to $\min(e_1,e_2)$ and $F(e_1,e_2)=\min(e_1,e_2)+1$. Thus $F(2,1)=F(3,1)=2$, which is $\operatorname{Fsf}(12)$ and $\operatorname{Fsf}(54)$ again, while $F(3,3)=4$.
A squarefree number. Every $e_i=1$, so every $x_J$ is $0$ or $1$ and the subsets with $x_J=1$ cover each index once: they form a partition of $\{1,\dots,r\}$. So $F(1,\dots,1)$ is the Bell number $B_r$, the number of partitions of an $r$-element set: $1,2,5,15,52,203$ for $r=1,\dots,6$. ^[6] Hence $\operatorname{Fsf}(30)=5$, namely $30$, $2\cdot15$, $3\cdot10$, $5\cdot6$, $2\cdot3\cdot5$; at $r=10$ it is $B_{10}=115975$.

5. Counting the solutions of a general signature

Fix any order $J_1,\dots,J_{2^r-1}$ of the non-empty subsets and choose $x_{J_1}$, then $x_{J_2}$, and so on. After the first $t-1$ choices let $\mathrm{rem}_i$ be $e_i$ minus what has been used. Then $x_{J_t}$ cannot exceed

$$\min_{i\in J_t}\mathrm{rem}_i,$$

since a larger value would drive some coordinate negative, and no later choice can raise it again. Once all subsets are assigned, the vector is a solution exactly when every $\mathrm{rem}_i$ is $0$. Every $x_J$ is chosen exactly once, so each solution appears exactly once: the enumeration is a bijection onto the solution set, whatever order we pick.

Two observations make it fast. After step $t$ only the pair $\bigl(t,\mathrm{rem}\bigr)$ matters, so equal states merge and their sub-count is reused. And taking large subsets first puts the strong constraints in front: a subset of size $k$ lowers $k$ coordinates at once, so the branching shrinks early. The hardest pattern below our limit is ten distinct primes each to the first power, with $1023$ subsets.

6. Regrouping the sum by signature

Every $n\ge2$ has one prime signature, so the classes are disjoint and cover $\{2,\dots,N\}$. Sorting the sum by class gives

$$S(N)=\sum_{\sigma}F(\sigma)\,A(\sigma,N),$$

where $A(\sigma,N)$ counts the integers $2\le n\le N$ of signature $\sigma$. The two factors are now independent problems: $F(\sigma)$ is combinatorics on exponents, $A(\sigma,N)$ is counting integers.

7. Which signatures occur below the limit

The smallest integer of signature $\sigma$ is

$$m(\sigma)=2^{e_1}3^{e_2}\cdots p_r^{e_r},$$

the largest exponent on the smallest prime. Two exchange arguments prove it: replacing any prime by a smaller unused one lowers the product, so the primes must be $2,3,\dots,p_r$; and if $p\lt q$ carry exponents $a\gt b$, then $p^{a-b}\le q^{a-b}$ gives $p^aq^b\le p^bq^a$, so sorting the exponents non-increasingly never increases the product. These representatives are exactly the products of primorials, a primorial being a product $2\cdot3\cdots p_j$ of the first $j$ primes. ^[5]

Since $m(\sigma)$ itself has signature $\sigma$, a signature occurs below $N$ exactly when $m(\sigma)\le N$. The patterns can therefore be grown prime by prime: extend a partial pattern by an exponent no larger than the previous one while the product stays $\le N$, emitting every prefix. At $N=10^{10}$ this yields $1961$ patterns, constrained by two products: $2\cdot3\cdot5\cdots29=6469693230\le10^{10}$ but $2\cdot3\cdots31=200560490130\gt10^{10}$, so $r\le10$; and $2^{33}=8589934592\le10^{10}\lt2^{34}$, so $e_1\le33$.

8. Counting the integers of one signature

A signature is unordered, but an integer assigns its exponents to its primes in increasing order. Splitting $A(\sigma,N)$ over the distinct orderings $a=(a_1,\dots,a_r)$ of the multiset $\sigma$,

$$A(\sigma,N)=\sum_{a}C(a,N),$$

$$\begin{aligned}C(a,N)=\#\{q_1\lt\dots\lt q_r\ :\\ \prod_{i=1}^{r}q_i^{a_i}\le N\},\end{aligned}$$

the $q_i$ being primes. Distinct orderings rather than all $r!$ of them is what keeps this exact: a repeated exponent would otherwise let the same integer appear twice. At $N=10^{10}$ the $1961$ signatures give $53558$ ordered tuples.

Each $C(a,N)$ is computed by choosing the primes from the left. With $q_1,\dots,q_{k-1}$ fixed, let $R=\lfloor N/\prod_{i\lt k}q_i^{a_i}\rfloor$; the floor costs nothing, the remaining product being an integer. Every later prime is at least $q_k$, hence

$$\begin{aligned}q_k^{T_k}\le\prod_{i\ge k}q_i^{a_i}\le R,\\ T_k=a_k+\dots+a_r,\end{aligned}$$

so $q_k\le\lfloor R^{1/T_k}\rfloor$. The bound is safe but not tight, since it lets the later primes equal $q_k$: with $R=25$ and two remaining exponents $1$ it still admits $q_k=5$, although the next prime is already $7$ and $5\cdot7\gt25$. Such an over-admitted candidate simply contributes $0$ one level down. At the last coordinate nothing branches: the admissible $q_r$ are the primes with $q_r^{a_r}\le R$ above $q_{r-1}$, and their number is $\pi\bigl(\lfloor R^{1/a_r}\rfloor\bigr)$ minus the index already reached, so one evaluation of $\pi$ settles the whole level. Different prefixes often leave the same level, first admissible prime index and $R$, so these sub-counts are reused too.

9. Evaluating the prime-counting function

Before the last coordinate at least two exponents remain, so $T_k\ge2$ and $q_k\le\lfloor\sqrt{R}\rfloor\le10^{5}$. An ordinary sieve up to $\lfloor\sqrt N\rfloor=10^{5}$ therefore supplies every prime the recursion can pick, all $9592$ of them, and $\pi(x)$ for $x\le10^{5}$. Only the last level asks for $\pi$ at larger arguments, up to $N$ itself.

Those are answered by Lehmer's exact formula. ^[7] Let $\varphi(x,a)$ count the integers up to $x$ divisible by none of the first $a$ primes; for instance $\varphi(10,2)=3$, counting $1,5,7$. Among the survivors of the first $a-1$ primes, those divisible by $p_a$ become, after division by $p_a$, the survivors up to $x/p_a$, so

$$\varphi(x,a)=\varphi(x,a-1)-\varphi(x/p_a,a-1).$$

For $a\le7$ the pattern repeats modulo $W=2\cdot3\cdot5\cdot7\cdot11\cdot13\cdot17$${}=510510$: writing $x=tW+u$ with $0\le u\lt W$ gives $t\varphi(W,a)+\varphi(u,a)$, both from a table. Now take $a=\pi(x^{1/4})$. A survivor's prime factors all exceed $x^{1/4}$, so it has at most three of them: it is $1$, a prime, or a product of two or three primes. With $P_2$ and $P_3$ counting the last two classes,

$$\pi(x)=\varphi(x,a)+a-1-P_2-P_3,$$

and $P_2$, $P_3$ are sums of $\pi$ at smaller arguments, which is what the two loops of pi subtract. No sieve beyond $\sqrt x$ is needed. The later Lagarias–Miller–Odlyzko refinement, with its $O(x^{2/3+\varepsilon})$ bound, is a different algorithm and does not describe this code. ^[8] For reference $\pi(10^{10})=455052511$. ^[7]

10. The worked example: $S(100)=193$

Fifteen signatures occur below $100$. Each line gives $\sigma$, then $F(\sigma)$, then the integers of that signature, then the product.

$(1)$, $F=1$: the $25$ primes, giving $25$.
$(1,1)$, $F=2$: $30$ products of two distinct primes, giving $60$.
$(1,1,1)$, $F=5$: $30,42,66,70,78$, giving $25$.
$(2)$, $F=1$: $4,9,25,49$, giving $4$.
$(2,1)$, $F=2$: $15$ numbers, giving $30$.
$(2,1,1)$, $F=6$: $60,84,90$, giving $18$. For $60$ the six are $2\cdot2\cdot3\cdot5$, $2\cdot2\cdot15$, $2\cdot3\cdot10$, $2\cdot5\cdot6$, $2\cdot30$ and $6\cdot10$.
$(2,2)$, $F=3$: $36,100$, giving $6$.
$(3)$, $F=1$: $8,27$, giving $2$.
$(3,1)$, $F=2$: $24,40,54,56,88$, giving $10$.
$(3,2)$, $F=3$: $72$, giving $3$.
$(4)$, $F=1$: $16,81$, giving $2$.
$(4,1)$, $F=2$: $48,80$, giving $4$.
$(5)$, $F=1$: $32$, giving $1$.
$(5,1)$, $F=2$: $96$, giving $2$.
$(6)$, $F=1$: $64$, giving $1$.

The products add to $193$ and the counts to $99$, the number of integers from $2$ to $100$, confirming that the classes partition the range. The line $(2,1)$ shows Section 8 in miniature: the ordering $(2,1)$ counts the $11$ numbers $p^2q$ with $p\lt q$, namely $12,20,28,44,45,52,$${}63,68,76,92,99$, and $(1,2)$ the $4$ numbers $pq^2$, namely $18,50,75,98$. And $54$ sits in the $(3,1)$ line through the ordering $(1,3)$.

11. Why the total is exact

Every link is an exact identity: Claim 1, the invariance of Section 3, the bijective sweep of Section 5, the partition of Section 6, and the disjoint orderings and once-only prefix recursion of Section 8. Nothing is rounded, and $\pi(x)$ is exact.

How the Code Works

The C++ file holds the full solver. ^[9] isqrt_u64, icbrt_u64, iroot4_u64 and kth_root_floor start from a floating-point root and walk up or down to the exact integer root, and pow_u64_capped returns cap + 1 instead of overflowing: the bounds of Section 8 are floors of real roots, so an off-by-one would silently change a count. PrimeTable runs the linear sieve up to $\max(\lfloor\sqrt N\rfloor,100)$, keeping the prime list and the table pi_small. PrimeCounter implements Section 9: build_phi_table tabulates $\varphi(x,a)$ for $a\le$ kSmallPhiPrimes${}=7$ over one period kSmallPhiMod${}=510510$; phi applies the recurrence, falls back to $\pi(x)-a+1$ once $p_a^2\gt x$, and caches for $a\le128$, $x\le5\cdot10^{10}$; pi assembles and caches the result.

fsf_from_signature is Section 5: it builds the $2^r-1$ masks, sorts them by decreasing population count, precomputes their bit lists in mask_bits, holds the remaining exponents in a ten-byte array and recurses over mask positions, subtracting and restoring as it goes. pack_small_vector squeezes those exponents into six bits each and FsfMemoKey pairs that word with the mask position; FsfMemoKeyHash only hashes the key while equality compares it in full, so distinct states cannot be confused when their hash codes collide.

generate_signatures_dfs is Section 7: it extends the pattern along the sieved primes, holding each exponent at most last_exp and the product within limit, and emits every prefix. generate_unique_permutations is the ordering step of Section 8, driven by a frequency array so each distinct ordering appears once. count_for_exponents is the recursion: tail_sum holds the $T_k$, the last level returns pi(root) - start_index, the others loop over primes from start_index while $p^{a_k}$ fits, and CountMemoKey is the remaining limit, start index and level.

solve_fast ties it together: one Task per ordered tuple, carrying the tuple and the $F$ value of its signature, with $\sum F\cdot C$ accumulated in an unsigned 128-bit integer and printed by to_string_u128. The checks in run_checkpoints rest on a separate brute force — squarefree_numbers sieves the squarefree numbers and brute_fsf_single enumerates non-decreasing sequences of squarefree factors directly, memoised on the remaining number and the first admissible factor index. The program demands $S(100)=193$ from both methods, their agreement at the limit $500$, and that the fast method reproduce its own value at $200000$ under two different execution settings. --skip-checkpoints omits them, --limit= changes the bound.

The Java file implements the same stages with camel-case names, but is not a line-for-line port. It accumulates in long, justified by the bound below, and corrects Math.sqrt, Math.cbrt and Math.pow estimates to exact integer roots. solveFast walks each signature and its orderings directly rather than building a task list, and main runs no checkpoints and parses no option: solve always takes $10^{10}$. Its memo tables are HashMaps keyed by the mixed 64-bit word itself — packRemaining(remaining, r) ^ pos * 0x9E3779B97F4A7C15L in fsfDfs, and the analogous mix in countDfs — not by a structural key, so distinct states mixing alike would share an entry.

The Python file is a bridge, not a third implementation: solve compiles the C++ source from the sibling solutionsCpp directory with clang++ or g++ at -std=c++17 -O2 when the cached binary is missing or stale, then runs it at the default limit. Because that source advertises --skip-checkpoints, run_cpp passes the option, so the packaged checks never run through this tab. parse_output prefers the last answer: line, then the last =, then the last non-empty line; here the second case applies.

Complexity and Verification

The sieve is linear, each composite being produced once from its least prime factor, so it costs $O(\sqrt N)$ time and memory, $10^{5}$ entries here; the wheel table is the largest single allocation, $510511\times8$ entries of 32 bits, about $16.3$ MB. The enumeration gives $1961$ signatures and $53558$ ordered tuples. One signature's sweep has at most $(2^r-1)\prod_i(e_i+1)$ states, each with at most $e_1+1$ branches; summed over all of them the $F$ sweeps create $16362770$ memo entries and the tuple recursion $5909014$. Those totals are cumulative — each table is discarded when its signature or ordering is done, the largest holding $523776$ and $556378$ entries. The run asks $\pi$ for $198458$ distinct arguments, $99999$ above $10^{5}$ and the largest $10^{10}$; those inside the sieve and the cached repeats are lookups, the rest expand the recurrence of Section 9. No closed form covers the whole cost, which follows the exponent patterns, so the state counts are the reproducible description; the measured $1.4$ seconds, and $0.02$ for the checkpoints, are one machine at -O2.

The accumulator width settles exactly. Over the $1961$ feasible signatures the largest $F$ is $360909$, at $(3,3,1,1,1,1,1,1,1)$, whose representative is $2^3\cdot3^3\cdot5\cdot7\cdot11\cdot13\cdot17\cdot19\cdot23$${}=8031343320\le10^{10}$. Hence

$$\begin{aligned}S(10^{10})&\le(10^{10}-1)\cdot360909\\&\lt3.61\cdot10^{15},\end{aligned}$$

some $2555$ times smaller than $2^{63}-1$, so Java's signed 64-bit accumulator is safe and the 128-bit one in C++ is conservative rather than necessary.

Beyond the packaged checkpoints, the following were verified here.

Section 3: for every $2\le n\le20000$, $F$ of the signature of $n$ equals $\operatorname{Fsf}(n)$ found by enumerating the factorisations of $n$ directly. No exception.
An independent program using neither signatures nor prime counting, which gets every $\operatorname{Fsf}(n)$ below a bound from the identity that the product of $(1-q^{-s})^{-1}$ over squarefree $q\ge2$ has constant term $1$, for the empty factorisation, and coefficient $\operatorname{Fsf}(n)$ at $n^{-s}$ for every $n\ge2$, ^[3] one pass per $q$. It gives $S(100)=193$, $S(1000)=3225$, $S(10^4)=51272$, $S(10^6)=11626006$, $S(10^7)=168329190$, $S(10^8)=2388394374$ — all six reproduced by the method of this article and by the C++ program with --limit=.
The prime counter, re-implemented from the same recurrences, gives $\pi(10^{5})=9592$, $\pi(10^{6})=78498$, $\pi(10^{7})=664579$, $\pi(10^{8})=5761455$, $\pi(10^{9})=50847534$ and $\pi(10^{10})=455052511$, matching the published table ^[7] and a direct sieve at every integer up to $2\cdot10^{6}$.
Record values, as a check on the large signatures: $\operatorname{Fsf}(900900)=1639$ is the largest below $10^{6}$, $\operatorname{Fsf}(6126120)=5443$ below $10^{7}$ and $\operatorname{Fsf}(91891800)=23992$ below $10^{8}$.
The Java key question, measured rather than assumed: replaying both Java memo tables over all $1961$ signatures and all $53558$ tuples at $N=10^{10}$ found no colliding key, every $F$ value was right, and the Java output matches the C++ output. The safety is empirical, a hash being non-injective by construction.
One limitation of the C++ file: its fixed ten-element arrays for the remaining exponents and tail sums fit the target exactly, since $r\le10$ at $10^{10}$. From the limit $200560490130$ upward, length-$11$ signatures become feasible and exceed them. Java sizes both arrays by the signature length.

Footnotes and References

Claim 1, the signature invariance, the minimality of $m(\sigma)$, the root bound of Section 8, Lehmer's formula as used in Section 9, and all the measured counts are derived or computed here; the references supply the rest.

Project Euler 362 — Squarefree Factors. The official statement lists the seven factorisations of 54, keeps only 3 × 3 × 6 and 2 × 3 × 3 × 3 as squarefree throughout, defines Fsf(n) and S(n) as the sum of Fsf(k) for k from 2 to n, gives S(100) = 193 and asks for S at 10 000 000 000.
Wikipedia — Square-free integer. A square-free integer is divisible by no square other than 1, equivalently no prime factor appears with an exponent greater than one; the asymptotic density of such numbers is 6/π².
OEIS A050320 — Number of ways n is a product of squarefree numbers larger than 1. For n ≥ 2 this is Fsf(n); its a(1) = 1 counts the empty factorisation. The entry gives the Dirichlet generating function as the product of 1/(1 − 1/n^s) over squarefree n larger than 1, the identity behind the independent check above; a(p^k) = 1 and a(primorial) = Bell number, the first and third cases of Section 4; and a(n) = A050325(A101296(n)), the invariance of Section 3. The companion entry OEIS A050325 — Number of factorizations into squarefree factors indexed by prime signatures is A050320 evaluated at A025487, that is F(σ) itself.
OEIS A001055 — The multiplicative partition function, the number of unordered factorisations of n into factors greater than 1. Its a(54) = 7 is the statement's count before the squarefree restriction, so Fsf never exceeds it.
OEIS A025487 — Least integer of each prime signature; also products of primorial numbers. A comment describes its terms as the numbers 2^k₁ · 3^k₂ ⋯ pₙ^kₙ with k₁ ≥ k₂ ≥ ⋯ ≥ kₙ, which is m(σ).
OEIS A000110 — Bell numbers, the number of partitions of a set of n labelled elements: 1, 1, 2, 5, 15, 52, 203, 877, 4140, 21147, 115975 for n = 0 to 10, the values F(1, …, 1) of Section 4. Wikipedia — Bell number gives the same definition and the recurrence B_(n+1) = Σ C(n, k) B_k.
Wikipedia — Prime-counting function. π(x) counts the primes not exceeding x; its table at powers of ten gives π(10^10) = 455052511, and it describes Meissel's method, Lehmer's extension and the partial sieve function used in Section 9.
Wikipedia — Meissel–Lehmer algorithm. It gives the same recurrence φ(x, 0) = ⌊x⌋, φ(x, a) = φ(x, a−1) − φ(x/p_a, a−1), and states that the Lagarias–Miller–Odlyzko refinement runs in O(x^(2/3+ε)) time and O(x^(1/3+ε)) space — the bound Section 9 declines to attribute to Lehmer's formula.
C++ — Euler362.cpp. The linked immutable revision holds the functions walked through above and the brute-force checkpoints, and is the implementation record for this article; the independent cross-checks are not part of it.

Problem 362 source code

C++

#include <algorithm>
#include <array>
#include <atomic>
#include <chrono>
#include <cmath>
#include <cstdint>
#include <functional>
#include <iostream>
#include <limits>
#include <stdexcept>
#include <string>
#include <thread>
#include <unordered_map>
#include <vector>

namespace {

using u64 = std::uint64_t;
using u32 = std::uint32_t;
using u16 = std::uint16_t;
using u8 = std::uint8_t;
using u128 = unsigned __int128;

constexpr u64 kDefaultLimit = 10'000'000'000ULL;
constexpr u64 kCheckpointKnownLimit = 100ULL;
constexpr u64 kCheckpointKnownExpected = 193ULL;
constexpr u64 kCheckpointBruteLimit = 500ULL;
constexpr u64 kThreadConsistencyLimit = 200'000ULL;

struct Options {
    u64 limit = kDefaultLimit;
    bool run_checkpoints = true;
    bool allow_multithreading = true;
    unsigned requested_threads = 0U;
};

bool parse_u64_after_prefix(const std::string& arg, const char* prefix, u64& value) {
    const std::string p(prefix);
    if (arg.rfind(p, 0) != 0) {
        return false;
    }

    const std::string tail = arg.substr(p.size());
    if (tail.empty()) {
        return false;
    }

    u64 parsed = 0ULL;
    for (const char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        const u64 digit = static_cast<u64>(c - '0');
        if (parsed > (std::numeric_limits<u64>::max() - digit) / 10ULL) {
            return false;
        }
        parsed = parsed * 10ULL + digit;
    }

    value = parsed;
    return true;
}

bool parse_unsigned_after_prefix(const std::string& arg,
                                 const char* prefix,
                                 unsigned& value) {
    u64 parsed = 0ULL;
    if (!parse_u64_after_prefix(arg, prefix, parsed)) {
        return false;
    }
    if (parsed > static_cast<u64>(std::numeric_limits<unsigned>::max())) {
        return false;
    }
    value = static_cast<unsigned>(parsed);
    return true;
}

bool parse_arguments(int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);

        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        if (arg == "--single-thread") {
            options.allow_multithreading = false;
            continue;
        }

        u64 parsed_u64 = 0ULL;
        if (parse_u64_after_prefix(arg, "--limit=", parsed_u64)) {
            options.limit = parsed_u64;
            continue;
        }

        unsigned parsed_unsigned = 0U;
        if (parse_unsigned_after_prefix(arg, "--threads=", parsed_unsigned)) {
            options.requested_threads = parsed_unsigned;
            continue;
        }

        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }

    if (options.limit < 2ULL) {
        std::cerr << "--limit must be at least 2.\n";
        return false;
    }

    return true;
}

unsigned choose_thread_count(bool allow_multithreading,
                             unsigned requested_threads,
                             std::size_t workload) {
    if (!allow_multithreading || workload < 2ULL) {
        return 1U;
    }

    unsigned threads = requested_threads;
    if (threads == 0U) {
        threads = std::thread::hardware_concurrency();
        if (threads == 0U) {
            threads = 1U;
        }
    }

    return std::max(1U, std::min<unsigned>(threads, static_cast<unsigned>(workload)));
}

u64 isqrt_u64(u64 x) {
    u64 r = static_cast<u64>(std::sqrt(static_cast<long double>(x)));
    while ((r + 1ULL) <= x / (r + 1ULL)) {
        ++r;
    }
    while (r > x / r) {
        --r;
    }
    return r;
}

u64 pow_u64_capped(u64 base, int exp, u64 cap) {
    u64 result = 1ULL;
    for (int i = 0; i < exp; ++i) {
        if (base != 0ULL && result > cap / base) {
            return cap + 1ULL;
        }
        result *= base;
    }
    return result;
}

u64 kth_root_floor(u64 n, int k) {
    if (k <= 1 || n <= 1ULL) {
        return n;
    }

    const long double approx =
        std::pow(static_cast<long double>(n), 1.0L / static_cast<long double>(k));
    u64 r = static_cast<u64>(approx);
    if (r == 0ULL) {
        r = 1ULL;
    }

    while (pow_u64_capped(r + 1ULL, k, n) <= n) {
        ++r;
    }
    while (pow_u64_capped(r, k, n) > n) {
        --r;
    }
    return r;
}

std::string to_string_u128(u128 value) {
    if (value == 0) {
        return "0";
    }

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

struct PrimeTable {
    int limit = 0;
    std::vector<int> primes;
    std::vector<int> pi_small;

    explicit PrimeTable(int sieve_limit) { build(sieve_limit); }

    void build(int sieve_limit) {
        limit = sieve_limit;
        std::vector<int> lp(static_cast<std::size_t>(limit) + 1ULL, 0);
        primes.reserve(static_cast<std::size_t>(limit) / 10ULL + 32ULL);

        for (int i = 2; i <= limit; ++i) {
            if (lp[static_cast<std::size_t>(i)] == 0) {
                lp[static_cast<std::size_t>(i)] = i;
                primes.push_back(i);
            }

            for (const int p : primes) {
                const long long v = 1LL * p * i;
                if (p > lp[static_cast<std::size_t>(i)] || v > limit) {
                    break;
                }
                lp[static_cast<std::size_t>(v)] = p;
            }
        }

        pi_small.assign(static_cast<std::size_t>(limit) + 1ULL, 0);
        int ptr = 0;
        for (int i = 1; i <= limit; ++i) {
            pi_small[static_cast<std::size_t>(i)] =
                pi_small[static_cast<std::size_t>(i - 1)];
            if (ptr < static_cast<int>(primes.size()) && primes[static_cast<std::size_t>(ptr)] == i) {
                ++pi_small[static_cast<std::size_t>(i)];
                ++ptr;
            }
        }
    }
};

u64 icbrt_u64(u64 x) {
    u64 r = static_cast<u64>(std::cbrt(static_cast<long double>(x)));
    while (static_cast<u128>(r + 1ULL) * (r + 1ULL) * (r + 1ULL) <= x) {
        ++r;
    }
    while (static_cast<u128>(r) * r * r > x) {
        --r;
    }
    return r;
}

u64 iroot4_u64(u64 x) {
    u64 r = static_cast<u64>(std::sqrt(std::sqrt(static_cast<long double>(x))));
    while (static_cast<u128>(r + 1ULL) * (r + 1ULL) * (r + 1ULL) * (r + 1ULL) <= x) {
        ++r;
    }
    while (static_cast<u128>(r) * r * r * r > x) {
        --r;
    }
    return r;
}

class PrimeCounter {
public:
    explicit PrimeCounter(const PrimeTable& table) : table_(table) { build_phi_table(); }

    u64 pi(u64 x) const {
        if (x <= static_cast<u64>(table_.limit)) {
            return static_cast<u64>(table_.pi_small[static_cast<std::size_t>(x)]);
        }

        auto& cache = pi_cache();
        const auto it = cache.find(x);
        if (it != cache.end()) {
            return it->second;
        }

        const u64 a = pi(iroot4_u64(x));
        const u64 b = pi(isqrt_u64(x));
        const u64 c = pi(icbrt_u64(x));

        i64 sum = static_cast<i64>(phi(x, static_cast<int>(a)));
        sum += static_cast<i64>((b + a - 2ULL) * (b - a + 1ULL) / 2ULL);

        for (u64 i = a + 1ULL; i <= b; ++i) {
            const u64 w = x / static_cast<u64>(table_.primes[static_cast<std::size_t>(i - 1ULL)]);
            sum -= static_cast<i64>(pi(w));
            if (i <= c) {
                const u64 bi = pi(isqrt_u64(w));
                for (u64 j = i; j <= bi; ++j) {
                    const u64 q =
                        w / static_cast<u64>(table_.primes[static_cast<std::size_t>(j - 1ULL)]);
                    sum -= static_cast<i64>(pi(q)) - static_cast<i64>(j - 1ULL);
                }
            }
        }

        const u64 ans = static_cast<u64>(sum);
        cache.emplace(x, ans);
        return ans;
    }

private:
    using i64 = std::int64_t;

    static constexpr int kSmallPhiPrimes = 7;   // 2,3,5,7,11,13,17
    static constexpr int kSmallPhiMod = 510510;
    static constexpr int kPhiCacheSMax = 128;
    static constexpr u64 kPhiCacheXMax = 50'000'000'000ULL;

    const PrimeTable& table_;
    std::vector<std::array<u32, kSmallPhiPrimes + 1>> phi_table_;

    void build_phi_table() {
        phi_table_.assign(static_cast<std::size_t>(kSmallPhiMod) + 1ULL, {});
        for (int n = 0; n <= kSmallPhiMod; ++n) {
            phi_table_[static_cast<std::size_t>(n)][0] = static_cast<u32>(n);
        }

        for (int s = 1; s <= kSmallPhiPrimes; ++s) {
            const int p = table_.primes[static_cast<std::size_t>(s - 1)];
            for (int n = 0; n <= kSmallPhiMod; ++n) {
                phi_table_[static_cast<std::size_t>(n)][static_cast<std::size_t>(s)] =
                    phi_table_[static_cast<std::size_t>(n)][static_cast<std::size_t>(s - 1)] -
                    phi_table_[static_cast<std::size_t>(n / p)][static_cast<std::size_t>(s - 1)];
            }
        }
    }

    static std::unordered_map<u64, u64>& pi_cache() {
        static thread_local std::unordered_map<u64, u64> cache;
        return cache;
    }

    static std::unordered_map<u64, u64>& phi_cache() {
        static thread_local std::unordered_map<u64, u64> cache;
        return cache;
    }

    u64 phi(u64 x, int s) const {
        if (s == 0) {
            return x;
        }

        if (s <= kSmallPhiPrimes) {
            const u64 block = x / static_cast<u64>(kSmallPhiMod);
            const u64 rem = x % static_cast<u64>(kSmallPhiMod);
            return block * static_cast<u64>(phi_table_[kSmallPhiMod][static_cast<std::size_t>(s)]) +
                   static_cast<u64>(
                       phi_table_[static_cast<std::size_t>(rem)][static_cast<std::size_t>(s)]);
        }

        if (x <= static_cast<u64>(table_.limit) &&
            static_cast<u64>(table_.primes[static_cast<std::size_t>(s - 1)]) *
                    static_cast<u64>(table_.primes[static_cast<std::size_t>(s - 1)]) >
                x) {
            return static_cast<u64>(table_.pi_small[static_cast<std::size_t>(x)] - s + 1);
        }

        const u64 ps = static_cast<u64>(table_.primes[static_cast<std::size_t>(s - 1)]);
        if (ps * ps > x) {
            return pi(x) - static_cast<u64>(s) + 1ULL;
        }

        if (s <= kPhiCacheSMax && x <= kPhiCacheXMax) {
            const u64 key = (x << 8U) ^ static_cast<u64>(s);
            auto& cache = phi_cache();
            const auto it = cache.find(key);
            if (it != cache.end()) {
                return it->second;
            }
            const u64 ans = phi(x, s - 1) - phi(x / ps, s - 1);
            cache.emplace(key, ans);
            return ans;
        }

        return phi(x, s - 1) - phi(x / ps, s - 1);
    }
};

struct FsfMemoKey {
    u16 mask_pos = 0U;
    u64 packed_state = 0ULL;

    bool operator==(const FsfMemoKey& other) const {
        return mask_pos == other.mask_pos && packed_state == other.packed_state;
    }
};

struct FsfMemoKeyHash {
    std::size_t operator()(const FsfMemoKey& key) const {
        const u64 mixed = key.packed_state ^ (static_cast<u64>(key.mask_pos) * 0x9E3779B97F4A7C15ULL);
        return static_cast<std::size_t>(mixed ^ (mixed >> 33U));
    }
};

u64 pack_small_vector(const std::array<u8, 10>& values, int size) {
    u64 packed = 0ULL;
    for (int i = 0; i < size; ++i) {
        packed |= static_cast<u64>(values[static_cast<std::size_t>(i)]) << (6U * i);
    }
    return packed;
}

u64 fsf_from_signature(const std::vector<u8>& signature) {
    const int r = static_cast<int>(signature.size());
    if (r == 0) {
        return 0ULL;
    }

    std::vector<int> masks;
    masks.reserve((1 << r) - 1);
    for (int mask = 1; mask < (1 << r); ++mask) {
        masks.push_back(mask);
    }
    std::sort(masks.begin(), masks.end(), [](int a, int b) {
        const int pa = __builtin_popcount(static_cast<unsigned>(a));
        const int pb = __builtin_popcount(static_cast<unsigned>(b));
        if (pa != pb) {
            return pa > pb;
        }
        return a < b;
    });

    std::vector<std::vector<u8>> mask_bits;
    mask_bits.reserve(masks.size());
    for (const int mask : masks) {
        std::vector<u8> bits;
        for (int i = 0; i < r; ++i) {
            if ((mask >> i) & 1) {
                bits.push_back(static_cast<u8>(i));
            }
        }
        mask_bits.push_back(std::move(bits));
    }

    std::array<u8, 10> remaining{};
    for (int i = 0; i < r; ++i) {
        remaining[static_cast<std::size_t>(i)] = signature[static_cast<std::size_t>(i)];
    }

    std::unordered_map<FsfMemoKey, u64, FsfMemoKeyHash> memo;
    memo.reserve(4096);

    std::function<u64(u16)> dfs = [&](u16 pos) -> u64 {
        if (pos == static_cast<u16>(mask_bits.size())) {
            for (int i = 0; i < r; ++i) {
                if (remaining[static_cast<std::size_t>(i)] != 0U) {
                    return 0ULL;
                }
            }
            return 1ULL;
        }

        const FsfMemoKey key{pos, pack_small_vector(remaining, r)};
        const auto it = memo.find(key);
        if (it != memo.end()) {
            return it->second;
        }

        const std::vector<u8>& bits = mask_bits[static_cast<std::size_t>(pos)];
        u8 max_take = std::numeric_limits<u8>::max();
        for (const u8 bit : bits) {
            max_take = std::min(max_take, remaining[static_cast<std::size_t>(bit)]);
        }

        u64 total = 0ULL;
        for (u32 take = 0U;; ++take) {
            total += dfs(static_cast<u16>(pos + 1U));
            if (take == max_take) {
                break;
            }
            for (const u8 bit : bits) {
                --remaining[static_cast<std::size_t>(bit)];
            }
        }

        for (const u8 bit : bits) {
            remaining[static_cast<std::size_t>(bit)] =
                static_cast<u8>(remaining[static_cast<std::size_t>(bit)] + max_take);
        }

        memo.emplace(key, total);
        return total;
    };

    return dfs(0U);
}

void generate_signatures_dfs(const std::vector<int>& seed_primes,
                             int prime_index,
                             int last_exp,
                             u64 current,
                             u64 limit,
                             std::vector<u8>& cur,
                             std::vector<std::vector<u8>>& out) {
    if (prime_index >= static_cast<int>(seed_primes.size())) {
        return;
    }

    const u64 p = static_cast<u64>(seed_primes[static_cast<std::size_t>(prime_index)]);
    u64 p_pow = 1ULL;

    for (int e = 1; e <= last_exp; ++e) {
        if (p_pow > limit / p) {
            break;
        }
        p_pow *= p;
        if (current > limit / p_pow) {
            break;
        }

        const u64 next_value = current * p_pow;
        cur.push_back(static_cast<u8>(e));
        out.push_back(cur);
        generate_signatures_dfs(seed_primes, prime_index + 1, e, next_value, limit, cur, out);
        cur.pop_back();
    }
}

std::vector<std::vector<u8>> generate_signatures(u64 limit, const PrimeTable& table) {
    std::vector<std::vector<u8>> signatures;
    std::vector<u8> cur;
    generate_signatures_dfs(table.primes, 0, 63, 1ULL, limit, cur, signatures);
    return signatures;
}

void generate_unique_permutations_rec(std::array<int, 64>& freq,
                                      std::vector<u8>& cur,
                                      std::size_t pos,
                                      std::vector<std::vector<u8>>& out) {
    if (pos == cur.size()) {
        out.push_back(cur);
        return;
    }

    for (int value = 1; value < static_cast<int>(freq.size()); ++value) {
        if (freq[static_cast<std::size_t>(value)] == 0) {
            continue;
        }
        --freq[static_cast<std::size_t>(value)];
        cur[pos] = static_cast<u8>(value);
        generate_unique_permutations_rec(freq, cur, pos + 1ULL, out);
        ++freq[static_cast<std::size_t>(value)];
    }
}

std::vector<std::vector<u8>> generate_unique_permutations(const std::vector<u8>& signature) {
    std::array<int, 64> freq{};
    for (const u8 e : signature) {
        ++freq[static_cast<std::size_t>(e)];
    }

    std::vector<std::vector<u8>> out;
    std::vector<u8> cur(signature.size(), 0U);
    generate_unique_permutations_rec(freq, cur, 0ULL, out);
    return out;
}

struct CountMemoKey {
    u64 limit = 0ULL;
    u32 start_index = 0U;
    u8 pos = 0U;

    bool operator==(const CountMemoKey& other) const {
        return limit == other.limit && start_index == other.start_index && pos == other.pos;
    }
};

struct CountMemoKeyHash {
    std::size_t operator()(const CountMemoKey& key) const {
        u64 h = key.limit;
        h ^= static_cast<u64>(key.start_index) * 0x9E3779B97F4A7C15ULL;
        h ^= static_cast<u64>(key.pos) * 0xC2B2AE3D27D4EB4FULL;
        h ^= (h >> 31U);
        return static_cast<std::size_t>(h);
    }
};

class OrderedPrimeTupleCounter {
public:
    OrderedPrimeTupleCounter(const PrimeTable& table, const PrimeCounter& prime_counter)
        : table_(table), prime_counter_(prime_counter) {}

    u64 count_for_exponents(const std::vector<u8>& exponents, u64 limit) const {
        const int r = static_cast<int>(exponents.size());
        if (r == 0) {
            return 0ULL;
        }
        if (r == 1) {
            const u64 root = kth_root_floor(limit, static_cast<int>(exponents[0]));
            return prime_counter_.pi(root);
        }

        std::array<int, 10> tail_sum{};
        int running = 0;
        for (int i = r - 1; i >= 0; --i) {
            running += static_cast<int>(exponents[static_cast<std::size_t>(i)]);
            tail_sum[static_cast<std::size_t>(i)] = running;
        }

        std::unordered_map<CountMemoKey, u64, CountMemoKeyHash> memo;
        memo.reserve(4096);

        std::function<u64(int, u32, u64)> dfs = [&](int pos, u32 start_index, u64 rem_limit) -> u64 {
            if (pos == r) {
                return 1ULL;
            }

            const CountMemoKey key{rem_limit, start_index, static_cast<u8>(pos)};
            const auto it = memo.find(key);
            if (it != memo.end()) {
                return it->second;
            }

            u64 result = 0ULL;

            if (pos == r - 1) {
                const int exp = static_cast<int>(exponents[static_cast<std::size_t>(pos)]);
                const u64 root = kth_root_floor(rem_limit, exp);
                const u64 upto = prime_counter_.pi(root);
                result = (upto > static_cast<u64>(start_index))
                             ? (upto - static_cast<u64>(start_index))
                             : 0ULL;
                memo.emplace(key, result);
                return result;
            }

            const int tail = tail_sum[static_cast<std::size_t>(pos)];
            const u64 max_p = kth_root_floor(rem_limit, tail);
            u64 upto = prime_counter_.pi(max_p);
            if (upto > table_.primes.size()) {
                upto = static_cast<u64>(table_.primes.size());
            }

            const int exp = static_cast<int>(exponents[static_cast<std::size_t>(pos)]);
            for (u64 idx = static_cast<u64>(start_index); idx < upto; ++idx) {
                const u64 p = static_cast<u64>(table_.primes[static_cast<std::size_t>(idx)]);
                const u64 p_pow = pow_u64_capped(p, exp, rem_limit);
                if (p_pow > rem_limit) {
                    break;
                }
                result += dfs(pos + 1, static_cast<u32>(idx + 1ULL), rem_limit / p_pow);
            }

            memo.emplace(key, result);
            return result;
        };

        return dfs(0, 0U, limit);
    }

private:
    const PrimeTable& table_;
    const PrimeCounter& prime_counter_;
};

struct Task {
    std::vector<u8> exponents;
    u64 fsf_value = 0ULL;
};

u128 solve_fast(u64 limit, bool allow_multithreading, unsigned requested_threads) {
    const u64 root_limit = std::max<u64>(isqrt_u64(limit), 100ULL);
    if (root_limit > static_cast<u64>(std::numeric_limits<int>::max())) {
        throw std::runtime_error("Sieve limit too large for this implementation.");
    }

    const PrimeTable table(static_cast<int>(root_limit));
    const PrimeCounter prime_counter(table);
    const OrderedPrimeTupleCounter tuple_counter(table, prime_counter);

    const std::vector<std::vector<u8>> signatures = generate_signatures(limit, table);

    std::vector<Task> tasks;
    tasks.reserve(16'000ULL);
    for (const auto& signature : signatures) {
        const u64 fsf = fsf_from_signature(signature);
        const std::vector<std::vector<u8>> permutations = generate_unique_permutations(signature);
        for (const auto& p : permutations) {
            tasks.push_back(Task{p, fsf});
        }
    }

    const unsigned threads =
        choose_thread_count(allow_multithreading, requested_threads, tasks.size());

    if (threads == 1U) {
        u128 total = 0;
        for (const Task& task : tasks) {
            const u64 count = tuple_counter.count_for_exponents(task.exponents, limit);
            total += static_cast<u128>(task.fsf_value) * static_cast<u128>(count);
        }
        return total;
    }

    std::atomic<std::size_t> next_task{0ULL};
    std::vector<u128> partials(threads, 0);
    std::vector<std::thread> workers;
    workers.reserve(threads);

    for (unsigned tid = 0; tid < threads; ++tid) {
        workers.emplace_back([&, tid]() {
            const OrderedPrimeTupleCounter local_counter(table, prime_counter);
            u128 local = 0;
            while (true) {
                const std::size_t idx = next_task.fetch_add(1ULL, std::memory_order_relaxed);
                if (idx >= tasks.size()) {
                    break;
                }
                const Task& task = tasks[idx];
                const u64 count = local_counter.count_for_exponents(task.exponents, limit);
                local += static_cast<u128>(task.fsf_value) * static_cast<u128>(count);
            }
            partials[static_cast<std::size_t>(tid)] = local;
        });
    }

    for (auto& th : workers) {
        th.join();
    }

    u128 total = 0;
    for (const u128 value : partials) {
        total += value;
    }
    return total;
}

std::vector<int> squarefree_numbers(int limit) {
    std::vector<int> sqf;
    std::vector<int> square_div(static_cast<std::size_t>(limit) + 1ULL, 0);
    for (int p = 2; p * p <= limit; ++p) {
        const int pp = p * p;
        for (int n = pp; n <= limit; n += pp) {
            square_div[static_cast<std::size_t>(n)] = 1;
        }
    }
    for (int n = 2; n <= limit; ++n) {
        if (square_div[static_cast<std::size_t>(n)] == 0) {
            sqf.push_back(n);
        }
    }
    return sqf;
}

u64 brute_fsf_single(u64 n,
                     std::size_t start_idx,
                     const std::vector<int>& sqf,
                     std::unordered_map<u64, u64>& memo) {
    const u64 key = (n << 16U) ^ static_cast<u64>(start_idx);
    const auto it = memo.find(key);
    if (it != memo.end()) {
        return it->second;
    }

    u64 total = 0ULL;
    for (std::size_t i = start_idx; i < sqf.size(); ++i) {
        const u64 d = static_cast<u64>(sqf[i]);
        if (d > n) {
            break;
        }
        if (n % d != 0ULL) {
            continue;
        }
        if (d == n) {
            ++total;
        } else {
            total += brute_fsf_single(n / d, i, sqf, memo);
        }
    }

    memo.emplace(key, total);
    return total;
}

u64 brute_s(u64 limit) {
    if (limit < 2ULL) {
        return 0ULL;
    }
    if (limit > 100'000ULL) {
        throw std::runtime_error("brute_s is intended only for small limits.");
    }

    const std::vector<int> sqf = squarefree_numbers(static_cast<int>(limit));
    u64 total = 0ULL;
    for (u64 n = 2ULL; n <= limit; ++n) {
        std::unordered_map<u64, u64> memo;
        total += brute_fsf_single(n, 0ULL, sqf, memo);
    }
    return total;
}

void run_checkpoints() {
    const u64 brute_100 = brute_s(kCheckpointKnownLimit);
    if (brute_100 != kCheckpointKnownExpected) {
        throw std::runtime_error("Brute-force checkpoint S(100)=193 failed.");
    }

    const u64 fast_100 = static_cast<u64>(solve_fast(kCheckpointKnownLimit, false, 1U));
    if (fast_100 != kCheckpointKnownExpected) {
        throw std::runtime_error("Fast checkpoint S(100)=193 failed.");
    }

    const u64 brute_small = brute_s(kCheckpointBruteLimit);
    const u64 fast_small = static_cast<u64>(solve_fast(kCheckpointBruteLimit, false, 1U));
    if (fast_small != brute_small) {
        throw std::runtime_error("Brute-force consistency checkpoint failed.");
    }

    const u64 single = static_cast<u64>(solve_fast(kThreadConsistencyLimit, false, 1U));
    const u64 threaded = static_cast<u64>(solve_fast(kThreadConsistencyLimit, true, 2U));
    if (single != threaded) {
        throw std::runtime_error("Thread consistency checkpoint failed.");
    }
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }

    try {
        if (options.run_checkpoints) {
            run_checkpoints();
        }

        const auto start = std::chrono::steady_clock::now();
        const u128 answer = solve_fast(
            options.limit, options.allow_multithreading, options.requested_threads);
        const auto end = std::chrono::steady_clock::now();
        const std::chrono::duration<double> elapsed = end - start;

        std::cout << "S(" << options.limit << ") = " << to_string_u128(answer) << '\n';
        std::cout << "Elapsed: " << elapsed.count() << " s\n";
    } catch (const std::exception& ex) {
        std::cerr << "Error: " << ex.what() << '\n';
        return 1;
    }

    return 0;
}

Python

from __future__ import annotations

import re
import shutil
import subprocess
from pathlib import Path

ANSWER_RE = re.compile(r"answer\s*:\s*(.+)$", re.IGNORECASE)
EQUAL_RE = re.compile(r"=\s*(.+)$")


def parse_output(stdout: str) -> str:
    lines = [line.strip() for line in stdout.splitlines() if line.strip()]
    if not lines:
        return ""
    answers = []
    equals = []
    for line in lines:
        m1 = ANSWER_RE.search(line)
        if m1:
            answers.append(m1.group(1).strip())
        m2 = EQUAL_RE.search(line)
        if m2:
            equals.append(m2.group(1).strip())
    if answers:
        return answers[-1]
    if equals:
        return equals[-1]
    return lines[-1]


def should_skip_cpp_checkpoints(src: Path) -> bool:
    try:
        text = src.read_text(encoding="utf-8", errors="ignore")
    except OSError:
        return False
    return "--skip-checkpoints" in text


def run_cpp(binary: Path, src: Path, root: Path) -> str:
    cmd = [str(binary)]
    if should_skip_cpp_checkpoints(src):
        cmd.append("--skip-checkpoints")

    try:
        return subprocess.check_output(cmd, text=True, cwd=root)
    except subprocess.CalledProcessError:
        return subprocess.check_output(cmd, text=True, cwd=src.parent)


def solve() -> str:
    problem_id = __file__.split("Euler")[-1].split(".")[0]
    root = Path(__file__).resolve().parent.parent
    src = root / "solutionsCpp" / f"Euler{problem_id}.cpp"
    binary = root / "solutionsCpp" / f".euler{problem_id}_py_bridge"

    if not binary.exists() or src.stat().st_mtime > binary.stat().st_mtime:
        compiler = shutil.which("clang++") or shutil.which("g++")
        if not compiler:
            raise RuntimeError("No C++ compiler found (clang++/g++).")
        subprocess.check_call([compiler, "-std=c++17", "-O2", str(src), "-o", str(binary)])

    output = run_cpp(binary=binary, src=src, root=root)
    parsed = parse_output(output)
    if not parsed:
        raise RuntimeError(f"Euler{problem_id} bridge produced empty output.")
    return parsed


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

Java

import java.util.*;

public class Euler362 {

    static long isqrt(long n) {
        if (n < 0)
            return 0;
        long r = (long) Math.sqrt(n);
        while ((r + 1) <= n / (r + 1))
            r++;
        while (r > n / r)
            r--;
        return r;
    }

    static long icbrt(long n) {
        if (n < 0)
            return 0;
        long r = (long) Math.cbrt(n);
        while ((r + 1) * (r + 1) * (r + 1) <= n)
            r++;
        while (r * r * r > n)
            r--;
        return r;
    }

    static long iroot4(long n) {
        if (n < 0)
            return 0;
        long r = (long) Math.sqrt(Math.sqrt(n));
        while ((r + 1) * (r + 1) * (r + 1) * (r + 1) <= n)
            r++;
        while (r * r * r * r > n)
            r--;
        return r;
    }

    static long powCapped(long base, int exp, long cap) {
        long res = 1;
        for (int i = 0; i < exp; i++) {
            if (base != 0 && res > cap / base)
                return cap + 1;
            res *= base;
        }
        return res;
    }

    static long kthRootFloor(long n, int k) {
        if (k <= 1 || n <= 1)
            return n;
        long r = (long) Math.pow(n, 1.0 / k);
        if (r == 0)
            r = 1;
        while (powCapped(r + 1, k, n) <= n)
            r++;
        while (powCapped(r, k, n) > n)
            r--;
        return r;
    }

    static class PrimeTable {
        int limit;
        int[] primes;
        int[] piSmall;

        PrimeTable(int limit) {
            this.limit = limit;
            int[] lp = new int[limit + 1];
            List<Integer> prs = new ArrayList<>();

            for (int i = 2; i <= limit; i++) {
                if (lp[i] == 0) {
                    lp[i] = i;
                    prs.add(i);
                }
                for (int p : prs) {
                    long v = (long) p * i;
                    if (p > lp[i] || v > limit)
                        break;
                    lp[(int) v] = p;
                }
            }

            primes = new int[prs.size()];
            for (int i = 0; i < prs.size(); i++)
                primes[i] = prs.get(i);

            piSmall = new int[limit + 1];
            int ptr = 0;
            for (int i = 1; i <= limit; i++) {
                piSmall[i] = piSmall[i - 1];
                if (ptr < primes.length && primes[ptr] == i) {
                    piSmall[i]++;
                    ptr++;
                }
            }
        }
    }

    static class PrimeCounter {
        PrimeTable table;
        Map<Long, Long> piCache = new HashMap<>();
        Map<Long, Long> phiCache = new HashMap<>();

        static final int K_SMALL_PHI_PRIMES = 7;
        static final int K_SMALL_PHI_MOD = 510510;
        static final int K_PHI_CACHE_S_MAX = 128;
        static final long K_PHI_CACHE_X_MAX = 50000000000L;

        int[][] phiTable;

        PrimeCounter(PrimeTable table) {
            this.table = table;
            phiTable = new int[K_SMALL_PHI_MOD + 1][K_SMALL_PHI_PRIMES + 1];
            for (int n = 0; n <= K_SMALL_PHI_MOD; n++)
                phiTable[n][0] = n;
            for (int s = 1; s <= K_SMALL_PHI_PRIMES; s++) {
                int p = table.primes[s - 1];
                for (int n = 0; n <= K_SMALL_PHI_MOD; n++) {
                    phiTable[n][s] = phiTable[n][s - 1] - phiTable[n / p][s - 1];
                }
            }
        }

        long phi(long x, int s) {
            if (s == 0)
                return x;
            if (s <= K_SMALL_PHI_PRIMES) {
                long block = x / K_SMALL_PHI_MOD;
                long rem = x % K_SMALL_PHI_MOD;
                return block * phiTable[K_SMALL_PHI_MOD][s] + phiTable[(int) rem][s];
            }

            if (x <= table.limit && (long) table.primes[s - 1] * table.primes[s - 1] > x) {
                return table.piSmall[(int) x] - s + 1;
            }

            long ps = table.primes[s - 1];
            if (ps * ps > x)
                return pi(x) - s + 1;

            if (s <= K_PHI_CACHE_S_MAX && x <= K_PHI_CACHE_X_MAX) {
                long key = (x << 8) ^ s;
                if (phiCache.containsKey(key))
                    return phiCache.get(key);
                long ans = phi(x, s - 1) - phi(x / ps, s - 1);
                phiCache.put(key, ans);
                return ans;
            }

            return phi(x, s - 1) - phi(x / ps, s - 1);
        }

        long pi(long x) {
            if (x <= table.limit)
                return table.piSmall[(int) x];
            if (piCache.containsKey(x))
                return piCache.get(x);

            long a = pi(iroot4(x));
            long b = pi(isqrt(x));
            long c = pi(icbrt(x));

            long ans = phi(x, (int) a) + (b + a - 2) * (b - a + 1) / 2;
            for (long i = a + 1; i <= b; i++) {
                long w = x / table.primes[(int) (i - 1)];
                ans -= pi(w);
                if (i <= c) {
                    long bi = pi(isqrt(w));
                    for (long j = i; j <= bi; j++) {
                        long q = w / table.primes[(int) (j - 1)];
                        ans -= pi(q) - (j - 1);
                    }
                }
            }
            piCache.put(x, ans);
            return ans;
        }
    }

    static long fsfFromSignature(List<Integer> signature) {
        int r = signature.size();
        if (r == 0)
            return 0;

        List<Integer> masks = new ArrayList<>();
        for (int m = 1; m < (1 << r); m++)
            masks.add(m);
        masks.sort((a, b) -> {
            int pa = Integer.bitCount(a);
            int pb = Integer.bitCount(b);
            if (pa != pb)
                return Integer.compare(pb, pa);
            return Integer.compare(a, b);
        });

        int[][] maskBits = new int[masks.size()][];
        for (int i = 0; i < masks.size(); i++) {
            int mask = masks.get(i);
            int[] bits = new int[Integer.bitCount(mask)];
            int idx = 0;
            for (int k = 0; k < r; k++) {
                if (((mask >> k) & 1) == 1)
                    bits[idx++] = k;
            }
            maskBits[i] = bits;
        }

        byte[] remaining = new byte[r];
        for (int i = 0; i < r; i++)
            remaining[i] = signature.get(i).byteValue();
        Map<Long, Long> memo = new HashMap<>();

        return fsfDfs(0, maskBits, remaining, r, memo);
    }

    static long packRemaining(byte[] remaining, int r) {
        long packed = 0;
        for (int i = 0; i < r; i++)
            packed |= (long) remaining[i] << (6 * i);
        return packed;
    }

    static long fsfDfs(int pos, int[][] maskBits, byte[] remaining, int r, Map<Long, Long> memo) {
        if (pos == maskBits.length) {
            for (int i = 0; i < r; i++)
                if (remaining[i] != 0)
                    return 0;
            return 1;
        }

        long key = packRemaining(remaining, r) ^ ((long) pos * 0x9E3779B97F4A7C15L);
        if (memo.containsKey(key))
            return memo.get(key);

        int[] bits = maskBits[pos];
        int maxTake = Integer.MAX_VALUE;
        for (int bit : bits)
            maxTake = Math.min(maxTake, remaining[bit]);

        long total = 0;
        for (int take = 0;; take++) {
            total += fsfDfs(pos + 1, maskBits, remaining, r, memo);
            if (take == maxTake)
                break;
            for (int bit : bits)
                remaining[bit]--;
        }

        for (int bit : bits)
            remaining[bit] += maxTake;

        memo.put(key, total);
        return total;
    }

    static void generateSignaturesDfs(int[] primes, int pIdx, int lastExp, long current, long limit, List<Integer> cur,
            List<List<Integer>> out) {
        if (pIdx >= primes.length)
            return;
        long p = primes[pIdx];
        long pPow = 1;

        for (int e = 1; e <= lastExp; e++) {
            if (pPow > limit / p)
                break;
            pPow *= p;
            if (current > limit / pPow)
                break;

            long nextVal = current * pPow;
            cur.add(e);
            out.add(new ArrayList<>(cur));
            generateSignaturesDfs(primes, pIdx + 1, e, nextVal, limit, cur, out);
            cur.remove(cur.size() - 1);
        }
    }

    static List<List<Integer>> generateSignatures(long limit, PrimeTable table) {
        List<List<Integer>> signatures = new ArrayList<>();
        generateSignaturesDfs(table.primes, 0, 63, 1, limit, new ArrayList<>(), signatures);
        return signatures;
    }

    static void generateUniquePermutationsRec(Map<Integer, Integer> freq, int[] cur, int pos, List<List<Integer>> out) {
        if (pos == cur.length) {
            List<Integer> list = new ArrayList<>();
            for (int v : cur)
                list.add(v);
            out.add(list);
            return;
        }
        for (int val : new ArrayList<>(freq.keySet())) {
            int f = freq.get(val);
            if (f > 0) {
                freq.put(val, f - 1);
                cur[pos] = val;
                generateUniquePermutationsRec(freq, cur, pos + 1, out);
                freq.put(val, f);
            }
        }
    }

    static List<List<Integer>> generateUniquePermutations(List<Integer> signature) {
        Map<Integer, Integer> freq = new HashMap<>();
        for (int e : signature)
            freq.put(e, freq.getOrDefault(e, 0) + 1);
        List<List<Integer>> out = new ArrayList<>();
        generateUniquePermutationsRec(freq, new int[signature.size()], 0, out);
        return out;
    }

    static class OrderedPrimeTupleCounter {
        PrimeTable table;
        PrimeCounter primeCounter;

        OrderedPrimeTupleCounter(PrimeTable table, PrimeCounter primeCounter) {
            this.table = table;
            this.primeCounter = primeCounter;
        }

        long countForExponents(List<Integer> exponents, long limit) {
            int r = exponents.size();
            if (r == 0)
                return 0;
            if (r == 1) {
                long root = kthRootFloor(limit, exponents.get(0));
                return primeCounter.pi(root);
            }

            int[] tailSum = new int[r];
            int running = 0;
            for (int i = r - 1; i >= 0; i--) {
                running += exponents.get(i);
                tailSum[i] = running;
            }

            Map<Long, Long> memo = new HashMap<>();
            return countDfs(0, 0, limit, r, exponents, tailSum, memo);
        }

        long countDfs(int pos, int startIndex, long remLimit, int r, List<Integer> exponents, int[] tailSum,
                Map<Long, Long> memo) {
            if (pos == r)
                return 1;

            long key = remLimit ^ ((long) startIndex * 0x9E3779B97F4A7C15L) ^ ((long) pos * 0xC2B2AE3D27D4EB4FL);
            if (memo.containsKey(key))
                return memo.get(key);

            long result = 0;
            if (pos == r - 1) {
                int exp = exponents.get(pos);
                long root = kthRootFloor(remLimit, exp);
                long upto = primeCounter.pi(root);
                result = (upto > startIndex) ? (upto - startIndex) : 0;
                memo.put(key, result);
                return result;
            }

            int tail = tailSum[pos];
            long maxP = kthRootFloor(remLimit, tail);
            long upto = primeCounter.pi(maxP);
            if (upto > table.primes.length)
                upto = table.primes.length;

            int exp = exponents.get(pos);
            for (long idx = startIndex; idx < upto; idx++) {
                long p = table.primes[(int) idx];
                long pPow = powCapped(p, exp, remLimit);
                if (pPow > remLimit)
                    break;
                result += countDfs(pos + 1, (int) (idx + 1), remLimit / pPow, r, exponents, tailSum, memo);
            }

            memo.put(key, result);
            return result;
        }
    }

    static long solveFast(long limit) {
        long rootLimit = Math.max(isqrt(limit), 100);
        PrimeTable table = new PrimeTable((int) rootLimit);
        PrimeCounter primeCounter = new PrimeCounter(table);
        OrderedPrimeTupleCounter tupleCounter = new OrderedPrimeTupleCounter(table, primeCounter);

        List<List<Integer>> signatures = generateSignatures(limit, table);
        long total = 0;

        for (List<Integer> signature : signatures) {
            long fsf = fsfFromSignature(signature);
            List<List<Integer>> perms = generateUniquePermutations(signature);
            for (List<Integer> p : perms) {
                long cnt = tupleCounter.countForExponents(p, limit);
                total += fsf * cnt;
            }
        }

        return total;
    }

    static String solve() {
        return Long.toString(solveFast(10000000000L));
    }

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