Problem 272: Modular Cubes, Part 2

Project Euler Problem 272 Solution

EulerSolve provides an optimized solution for Project Euler Problem 272, Modular Cubes, Part 2, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Let $\rho(n)$ be the number of residues $x \pmod n$ satisfying $$x^3\equiv1\pmod n.$$ The problem asks for the sum of all integers $n\le10^{11}$ having exactly $242$ nontrivial cube roots of unity. Since the trivial root $x=1$ is always present, this is the same as asking for $$\rho(n)=243=3^5.$$ Mathematical Approach 1. Root Counts for Prime Powers The whole method starts from $\rho(p^a)$. For a prime $p\neq3$, the unit group modulo $p^a$ has order $$\varphi(p^a)=p^{a-1}(p-1).$$ The equation $u^3=1$ in a cyclic group has exactly $\gcd(3,\varphi(p^a))$ solutions. Therefore: 1. if $p\equiv2\pmod3$, then $\gcd(3,\varphi(p^a))=1$, so $\rho(p^a)=1$, 2. if $p\equiv1\pmod3$, then $\gcd(3,\varphi(p^a))=3$, so $\rho(p^a)=3$. The prime $p=3$ is special. The code uses the standard fact $$\rho(3)=1,\qquad \rho(3^a)=3\quad(a\ge2).$$ For example, modulo $9$, the three roots are $1,4,7$. 2. The Entire Root Count Depends Only on “Special” Factors By the Chinese Remainder Theorem, $\rho(n)$ is multiplicative across coprime factors. So if $$n=\prod p_i^{a_i},$$ then $$\rho(n)=\prod \rho(p_i^{a_i}).$$ That means the only factors that matter are: 1. distinct primes $p\equiv1\pmod3$, each contributing one factor $3$, regardless of exponent, 2. the prime $3$, but only if its exponent is at least $2$, contributing one factor $3$....

Detailed mathematical approach

Problem Summary

For a positive integer $n$, $C(n)$ is the number of integers $x$ with $1\lt x\lt n$ and $x^3\equiv1\pmod n$, that is, the number of cube roots of $1$ modulo $n$ other than $x=1$ itself. For $n=91$ there are eight of them, namely $9$, $16$, $22$, $29$, $53$, $74$, $79$ and $81$, so $C(91)=8$. We must find the sum of all positive integers $n\le10^{11}$ with $C(n)=242$. ^[1]

The whole solution in one sentence: modulo $n$ the number $1$ has exactly $3^{s}$ cube roots, where $s$ counts the distinct prime factors of $n$ that leave remainder $1$ when divided by $3$, plus one more if $9$ divides $n$, so we list every possible "special part" of $n$ with $s=5$ and add up all of its admissible multiples at once from a table of prefix sums.

Mathematical Approach

Sections 1–6 show that $C(n)$ depends only on which primes divide $n$, and for the prime $3$ on whether $9$ divides $n$; Sections 7–11 turn the sum over all $n\le10^{11}$ into a sum over special parts.

1. Small moduli

Let $\rho(n)$ be the number of residues $x$ modulo $n$ with $x^3\equiv1\pmod n$. For $n\ge2$ the residue $1$ is always a root and $0$ never is, so $C(n)=\rho(n)-1$; for $n=1$ the range $1\lt x\lt1$ is empty and $C(1)=0$. Trying every residue for a few small moduli gives:

modulo $2$, $3$, $4$, $5$, $8$, $16$ and $25$ the only root is $1$;
modulo $7$ the roots are $1,2,4$, modulo $13$ they are $1,3,9$, modulo $19$ they are $1,7,11$, and modulo $49=7^2$ they are $1,18,30$;
modulo $9$ the roots are $1,4,7$ and modulo $27$ they are $1,10,19$, although modulo $3$ only $1$ is a root;
modulo $91=7\cdot13$ there are nine roots: $1$ and the eight numbers of the statement.

Every count is a power of $3$: the primes $7$, $13$ and $19$, which leave remainder $1$ when divided by $3$, give three roots, $2$ and $5$ give one, squaring $7$ changes nothing, and $3$ needs the exponent $2$. Sections 2–5 explain why.

2. Every root is a unit

If $x^3\equiv1\pmod n$, then $x\cdot x^2\equiv1$, so $x^2$ is an inverse of $x$ modulo $n$ and $\gcd(x,n)=1$. The residues coprime to $n$ form a group under multiplication, the group of units $(\mathbb Z/n\mathbb Z)^{\times}$, which has $\varphi(n)$ elements, where $\varphi$ is Euler's totient function. ^[2] Hence $\rho(n)$ is the number of elements $u$ of this group with $u^3=1$.

3. Cube roots of $1$ in a group

Suppose first that a group $G$ with $m$ elements is cyclic, that is, generated by one element $g$: every element equals $g^k$ for exactly one $k$ with $0\le k\lt m$, and $g^j=1$ exactly when $m$ divides $j$. Then $(g^k)^3=1$ means $m\mid3k$. If $3\nmid m$, this forces $m\mid k$, so $k=0$. If $3\mid m$, it means $(m/3)\mid k$, so $k\in\{0,m/3,2m/3\}$. Hence $u^3=1$ has exactly $\gcd(3,m)$ solutions in $G$. ^[3] Modulo $7$, for example, the powers $3^1,\dots,3^6$ are $3,2,6,4,5,1$, so $3$ generates the group of order $6$, and $k=0,2,4$ give the three roots $1,2,4$ of Section 1.

The order of an element $u$ is the least $d\ge1$ with $u^d=1$; modulo $7$, for instance, $2$ has order $3$, since its powers are $2,4,1$. If $u^3=1$, then $d\mid3$, because otherwise the remainder $r$ of $3$ divided by $d$ would satisfy $0\lt r\lt d$ and $u^r=1$. If $G$ is any finite group with $3\nmid|G|$, cyclic or not, then $d$ also divides $|G|$ by Lagrange's theorem. ^[4] So $d=1$, and $u=1$ is the only solution.

4. Prime powers

For an odd prime $p$, the group of units modulo $p^a$ is cyclic of order $\varphi(p^a)=p^{a-1}(p-1)$. ^[2] Section 3 therefore gives:

if $p\equiv1\pmod3$, then $3\mid p-1$ and $\rho(p^a)=3$ for every $a\ge1$;
if $p\equiv2\pmod3$, including $p=2$, the group has $p^{a-1}(p-1)$ elements, a number not divisible by $3$, so $\rho(p^a)=1$ by the second paragraph of Section 3, which also covers the non-cyclic groups modulo $2^a$ with $a\ge3$;
for $p=3$ the order is $2\cdot3^{a-1}$, so $\rho(3)=1$ and $\rho(3^a)=3$ for $a\ge2$.

For $a\ge2$ the three roots modulo $3^a$ are $1+t\cdot3^{a-1}$ with $t=0,1,2$, because

