Problem 152: Writing 1/2 as a Sum of Inverse Squares

Project Euler Problem 152 Solution

EulerSolve provides an optimized solution for Project Euler Problem 152, Writing 1/2 as a Sum of Inverse Squares, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Project Euler 152 asks for the number of subsets $S \subseteq \{2,3,\dots,80\}$ such that $$\sum_{n \in S}\frac{1}{n^2}=\frac12.$$ A brute-force search would inspect $2^{79}$ subsets, so the real task is to discover arithmetic structure that forbids almost all candidates before any global enumeration begins. The implementations do exactly that: they use local congruence conditions coming from large odd prime powers, turn the surviving rational equation into an integer subset-sum problem, combine the constrained prime blocks by convolution, and finish the remaining search with meet-in-the-middle. Mathematical Approach The key point is that the identity $\sum 1/n^2 = 1/2$ is rigid modulo $p^2$ for carefully chosen odd primes. Once those local conditions are enforced, the set of viable denominators becomes small enough that the exact count can be obtained with a structured subset-sum computation. Congruence Conditions from Maximal Odd Prime Powers Fix an odd prime $p$, and let $P_p=p^e$ be the largest power of $p$ not exceeding 80. For the actual problem, examples are $$P_3=27,\qquad P_5=25,\qquad P_7=49,\qquad P_{13}=13.$$ Now define the top group $$G_p=\{n\in\{2,\dots,80\}: P_p \mid n\}.$$ Let $M=\operatorname{lcm}(2,3,\dots,80)$. Since $v_p(M)=e$, we may write $M=P_pM_p$ with $p\nmid M_p$....

Detailed mathematical approach

Problem Summary

Project Euler 152 asks for the number of subsets $S \subseteq \{2,3,\dots,80\}$ such that

$$\sum_{n \in S}\frac{1}{n^2}=\frac12.$$

A brute-force search would inspect $2^{79}$ subsets, so the real task is to discover arithmetic structure that forbids almost all candidates before any global enumeration begins. The implementations do exactly that: they use local congruence conditions coming from large odd prime powers, turn the surviving rational equation into an integer subset-sum problem, combine the constrained prime blocks by convolution, and finish the remaining search with meet-in-the-middle.

Mathematical Approach

The key point is that the identity $\sum 1/n^2 = 1/2$ is rigid modulo $p^2$ for carefully chosen odd primes. Once those local conditions are enforced, the set of viable denominators becomes small enough that the exact count can be obtained with a structured subset-sum computation.

Congruence Conditions from Maximal Odd Prime Powers

Fix an odd prime $p$, and let $P_p=p^e$ be the largest power of $p$ not exceeding 80. For the actual problem, examples are

$$P_3=27,\qquad P_5=25,\qquad P_7=49,\qquad P_{13}=13.$$

Now define the top group

$$G_p=\{n\in\{2,\dots,80\}: P_p \mid n\}.$$

Let $M=\operatorname{lcm}(2,3,\dots,80)$. Since $v_p(M)=e$, we may write $M=P_pM_p$ with $p\nmid M_p$. Multiplying the target equation by $2M^2$ gives

$$\sum_{n\in S} 2\left(\frac{M}{n}\right)^2=M^2.$$

Reduce this congruence modulo $p^2$. If $n\notin G_p$, then $n$ contains fewer than $e$ copies of $p$, so $(M/n)^2$ is divisible by $p^2$, and that term vanishes modulo $p^2$. If $n=P_pm\in G_p$, then $p\nmid m$ and the surviving term is $2(M_p/m)^2$. Because $2M_p^2$ is invertible modulo $p^2$, the global identity forces the local condition

$$\sum_{n\in S\cap G_p}\left(\frac{n}{P_p}\right)^{-2}\equiv 0 \pmod{p^2},$$

where $x^{-2}$ means the inverse of $x^2$ modulo $p^2$. This explains why the implementations inspect only the numbers in the top group for each odd prime: every other denominator is invisible mod $p^2$ at that prime.

The prime $2$ is deliberately excluded from this filter. After multiplying by $2M^2$, the factor 2 is not invertible modulo $4$, so the same clean argument does not apply. Powers of 2 are therefore left to the later global subset-sum stage.

What the Local Filters Do at $L=80$

For each odd prime, the implementations enumerate every subset of $G_p$ and keep only the masks satisfying the congruence above. These groups are tiny, so exhaustive local search is cheap. For $L=80$, the pruning is extremely strong.

Most odd-prime top groups allow only the empty choice. For example, $G_3=\{27,54\}$, $G_5=\{25,50,75\}$, $G_7=\{49\}$, and $G_{11}=\{11,22,33,44,55,66,77\}$ all collapse immediately: no non-empty subset satisfies the required congruence modulo $p^2$.

The only nontrivial group is

$$G_{13}=\{13,26,39,52,65,78\}.$$

Here the valid local choices are exactly $\varnothing$ and $\{13,39,52\}$. Indeed, after normalization by 13 the selected residues are $1,3,4$, and

$$1^{-2}+3^{-2}+4^{-2}\equiv 1+94+74=169\equiv 0 \pmod{169}.$$

So the entire 13-part of any global solution is a binary decision: either choose none of these numbers, or choose the three-number block $\{13,39,52\}$. Numbers such as 26, 65, and 78 never appear because they belong to no valid local mask.

After all odd-prime filters are applied, the candidate set shrinks from 79 numbers down to 36. Of those 36, three belong to the single nontrivial 13-group, and the other 33 are unconstrained by any odd-prime top group and become the free part of the search.

Turning the Rational Equation into Integer Weights

