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

We must count the sets $S$ of distinct integers between $2$ and $80$ inclusive for which

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

The statement gives the example {2, 3, 4, 5, 7, 12, 15, 20, 28, 35} and says that the integers from $2$ to $45$ allow exactly three such sets, the other two being {2, 3, 4, 6, 7, 9, 10, 20, 28, 35, 36, 45} and {2, 3, 4, 6, 7, 9, 12, 15, 28, 30, 35, 36, 45}. ^[1] Sums of distinct unit fractions are called Egyptian fractions, and by a theorem of R. L. Graham the finite sums of reciprocals of distinct squares are exactly the rational numbers in $[0,\pi^2/6-1)$ and $[1,\pi^2/6)$. ^[2]

The whole solution in one sentence: a prime factor that occurs in only a few numbers up to $80$ must cancel out of the denominators, a test modulo the square of the prime shows that this leaves only $33$ candidates and the optional block $13,39,52$, and the subsets with the right sum are counted by splitting the $33$ numbers into two halves and matching the sums of the halves.

Mathematical Approach

Trying all $2^{79}\approx6\cdot10^{23}$ subsets is out of the question. We first find out, prime by prime, which numbers can occur at all, and then count the subsets of the few remaining numbers exactly.

1. A prime that occurs only once

Suppose that $79\in S$. Then

$$\frac1{79^2}=\frac12-\sum_{n\in S,\ n\ne79}\frac1{n^2}.$$

No other number up to $80$ is a multiple of $79$, so no fraction on the right has the factor $79$ in its denominator, and neither does their sum, while the left side has the denominator $79^2$. This is impossible, so $79\notin S$. In the same way no prime from $41$ to $79$ occurs, since each of them is the only multiple of itself up to $80$. In terms of the $p$-adic valuation $v_p$, the exponent of $p$ in the numerator of a nonzero fraction minus that in its denominator, the argument uses $v_p(x+y)=\min(v_p(x),v_p(y))$ for $v_p(x)\ne v_p(y)$. ^[3]

2. Several multiples of one prime

For a prime $p\ge11$ we have $p^2\gt80$, so the multiples of $p$ up to $80$ are the numbers $pm$ with $m\lt p$. Suppose $S$ contains $pm$ exactly for the $m$ in a set $A$. They contribute

$$\frac1{p^2}\sum_{m\in A}\frac1{m^2}=\frac1{p^2}\cdot\frac ab,$$

with $a/b$ in lowest terms and $p\nmid b$. All other members of $S$ are free of the factor $p$, so by Section 1 this contribution must have a denominator without $p$: the numerator $a$ must be divisible by $p^2$. For the primes $17$ to $37$ we have $m\le4$, and over the common denominator $144$, which $p$ does not divide, the numerator is at most $144+36+16+9=205\lt17^2$, so only $A=\varnothing$ is possible. For $p=11$ a check of all $127$ nonempty sets $A\subseteq\{1,\dots,7\}$ finds none that works. For $p=13$ exactly one of the $63$ nonempty sets $A\subseteq\{1,\dots,6\}$ works, namely $A=\{1,3,4\}$:

$$\begin{aligned}1+\frac1{3^2}+\frac1{4^2}&=\frac{144+16+9}{144}\\&=\frac{13^2}{12^2}.\end{aligned}$$

Divided by $13^2$ this reads $\frac1{13^2}+\frac1{39^2}+\frac1{52^2}=\frac1{12^2}$. So $13$, $39$ and $52$ can only occur together, as one block worth $1/144$, and $26$, $65$, $78$ never occur.

3. Testing modulo $p^2$

The programs test this condition with remainders. Two integers are congruent modulo $q$ when their difference is a multiple of $q$. If $p\nmid x$, there is an integer $x^{-1}$, the inverse of $x$ modulo $p^2$, with $x\,x^{-1}\equiv1\pmod{p^2}$; it exists because $\gcd(x,p^2)=1$, and the extended Euclidean algorithm finds it. ^[4] ^[5] Write $m^{-2}$ for the inverse of $m^2$. If $D$ is a common denominator of the fractions $1/m^2$ with $p\nmid D$, then $m^{-2}\equiv(D/m^2)\,D^{-1}$, so

$$\begin{aligned}\sum_{m\in A}m^{-2}&\equiv N\cdot D^{-1}\pmod{p^2},\\ N&=\sum_{m\in A}\frac{D}{m^2},\end{aligned}$$

and the left side is $0$ exactly when $p^2$ divides the numerator $N$, the condition of Section 2. For $p=13$: $9\cdot94=846=5\cdot169+1$ and $16\cdot74=1184=7\cdot169+1$, so the inverses of $1$, $9$, $16$ modulo $169$ are $1$, $94$, $74$, and indeed $1+94+74=169\equiv0$.