$$\begin{aligned}\bigl(1+t\cdot3^{a-1}\bigr)^3&=1+t\cdot3^{a}\\&\quad+t^2\cdot3^{2a-1}\\&\quad+t^3\cdot3^{3a-3},\end{aligned}$$

and $2a-1\ge a$ and $3a-3\ge a$ when $a\ge2$. For $a=2$ and $a=3$ this gives $1,4,7$ and $1,10,19$, as found in Section 1. So the exponent of a prime $p\ne3$ never matters, and for $p=3$ only whether $9$ divides $n$ matters.

5. Gluing the prime powers together

Write $n=p_1^{a_1}\cdots p_r^{a_r}$. By the Chinese remainder theorem, sending $x\bmod n$ to the tuple of its residues modulo the $p_i^{a_i}$ is a bijection that respects multiplication, ^[5] so $x^3\equiv1\pmod n$ holds exactly when it holds modulo every $p_i^{a_i}$. The roots modulo $n$ correspond to the tuples of roots modulo the prime powers, so $\rho(n)=\prod_i\rho(p_i^{a_i})$, and Section 4 gives

$$\begin{aligned}\rho(n)&=3^{s(n)},\\s(n)&=\omega_1(n)+[\,9\mid n\,],\end{aligned}$$

where $\omega_1(n)$ is the number of distinct primes $p\equiv1\pmod3$ that divide $n$, and $[\,9\mid n\,]$ is $1$ if $9$ divides $n$ and $0$ otherwise. The OEIS entry for the number of solutions of $x^3\equiv1\pmod n$ records the same formula. ^[6] For $91=7\cdot13$, the roots $1,2,4$ modulo $7$ and $1,3,9$ modulo $13$ form $3\cdot3=9$ pairs, and each pair belongs to exactly one root modulo $91$: $1\leftrightarrow(1,1)$, $9\leftrightarrow(2,9)$, $16\leftrightarrow(2,3)$, $22\leftrightarrow(1,9)$, $29\leftrightarrow(1,3)$, $53\leftrightarrow(4,1)$, $74\leftrightarrow(4,9)$, $79\leftrightarrow(2,1)$ and $81\leftrightarrow(4,3)$. Hence $C(91)=9-1=8$, as in the statement.

6. The condition $C(n)=242$

Since $C(1)=0$, the condition $C(n)=242$ means $3^{s(n)}=243=3^5$, that is, $s(n)=5$. This splits into two disjoint cases:

Case A: $9\nmid n$, and $n$ has exactly five distinct prime factors $p\equiv1\pmod3$;
Case B: $9\mid n$, and $n$ has exactly four of them.

All other prime factors of $n$ satisfy $p\equiv2\pmod3$, the prime $2$ included, with any exponents, and in Case A the prime $3$ may also divide $n$ once. For example, $7\cdot13\cdot19\cdot31\cdot37=1983163$ and its doubles and triples belong to Case A, and $9\cdot7\cdot13\cdot19\cdot31=482391$ and $27\cdot7\cdot13\cdot19\cdot31$ belong to Case B, but $9\cdot1983163$ has $s=6$ and therefore $C=728$.

7. Special part and plain cofactor

Call a positive integer plain if all of its prime factors are $\equiv2\pmod3$, and let $\mathcal P$ be the set of plain numbers; the plain numbers up to $20$ are $1,2,4,5,8,10,11,16,17,20$. The full power of a prime $p$ in $n$ is the largest power of $p$ that divides $n$, such as $7^2$ in $98$. Every $n$ with $s(n)=5$ splits in exactly one way as $n=P\cdot m$:

in Case A, the special part $P$ is the product of the full powers of the five primes $p\equiv1\pmod3$ that divide $n$, and the cofactor $m=n/P$ is plain or $3$ times a plain number;
in Case B, $P$ is the full power $3^a$ of $3$ in $n$, with $a\ge2$, times the full powers of the four primes $p\equiv1\pmod3$, and the cofactor $m=n/P$ is plain.

Conversely, multiplying such a $P$ by such an $m$ always gives a number with $s=5$: the cofactor adds no prime $p\equiv1\pmod3$, and it adds at most one factor $3$ in Case A and none in Case B. Let $\mathcal A$ and $\mathcal B$ be the sets of special parts $P\le L$ of the two cases, with $L=10^{11}$, and define

$$F(M)=\sum_{m\in\mathcal P,\ m\le M}m .$$

For a fixed $P$, the plain cofactors are the plain $m\le M_P=\lfloor L/P\rfloor$, with sum $F(M_P)$. In Case A, the cofactors $3m'$ with $m'$ plain and $3m'\le M_P$, that is $m'\le\lfloor M_P/3\rfloor$, add $3F(\lfloor M_P/3\rfloor)$, so there the cofactors sum to $G(M_P)$ with $G(M)=F(M)+3F(\lfloor M/3\rfloor)$. Adding $n=P\cdot m$ over all admissible pairs gives the requested sum:

$$\begin{aligned}S(L)={}&\sum_{P\in\mathcal A}P\,G(M_P)\\&+\sum_{P\in\mathcal B}P\,F(M_P).\end{aligned}$$

8. A worked example

The same method answers the smaller question "which $n\le300$ have $C(n)=8$?", where $s(n)=2$: in Case A the special part has two primes $p\equiv1\pmod3$, and in Case B it is $3^a$ with $a\ge2$ times a single such prime power. The plain numbers up to $4$ are $1,2,4$, so $F(1)=1$, $F(2)=F(3)=3$ and $F(4)=7$.

Case A: $P=91$ has $M_P=3$ and contributes $91\,G(3)=91\bigl(F(3)+3F(1)\bigr)$${}=91\cdot6=546$, the sum of $91$, $182$ and $273$. $P=133=7\cdot19$ has $M_P=2$ and contributes $133\,G(2)=133\cdot3=399$, the sum of $133$ and $266$. $P=217=7\cdot31$, $247=13\cdot19$ and $259=7\cdot37$ have $M_P=1$ and contribute themselves. The total is $1668$.
Case B: $P=63=9\cdot7$ has $M_P=4$ and contributes $63\cdot F(4)=441$, the sum of $63$, $126$ and $252$. $P=117=9\cdot13$ contributes $117\cdot F(2)=351$, the sum of $117$ and $234$. $P=171=9\cdot19$, $189=27\cdot7$ and $279=9\cdot31$ contribute themselves. The total is $1431$.

So there are $16$ such numbers, with sum $1668+1431=3099$, and a direct search over all $n\le300$ finds exactly the same $16$ numbers. Note that $273=3\cdot91$ enters through the single factor $3$ allowed in Case A, while $189=3\cdot63$ is not a plain multiple of $63$: its own special part is $27\cdot7$.

9. How large the pieces can be