Once impossible denominators are removed, let $C$ be the least common multiple of all surviving candidates. For each remaining $n$, define

$$w_n=\frac{C^2}{n^2},\qquad T=\frac{C^2}{2}.$$

Then the original equation is equivalent to the exact integer subset-sum equation

$$\sum_{n\in S} w_n=T.$$

This conversion is crucial. It removes all rational arithmetic from the search itself: every valid local prime choice becomes an integer partial sum, and the final counting problem is purely additive.

Convolving the Constrained Prime Blocks

For each active odd-prime group $g$, let $A_g(t)$ denote the number of valid masks in that group whose total integer weight is $t$. Because the groups are disjoint for this $L=80$ instance, their contributions can be combined independently.

The implementations use the convolution recurrence

$$B_0(0)=1,\qquad B_{r+1}(u)=\sum_t B_r(u-t)\,A_{r+1}(t),$$

where $B_r(u)$ counts how many ways the first $r$ active groups can contribute total weight $u$. This is not yet the full answer; it is only the count of all legal choices forced by the odd-prime congruences.

For Problem 152, that recurrence is almost trivial after pruning, because only the 13-group survives as an active group. Its contribution map has just two entries: one for choosing nothing from $G_{13}$, and one for choosing the block $\{13,39,52\}$.

Meet-in-the-Middle on the Free Part

The remaining 33 denominators lie outside all odd-prime top groups. Let their integer weights be $f_1,\dots,f_{33}$. A direct search over these numbers would still cost $2^{33}$ subset checks, so the implementations split them into two halves, of sizes 16 and 17.

Define $L(x)$ as the number of left-half subsets with weight sum $x$, and $R(y)$ as the number of right-half subsets with weight sum $y$. Then the final number of solutions is

$$\sum_u B(u)\sum_x L(x)\,R(T-u-x).$$

In words: choose a legal contribution from the constrained odd-prime groups, choose a subset from the left free half, and ask whether the right free half can supply the exact missing weight. Hash tables make the final lookup fast.

Worked Examples

The $13$-group example above is the cleanest illustration of the local filter at the real limit $L=80$: it shows that the congruence step does not merely remove isolated denominators, but can force a whole block decision.

A smaller checkpoint instance also appears in the implementations. For $L=45$, one valid identity is

$$\frac12=\frac1{2^2}+\frac1{3^2}+\frac1{4^2}+\frac1{5^2}+\frac1{7^2}+\frac1{12^2}+\frac1{15^2}+\frac1{20^2}+\frac1{28^2}+\frac1{35^2}.$$

Inside the $7$-group for $L=45$, the chosen denominators are $\{7,28,35\}$. After dividing by 7 they become $\{1,4,5\}$, and

$$1^{-2}+4^{-2}+5^{-2}\equiv 1+46+2=49\equiv 0 \pmod{49}.$$

So the local congruence is genuinely visible in an explicit solution. The point of the full algorithm is that it enforces all such local conditions first and only then solves the remaining global sum exactly.

How the Code Works

Build the Odd-Prime Feasibility Data

The C++, Python, and Java implementations first sieve the primes up to the limit, determine the largest power $P_p\le 80$ for each odd prime, and form the corresponding top groups $G_p$. Every subset of each small group is tested modulo $p^2$ by using modular inverses of the normalized squares. From this, the implementation learns two things: which local masks are legal, and which denominators never occur in any legal mask and can therefore be discarded immediately.

Convert Legal Local Choices into Weighted Partial Sums

After pruning, the implementation computes a common denominator $C$, turns every surviving reciprocal square into the integer weight $C^2/n^2$, and sets the target to $C^2/2$. Each valid mask in each active odd-prime group is converted into a partial weight sum, and the constrained groups are combined by the convolution recurrence described above. The result is a map from constrained partial sum to the number of ways to obtain it.

Enumerate the Free Halves and Match Complements

The remaining free weights are split into two halves. Each implementation enumerates all subset sums of the left half and all subset sums of the right half, storing multiplicities rather than raw subset lists. The final answer is obtained by scanning the left sums and the constrained sums and looking up the unique complementary right sum needed to hit the target.

The C++ and Python implementations can parallelize this final matching pass when the table is large enough. The Java implementation follows the same arithmetic, but keeps the control flow single-threaded and compact. The C++ and Python versions also validate the method on smaller checkpoint cases, including three explicit identities for $L=45$, the fact that there are exactly three solutions up to 45, and a brute-force comparison at $L=18$.

Complexity Analysis

If the active odd-prime groups have sizes $g_1,g_2,\dots$, building the local legality tables costs

$$O\!\left(\sum_i 2^{g_i}g_i\right),$$

because each group is searched exhaustively. If $A_i$ is the number of distinct weighted options produced by group $i$, the constrained-group convolution costs roughly $O\!\left(\sum_i |B_{i-1}|\,|A_i|\right)$, where $|B_{i-1}|$ is the number of partial sums already accumulated.

The dominant phase is the meet-in-the-middle search on the free part. If $f$ free denominators remain, the subset-sum tables require about $2^{\lfloor f/2\rfloor}+2^{\lceil f/2\rceil}$ subset enumerations, plus hash lookups to match complementary sums. For the actual Euler 152 instance, the local filters reduce the problem to one nontrivial odd-prime block and $f=33$ free numbers, so the global search is closer to $2^{16}+2^{17}$ than to $2^{79}$.

Memory usage is dominated by the hash tables for left-half and right-half subset sums and by the map of constrained partial sums. That is exactly the right tradeoff here: a moderate amount of memory replaces an astronomically large brute-force search.

Footnotes and References

