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

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$.

The code encodes this with

$$\rho(n)=3^{s(n)},\qquad s(n)=\omega_1(n)+\mathbf 1_{v_3(n)\ge2},$$

where $\omega_1(n)$ counts distinct prime divisors congruent to $1\pmod3$.

3. Target Condition $\rho(n)=243$

Because $243=3^5$, we need

$$s(n)=5.$$

This splits the search into two disjoint cases:

1. Case A: exactly five distinct primes $p\equiv1\pmod3$, and $v_3(n)<2$,

2. Case B: exactly four such primes, and $v_3(n)\ge2$.

Those are the only possibilities, because every special contributor adds exactly one factor $3$ to $\rho(n)$.

4. What the Remaining Cofactor Is Allowed to Contain

Once the special contributors have been fixed, the rest of $n$ must not create any new factor $3$ in $\rho(n)$.

So the remaining cofactor may contain only primes $q\equiv2\pmod3$, because such primes contribute $\rho(q^a)=1$.

There is one extra detail:

1. in Case A, a single factor $3$ is allowed, since $v_3(n)=1$ still contributes nothing extra,

2. in Case B, no additional factor $3$ may appear in the cofactor, because the $3^a$ contribution is already handled explicitly by taking $a\ge2$.

5. Why the Search Uses Prime Powers, Not Just Primes

A prime $p\equiv1\pmod3$ contributes the same factor $3$ to $\rho(n)$ no matter whether it appears as $p$, $p^2$, or $p^{17}$. Therefore the DFS must enumerate prime powers, not merely distinct primes.

The same phenomenon occurs for $3$: once the exponent reaches $2$, every larger power $3^a$ still contributes only one factor $3$ to $\rho(n)$.

This is exactly why the recursion loops over

$$p,\;p^2,\;p^3,\dots$$

as long as the global product remains within the limit.

6. Prefix Sums for the “Ordinary” Cofactor

The code precomputes the function

$$F(M)=\sum_{\substack{m\le M\\ p\mid m\Rightarrow p\equiv2\pmod3}} m,$$

the sum of all integers up to $M$ whose prime factors are all $2\pmod3$.

This is stored as a prefix table prefix_allowed_sum_.

Then:

1. in Case B, the terminal contribution is just $F(M)$,

2. in Case A, the cofactor may also include one extra factor $3$, so the code uses

$$F(M)+3F(\lfloor M/3\rfloor).$$

That makes each leaf contribution $O(1)$.

7. DFS Structure and Pruning

The search explores only primes $p\equiv1\pmod3$. To avoid useless branches, the code precomputes min_prod_[r][i], the smallest possible product of $r$ future special primes starting from index $i$.

If even that minimal remaining product already exceeds the limit, the branch is abandoned immediately. This pruning is the reason the search is practical up to $10^{11}$.

8. Worked Structural Examples

A number like

$$n=7\cdot13\cdot19\cdot31\cdot37$$

has exactly five special $1\pmod3$ primes and no $3^2$, so $\rho(n)=3^5=243$ and it belongs to Case A.

A number like

$$n=9\cdot7\cdot13\cdot19\cdot31$$

has four special $1\pmod3$ primes plus $v_3(n)\ge2$, so again $\rho(n)=3^5=243$; this is Case B.

In contrast, multiplying either example by a new prime such as $43\equiv1\pmod3$ would push the root count to $3^6$, which is too large.

9. Validation Logic in the Code

The program checks three things before solving the full problem:

1. it verifies that modulo $91$ there are exactly $9$ cube roots of unity, so the nontrivial count is $8$,

2. for every $n\le2000$, it compares the formula-based root count against a brute-force count of solutions to $x^3\equiv1\pmod n$,

3. up to $5{,}000{,}000$, it compares the fast summation algorithm with a brute-force sum over all $n$ having $s(n)=5$.

These checkpoints validate both the arithmetic formula and the large-scale enumeration strategy.

How the Code Works

special_factor_count computes $s(n)$, the exponent of $3$ in $\rho(n)$.

build_special_primes generates all primes $p\equiv1\pmod3$ that could possibly appear in a valid solution under the $10^{11}$ limit.

build_cofactor_prefix precomputes the prefix sums for cofactors made only of primes $2\pmod3$.

solve_case_without_9 handles Case A, and solve_case_with_9 handles Case B.

Each case uses DFS over special prime powers and aggregates whole blocks of remaining cofactors in constant time via the prefix table.

Complexity Analysis

The algorithm is dominated by the number of reachable DFS nodes after pruning. That search space is vastly smaller than the naive set of all integers up to $10^{11}$.

Terminal aggregation is $O(1)$ because the cofactor sums are precomputed. So, in practice, the runtime is determined almost entirely by the number of viable special-factor branches.

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