The five smallest primes $p\equiv1\pmod3$ are $7$, $13$, $19$, $31$ and $37$. Every special part of Case A is therefore at least $1983163$, so $M_P\le\lfloor10^{11}/1983163\rfloor=50424$, and every special part of Case B is at least $482391$, so $M_P\le\lfloor10^{11}/482391\rfloor=207300$. A table of $F$ up to $207300$ covers every lookup, and the smallest $n$ with $C(n)=242$ is $482391$.

The largest possible prime $p\equiv1\pmod3$ is bounded similarly. In Case B the other three such primes are at least $7$, $13$ and $19$, and the power of $3$ is at least $9$, so $p\le\lfloor10^{11}/15561\rfloor=6426322$ with $15561=9\cdot7\cdot13\cdot19$. In Case A the bound $p\le\lfloor10^{11}/53599\rfloor=1865706$, with $53599=7\cdot13\cdot19\cdot31$, is smaller. Up to $6426322$ there are $219949$ primes $p\equiv1\pmod3$, about half of the $440086$ primes in that range, as the prime number theorem for arithmetic progressions predicts. ^[7]

10. Listing the special parts

A depth-first search builds every special part from its primes $p\equiv1\pmod3$ in increasing order, each with an exponent $e\ge1$: five prime powers with product at most $L'=L$ in Case A, and four with product at most $L'=\lfloor L/9\rfloor$ in Case B, where the power of $3$ is at least $9$. Let $q_0\lt q_1\lt\cdots$ be the primes $p\equiv1\pmod3$, and let $\mu_r(i)=q_iq_{i+1}\cdots q_{i+r-1}$ be the product of $r$ consecutive ones, the smallest product of $r$ distinct such primes that are all at least $q_i$; for $r=0$ it is the empty product $\mu_0(i)=1$.

Suppose that the prime powers chosen so far have product $\Pi$, that $r$ primes are still missing, and that the next one is $q_i$ with exponent $e$. The other $r-1$ primes are larger than $q_i$, so their product is at least $\mu_{r-1}(i+1)$, and a completion exists exactly when

$$\Pi\cdot q_i^{\,e}\cdot\mu_{r-1}(i+1)\le L'.$$

Indeed, if the inequality holds, the next $r-1$ consecutive primes with exponent $1$ complete the product, and if it fails, every completion is even larger. As $i$ grows, $q_i$ increases and $\mu_{r-1}(i+1)$ does not decrease (for $r=1$ it is the constant $1$), so once the test fails with $e=1$, it fails for every later prime, and the loop over $i$ stops. Every node that the search enters therefore has a special part below it, every special part has at most five nodes above it, and every loop step either enters a node or ends the loop, so the work is proportional to the number of terms added. In Case A every special part $P$ yields the term of Section 7; in Case B every product $Q$ of four prime powers is combined with $3^a=9,27,81,\dots$ as long as $Q\cdot3^a\le L$.

11. The table of plain numbers

A number $m\ge2$ is plain exactly when its smallest prime factor $p$ satisfies $p\equiv2\pmod3$ and $m/p$ is plain, because the prime factors of $m$ are $p$ and those of $m/p$. With the smallest prime factor of every number up to $207300$ from a linear sieve, ^[8] one pass from $2$ upwards decides every $m$, and running sums give all values $F(M)$. For instance, $10$ is plain because $2\equiv2\pmod3$ and $5$ is plain, while $14=2\cdot7$ is not, because $7$ is not plain.

How the Code Works

The C++ program ^[9] runs run_validation_checkpoints, builds an Euler272Solver for $L=10^{11}$ (kTargetLimit) and prints the value of its solve. build_spf is the linear sieve of Section 11, and special_factor_count factors $n$ with this table and computes $s(n)$ of Section 5, adding one for every prime $p\equiv1\pmod3$ and one for the prime $3$ when its exponent is at least $2$. cube_root_count_formula returns $3^{s(n)}$, or $0$ for $n=1$; cube_root_count_bruteforce counts the $x$ with $1\le x\lt n$ and $x^3\equiv1\pmod n$ directly, with 128-bit intermediate products; and brute_sum_for_limit adds every $n$ up to a limit with $s(n)=5$.

The constructor of Euler272Solver prepares three tables. build_special_primes runs a sieve of Eratosthenes ^[10] up to $\lfloor L/15561\rfloor$, the bound of Section 9 (constant kPrimeBoundDenominator), and keeps the primes $p\equiv1\pmod3$, starting with $7$, in special_primes_. build_min_product_table fills min_prod_[r][i] with $\mu_r(i)$ of Section 10 for $r=0,\dots,4$, from the end of the list backwards; a product that would exceed 64 bits, or that would need primes beyond the end of the list, is stored as kInfinity. build_cofactor_prefix sizes the plain table with the two bounds of Section 9, marks the plain numbers in the byte array valid by the rule of Section 11 and stores the running sums $F(M)$ in prefix_allowed_sum_. The C++ and Java constant for $7\cdot13\cdot19\cdot31\cdot37=1983163$ (kMinProd5Special, MIN_PROD_5_SPECIAL) reads $1983709$; nothing changes, because the table size $207300$ comes from the larger Case B bound.

sum_allowed_without_3(m) returns $F(m)$ and sum_allowed_with_optional_single_3(m) returns $G(m)=F(m)+3F(\lfloor m/3\rfloor)$; both first clamp $m$ to the table size, a step never needed at this limit. solve_case_without_9 handles Case A. Its outer loop runs over the index $i$ of the smallest prime while $q_i\cdot\mu_4(i+1)\le L$; for every power of $q_i$ up to $\lfloor L/\mu_4(i+1)\rfloor$ it calls dfs_case_without_9(i + 1, 4, power) and adds the results. The recursive function is the search of Section 10, with product as $\Pi$ and remaining as $r$. For each index it reads min_rest, which is $\mu_{r-1}(i+1)$, leaves the loop at the first index that fails, and tries the powers of $q_i$ up to max_prime_power, which is $\lfloor\lfloor L/\Pi\rfloor/\mu_{r-1}(i+1)\rfloor$. The tests are written as divisions, such as product > limit_ / min_rest and power > max_prime_power / p, so no product can overflow. When remaining reaches $0$, the function returns $\Pi\cdot G(\lfloor L/\Pi\rfloor)$.

solve_case_with_9 and dfs_case_with_9 do the same for four primes with the bound special_limit, which is $\lfloor L/9\rfloor$. At a leaf they run through $3^a=9,27,\dots$ while $\Pi\cdot3^a\le L$ and add $\Pi\cdot3^a\cdot F(\lfloor L/(\Pi\cdot3^a)\rfloor)$. Products and table entries are 64-bit integers, the totals are unsigned 128-bit integers, and to_string_u128 converts the final total to decimal.

Before the real computation, run_validation_checkpoints makes three checks: the direct count finds $9$ roots modulo $91$, that is $C(91)=8$; the formula agrees with the direct count for every $n\le2000$; and a solver built for $L=5\cdot10^6$ returns the same sum as brute_sum_for_limit, namely $635\,561\,855$ from $191$ numbers. If a check fails, the program reports it and stops without solving.