Problem page: https://projecteuler.net/problem=152
Subset-sum problem: Wikipedia - Subset sum problem
$p$-adic valuation: Wikipedia - p-adic valuation
Least common multiple: Wikipedia - Least common multiple
Modular multiplicative inverse: Wikipedia - Modular multiplicative inverse

Problem 152 source code

C++

#include <algorithm>
#include <atomic>
#include <cstdint>
#include <cstdlib>
#include <iostream>
#include <limits>
#include <numeric>
#include <stdexcept>
#include <string>
#include <thread>
#include <unordered_map>
#include <utility>
#include <vector>

namespace {

using i64 = std::int64_t;
using u64 = std::uint64_t;

struct Options {
    int limit = 80;
    bool run_checkpoints = true;
    unsigned requested_threads = 0U;
};

struct PrimeTopGroup {
    int prime = 0;
    int prime_power = 1;
    std::vector<int> numbers;
    std::vector<u64> valid_masks;
    u64 used_union_mask = 0ULL;
};

bool parse_unsigned_after_prefix(const std::string& arg,
                                 const std::string& prefix,
                                 unsigned& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }

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

    value = static_cast<unsigned>(parsed);
    return true;
}

bool parse_int_after_prefix(const std::string& arg,
                            const std::string& prefix,
                            int& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }

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

    value = static_cast<int>(parsed);
    return true;
}

bool parse_arguments(const int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        if (parse_unsigned_after_prefix(arg, "--threads=", options.requested_threads)) {
            continue;
        }
        if (parse_int_after_prefix(arg, "--limit=", options.limit)) {
            continue;
        }

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

    return true;
}

unsigned pick_thread_count(const unsigned requested) {
    if (requested > 0U) {
        return requested;
    }
    unsigned hw = std::thread::hardware_concurrency();
    if (hw == 0U) {
        hw = 4U;
    }
    return hw;
}

std::vector<int> sieve_primes(const int limit) {
    std::vector<bool> is_prime(static_cast<std::size_t>(limit) + 1ULL, true);
    if (limit >= 0) {
        is_prime[0] = false;
    }
    if (limit >= 1) {
        is_prime[1] = false;
    }

    for (int p = 2; p * p <= limit; ++p) {
        if (!is_prime[static_cast<std::size_t>(p)]) {
            continue;
        }
        for (int q = p * p; q <= limit; q += p) {
            is_prime[static_cast<std::size_t>(q)] = false;
        }
    }

    std::vector<int> primes;
    for (int p = 2; p <= limit; ++p) {
        if (is_prime[static_cast<std::size_t>(p)]) {
            primes.push_back(p);
        }
    }
    return primes;
}

i64 checked_lcm(const i64 a, const i64 b) {
    const i64 g = std::gcd(a, b);
    const __int128 scaled = static_cast<__int128>(a / g) * static_cast<__int128>(b);
    if (scaled > static_cast<__int128>(std::numeric_limits<i64>::max())) {
        throw std::overflow_error("LCM overflow");
    }
    return static_cast<i64>(scaled);
}

i64 mod_inverse(i64 a, i64 mod) {
    i64 t = 0;
    i64 new_t = 1;
    i64 r = mod;
    i64 new_r = a % mod;

    while (new_r != 0) {
        const i64 q = r / new_r;

        const i64 next_t = t - q * new_t;
        t = new_t;
        new_t = next_t;

        const i64 next_r = r - q * new_r;
        r = new_r;
        new_r = next_r;
    }

    if (r != 1) {
        throw std::runtime_error("Modular inverse does not exist");
    }

    t %= mod;
    if (t < 0) {
        t += mod;
    }
    return t;
}

std::vector<PrimeTopGroup> build_odd_prime_top_groups(const int limit) {
    const std::vector<int> primes = sieve_primes(limit);
    std::vector<PrimeTopGroup> groups;

    for (const int p : primes) {
        if (p == 2) {
            continue;
        }

        int prime_power = p;
        while (prime_power <= limit / p) {
            prime_power *= p;
        }

        PrimeTopGroup group;
        group.prime = p;
        group.prime_power = prime_power;

        const int next_power = prime_power * p;
        for (int n = 2; n <= limit; ++n) {
            if (n % prime_power != 0) {
                continue;
            }
            if (next_power <= limit && n % next_power == 0) {
                continue;
            }
            group.numbers.push_back(n);
        }

        if (group.numbers.empty()) {
            continue;
        }

        if (group.numbers.size() > 20U) {
            throw std::runtime_error("Unexpectedly large odd-prime top group");
        }

        // For odd prime p, modulo p^2 only numbers with maximal p-adic valuation survive.
        // Their normalized coefficients are inverses of (n / p^e)^2 modulo p^2.
        const i64 mod = static_cast<i64>(p) * static_cast<i64>(p);
        std::vector<i64> coeffs;
        coeffs.reserve(group.numbers.size());
        for (const int n : group.numbers) {
            const i64 m = static_cast<i64>(n / prime_power);
            const i64 sq = (m * m) % mod;
            coeffs.push_back(mod_inverse(sq, mod));
        }

        const u64 mask_count = 1ULL << group.numbers.size();
        group.valid_masks.reserve(static_cast<std::size_t>(mask_count));

        for (u64 mask = 0ULL; mask < mask_count; ++mask) {
            i64 residue = 0;
            for (std::size_t i = 0; i < coeffs.size(); ++i) {
                if ((mask >> i) & 1ULL) {
                    residue += coeffs[i];
                    residue %= mod;
                }
            }
            if (residue == 0) {
                group.valid_masks.push_back(mask);
                group.used_union_mask |= mask;
            }
        }

        groups.push_back(std::move(group));
    }

    return groups;
}