4. The test for every odd prime

For $p=3,5,7$ some numbers up to $80$ are divisible by $p^2$, and one more idea is needed: only the multiples of the highest power of $p$ matter. Let $L$ be the upper limit, $p$ an odd prime, $p^e$ the largest power of $p$ with $p^e\le L$, and let the top group $G_p$ be the set of multiples of $p^e$ up to $L$. They are the numbers $p^em$ with $m\lt p$, so $p\nmid m$. The key claim is that for every solution $S$, the multipliers $m$ of the members of $S\cap G_p$ satisfy

$$\sum_{p^em\in S}m^{-2}\equiv0\pmod{p^2}.$$

To prove it, let $M=\operatorname{lcm}(2,\dots,L)$ and write $M=p^eM'$ with $p\nmid M'$. Multiplying the equation by $2M^2$ gives the equation of integers $\sum_{n\in S}2(M/n)^2=M^2$, whose right side is divisible by $p^2$. If $n\notin G_p$, then $n$ has at most $e-1$ factors $p$, so $p$ divides $M/n$ and $p^2$ divides its square. If $n=p^em$, then $M/n=M'/m$, and $(M'/m)\,m=M'$ gives $M/n\equiv M'm^{-1}\pmod{p^2}$. Modulo $p^2$ only

$$2M'^2\sum_{p^em\in S}m^{-2}\equiv0\pmod{p^2}$$

remains, and since $p$ is odd and does not divide $M'$, the factor $2M'^2$ is invertible and can be cancelled. For $L=80$ and $p\ge11$, where $e=1$, this is the test of Section 3. Modulo $p^3$ the members with exactly $e-1$ factors $p$ would no longer drop out, which is why the test uses $p^2$. For $p=2$ the factor $2M'^2$ cannot be cancelled, and the programs test no top group of $2$.

5. The test for the limit 80

For $L=80$ the top groups of $3$, $5$ and $7$ are $\{27,54\}$, $\{25,50,75\}$ and $\{49\}$. Modulo $9$ the inverses of $1^2$ and $2^2$ are $1$ and $7$, with the subset sums $1$, $7$, $8$; modulo $25$ those of $1^2$, $2^2$, $3^2$ are $1$, $19$, $14$, whose subset sums leave the remainders $1,8,9,14,15,19,20$; modulo $49$ the only sum is $1$. None is $0$, so $25$, $27$, $49$, $50$, $54$, $75$ never occur, while smaller multiples such as $9$, $35$ or $63$ are not affected. With Sections 1 and 2, exactly $36$ numbers survive: the block $13,39,52$ and the $33$ free numbers

$$\begin{aligned}&2,3,4,5,6,7,8,9,10,12,14,15,16,\\&18,20,21,24,28,30,32,35,36,40,\\&42,45,48,56,60,63,64,70,72,80,\end{aligned}$$

which lie in no top group: they are the integers $n$ from $2$ to $80$ of the form $2^a3^b5^c7^d$ with $b\le2$ and $c,d\le1$.

6. Worked example: the limit 45

For $L=45$ the highest power of $7$ is $7$ itself, so the top group of $7$ is $\{7,14,\dots,42\}$ with $m=1,\dots,6$. The inverses of $m^2$ modulo $49$ are $1,37,11,46,2,15$, and exactly three nonempty subsets have a sum divisible by $49$: $1+37+11=49$, $1+46+2=49$ and $37+46+15=98$. They are $\{7,14,21\}$, $\{7,28,35\}$ and $\{14,28,42\}$, worth

$$\begin{aligned}\frac1{7^2}+\frac1{14^2}+\frac1{21^2}&=\frac1{36},\\ \frac1{7^2}+\frac1{28^2}+\frac1{35^2}&=\frac9{400},\\ \frac1{14^2}+\frac1{28^2}+\frac1{42^2}&=\frac1{144}.\end{aligned}$$

Every other top group for $L=45$ admits only the empty choice, and $19$ free numbers remain: $2,3,4,5,6,8,9,10,12$ and $15,16,18,20,24,30,32,36,40,45$.

7. Whole numbers instead of fractions

Let $C$ be the least common multiple of the surviving numbers. Then every $1/n^2$ becomes the whole number $w_n=C^2/n^2$, and $1/2$ becomes $T=C^2/2$. For $L=80$ we get $C=262080=2^6\cdot3^2\cdot5\cdot7\cdot13$, $C^2=68685926400$ and $T=34342963200$, and the block $13,39,52$ weighs $C^2/144=476985600$. We must count the subsets of the free numbers with weight $T$, and those with weight $T-476985600$.

8. Meet in the middle