The Python file implements the same tables and the two recursive searches as nested functions inside solve. Its integers never overflow, float('inf') plays the role of kInfinity, and the recursion is at most five levels deep; it defines special_factor_count but runs none of the checks. The Java class mirrors the C++ solver with long arithmetic, also without the checks; its top-level method collects the smallest-prime indices of both cases, evaluates each of them separately and adds the results. Every single term is below $3\cdot10^{15}$, and the Java total agrees with the 128-bit C++ total, so nothing overflows at this limit.

Complexity and Verification

Let $K$ be the number of terms that the search adds. The running time is $O(U\log\log U+M+K)$, where $U=6\,426\,322$ is the sieve bound and $M=207\,300$ the size of the plain table: the sieve keeps $219\,949$ primes, the linear sieve and the prefix sums are linear in $M$, and min_prod_ has $5\cdot219\,950$ entries. The search dominates, and by Section 10 its work is proportional to $K$: there are $19\,372\,229$ special parts in Case A, and $21\,438\,638$ products of four prime powers in Case B, which together with their powers of $3$ give $30\,638\,739$ terms. These $K\approx5\cdot10^7$ constant-time terms stand for all $157\,012\,165$ numbers $n\le10^{11}$ with $C(n)=242$, while the search has fewer than $2\cdot10^5$ inner nodes. On one desktop computer the C++ program, compiled with g++ -O2, finishes in about a quarter of a second including its checks, the Java version in under a second, and the Python version in about $15$ seconds. Memory is $O(U+M)$: the largest C++ arrays are the sieve bytes (about $6.4$ MB, released before min_prod_ is built), min_prod_ (about $8.8$ MB) and, during the checkpoint up to $5\cdot10^6$, a $20$ MB smallest-prime-factor table; the measured peaks were about $43$ MB for C++, $86$ MB for Java and $98$ MB for Python.

For this article we also checked the following. The formula $\rho(n)=3^{s(n)}$ agrees with a direct count for every $n\le5000$, and the roots modulo $91$ with their residue pairs, the worked example of Section 8 and the value $C(482391)=242$ were confirmed by direct search. Copies of the C++ solver compiled with other limits agree with a segmented sieve that computes $s(n)$ for every $n$: the sum is $1\,432\,459\,680\,651$ from $24\,038$ numbers for $L=10^8$, $316\,174\,709\,149\,021$ from $555\,292$ numbers for $L=10^9$, and $55\,516\,771\,054\,905\,258$ from $10\,049\,916$ numbers for $L=10^{10}$. At $L=10^{11}$ the three published programs, an independent implementation that builds the plain table by crossing out the multiples of all other primes, and a direct sieve over every $n\le10^{11}$ give the same total. These checks are not part of the published solutions.

Footnotes and References

The references below supply the background facts used above; the decomposition into special parts and plain cofactors, the pruning test and the operation counts are derived in this article.

Project Euler 272 — Modular Cubes, Part 2. The official statement defines C(n) as the number of integers x with 1 < x < n and x³ ≡ 1 (mod n), gives C(91) = 8 together with the eight values of x, and asks for the sum of all positive n ≤ 10¹¹ with C(n) = 242.
Wikipedia — Multiplicative group of integers modulo n. The residues coprime to n form a group of order φ(n). It is cyclic exactly when n is 1, 2, 4, pᵏ or 2pᵏ for an odd prime p; in particular it is cyclic modulo every power of an odd prime and not cyclic modulo 2ᵏ for k ≥ 3.
Wikipedia — Cyclic group. A cyclic group of order n has exactly one subgroup of order d for every divisor d of n, and exactly d elements whose order divides d; for d = 3 this is the count of three cube roots of 1 in Section 3.
Wikipedia — Lagrange's theorem (group theory). The order of a subgroup of a finite group divides the order of the group; in particular the order of every element divides the order of the group, which settles the case 3 ∤ |G| in Section 3.
Wikipedia — Chinese remainder theorem. For pairwise coprime moduli n₁, …, nₖ with product N, the map sending x mod N to (x mod n₁, …, x mod nₖ) is a ring isomorphism, so a congruence modulo N holds exactly when it holds modulo every factor, as used in Section 5.
OEIS A060839 — Number of solutions to x^3 == 1 (mod n). The entry states that the sequence is multiplicative with a(3) = 1, a(3ᵉ) = 3 for e ≥ 2, a(pᵉ) = 3 for primes p ≡ 1 (mod 3) and a(pᵉ) = 1 for primes p ≡ 2 (mod 3), which is the formula of Section 5.
Wikipedia — Prime number theorem, section "Prime number theorem for arithmetic progressions". De la Vallée Poussin proved, as Dirichlet and Legendre had conjectured, that for coprime a and d the number of primes p ≤ x with p ≡ a (mod d) is asymptotic to Li(x)/φ(d); for d = 3 each of the two classes receives about half of the primes.
cp-algorithms — Linear Sieve. The linear sieve computes the least prime factor of every number up to n in O(n) time, generating each composite exactly once from its least prime factor; build_spf in C++ and Python and buildSpf in Java follow this scheme.
C++ — Euler272.cpp. The linked immutable C++ revision contains the solver class, the two searches and the checks that run before the main computation. It is the implementation record for this article.
Wikipedia — Sieve of Eratosthenes. The sieve crosses out the multiples of each prime p, starting at p², and runs in O(n log log n) time; build_special_primes uses it for the primes up to ⌊L/15561⌋.

Problem 272 source code

C++

#include <algorithm>
#include <atomic>
#include <cstdint>
#include <iostream>
#include <limits>
#include <string>
#include <thread>
#include <vector>

namespace {

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

constexpr u64 kTargetLimit = 100'000'000'000ULL;
constexpr int kValidationRootLimit = 2'000;
constexpr int kValidationSumLimit = 5'000'000;

constexpr u64 kMinProd4Special = 53'599ULL;   // 7*13*19*31
constexpr u64 kMinProd5Special = 1'983'709ULL; // 7*13*19*31*37
constexpr u64 kPrimeBoundDenominator = 15'561ULL; // 9*7*13*19
constexpr u64 kInfinity = std::numeric_limits<u64>::max();

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

    std::string digits;
    while (value > 0) {
        const int digit = static_cast<int>(value % 10);
        digits.push_back(static_cast<char>('0' + digit));
        value /= 10;
    }

    std::reverse(digits.begin(), digits.end());
    return digits;
}