std::unordered_map<i64, u64> subset_sum_counts(const std::vector<i64>& weights) {
    if (weights.size() > 30U) {
        throw std::runtime_error("Subset half is unexpectedly large");
    }

    const u64 mask_count = 1ULL << weights.size();
    std::unordered_map<i64, u64> counts;
    counts.reserve(static_cast<std::size_t>(mask_count * 1.3));

    for (u64 mask = 0ULL; mask < mask_count; ++mask) {
        i64 sum = 0;
        for (std::size_t i = 0; i < weights.size(); ++i) {
            if ((mask >> i) & 1ULL) {
                sum += weights[i];
            }
        }
        ++counts[sum];
    }

    return counts;
}

u64 count_with_meet_in_middle(const std::vector<i64>& free_weights,
                              const std::vector<std::pair<i64, u64>>& base_states,
                              const i64 target,
                              const unsigned thread_count) {
    const std::size_t mid = free_weights.size() / 2U;
    const std::vector<i64> left_weights(free_weights.begin(), free_weights.begin() + static_cast<std::ptrdiff_t>(mid));
    const std::vector<i64> right_weights(free_weights.begin() + static_cast<std::ptrdiff_t>(mid), free_weights.end());

    const std::unordered_map<i64, u64> left_counts = subset_sum_counts(left_weights);
    const std::unordered_map<i64, u64> right_counts = subset_sum_counts(right_weights);

    std::vector<std::pair<i64, u64>> left_entries;
    left_entries.reserve(left_counts.size());
    for (const auto& [sum, ways] : left_counts) {
        left_entries.push_back({sum, ways});
    }

    auto count_chunk = [&](const std::size_t begin, const std::size_t end) {
        u64 local = 0ULL;
        for (std::size_t i = begin; i < end; ++i) {
            const i64 left_sum = left_entries[i].first;
            const u64 left_ways = left_entries[i].second;

            for (const auto& [base_sum, base_ways] : base_states) {
                const i64 need = target - base_sum - left_sum;
                const auto it = right_counts.find(need);
                if (it == right_counts.end()) {
                    continue;
                }
                local += left_ways * it->second * base_ways;
            }
        }
        return local;
    };

    if (thread_count <= 1U || left_entries.size() < 12000U) {
        return count_chunk(0U, left_entries.size());
    }

    const unsigned workers = std::min<unsigned>(thread_count, static_cast<unsigned>(left_entries.size()));
    std::atomic<std::size_t> next(0U);
    constexpr std::size_t kChunk = 256U;

    std::vector<u64> partial(workers, 0ULL);
    std::vector<std::thread> pool;
    pool.reserve(workers);

    for (unsigned t = 0U; t < workers; ++t) {
        pool.emplace_back([&, t]() {
            u64 local = 0ULL;
            while (true) {
                const std::size_t begin = next.fetch_add(kChunk, std::memory_order_relaxed);
                if (begin >= left_entries.size()) {
                    break;
                }
                const std::size_t end = std::min(begin + kChunk, left_entries.size());
                local += count_chunk(begin, end);
            }
            partial[t] = local;
        });
    }

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

    u64 total = 0ULL;
    for (const u64 v : partial) {
        total += v;
    }
    return total;
}