Instead of all $2^{33}\approx8.6\cdot10^9$ subsets of the free numbers, the meet-in-the-middle method of Horowitz and Sahni ^[6] works with two halves: the programs take the $16$ smallest free numbers as the left half $X$ and the other $17$ as the right half $Y$. Let $\ell(x)$ and $r(y)$ be the numbers of subsets of $X$ and of $Y$ with weights $x$ and $y$, and $\beta(u)$ the number of allowed choices from the top groups with weight $u$; for $L=80$, $\beta(0)=\beta(C^2/144)=1$. The number of solutions is

$$\sum_{u}\beta(u)\sum_{x}\ell(x)\,r(T-u-x).$$

The top groups never overlap: $q=p^e$ satisfies $qp\gt L$ and $q\ge p$, so $q^2\gt L$, and a common multiple of the top powers of two primes is at least their product, which exceeds $L$. Hence every solution splits in exactly one way into its parts in the top groups, in $X$ and in $Y$; by Section 4 the first part is an allowed choice, and the three weights add up to $T$. Conversely, every such triple of parts is a solution, so the double sum counts every solution exactly once. The tables count subsets because different subsets can share a weight: $\{12\}$ and $\{15,20\}$ both lie in $X$, and $1/12^2=1/15^2+1/20^2$. For $L=80$ the $2^{16}=65536$ subsets of $X$ have $52224$ different weights and the $2^{17}=131072$ subsets of $Y$ have $78544$. Horowitz and Sahni match sorted lists; the programs use hash tables.

9. Worked example continued

For $L=45$ the halves are $X=\{2,3,4,5,6,8,9,10,12\}$ and the other ten free numbers, and the choices from the top group of $7$ are worth $0$, $1/36$, $9/400$ and $1/144$. For a check by hand we measure in units of $1/3600$. The left subset $\{2,3,4,5,12\}$ is worth $900+400+225+144+25=1694$, the block $\{7,28,35\}$ is worth $81$ and $1/2$ is $1800$, so the right half has to supply $1800-81-1694=25$, that is $1/144$. Exactly one right subset is worth that, $\{15,20\}$, since $16+9=25$, and this gives the first solution of the statement. Over all left weights only the block $\{7,28,35\}$ finds partners, three in all, which are exactly the three solutions of the statement.

10. What the tests leave to the search

The programs test no top group of $2$, so some useless numbers stay. If $2^A$ with $A\ge2$ is the highest power of $2$ dividing a member of $S$, multiply the equation by $4^AK^2$, where $K$ is odd and divisible by the odd parts of all members. The right side becomes $2^{2A-1}K^2$, and modulo $4$ only the odd squares from the members divisible by $2^A$ remain, each with the remainder $1$; so their number is a positive multiple of $4$. Up to $80$ the numbers with exactly $6$, $5$ or $4$ factors $2$ are $\{64\}$, $\{32\}$ and $\{16,48,80\}$, fewer than four each time, so $16,32,48,64,80$ never occur; the tables simply find no solution with them.

How the Code Works

The C++ program ^[7] starts in main: parse_arguments accepts --limit=L, default $80$, and --skip-checkpoints, and a limit outside $[2,80]$ is refused. run_checkpoints runs first unless skipped; then the result of count_representations is printed.

count_representations calls build_odd_prime_top_groups, which takes the primes from sieve_primes, a sieve of Eratosthenes, ^[8] and for each odd prime multiplies prime_power by $p$ while it stays at most limit / p, which gives the $p^e$ of Section 4. For every multiple of prime_power it computes the inverse of $m^2$ modulo $p^2$ with mod_inverse, the extended Euclidean algorithm of ^[5]. Then all $2^{|G_p|}$ bit masks are tested, bit $i$ standing for the $i$-th number of the group; the masks whose coefficients add up to $0$ modulo $p^2$ go into valid_masks and are OR-ed into used_union_mask.

number_to_group and number_to_group_bit record where each number lies; an overlap of groups, impossible by Section 8, or a group of more than $20$ numbers would stop the program. A group is active unless its only valid mask is the empty one, and a number is kept if it lies in no group, which makes it free, or in some valid mask. checked_lcm forms $C$ from the kept numbers with a 128-bit overflow check, and $C^2$, the target and the weights are checked for overflow and exact division. The loop over the active groups builds base_states, the table $\beta$ of Section 8, by combining every old state with every entry of option_counts, which counts the valid masks of a group by weight; for $L=80$ it is $\{0\mapsto1,\ 476985600\mapsto1\}$.

count_with_meet_in_middle splits the free weights at mid $=\lfloor f/2\rfloor$ for $f$ free numbers. subset_sum_counts builds for one half an unordered_map from weight to number of subsets, running over all masks and adding the weights of the set bits. Finally, for every left entry and every base state, need $=T-u-x$ is looked up in the right table, and the product of the three counts is added to an unsigned 64-bit total.