std::vector<int> build_spf(int limit) {
    std::vector<int> spf(limit + 1, 0);
    std::vector<int> primes;
    if (limit >= 1) {
        spf[1] = 1;
    }

    for (int i = 2; i <= limit; ++i) {
        if (spf[i] == 0) {
            spf[i] = i;
            primes.push_back(i);
        }

        for (const int p : primes) {
            const std::int64_t composite = static_cast<std::int64_t>(i) * p;
            if (composite > limit || p > spf[i]) {
                break;
            }
            spf[static_cast<std::size_t>(composite)] = p;
        }
    }

    return spf;
}

int special_factor_count(int n, const std::vector<int>& spf) {
    int count = 0;
    while (n > 1) {
        const int p = spf[n];
        int exponent = 0;
        while (n % p == 0) {
            n /= p;
            ++exponent;
        }

        if (p == 3) {
            if (exponent >= 2) {
                ++count;
            }
        } else if (p % 3 == 1) {
            ++count;
        }
    }
    return count;
}

u64 cube_root_count_formula(int n, const std::vector<int>& spf) {
    if (n == 1) {
        return 0;
    }

    const int special = special_factor_count(n, spf);
    u64 roots = 1;
    for (int i = 0; i < special; ++i) {
        roots *= 3;
    }
    return roots;
}

u64 cube_root_count_bruteforce(u64 n) {
    u64 count = 0;
    for (u64 x = 1; x < n; ++x) {
        const u64 x2 = static_cast<u64>((static_cast<u128>(x) * x) % n);
        const u64 x3 = static_cast<u64>((static_cast<u128>(x2) * x) % n);
        if (x3 == 1) {
            ++count;
        }
    }
    return count;
}

u128 brute_sum_for_limit(int limit) {
    const std::vector<int> spf = build_spf(limit);
    u128 sum = 0;

    for (int n = 1; n <= limit; ++n) {
        if (special_factor_count(n, spf) == 5) {
            sum += static_cast<u64>(n);
        }
    }

    return sum;
}

class Euler272Solver {
public:
    explicit Euler272Solver(u64 limit) : limit_(limit) {
        build_special_primes();
        build_min_product_table();
        build_cofactor_prefix();
    }

    u128 solve(bool allow_multithreading, unsigned requested_threads = 0) const {
        unsigned threads = requested_threads;
        if (threads == 0) {
            threads = std::thread::hardware_concurrency();
            if (threads == 0) {
                threads = 1;
            }
        }

        if (!allow_multithreading) {
            threads = 1;
        }

        const u128 case_without_9 = solve_case_without_9(threads);
        const u128 case_with_9 = solve_case_with_9(threads);
        return case_without_9 + case_with_9;
    }

private:
    u64 limit_ = 0;
    u64 max_cofactor_limit_ = 0;

    std::vector<u64> special_primes_;          // primes p == 1 (mod 3)
    std::vector<std::vector<u64>> min_prod_;   // min_prod_[r][i] = product of r consecutive primes from i
    std::vector<u64> prefix_allowed_sum_;      // F1(M): sum of allowed cofactors <= M

    void build_special_primes() {
        special_primes_.clear();

        const u64 upper_u64 = limit_ / kPrimeBoundDenominator;
        if (upper_u64 < 7) {
            return;
        }

        const int upper = static_cast<int>(upper_u64);
        std::vector<std::uint8_t> is_prime(static_cast<std::size_t>(upper) + 1, 1);
        is_prime[0] = 0;
        is_prime[1] = 0;

        for (int p = 2; static_cast<std::int64_t>(p) * p <= upper; ++p) {
            if (!is_prime[p]) {
                continue;
            }
            for (int multiple = p * p; multiple <= upper; multiple += p) {
                is_prime[multiple] = 0;
            }
        }

        for (int p = 7; p <= upper; ++p) {
            if (is_prime[p] && p % 3 == 1) {
                special_primes_.push_back(static_cast<u64>(p));
            }
        }
    }

    void build_min_product_table() {
        const int n = static_cast<int>(special_primes_.size());
        min_prod_.assign(5, std::vector<u64>(static_cast<std::size_t>(n) + 1, kInfinity));

        for (int i = 0; i <= n; ++i) {
            min_prod_[0][i] = 1;
        }

        for (int r = 1; r <= 4; ++r) {
            min_prod_[r][n] = kInfinity;
        }

        for (int i = n - 1; i >= 0; --i) {
            for (int r = 1; r <= 4; ++r) {
                const u64 tail = min_prod_[r - 1][i + 1];
                if (tail == kInfinity || special_primes_[i] > kInfinity / tail) {
                    min_prod_[r][i] = kInfinity;
                } else {
                    min_prod_[r][i] = special_primes_[i] * tail;
                }
            }
        }
    }

    void build_cofactor_prefix() {
        const u64 max_case_without_9 = limit_ / kMinProd5Special;
        const u64 max_case_with_9 = limit_ / (9ULL * kMinProd4Special);
        max_cofactor_limit_ = std::max(max_case_without_9, max_case_with_9);

        prefix_allowed_sum_.assign(static_cast<std::size_t>(max_cofactor_limit_) + 1, 0);
        if (max_cofactor_limit_ == 0) {
            return;
        }

        const int max_limit = static_cast<int>(max_cofactor_limit_);
        const std::vector<int> spf = build_spf(max_limit);

        std::vector<std::uint8_t> valid(static_cast<std::size_t>(max_limit) + 1, 0);
        valid[1] = 1;

        for (int n = 2; n <= max_limit; ++n) {
            const int p = spf[n];
            valid[n] = static_cast<std::uint8_t>(valid[n / p] && (p % 3 == 2));
        }

        for (u64 n = 1; n <= max_cofactor_limit_; ++n) {
            prefix_allowed_sum_[n] = prefix_allowed_sum_[n - 1] + (valid[n] ? n : 0);
        }
    }

    u64 sum_allowed_without_3(u64 m) const {
        if (m > max_cofactor_limit_) {
            m = max_cofactor_limit_;
        }
        return prefix_allowed_sum_[m];
    }

    u128 sum_allowed_with_optional_single_3(u64 m) const {
        if (m > max_cofactor_limit_) {
            m = max_cofactor_limit_;
        }

        const u64 no_three = sum_allowed_without_3(m);
        const u64 one_three = 3ULL * sum_allowed_without_3(m / 3);
        return static_cast<u128>(no_three) + one_three;
    }

    u128 dfs_case_without_9(int start, int remaining, u64 product) const {
        if (remaining == 0) {
            const u64 cofactor_limit = limit_ / product;
            return static_cast<u128>(product) * sum_allowed_with_optional_single_3(cofactor_limit);
        }

        const int n = static_cast<int>(special_primes_.size());
        u128 total = 0;

        for (int i = start; i < n; ++i) {
            const u64 min_rest = min_prod_[remaining - 1][i + 1];
            if (min_rest == kInfinity) {
                break;
            }
            if (product > limit_ / min_rest) {
                break;
            }

            const u64 max_prime_power = (limit_ / product) / min_rest;
            const u64 p = special_primes_[i];
            if (p > max_prime_power) {
                break;
            }

            for (u64 power = p; power <= max_prime_power;) {
                total += dfs_case_without_9(i + 1, remaining - 1, product * power);
                if (power > max_prime_power / p) {
                    break;
                }
                power *= p;
            }
        }

        return total;
    }