u64 count_representations(const int limit, const unsigned thread_count) {
    if (limit < 2) {
        return 0ULL;
    }
    if (limit > 80) {
        throw std::runtime_error("This implementation targets limits up to 80.");
    }

    const std::vector<PrimeTopGroup> groups = build_odd_prime_top_groups(limit);

    std::vector<int> number_to_group(static_cast<std::size_t>(limit) + 1ULL, -1);
    std::vector<int> number_to_group_bit(static_cast<std::size_t>(limit) + 1ULL, -1);

    for (std::size_t g = 0; g < groups.size(); ++g) {
        for (std::size_t i = 0; i < groups[g].numbers.size(); ++i) {
            const int n = groups[g].numbers[i];
            if (number_to_group[static_cast<std::size_t>(n)] != -1) {
                throw std::runtime_error("Odd-prime top groups overlapped unexpectedly.");
            }
            number_to_group[static_cast<std::size_t>(n)] = static_cast<int>(g);
            number_to_group_bit[static_cast<std::size_t>(n)] = static_cast<int>(i);
        }
    }

    std::vector<const PrimeTopGroup*> active_groups;
    active_groups.reserve(groups.size());

    for (const PrimeTopGroup& group : groups) {
        if (group.valid_masks.size() == 1U && group.valid_masks.front() == 0ULL) {
            continue;
        }
        active_groups.push_back(&group);
    }

    std::vector<char> number_is_possible(static_cast<std::size_t>(limit) + 1ULL, false);
    std::vector<int> free_numbers;
    free_numbers.reserve(static_cast<std::size_t>(limit));

    for (int n = 2; n <= limit; ++n) {
        const int gid = number_to_group[static_cast<std::size_t>(n)];
        if (gid == -1) {
            number_is_possible[static_cast<std::size_t>(n)] = true;
            free_numbers.push_back(n);
            continue;
        }

        const PrimeTopGroup& group = groups[static_cast<std::size_t>(gid)];
        const int bit = number_to_group_bit[static_cast<std::size_t>(n)];
        if (((group.used_union_mask >> bit) & 1ULL) != 0ULL) {
            number_is_possible[static_cast<std::size_t>(n)] = true;
        }
    }

    i64 common = 1;
    for (int n = 2; n <= limit; ++n) {
        if (!number_is_possible[static_cast<std::size_t>(n)]) {
            continue;
        }
        common = checked_lcm(common, static_cast<i64>(n));
    }

    const __int128 sq = static_cast<__int128>(common) * static_cast<__int128>(common);
    if (sq > static_cast<__int128>(std::numeric_limits<i64>::max())) {
        throw std::overflow_error("Common denominator square overflow");
    }
    const i64 common_sq = static_cast<i64>(sq);
    if ((common_sq % 2LL) != 0LL) {
        throw std::runtime_error("Unexpected odd denominator square");
    }
    const i64 target = common_sq / 2LL;

    std::vector<i64> weight_by_number(static_cast<std::size_t>(limit) + 1ULL, 0LL);
    for (int n = 2; n <= limit; ++n) {
        if (!number_is_possible[static_cast<std::size_t>(n)]) {
            continue;
        }
        const i64 nn = static_cast<i64>(n) * static_cast<i64>(n);
        if ((common_sq % nn) != 0LL) {
            throw std::runtime_error("Weight is not integral");
        }
        weight_by_number[static_cast<std::size_t>(n)] = common_sq / nn;
    }

    std::unordered_map<i64, u64> base_states;
    base_states.reserve(8U);
    base_states[0LL] = 1ULL;

    for (const PrimeTopGroup* group_ptr : active_groups) {
        const PrimeTopGroup& group = *group_ptr;

        std::unordered_map<i64, u64> option_counts;
        option_counts.reserve(group.valid_masks.size() * 2ULL + 1ULL);

        for (const u64 mask : group.valid_masks) {
            i64 sum = 0LL;
            for (std::size_t i = 0; i < group.numbers.size(); ++i) {
                if (((mask >> i) & 1ULL) == 0ULL) {
                    continue;
                }
                const int n = group.numbers[i];
                sum += weight_by_number[static_cast<std::size_t>(n)];
            }
            ++option_counts[sum];
        }

        std::unordered_map<i64, u64> next_states;
        next_states.reserve(base_states.size() * option_counts.size() * 2ULL + 1ULL);

        for (const auto& [base_sum, base_ways] : base_states) {
            for (const auto& [opt_sum, opt_ways] : option_counts) {
                next_states[base_sum + opt_sum] += base_ways * opt_ways;
            }
        }

        base_states.swap(next_states);
    }

    std::vector<i64> free_weights;
    free_weights.reserve(free_numbers.size());
    for (const int n : free_numbers) {
        free_weights.push_back(weight_by_number[static_cast<std::size_t>(n)]);
    }

    std::vector<std::pair<i64, u64>> base_entries;
    base_entries.reserve(base_states.size());
    for (const auto& [sum, ways] : base_states) {
        base_entries.push_back({sum, ways});
    }

    return count_with_meet_in_middle(free_weights, base_entries, target, thread_count);
}

u64 brute_force_count(const int limit) {
    if (limit < 2) {
        return 0ULL;
    }

    i64 common = 1;
    for (int n = 2; n <= limit; ++n) {
        common = checked_lcm(common, static_cast<i64>(n));
    }

    const __int128 sq = static_cast<__int128>(common) * static_cast<__int128>(common);
    if (sq > static_cast<__int128>(std::numeric_limits<i64>::max())) {
        throw std::overflow_error("Brute-force denominator square overflow");
    }
    const i64 common_sq = static_cast<i64>(sq);
    if ((common_sq % 2LL) != 0LL) {
        throw std::runtime_error("Unexpected odd brute-force denominator square");
    }
    const i64 target = common_sq / 2LL;

    std::vector<i64> weights;
    weights.reserve(static_cast<std::size_t>(limit - 1));
    for (int n = 2; n <= limit; ++n) {
        weights.push_back(common_sq / (static_cast<i64>(n) * static_cast<i64>(n)));
    }

    if (weights.size() > 26U) {
        throw std::runtime_error("Brute-force checkpoint requested above safe size.");
    }

    const u64 mask_count = 1ULL << weights.size();
    u64 count = 0ULL;
    for (u64 mask = 0ULL; mask < mask_count; ++mask) {
        i64 sum = 0LL;
        for (std::size_t i = 0; i < weights.size(); ++i) {
            if ((mask >> i) & 1ULL) {
                sum += weights[i];
            }
        }
        if (sum == target) {
            ++count;
        }
    }

    return count;
}

bool subset_sums_to_half(const std::vector<int>& subset) {
    i64 common = 2;
    for (const int n : subset) {
        common = checked_lcm(common, static_cast<i64>(n));
    }

    const __int128 sq = static_cast<__int128>(common) * static_cast<__int128>(common);
    if (sq > static_cast<__int128>(std::numeric_limits<i64>::max())) {
        throw std::overflow_error("Checkpoint denominator square overflow");
    }
    const i64 common_sq = static_cast<i64>(sq);

    i64 lhs = 0LL;
    for (const int n : subset) {
        lhs += common_sq / (static_cast<i64>(n) * static_cast<i64>(n));
    }

    return lhs * 2LL == common_sq;
}