Python mirrors these functions with dictionaries and unbounded integers, without the overlap and group-size checks; its main returns $1$ after an error, but the script ignores this value, so the exit status is then $0$. Java fixes $L=80$ and has no options or checkpoints; apart from the helpers sievePrimes, modInverse, gcd and subsetSums, everything runs in main, with HashMap<Long, Long>, int masks and unchecked long arithmetic. This suffices, because all weight sums stay below $C^2\approx6.9\cdot10^{10}$, since $\sum_{n\ge2}1/n^2\lt1$, and all counts are at most $2^{34}$. Its last loop puts the base states outside and the left table inside, the same sum in another order.

Complexity and Verification

For $L=80$ the local tests examine the $278$ subsets of the $21$ top groups. The two tables cost $16\cdot2^{16}+17\cdot2^{17}=3276800$ bit tests and $196608$ hash updates, and the final loop checks $2\cdot52224=104448$ pairs of a left weight and a base state. The memory is two hash tables with $52224$ and $78544$ entries, plus a list copy of the left one in C++ and Python. For $f$ free numbers, few base states and hash operations of constant average cost this means time $O(f\,2^{f/2})$ and memory $O(2^{f/2})$; the factor $f$ comes from scanning each mask bit by bit, which the merging of ^[6] avoids. On a desktop computer each program, checkpoints included, takes well under a second.

Before the main count, run_checkpoints in C++ and Python checks with exact integers that the three sets of the statement sum to $1/2$ (subset_sums_to_half), that the count for $L=45$ is $3$, and that brute_force_count(18), a loop over all $2^{17}$ subsets of $\{2,\dots,18\}$, agrees with the fast method. No solution exists below $35$, so this comparison only shows that both find nothing; the real test is the checkpoint for $L=45$. The Java version has no checks. Separate computations for this article, not part of the published solutions, confirmed the following:

an exact meet in the middle over all numbers $2,\dots,L$, with 128-bit integers and no prime test, finds no solution for $L\le34$, only the first example for $35\le L\le44$ and three solutions for $L=45$ and $46$, agreeing with count_representations for every $L\le46$;
for every $L\le80$ and every odd prime, the modular test accepts exactly the subsets of the top group whose sum of reciprocal squares, computed with exact fractions, has at most $2e-2$ factors $p$ in its denominator;
for $L=80$, an exact meet in the middle over all $45$ numbers up to $80$ without prime factors other than $2,3,5,7,13$, with no test at all for these primes, gives the same count as the programs;
a listing of all solutions for $L=80$, each checked with exact fractions, shows that all of them contain $2$ and $3$, that they have $10$ to $24$ members, that the only one with $10$ members is the first example of the statement, and that $16,32,48,64,80$ never occur;
the C++ program, built with GCC and with Clang, and the Python program give the same counts for every $L$ from $2$ to $80$, and all three programs print the same number for $L=80$.

Footnotes and References

The references below supply the background facts; the local test, the $36$ surviving numbers, the worked examples and the remark on the prime $2$ are derived in this article.

Project Euler 152 — Sums of Square Reciprocals. The official statement writes 1/2 with the squares of 2, 3, 4, 5, 7, 12, 15, 20, 28 and 35, says that the integers from 2 to 45 give exactly three ways and lists the other two, and asks for the number of ways with distinct integers between 2 and 80 inclusive.
Wikipedia — Egyptian fraction. Graham (1964) characterised the sums of distinct unit fractions with nth-power denominators; for square denominators they are exactly the rational numbers in [0, π²/6 − 1) ∪ [1, π²/6).
Wikipedia — p-adic valuation. The p-adic valuation of a nonzero integer is the exponent of the highest power of p dividing it, extended to rational numbers by ν_p(r/s) = ν_p(r) − ν_p(s); ν_p(r + s) = min(ν_p(r), ν_p(s)) whenever ν_p(r) ≠ ν_p(s).
Wikipedia — Modular multiplicative inverse. An integer a has an inverse modulo m if and only if gcd(a, m) = 1, and the inverse is then unique modulo m; it can be computed with the extended Euclidean algorithm.
Wikipedia — Extended Euclidean algorithm. Its section on modular integers computes the inverse of a modulo n by the loop over t, newt, r, newr that mod_inverse follows, including the final correction of a negative t.
Wikipedia — Subset sum problem. The algorithm of Horowitz and Sahni (1974) splits the elements into two halves, builds the sorted lists of the subset sums of both halves by repeated merging and matches them, in time and space O(2^(n/2)).
C++ — Euler152.cpp. The linked immutable revision, with build_odd_prime_top_groups, count_with_meet_in_middle and run_checkpoints, is the implementation record for this article.
Wikipedia — Sieve of Eratosthenes. The sieve marks the multiples of each prime as composite, starting from its square, and takes time O(n log log n).

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