    u128 dfs_case_with_9(int start, int remaining, u64 product, u64 special_limit) const {
        if (remaining == 0) {
            u128 total = 0;
            const u64 max_three_power = limit_ / product;

            for (u64 three_power = 9; three_power <= max_three_power;) {
                const u64 cofactor_limit = limit_ / (product * three_power);
                total += static_cast<u128>(product) * three_power *
                         sum_allowed_without_3(cofactor_limit);
                if (three_power > max_three_power / 3) {
                    break;
                }
                three_power *= 3;
            }
            return total;
        }

        const int n = static_cast<int>(special_primes_.size());
        u128 total = 0;

        for (int i = start; i < n; ++i) {
            const u64 min_rest = min_prod_[remaining - 1][i + 1];
            if (min_rest == kInfinity) {
                break;
            }
            if (product > special_limit / min_rest) {
                break;
            }

            const u64 max_prime_power = (special_limit / product) / min_rest;
            const u64 p = special_primes_[i];
            if (p > max_prime_power) {
                break;
            }

            for (u64 power = p; power <= max_prime_power;) {
                total += dfs_case_with_9(i + 1, remaining - 1, product * power, special_limit);
                if (power > max_prime_power / p) {
                    break;
                }
                power *= p;
            }
        }

        return total;
    }