bool run_checkpoints(const unsigned thread_count) {
    const std::vector<int> known_one = {2, 3, 4, 5, 7, 12, 15, 20, 28, 35};
    const std::vector<int> known_two = {2, 3, 4, 6, 7, 9, 10, 20, 28, 35, 36, 45};
    const std::vector<int> known_three = {2, 3, 4, 6, 7, 9, 12, 15, 28, 30, 35, 36, 45};

    if (!subset_sums_to_half(known_one)) {
        std::cerr << "Checkpoint failed: first known decomposition is invalid.\n";
        return false;
    }
    if (!subset_sums_to_half(known_two)) {
        std::cerr << "Checkpoint failed: second known decomposition is invalid.\n";
        return false;
    }
    if (!subset_sums_to_half(known_three)) {
        std::cerr << "Checkpoint failed: third known decomposition is invalid.\n";
        return false;
    }

    const u64 count_45 = count_representations(45, thread_count);
    if (count_45 != 3ULL) {
        std::cerr << "Checkpoint failed: expected 3 solutions for 2 <= n <= 45, got "
                  << count_45 << ".\n";
        return false;
    }

    const int brute_limit = 18;
    const u64 brute = brute_force_count(brute_limit);
    const u64 fast = count_representations(brute_limit, 1U);
    if (brute != fast) {
        std::cerr << "Checkpoint failed: fast/brute mismatch at limit " << brute_limit
                  << " (fast=" << fast << ", brute=" << brute << ").\n";
        return false;
    }

    return true;
}

}  // namespace

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

    if (options.limit < 2 || options.limit > 80) {
        std::cerr << "Please use --limit in [2, 80].\n";
        return 1;
    }

    const unsigned thread_count = pick_thread_count(options.requested_threads);

    try {
        if (options.run_checkpoints && !run_checkpoints(thread_count)) {
            return 1;
        }

        const u64 answer = count_representations(options.limit, thread_count);
        std::cout << answer << '\n';
    } catch (const std::exception& ex) {
        std::cerr << "Error: " << ex.what() << '\n';
        return 1;
    }

    return 0;
}

Python

import math
from collections import defaultdict
from concurrent.futures import ThreadPoolExecutor
import sys

def sieve_primes(limit):
    if limit < 2:
        return []
    is_prime = [True] * (limit + 1)
    is_prime[0] = is_prime[1] = False
    for p in range(2, int(limit ** 0.5) + 1):
        if is_prime[p]:
            for q in range(p * p, limit + 1, p):
                is_prime[q] = False
    return [p for p in range(2, limit + 1) if is_prime[p]]

def mod_inverse(a, mod):
    t, new_t = 0, 1
    r, new_r = mod, a % mod
    while new_r != 0:
        q = r // new_r
        t, new_t = new_t, t - q * new_t
        r, new_r = new_r, r - q * new_r
    if r != 1:
        raise ValueError("Modular inverse does not exist")
    t %= mod
    if t < 0:
        t += mod
    return t

def build_odd_prime_top_groups(limit):
    primes = sieve_primes(limit)
    groups = []

    for p in primes:
        if p == 2:
            continue

        prime_power = p
        while prime_power <= limit // p:
            prime_power *= p

        numbers = []
        next_power = prime_power * p
        for n in range(2, limit + 1):
            if n % prime_power != 0:
                continue
            if next_power <= limit and n % next_power == 0:
                continue
            numbers.append(n)

        if not numbers:
            continue

        mod = p * p
        coeffs = []
        for n in numbers:
            m = n // prime_power
            sq = (m * m) % mod
            coeffs.append(mod_inverse(sq, mod))

        valid_masks = []
        used_union_mask = 0
        mask_count = 1 << len(numbers)

        for mask in range(mask_count):
            residue = 0
            for i, coeff in enumerate(coeffs):
                if (mask >> i) & 1:
                    residue += coeff
                    residue %= mod
            if residue == 0:
                valid_masks.append(mask)
                used_union_mask |= mask

        groups.append({
            'prime': p,
            'prime_power': prime_power,
            'numbers': numbers,
            'valid_masks': valid_masks,
            'used_union_mask': used_union_mask
        })

    return groups

def subset_sum_counts(weights):
    mask_count = 1 << len(weights)
    counts = defaultdict(int)

    for mask in range(mask_count):
        total = 0
        for i, w in enumerate(weights):
            if (mask >> i) & 1:
                total += w
        counts[total] += 1

    return dict(counts)

def count_with_meet_in_middle(free_weights, base_states, target, thread_count):
    mid = len(free_weights) // 2
    left_weights = free_weights[:mid]
    right_weights = free_weights[mid:]

    left_counts = subset_sum_counts(left_weights)
    right_counts = subset_sum_counts(right_weights)

    left_entries = list(left_counts.items())

    def count_chunk(begin, end):
        local_count = 0
        for i in range(begin, end):
            left_sum, left_ways = left_entries[i]
            for base_sum, base_ways in base_states:
                need = target - base_sum - left_sum
                if need in right_counts:
                    local_count += left_ways * right_counts[need] * base_ways
        return local_count

    if thread_count <= 1 or len(left_entries) < 12000:
        return count_chunk(0, len(left_entries))

    workers = min(thread_count, len(left_entries))
    chunk_size = 256
    next_idx = [0]
    lock = __import__('threading').Lock()

    def worker(t):
        local_count = 0
        while True:
            with lock:
                begin = next_idx[0]
                next_idx[0] += chunk_size
            if begin >= len(left_entries):
                break
            end = min(begin + chunk_size, len(left_entries))
            local_count += count_chunk(begin, end)
        return local_count

    with ThreadPoolExecutor(max_workers=workers) as executor:
        results = list(executor.map(worker, range(workers)))

    return sum(results)

def count_representations(limit, thread_count):
    if limit < 2:
        return 0

    groups = build_odd_prime_top_groups(limit)

    number_to_group = [-1] * (limit + 1)
    number_to_group_bit = [-1] * (limit + 1)

    for g_idx, group in enumerate(groups):
        for i, n in enumerate(group['numbers']):
            number_to_group[n] = g_idx
            number_to_group_bit[n] = i

    active_groups = []
    for group in groups:
        if len(group['valid_masks']) == 1 and group['valid_masks'][0] == 0:
            continue
        active_groups.append(group)

    number_is_possible = [False] * (limit + 1)
    free_numbers = []

    for n in range(2, limit + 1):
        gid = number_to_group[n]
        if gid == -1:
            number_is_possible[n] = True
            free_numbers.append(n)
            continue

        group = groups[gid]
        bit = number_to_group_bit[n]
        if (group['used_union_mask'] >> bit) & 1:
            number_is_possible[n] = True

    common = 1
    for n in range(2, limit + 1):
        if not number_is_possible[n]:
            continue
        common = common * n // math.gcd(common, n)

    common_sq = common * common
    if common_sq % 2 != 0:
        raise ValueError("Unexpected odd denominator square")

    target = common_sq // 2

    weight_by_number = [0] * (limit + 1)
    for n in range(2, limit + 1):
        if not number_is_possible[n]:
            continue
        nn = n * n
        if common_sq % nn != 0:
            raise ValueError("Weight is not integral")
        weight_by_number[n] = common_sq // nn

    base_states = {0: 1}

    for group in active_groups:
        option_counts = defaultdict(int)

        for mask in group['valid_masks']:
            total = 0
            for i, n in enumerate(group['numbers']):
                if (mask >> i) & 1:
                    total += weight_by_number[n]
            option_counts[total] += 1

        next_states = defaultdict(int)
        for base_sum, base_ways in base_states.items():
            for opt_sum, opt_ways in option_counts.items():
                next_states[base_sum + opt_sum] += base_ways * opt_ways

        base_states = next_states

    free_weights = [weight_by_number[n] for n in free_numbers]

    base_entries = list(base_states.items())

    return count_with_meet_in_middle(free_weights, base_entries, target, thread_count)