    u128 solve_case_without_9(unsigned threads) const {
        if (special_primes_.empty()) {
            return 0;
        }

        std::vector<int> roots;
        roots.reserve(special_primes_.size());

        for (int i = 0; i < static_cast<int>(special_primes_.size()); ++i) {
            const u64 min_rest = min_prod_[4][i + 1];
            if (min_rest == kInfinity) {
                break;
            }
            if (special_primes_[i] > limit_ / min_rest) {
                break;
            }
            roots.push_back(i);
        }

        if (roots.empty()) {
            return 0;
        }

        auto solve_root = [&](int root_index) {
            const u64 p = special_primes_[root_index];
            const u64 min_rest = min_prod_[4][root_index + 1];
            const u64 max_prime_power = limit_ / min_rest;

            u128 local = 0;
            for (u64 power = p; power <= max_prime_power;) {
                local += dfs_case_without_9(root_index + 1, 4, power);
                if (power > max_prime_power / p) {
                    break;
                }
                power *= p;
            }
            return local;
        };

        if (threads <= 1 || roots.size() < 2) {
            u128 total = 0;
            for (const int root : roots) {
                total += solve_root(root);
            }
            return total;
        }

        threads = std::min<unsigned>(threads, static_cast<unsigned>(roots.size()));

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

        for (unsigned tid = 0; tid < threads; ++tid) {
            workers.emplace_back([&, tid]() {
                u128 local = 0;
                while (true) {
                    const std::size_t index = next.fetch_add(1, std::memory_order_relaxed);
                    if (index >= roots.size()) {
                        break;
                    }
                    local += solve_root(roots[index]);
                }
                partial[tid] = local;
            });
        }

        for (std::thread& worker : workers) {
            worker.join();
        }

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

    u128 solve_case_with_9(unsigned threads) const {
        if (limit_ < 9 || special_primes_.empty()) {
            return 0;
        }

        const u64 special_limit = limit_ / 9;

        std::vector<int> roots;
        roots.reserve(special_primes_.size());

        for (int i = 0; i < static_cast<int>(special_primes_.size()); ++i) {
            const u64 min_rest = min_prod_[3][i + 1];
            if (min_rest == kInfinity) {
                break;
            }
            if (special_primes_[i] > special_limit / min_rest) {
                break;
            }
            roots.push_back(i);
        }

        if (roots.empty()) {
            return 0;
        }

        auto solve_root = [&](int root_index) {
            const u64 p = special_primes_[root_index];
            const u64 min_rest = min_prod_[3][root_index + 1];
            const u64 max_prime_power = special_limit / min_rest;

            u128 local = 0;
            for (u64 power = p; power <= max_prime_power;) {
                local += dfs_case_with_9(root_index + 1, 3, power, special_limit);
                if (power > max_prime_power / p) {
                    break;
                }
                power *= p;
            }
            return local;
        };

        if (threads <= 1 || roots.size() < 2) {
            u128 total = 0;
            for (const int root : roots) {
                total += solve_root(root);
            }
            return total;
        }

        threads = std::min<unsigned>(threads, static_cast<unsigned>(roots.size()));

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

        for (unsigned tid = 0; tid < threads; ++tid) {
            workers.emplace_back([&, tid]() {
                u128 local = 0;
                while (true) {
                    const std::size_t index = next.fetch_add(1, std::memory_order_relaxed);
                    if (index >= roots.size()) {
                        break;
                    }
                    local += solve_root(roots[index]);
                }
                partial[tid] = local;
            });
        }

        for (std::thread& worker : workers) {
            worker.join();
        }

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

bool run_validation_checkpoints() {
    {
        const u64 roots_mod_91 = cube_root_count_bruteforce(91);
        if (roots_mod_91 != 9) {
            std::cerr << "Checkpoint failed for n=91: got C(91)=" << (roots_mod_91 - 1)
                      << ", expected 8\n";
            return false;
        }
    }

    {
        const std::vector<int> spf = build_spf(kValidationRootLimit);
        for (int n = 1; n <= kValidationRootLimit; ++n) {
            const u64 formula = cube_root_count_formula(n, spf);
            const u64 brute = cube_root_count_bruteforce(static_cast<u64>(n));
            if (formula != brute) {
                std::cerr << "Root-count checkpoint failed for n=" << n << ": got "
                          << formula << ", expected " << brute << "\n";
                return false;
            }
        }
    }

    {
        const u128 brute = brute_sum_for_limit(kValidationSumLimit);
        const Euler272Solver solver(kValidationSumLimit);
        const u128 fast = solver.solve(false, 1);

        if (fast != brute) {
            std::cerr << "Sum checkpoint failed for n<=" << kValidationSumLimit << ": got "
                      << to_string_u128(fast) << ", expected " << to_string_u128(brute)
                      << "\n";
            return false;
        }
    }

    return true;
}

} // namespace

int main() {
    if (!run_validation_checkpoints()) {
        return 1;
    }

    unsigned threads = std::thread::hardware_concurrency();
    if (threads == 0) {
        threads = 4;
    }

    const Euler272Solver solver(kTargetLimit);
    const u128 answer = solver.solve(true, threads);

    std::cout << to_string_u128(answer) << '\n';
    return 0;
}

Python

import math

def solve():
    LIMIT = 100_000_000_000

    def build_spf(limit):
        spf = list(range(limit + 1))
        primes = []
        spf[0] = 0
        if limit >= 1:
            spf[1] = 1
        for i in range(2, limit + 1):
            if spf[i] == i:
                primes.append(i)
            for p in primes:
                if p * i > limit or p > spf[i]:
                    break
                spf[p * i] = p
        return spf

    def special_factor_count(n, spf):
        count = 0
        while n > 1:
            p = spf[n]
            e = 0
            while n % p == 0:
                n //= p
                e += 1
            if p == 3:
                if e >= 2:
                    count += 1
            elif p % 3 == 1:
                count += 1
        return count

    # Sieve special (p ≡ 1 mod 3) primes
    PRIME_BOUND_DENOM = 15561  # 9*7*13*19
    upper = LIMIT // PRIME_BOUND_DENOM
    sp_sieve = bytearray(b'\x01' * (upper + 1))
    sp_sieve[0] = 0
    sp_sieve[1] = 0
    for p in range(2, math.isqrt(upper) + 1):
        if sp_sieve[p]:
            sp_sieve[p*p::p] = bytearray(len(sp_sieve[p*p::p]))
    special_primes = [p for p in range(7, upper + 1) if sp_sieve[p] and p % 3 == 1]

    INF = float('inf')
    n_sp = len(special_primes)
    min_prod = [[INF]*(n_sp + 1) for _ in range(5)]
    for i in range(n_sp + 1):
        min_prod[0][i] = 1
    for i in range(n_sp - 1, -1, -1):
        for r in range(1, 5):
            tail = min_prod[r-1][i+1]
            if tail == INF or special_primes[i] > LIMIT:
                min_prod[r][i] = INF
            else:
                min_prod[r][i] = special_primes[i] * tail

    # Build cofactor prefix sums
    MIN_PROD5 = 7 * 13 * 19 * 31 * 37
    MIN_PROD4 = 7 * 13 * 19 * 31
    max_cf1 = LIMIT // MIN_PROD5 if MIN_PROD5 <= LIMIT else 0
    max_cf2 = LIMIT // (9 * MIN_PROD4) if 9 * MIN_PROD4 <= LIMIT else 0
    max_cofactor = max(max_cf1, max_cf2)

    spf = build_spf(max_cofactor)
    valid = bytearray(max_cofactor + 1)
    valid[1] = 1
    for nn in range(2, max_cofactor + 1):
        p = spf[nn]
        valid[nn] = 1 if (valid[nn // p] and p % 3 == 2) else 0

    prefix = [0] * (max_cofactor + 2)
    for nn in range(1, max_cofactor + 1):
        prefix[nn] = prefix[nn - 1] + (nn if valid[nn] else 0)

    def sum_allowed(m):
        if m > max_cofactor:
            m = max_cofactor
        return prefix[m]

    def sum_allowed_opt3(m):
        if m > max_cofactor:
            m = max_cofactor
        return sum_allowed(m) + 3 * sum_allowed(m // 3)

    def dfs_without_9(start, remaining, product):
        if remaining == 0:
            cf_limit = LIMIT // product
            return product * sum_allowed_opt3(cf_limit)
        total = 0
        for i in range(start, n_sp):
            mr = min_prod[remaining - 1][i + 1]
            if mr == INF:
                break
            if product > LIMIT // mr:
                break
            mp = (LIMIT // product) // mr
            p = special_primes[i]
            if p > mp:
                break
            power = p
            while power <= mp:
                total += dfs_without_9(i + 1, remaining - 1, product * power)
                if power > mp // p:
                    break
                power *= p
        return total

    def dfs_with_9(start, remaining, product, sl):
        if remaining == 0:
            total = 0
            mtp = LIMIT // product
            tp = 9
            while tp <= mtp:
                cf_limit = LIMIT // (product * tp)
                total += product * tp * sum_allowed(cf_limit)
                if tp > mtp // 3:
                    break
                tp *= 3
            return total
        total = 0
        for i in range(start, n_sp):
            mr = min_prod[remaining - 1][i + 1]
            if mr == INF:
                break
            if product > sl // mr:
                break
            mp = (sl // product) // mr
            p = special_primes[i]
            if p > mp:
                break
            power = p
            while power <= mp:
                total += dfs_with_9(i + 1, remaining - 1, product * power, sl)
                if power > mp // p:
                    break
                power *= p
        return total

    case1 = 0
    for i in range(n_sp):
        mr = min_prod[4][i + 1]
        if mr == INF:
            break
        if special_primes[i] > LIMIT // mr:
            break
        p = special_primes[i]
        mp = LIMIT // mr
        power = p
        while power <= mp:
            case1 += dfs_without_9(i + 1, 4, power)
            if power > mp // p:
                break
            power *= p

    sl = LIMIT // 9
    case2 = 0
    for i in range(n_sp):
        mr = min_prod[3][i + 1]
        if mr == INF:
            break
        if special_primes[i] > sl // mr:
            break
        p = special_primes[i]
        mp = sl // mr
        power = p
        while power <= mp:
            case2 += dfs_with_9(i + 1, 3, power, sl)
            if power > mp // p:
                break
            power *= p

    return str(case1 + case2)

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

Java

import java.util.*;
import java.util.concurrent.*;

public class Euler272 {
    static final long TARGET_LIMIT = 100_000_000_000L;
    static final long MIN_PROD_4_SPECIAL = 53_599L;
    static final long MIN_PROD_5_SPECIAL = 1_983_709L;
    static final long PRIME_BOUND_DENOMINATOR = 15_561L;
    static final long INFINITY = Long.MAX_VALUE;

    static int[] buildSpf(int limit) {
        int[] spf = new int[limit + 1];
        if (limit >= 1)
            spf[1] = 1;
        List<Integer> primes = new ArrayList<>();

        for (int i = 2; i <= limit; ++i) {
            if (spf[i] == 0) {
                spf[i] = i;
                primes.add(i);
            }
            for (int p : primes) {
                long composite = (long) i * p;
                if (composite > limit || p > spf[i])
                    break;
                spf[(int) composite] = p;
            }
        }
        return spf;
    }

    static class Solver {
        long limit;
        long maxCofactorLimit;
        List<Long> specialPrimes = new ArrayList<>();
        long[][] minProd;
        long[] prefixAllowedSum;

        Solver(long limit) {
            this.limit = limit;
            buildSpecialPrimes();
            buildMinProductTable();
            buildCofactorPrefix();
        }

        void buildSpecialPrimes() {
            long upperU64 = limit / PRIME_BOUND_DENOMINATOR;
            if (upperU64 < 7)
                return;

            int upper = (int) upperU64;
            boolean[] isPrime = new boolean[upper + 1];
            Arrays.fill(isPrime, true);
            isPrime[0] = false;
            isPrime[1] = false;

            for (int p = 2; (long) p * p <= upper; ++p) {
                if (isPrime[p]) {
                    for (int multiple = p * p; multiple <= upper; multiple += p) {
                        isPrime[multiple] = false;
                    }
                }
            }

            for (int p = 7; p <= upper; ++p) {
                if (isPrime[p] && p % 3 == 1) {
                    specialPrimes.add((long) p);
                }
            }
        }

        void buildMinProductTable() {
            int n = specialPrimes.size();
            minProd = new long[5][n + 1];
            for (int r = 0; r < 5; ++r) {
                Arrays.fill(minProd[r], INFINITY);
            }
            for (int i = 0; i <= n; ++i) {
                minProd[0][i] = 1;
            }
            for (int i = n - 1; i >= 0; --i) {
                for (int r = 1; r <= 4; ++r) {
                    long tail = minProd[r - 1][i + 1];
                    if (tail == INFINITY || specialPrimes.get(i) > INFINITY / tail) {
                        minProd[r][i] = INFINITY;
                    } else {
                        minProd[r][i] = specialPrimes.get(i) * tail;
                    }
                }
            }
        }

        void buildCofactorPrefix() {
            long maxCaseWithout9 = limit / MIN_PROD_5_SPECIAL;
            long maxCaseWith9 = limit / (9L * MIN_PROD_4_SPECIAL);
            maxCofactorLimit = Math.max(maxCaseWithout9, maxCaseWith9);

            if (maxCofactorLimit == 0) {
                prefixAllowedSum = new long[1];
                return;
            }

            int maxLimit = (int) maxCofactorLimit;
            int[] spf = buildSpf(maxLimit);
            boolean[] valid = new boolean[maxLimit + 1];
            valid[1] = true;

            for (int n = 2; n <= maxLimit; ++n) {
                int p = spf[n];
                valid[n] = valid[n / p] && (p % 3 == 2);
            }

            prefixAllowedSum = new long[maxLimit + 1];
            for (int n = 1; n <= maxLimit; ++n) {
                prefixAllowedSum[n] = prefixAllowedSum[n - 1] + (valid[n] ? n : 0);
            }
        }

        long sumAllowedWithout3(long m) {
            if (m > maxCofactorLimit)
                m = maxCofactorLimit;
            return prefixAllowedSum[(int) m];
        }

        long sumAllowedWithOptionalSingle3(long m) {
            if (m > maxCofactorLimit)
                m = maxCofactorLimit;
            long noThree = sumAllowedWithout3(m);
            long oneThree = 3L * sumAllowedWithout3(m / 3);
            return noThree + oneThree;
        }

        long dfsCaseWithout9(int start, int remaining, long product) {
            if (remaining == 0) {
                long cofactorLimit = limit / product;
                return product * sumAllowedWithOptionalSingle3(cofactorLimit);
            }

            long total = 0;
            int n = specialPrimes.size();

            for (int i = start; i < n; ++i) {
                long minRest = minProd[remaining - 1][i + 1];
                if (minRest == INFINITY)
                    break;
                if (product > limit / minRest)
                    break;

                long maxPrimePower = (limit / product) / minRest;
                long p = specialPrimes.get(i);
                if (p > maxPrimePower)
                    break;

                for (long power = p; power <= maxPrimePower;) {
                    total += dfsCaseWithout9(i + 1, remaining - 1, product * power);
                    if (power > maxPrimePower / p)
                        break;
                    power *= p;
                }
            }

            return total;
        }

        long dfsCaseWith9(int start, int remaining, long product, long specialLimit) {
            if (remaining == 0) {
                long total = 0;
                long maxThreePower = limit / product;

                for (long threePower = 9; threePower <= maxThreePower;) {
                    long cofactorLimit = limit / (product * threePower);
                    total += product * threePower * sumAllowedWithout3(cofactorLimit);
                    if (threePower > maxThreePower / 3)
                        break;
                    threePower *= 3;
                }
                return total;
            }

            long total = 0;
            int n = specialPrimes.size();

            for (int i = start; i < n; ++i) {
                long minRest = minProd[remaining - 1][i + 1];
                if (minRest == INFINITY)
                    break;
                if (product > specialLimit / minRest)
                    break;

                long maxPrimePower = (specialLimit / product) / minRest;
                long p = specialPrimes.get(i);
                if (p > maxPrimePower)
                    break;

                for (long power = p; power <= maxPrimePower;) {
                    total += dfsCaseWith9(i + 1, remaining - 1, product * power, specialLimit);
                    if (power > maxPrimePower / p)
                        break;
                    power *= p;
                }
            }

            return total;
        }

        long solveParallel() {
            List<Integer> rootsWithout = new ArrayList<>();
            for (int i = 0; i < specialPrimes.size(); ++i) {
                long minRest = minProd[4][i + 1];
                if (minRest == INFINITY)
                    break;
                if (specialPrimes.get(i) > limit / minRest)
                    break;
                rootsWithout.add(i);
            }

            long specialLimit = limit / 9;
            List<Integer> rootsWith = new ArrayList<>();
            for (int i = 0; i < specialPrimes.size(); ++i) {
                long minRest = minProd[3][i + 1];
                if (minRest == INFINITY)
                    break;
                if (specialPrimes.get(i) > specialLimit / minRest)
                    break;
                rootsWith.add(i);
            }

            int threads = Math.max(1, Runtime.getRuntime().availableProcessors());
            ExecutorService executor = Executors.newFixedThreadPool(threads);
            List<Future<Long>> futuresWithout = new ArrayList<>();
            List<Future<Long>> futuresWith = new ArrayList<>();

            for (int root : rootsWithout) {
                futuresWithout.add(executor.submit(() -> {
                    long p = specialPrimes.get(root);
                    long minRest = minProd[4][root + 1];
                    long maxPrimePower = limit / minRest;
                    long local = 0;
                    for (long power = p; power <= maxPrimePower;) {
                        local += dfsCaseWithout9(root + 1, 4, power);
                        if (power > maxPrimePower / p)
                            break;
                        power *= p;
                    }
                    return local;
                }));
            }

            for (int root : rootsWith) {
                futuresWith.add(executor.submit(() -> {
                    long p = specialPrimes.get(root);
                    long minRest = minProd[3][root + 1];
                    long maxPrimePower = specialLimit / minRest;
                    long local = 0;
                    for (long power = p; power <= maxPrimePower;) {
                        local += dfsCaseWith9(root + 1, 3, power, specialLimit);
                        if (power > maxPrimePower / p)
                            break;
                        power *= p;
                    }
                    return local;
                }));
            }

            long total = 0;
            try {
                for (Future<Long> f : futuresWithout)
                    total += f.get();
                for (Future<Long> f : futuresWith)
                    total += f.get();
            } catch (Exception e) {
                e.printStackTrace();
            }

            executor.shutdown();
            return total;
        }
    }

    public static String solve() {
        Solver solver = new Solver(TARGET_LIMIT);
        return String.valueOf(solver.solveParallel());
    }

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