def brute_force_count(limit):
    if limit < 2:
        return 0

    common = 1
    for n in range(2, limit + 1):
        common = common * n // math.gcd(common, n)

    common_sq = common * common
    if common_sq % 2 != 0:
        raise ValueError("Unexpected odd brute-force denominator square")

    target = common_sq // 2

    weights = []
    for n in range(2, limit + 1):
        weights.append(common_sq // (n * n))

    if len(weights) > 26:
        raise ValueError("Brute-force checkpoint requested above safe size.")

    count = 0
    mask_count = 1 << len(weights)
    for mask in range(mask_count):
        total = 0
        for i, w in enumerate(weights):
            if (mask >> i) & 1:
                total += w
        if total == target:
            count += 1

    return count

def subset_sums_to_half(subset):
    common = 2
    for n in subset:
        common = common * n // math.gcd(common, n)

    common_sq = common * common

    lhs = 0
    for n in subset:
        lhs += common_sq // (n * n)

    return lhs * 2 == common_sq

def run_checkpoints(thread_count):
    known_one = [2, 3, 4, 5, 7, 12, 15, 20, 28, 35]
    known_two = [2, 3, 4, 6, 7, 9, 10, 20, 28, 35, 36, 45]
    known_three = [2, 3, 4, 6, 7, 9, 12, 15, 28, 30, 35, 36, 45]

    if not subset_sums_to_half(known_one):
        sys.stderr.write("Checkpoint failed: first known decomposition is invalid.\n")
        return False
    if not subset_sums_to_half(known_two):
        sys.stderr.write("Checkpoint failed: second known decomposition is invalid.\n")
        return False
    if not subset_sums_to_half(known_three):
        sys.stderr.write("Checkpoint failed: third known decomposition is invalid.\n")
        return False

    count_45 = count_representations(45, thread_count)
    if count_45 != 3:
        sys.stderr.write(f"Checkpoint failed: expected 3 solutions for 2 <= n <= 45, got {count_45}.\n")
        return False

    brute_limit = 18
    brute = brute_force_count(brute_limit)
    fast = count_representations(brute_limit, 1)
    if brute != fast:
        sys.stderr.write(f"Checkpoint failed: fast/brute mismatch at limit {brute_limit} (fast={fast}, brute={brute}).\n")
        return False

    return True

def main():
    args = sys.argv[1:]

    limit = 80
    run_checkpoints_flag = True
    requested_threads = 0

    for arg in args:
        if arg == "--skip-checkpoints":
            run_checkpoints_flag = False
        elif arg.startswith("--threads="):
            try:
                requested_threads = int(arg[10:])
            except ValueError:
                sys.stderr.write(f"Invalid thread count: {arg}\n")
                return 1
        elif arg.startswith("--limit="):
            try:
                limit = int(arg[8:])
            except ValueError:
                sys.stderr.write(f"Invalid limit: {arg}\n")
                return 1
        else:
            sys.stderr.write(f"Unknown argument: {arg}\n")
            return 1

    if limit < 2 or limit > 80:
        sys.stderr.write("Please use --limit in [2, 80].\n")
        return 1

    if requested_threads == 0:
        import os
        hw = os.cpu_count()
        requested_threads = hw if hw else 4

    try:
        if run_checkpoints_flag and not run_checkpoints(requested_threads):
            return 1

        answer = count_representations(limit, requested_threads)
        print(answer)
    except Exception as ex:
        sys.stderr.write(f"Error: {ex}\n")
        return 1

if __name__ == "__main__":
    main()

Java

import java.util.*;

public class Euler152 {
    static int[] sievePrimes(int limit) {
        boolean[] np = new boolean[limit + 1];
        np[0] = np[1] = true;
        for (int p = 2; p * p <= limit; p++)
            if (!np[p])
                for (int q = p * p; q <= limit; q += p)
                    np[q] = true;
        int c = 0;
        for (int i = 2; i <= limit; i++)
            if (!np[i])
                c++;
        int[] pr = new int[c];
        c = 0;
        for (int i = 2; i <= limit; i++)
            if (!np[i])
                pr[c++] = i;
        return pr;
    }

    static long modInverse(long a, long mod) {
        long t = 0, nt = 1, r = mod, nr = a % mod;
        while (nr != 0) {
            long q = r / nr;
            long tmp = nt;
            nt = t - q * nt;
            t = tmp;
            tmp = nr;
            nr = r - q * nr;
            r = tmp;
        }
        return ((t % mod) + mod) % mod;
    }

    public static void main(String[] args) {
        int limit = 80;
        int[] primes = sievePrimes(limit);
        // Build odd prime top groups
        List<int[]> groupNumbers = new ArrayList<>();
        List<int[][]> groupValid = new ArrayList<>();
        List<Integer> groupUsedUnion = new ArrayList<>();

        for (int p : primes) {
            if (p == 2)
                continue;
            int pp = p;
            while (pp <= limit / p)
                pp *= p;
            int npp = pp * p;
            List<Integer> nums = new ArrayList<>();
            for (int n = 2; n <= limit; n++) {
                if (n % pp != 0)
                    continue;
                if (npp <= limit && n % npp == 0)
                    continue;
                nums.add(n);
            }
            if (nums.isEmpty())
                continue;
            long mod = (long) p * p;
            long[] coeffs = new long[nums.size()];
            for (int i = 0; i < nums.size(); i++) {
                long m = nums.get(i) / pp;
                long sq = (m * m) % mod;
                coeffs[i] = modInverse(sq, mod);
            }
            List<Integer> validMasks = new ArrayList<>();
            int unionMask = 0;
            for (int mask = 0; mask < (1 << nums.size()); mask++) {
                long res = 0;
                for (int i = 0; i < nums.size(); i++)
                    if ((mask >> i & 1) == 1)
                        res = (res + coeffs[i]) % mod;
                if (res == 0) {
                    validMasks.add(mask);
                    unionMask |= mask;
                }
            }
            int[] na = nums.stream().mapToInt(Integer::intValue).toArray();
            int[][] va = { validMasks.stream().mapToInt(Integer::intValue).toArray() };
            groupNumbers.add(na);
            groupValid.add(va);
            groupUsedUnion.add(unionMask);
        }

        int[] numToGroup = new int[limit + 1], numToBit = new int[limit + 1];
        Arrays.fill(numToGroup, -1);
        Arrays.fill(numToBit, -1);
        for (int g = 0; g < groupNumbers.size(); g++) {
            int[] gn = groupNumbers.get(g);
            for (int i = 0; i < gn.length; i++) {
                numToGroup[gn[i]] = g;
                numToBit[gn[i]] = i;
            }
        }

        boolean[] possible = new boolean[limit + 1];
        List<Integer> freeNums = new ArrayList<>();
        for (int n = 2; n <= limit; n++) {
            int g = numToGroup[n];
            if (g == -1) {
                possible[n] = true;
                freeNums.add(n);
            } else if ((groupUsedUnion.get(g) >> numToBit[n] & 1) == 1)
                possible[n] = true;
        }

        long common = 1;
        for (int n = 2; n <= limit; n++)
            if (possible[n]) {
                common = common / gcd(common, n) * n;
            }
        long commonSq = common * common;
        long target = commonSq / 2;

        long[] weight = new long[limit + 1];
        for (int n = 2; n <= limit; n++)
            if (possible[n])
                weight[n] = commonSq / ((long) n * n);

        // Process groups
        Map<Long, Long> baseStates = new HashMap<>();
        baseStates.put(0L, 1L);
        for (int g = 0; g < groupNumbers.size(); g++) {
            int[] gn = groupNumbers.get(g);
            int[] vm = groupValid.get(g)[0];
            if (vm.length == 1 && vm[0] == 0)
                continue;
            Map<Long, Long> optCounts = new HashMap<>();
            for (int mask : vm) {
                long tot = 0;
                for (int i = 0; i < gn.length; i++)
                    if ((mask >> i & 1) == 1)
                        tot += weight[gn[i]];
                optCounts.merge(tot, 1L, Long::sum);
            }
            Map<Long, Long> next = new HashMap<>();
            for (var e1 : baseStates.entrySet())
                for (var e2 : optCounts.entrySet())
                    next.merge(e1.getKey() + e2.getKey(), e1.getValue() * e2.getValue(), Long::sum);
            baseStates = next;
        }

        // Meet in middle on free numbers
        long[] fw = freeNums.stream().mapToLong(n -> weight[n]).toArray();
        int mid = fw.length / 2;
        Map<Long, Long> leftCounts = subsetSums(fw, 0, mid);
        Map<Long, Long> rightCounts = subsetSums(fw, mid, fw.length);

        long count = 0;
        for (var be : baseStates.entrySet())
            for (var le : leftCounts.entrySet()) {
                long need = target - be.getKey() - le.getKey();
                Long rc = rightCounts.get(need);
                if (rc != null)
                    count += be.getValue() * le.getValue() * rc;
            }
        System.out.println(count);
    }

    static Map<Long, Long> subsetSums(long[] w, int from, int to) {
        Map<Long, Long> counts = new HashMap<>();
        for (int mask = 0; mask < (1 << (to - from)); mask++) {
            long s = 0;
            for (int i = 0; i < to - from; i++)
                if ((mask >> i & 1) == 1)
                    s += w[from + i];
            counts.merge(s, 1L, Long::sum);
        }
        return counts;
    }

    static long gcd(long a, long b) {
        while (b != 0) {
            long t = b;
            b = a % b;
            a = t;
        }
        return a;
    }
}