Problem 241: Perfection Quotients

Project Euler Problem 241 Solution

EulerSolve provides an optimized solution for Project Euler Problem 241, Perfection Quotients, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For a positive integer $n$, the perfection quotient is $\sigma(n)/n$, where $\sigma(n)$ is the sum of the positive divisors of $n$. The task is to sum every $n \le 10^{18}$ whose perfection quotient is a half-integer. For the bound used here, the implementations split the search into the six target values $$\frac{3}{2},\ \frac{5}{2},\ \frac{7}{2},\ \frac{9}{2},\ \frac{11}{2},\ \frac{13}{2}.$$ A direct scan up to $10^{18}$ is hopeless. The successful strategy is to build $n$ prime-power by prime-power, while tracking an exact residual quotient that measures how far the current partial factorization is from the chosen target. Mathematical Approach Fix one target value $T=\frac{A}{2}$. We want all integers $n$ satisfying $$\frac{\sigma(n)}{n}=T.$$ The quotient factors cleanly over prime powers If $$n=\prod_i p_i^{e_i},$$ then multiplicativity gives $$\sigma(n)=\prod_i \sigma(p_i^{e_i}),\qquad \frac{\sigma(n)}{n}=\prod_i \frac{\sigma(p_i^{e_i})}{p_i^{e_i}}.$$ For a single prime power, $$\sigma(p^e)=1+p+p^2+\cdots+p^e=\frac{p^{e+1}-1}{p-1},$$ so its contribution to the perfection quotient is $$\frac{\sigma(p^e)}{p^e}=1+\frac1p+\frac1{p^2}+\cdots+\frac1{p^e}.$$ This is always greater than 1. Therefore every time we append a new prime power to $n$, the quotient $\sigma(n)/n$ increases, and the reciprocal quantity $n/\sigma(n)$ decreases....

Detailed mathematical approach

Problem Summary

For a positive integer $n$, let $\sigma(n)$ be the sum of all positive divisors of $n$; for example $\sigma(6)=1+2+3+6=12$. A perfect number satisfies $\sigma(n)=2n$, and the perfection quotient of $n$ is $p(n)=\sigma(n)/n$. We must find the sum of all positive integers $n\le10^{18}$ for which $p(n)$ has the form $k+\tfrac12$ with an integer $k$. ^[1] The statement gives no checkpoint for the requested sum. From now on we write $\sigma(n)/n$ instead of $p(n)$ and keep the letter $p$ for primes.

The whole solution in one sentence: we build $n$ one prime power at a time and keep track of the fraction that the rest of $n$ still has to supply, whose denominator must divide that rest and therefore names the next prime, and we abandon a branch as soon as the bound $10^{18}$ or this fraction shows that it cannot be completed.

Mathematical Approach

We first limit the possible values of $\sigma(n)/n$ and then build $n$ prime power by prime power, every step forced by divisibility.

1. Small cases

A sieve that adds every $d$ to $\sigma$ of each multiple of $d$ finds five numbers with a half-integral quotient up to $10^5$:

$\sigma(2)=3$, so $\sigma(2)/2=\tfrac32$;
$\sigma(24)=1+2+3+4$${}+6+8+12+24=60$, so $\sigma(24)/24=\tfrac52$;
$4320=2^5\cdot3^3\cdot5$, $4680=2^3\cdot3^2\cdot5\cdot13$ and $26208=2^5\cdot3^2\cdot7\cdot13$ all have the quotient $\tfrac72$.

Up to $10^{10}$ there are only $11$ such numbers, the first terms of the OEIS sequence A159907, ^[2] so a scan of all numbers up to $10^{18}$ is hopeless. All of them are even, and this holds in general: $\sigma(n)/n=k+\tfrac12$ means $2\sigma(n)=(2k+1)\,n$; the left side is even and $2k+1$ is odd, so $n$ must be even.

2. The quotient is a product over prime powers

The function $\sigma$ is multiplicative: $\sigma(ab)=\sigma(a)\,\sigma(b)$ whenever $\gcd(a,b)=1$, because every divisor of $ab$ is, in exactly one way, a divisor of $a$ times a divisor of $b$. ^[3] Dividing by $n$ therefore gives

$$\frac{\sigma(n)}{n}=\prod_{p\mid n}\frac{\sigma(p^{e_p})}{p^{e_p}},$$

where $e_p$ is the exponent of $p$ in $n$, and for a single prime power

$$\frac{\sigma(p^e)}{p^e}=1+\frac1p+\frac1{p^2}+\dots+\frac1{p^e}.$$

For example, $24=2^3\cdot3$ gives $\tfrac{15}8\cdot\tfrac43=\tfrac52$, and $4320=2^5\cdot3^3\cdot5$ gives $\tfrac{63}{32}\cdot\tfrac{40}{27}\cdot\tfrac65=\tfrac72$. Three facts about one such factor are used throughout: it is larger than $1$ as soon as $e\ge1$, so $\sigma(m)/m\gt1$ for every $m\gt1$; it grows with $e$; and it stays below the sum of the whole geometric series, $1/(1-1/p)=p/(p-1)$.

3. Only six values are possible

The product of the first $15$ primes, $2\cdot3\cdot5\cdots47$${}=614\,889\,782\,588\,491\,410$, is below $10^{18}$, but multiplying it by the next prime $53$ exceeds $10^{18}$. So a number $n\le10^{18}$ has at most $15$ distinct primes. For $n\gt1$, Section 2 shows that $\sigma(n)/n$ is smaller than the product of $p/(p-1)$ over the primes of $n$. Since $p/(p-1)$ decreases as $p$ grows, and the $i$-th smallest prime of $n$ is at least the $i$-th prime, this product is at most

$$\prod_{p\le47}\frac{p}{p-1}\approx7.2096.$$

Hence $k+\tfrac12\lt7.21$, and $k\ge1$ because $\sigma(n)/n\ge1$. The quotient must therefore be one of the six targets

$$t=\frac A2,\qquad A\in\{3,5,7,9,11,13\}.$$

In fact even $\tfrac{13}2$ never occurs. The smallest $n\le10^{18}$ with the largest quotient is superabundant, that is, its quotient beats that of every smaller number, and by a theorem of Alaoglu and Erdős such numbers have the form $2^{e_1}3^{e_2}5^{e_3}\cdots$ with $e_1\ge e_2\ge\dots$; ^[4] a search over these gives $\sigma(n)/n\lt6.4452$ for all $n\le10^{18}$. The programs search $\tfrac{13}2$ anyway and find nothing.

4. The fraction that is still missing

Fix a target $t=A/2$ and a solution $N\le10^{18}$, so $\sigma(N)/N=t$. Call $n$ a full part of $N$ if $n$ divides $N$ and contains each of its primes with the same exponent as $N$ (a unitary divisor); then $m=N/n$ is coprime to $n$. For example, $8$ is a full part of $24$, but $4$ is not. By multiplicativity $\sigma(N)/N$ is the product of $\sigma(n)/n$ and $\sigma(m)/m$, so the rest $m$ has to supply exactly

$$\frac{\sigma(m)}{m}=t\cdot\frac{n}{\sigma(n)}=\frac uv,\qquad \gcd(u,v)=1.$$

This reduced fraction, the missing fraction, depends only on $t$ and $n$, not on $N$; the programs store it as the pair rn, rs. Cross-multiplying gives the key equation

$$u\,m=v\,\sigma(m).$$

5. The denominator divides the rest

The number $v$ divides $u\,m$ and shares no prime with $u$, so it divides $m$. Four consequences drive the whole search:

Size. $m\ge v$, so $N=n\,m\ge n\,v$. If $n\,v\gt10^{18}$, no solution has $n$ as a full part.
New primes. Every prime $p$ of $v$ divides $m$, which is coprime to $n$. If such a $p$ already divides $n$, there is no solution.
Exponents. If $p^e$ is the exact power of $p$ in $v$, then $p^e$ divides $m$, so $N$ contains $p$ at least $e$ times.
The end. $u/v=\sigma(m)/m$ is at least $1$, with equality only for $m=1$. So $u\lt v$ is impossible, and $u=v=1$ means $N=n$.

For $t=\tfrac52$ and $n=8$ the missing fraction is $\tfrac52\cdot\tfrac8{15}=\tfrac43$: every solution with the full part $8$ is a multiple of $3$ and at least $8\cdot3=24$, and $24$ itself is one.

6. One step of the search

The search starts from $n=1$, whose missing fraction is $t=A/2$ itself, already in lowest terms because $A$ is odd. At a partial number $n$, the product of the prime powers chosen so far, with $v\gt1$, let $p$ be the smallest prime of $v$ and $e$ its exponent in $v$. By Section 5, every solution with the full part $n$ contains $p$ with some exponent $f\ge e$, so the search tries $f=e,e+1,\dots$ and moves to

$$n'=n\,p^f,\qquad \frac{u'}{v'}=\frac uv\cdot\frac{p^f}{\sigma(p^f)},$$

where the new fraction is again reduced by its greatest common divisor. The prime $p$ disappears from the denominator: $\sigma(p^f)=1+p+\dots+p^f$ leaves the remainder $1$ on division by $p$, and $p$ does not divide $u$, so the numerator $u\,p^f$ contains $p$ exactly $f$ times and the denominator $v\,\sigma(p^f)$ exactly $e\le f$ times. Three tests limit the exponents:

if $n\,p^f\gt10^{18}$, every larger exponent is too large as well, and the loop stops;
if $u'\lt v'$, the loop also stops: the factor $p^f/\sigma(p^f)$ shrinks as $f$ grows (Section 2), so larger exponents only push the fraction further below $1$;
if $n'\,v'\gt10^{18}$, this exponent is skipped, but the loop goes on, because $v'$ does not change monotonically with $f$.

Any prime of $v$ would be correct here; always taking the smallest gives every solution a single path through the search (Section 8).

7. Worked example: the target 5/2 up to 100

Here is the complete search for $t=\tfrac52$ with the bound $100$ in place of $10^{18}$. The root $n=1$ has the missing fraction $\tfrac52$, so $p=2$ and $e=1$.

$f=1$: $n=2$ with the fraction $\tfrac52\cdot\tfrac23=\tfrac53$, so the next prime is $3$. With $3^1$ we reach $n=6$ and $\tfrac53\cdot\tfrac34=\tfrac54$; the denominator $4$ asks for the prime $2$, which already divides $6$, so this branch is dead. With $3^2$ and $3^3$ we get $n=18$ with $\tfrac{15}{13}$ and $n=54$ with $\tfrac98$; the products $n\,v=234$ and $432$ exceed $100$, so both are skipped, and $3^4$ gives $n=162\gt100$, which ends the loop.
$f=2$: $n=4$ with $\tfrac52\cdot\tfrac47=\tfrac{10}7$, next prime $7$. Then $n=28$ has the fraction $\tfrac{10}7\cdot\tfrac78=\tfrac54$ and $n\,v=112\gt100$, so it is skipped, and $7^2$ gives $n=196\gt100$.
$f=3$: $n=8$ with $\tfrac43$. With $3^1$ we reach $n=24$ and $\tfrac43\cdot\tfrac34=1$: a solution. With $3^2$ the fraction $\tfrac43\cdot\tfrac9{13}=\tfrac{12}{13}$ is below $1$, so the loop stops.
$f=4,5,6$: $n=16,32,64$ with the fractions $\tfrac{40}{31}$, $\tfrac{80}{63}$ and $\tfrac{160}{127}$; the products $n\,v=496$, $2016$ and $8128$ all exceed $100$. Finally $2^7=128\gt100$.

The search visits only the six partial numbers $1,2,6,4,8,24$ and finds the single solution $24$. A longer path shows the exponent rule at work. For $t=\tfrac72$, the partial number $32$ has the fraction $\tfrac72\cdot\tfrac{32}{63}=\tfrac{16}9$, so $3$ must appear at least twice. With $3^3$ the fraction becomes $\tfrac{16}9\cdot\tfrac{27}{40}=\tfrac65$ at $n=864$, and $5^1$ completes $4320=2^5\cdot3^3\cdot5$. With $3^2$ instead, $n=288$ has the fraction $\tfrac{16}{13}$; the prime $13$ leads to $\tfrac87$ at $n=3744$, and $7^1$ completes $26208$.

8. Every solution is found exactly once

Let $N\le10^{18}$ be a solution for the target $t$. The search walks along a chain of full parts of $N$, starting from the root $n=1$, and counts $N$ at its end. Suppose the search visits a full part $n\ne N$ of $N$ with missing fraction $u/v$ and $v\gt1$. By Section 5 the smallest prime $p$ of $v$ divides $m=N/n$ but not $n$, so the branch is not declared dead, and $N$ contains $p$ with some exponent $f\ge e$. For every exponent $g\le f$ we have $n\,p^g\le N\le10^{18}$, and the fraction after $p^g$ is at least the fraction after $p^f$, which equals $\sigma(m')/m'\ge1$ with $m'=N/(n\,p^f)$. So no stop rule ends the loop before $f$. At $f$ the size test passes, because $n\,p^f$ is again a full part of $N$ and Section 5 gives $n\,p^f\,v'\le N$. The search therefore moves on to the full part $n\,p^f$, which has one prime more than $n$, and after at most as many steps as $N$ has distinct primes it reaches $n=N$, where $u=v=1$, and adds $N$. A full part $n\ne N$ with $v=1$ would leave an odd rest $m$ with $\sigma(m)/m=u\ge2$; Section 9 shows that no such $m$ exists in our range.

Conversely, the search adds a number only when its missing fraction equals $1$, that is, when $\sigma(n)/n=t$, and only after the test $n\,v=n\le10^{18}$. No number is added twice: along a path the prime chosen at each step is determined by the fraction, two branches of one loop differ in the exponent of that prime, which no later step changes because a prime of $n$ is never added again, and different targets have different quotients. The sum of the six searches is therefore the requested sum.

9. When the denominator is 1

One state is not handled by these rules: a missing fraction $u/1$ with $u\ge2$. Every partial number after the root is even, because the denominator $2$ of the root forces the prime $2$ first, so the rest $m$ would have to be an odd number with $\sigma(m)=u\,m$. Such numbers are called odd multiply perfect numbers, and apart from $1$ none is known. ^[5] The programs simply stop at such a state, and for the bound $10^{18}$ this is safe:

The rest $m\le10^{18}/2$ is odd, so it has at most $14$ distinct primes, because the first $15$ odd primes multiply to about $1.6\cdot10^{19}$. As in Section 3, $\sigma(m)/m$ is below the product of $p/(p-1)$ over the odd primes $p\le47$, about $3.605$, so $u$ is $2$ or $3$.
For $u=2$, $m$ would be an odd perfect number, and an odd perfect number would have to exceed $10^{1500}$. ^[6]
For $u=3$, $\sigma(m)=3m$ is odd. For odd $m$, every factor $\sigma(p^e)=1+p+\dots+p^e$ is a sum of $e+1$ odd numbers, so $\sigma(m)$ is odd only if every exponent is even, that is, if $m$ is a square. ^[3] An odd triperfect number would have to be a square larger than $10^{70}$ with at least $12$ distinct primes. ^[5]

The search reaches such a state $14$ times at the bound $10^{18}$, twelve times with $u=2$ and twice with $u=3$. For example, for $t=\tfrac72$ the partial number $4$ has the fraction $\tfrac72\cdot\tfrac47=2$, and for $t=\tfrac92$ the partial number $2$ has $\tfrac92\cdot\tfrac23=3$. For this article the case $u=3$ was also checked directly: none of the $353\,553\,391$ odd numbers $s\le707\,106\,781$ satisfies $\sigma(s^2)=3s^2$, so no odd square up to $5\cdot10^{17}$ is triperfect.

10. Smallest prime factors of the denominators

Every step needs the smallest prime of $v$, where $v\le5\cdot10^{17}$ because $n$ is even after the root and $n\,v\le10^{18}$. A primality test first handles the twelve primes $2,3,\dots,37$ and their multiples and then applies the Miller–Rabin test: write $v-1=2^s d$ with $d$ odd; a base $a$ that is not a multiple of $v$ proves $v$ composite if $a^d\not\equiv1$ and $a^{2^r d}\not\equiv-1\pmod v$ for all $0\le r\lt s$. ^[7] Such a base never proves a prime composite, and a base that is a multiple of $v$ is skipped; for instance $28178=2\cdot73\cdot193$ is skipped for the prime $73$. For every composite number below $2^{64}$, one of the seven bases $2$, $325$, $9375$, $28178$, $450775$, $9780504$ and $1795265022$ proves it composite, so the test is exact in this range. ^[8] A composite $v$ is split by Pollard's rho method, which iterates $x\mapsto x^2+c$ modulo $v$ from a random start and compares the sequence with a copy that runs twice as fast; the greatest common divisor of their difference and $v$ usually reveals a factor after about $\sqrt q$ steps, where $q$ is the smallest prime factor, and a new random $c$ is tried when it only finds $v$ itself. ^[9] Splitting recursively until only primes remain and taking the smallest one gives the smallest prime factor, whatever random numbers were used.

How the Code Works

The C++ program keeps one Solver object per target. Its fields are the bound limit, the odd numerator target_num $=A$, the running sum and the cache min_factor_cache, a std::map from a denominator to its smallest prime. run calls dfs(1, A, 2), the root of Section 6. dfs(n, rn, rs) receives a partial number with its missing fraction $u/v$ and applies the tests of Sections 5 and 9 in this order: it returns if $n\,v$ exceeds the bound, adds $n$ to sum if $u=v$, and returns if $u\lt v$ or $v=1$. Otherwise get_min_factor supplies the smallest prime $p$ of $v$, from the cache or from min_factor_u64; the call returns if $p$ divides $n$, and it counts the exponent $e$ of $p$ in $v$. The loop over $i=1,2,\dots$ updates p_pow $=p^i$ and sigma $=\sigma(p^i)$ together. It stops when the next power would exceed the bound, skips $i\lt e$, stops when $n\,p^i$ exceeds the bound, reduces the new fraction with gcd_u128, stops when the fraction falls below $1$, skips the exponent when $n'\,v'$ exceeds the bound, and otherwise calls dfs for $n'$: these are the tests of Section 6.

The factoring functions follow Section 10. is_prime_u64 divides by the twelve primes up to $37$ and then runs the witness test for the seven bases, with modular products formed in 128 bits by mul_mod_u64 and powers by pow_mod_u64. pollard_rho returns $2$ or $3$ at once for multiples of these primes and otherwise takes its constant $c$ and its start value from rand_u64, a SplitMix64 generator seeded from the clock. min_factor_u64 returns $2$ for even numbers and otherwise calls min_factor_rec, which splits recursively and returns the smaller of the smallest prime factors of the two parts. solve runs the six targets $A=3,5,\dots,13$ and adds their sums, and main accepts -n LIMIT to replace the default bound kDefaultLimit $=10^{18}$ and prints the 128-bit total with to_string_u128.

In a visited state $n\le10^{18}$ and $n\,v\le10^{18}$, so $n$ and $v$ fit in 64 bits, and so does $u\le\tfrac{13}2\,v$. The products $u\,p^i$ and $v\,\sigma(p^i)$ before the reduction and the product $n'\,v'$ of the size test are formed in unsigned 128-bit arithmetic. An instrumented copy of the search at the bound $10^{18}$ found the unreduced products always below $2^{63}$ and $n'\,v'$ always below $2^{120}$, so nothing overflows, and the numerators and denominators passed to the next call are below $2^{32}$. The Python version mirrors the C++ program function by function, with unbounded integers, a dictionary as cache and solve(limit), which runs the six targets through the helper worker. The Java version keeps n, rn and rs in long variables but uses BigInteger for the products inside dfs, for the running sum and in mulMod; its cache is a HashMap filled through computeIfAbsent, and its bound LIMIT is fixed. Its solve ignores any exception raised while the result of a target is collected, so a failed target would silently be missing from the total; none fails in the runs below.

Complexity and Verification

No formula predicts the size of the search tree, so we measured it. At the bound $10^{18}$, dfs is called $12\,979$ times in total: $2$, $178$, $934$, $2603$, $4852$ and $4410$ times for $A=3,5,\dots,13$. These searches find $1$, $3$, $9$, $7$, $2$ and $0$ solutions, $22$ in all. No path chooses more than $9$ prime powers, and $4295$ branches end because the next prime already divides $n$. Within one call the exponent loop runs at most $59$ times, because $p^f\le10^{18}\lt2^{60}$; in total the six searches reduce $51\,295$ fractions with a greatest common divisor. The six caches hold $7453$ entries for $4063$ different denominators, the largest being $2^{30}-1$. For numbers of this size a Miller–Rabin test costs a few hundred modular multiplications, and a composite one has a prime factor below $2^{16}$, which Pollard's rho finds after a few hundred steps on average. Altogether this is a few million operations: on a desktop computer the C++ program finishes in a few hundredths of a second, and the Python and Java versions in well under a second. The memory is dominated by the caches, a few thousand entries.

The programs contain no self-test, so the following checks were made for this article. An independently written segmented sieve of $\sigma(n)$ for every $n\le10^{10}$ finds exactly the $11$ numbers of Section 1, with the sums $26$, $9026$, $35234$, $335830146$ and $8919474306$ for $n\le10^2$, $10^4$, $10^5$, $10^9$ and $10^{10}$. The C++ program with -n, the Python function solve(limit) and a copy of the Java program with a changed bound give the same sums for these bounds, and at $10^{12}$ all three give $362349699202$. An exact-integer copy of the search, factorising with the computer algebra system SymPy, visits the same states and confirms the counts above, the integer sizes of the previous section and the $14$ stops of Section 9. The searches described in Sections 3 and 9 were also made for this article. At the full bound the three programs print the same value.

Footnotes and References

The references below supply the background facts used above. The missing fraction, the forced next prime, the size bound $n\,v\le10^{18}$ and the proof that every solution is found exactly once are derived in this article.

Project Euler 241 — Perfection Quotients. The official statement defines σ(n) as the sum of the divisors of n, with σ(6) = 12, calls p(n) = σ(n)/n the perfection quotient and asks for the sum of all positive integers n ≤ 10¹⁸ for which p(n) = k + 1/2 with an integer k.
OEIS A159907 — Numbers with half-integral abundancy index. It lists the numbers m with σ(m)/m = k + 1/2, starting 2, 24, 4320, 4680, 26208, 8910720, as in Section 1; a comment notes that all terms are even.
Wikipedia — Divisor function. The function σ is multiplicative, σ(p^a) = 1 + p + … + p^a for a prime power, and σ(n) is odd if and only if n is a square or twice a square.
Wikipedia — Superabundant number. A number n is superabundant if σ(m)/m < σ(n)/n for every m < n. Alaoglu and Erdős (1944) proved that a superabundant number is a product of the first k primes with non-increasing exponents a₁ ≥ a₂ ≥ … ≥ aₖ.
Wikipedia — Multiply perfect number. A number n is k-perfect if σ(n) = kn. It is unknown whether there are odd multiply perfect numbers other than 1, and an odd triperfect number would have to be a square exceeding 10⁷⁰ with at least 12 distinct prime factors.
Wikipedia — Perfect number. It is unknown whether any odd perfect number exists; Ochem and Rao (2012) proved that an odd perfect number would have to be greater than 10¹⁵⁰⁰.
Wikipedia — Miller–Rabin primality test. For odd n and a base a coprime to n, write n − 1 = 2ˢ·d with d odd; then n is a strong probable prime to the base a if a^d ≡ 1 or a^(2ʳ·d) ≡ −1 (mod n) for some 0 ≤ r < s. Every prime passes for every such base, and a base that a composite number fails is a witness of its compositeness.
Deterministic variants of the Miller–Rabin primality test. The table of best known base sets lists the seven bases 2, 325, 9375, 28178, 450775, 9780504 and 1795265022, found by Jim Sinclair, as correct for every number below 2⁶⁴. Its remarks add that a base equal to 0 after reduction modulo n must count as a pass.
Wikipedia — Pollard's rho algorithm. The method iterates x² + c modulo n, compares the sequence with one that moves twice as fast through gcd(|x − y|, n), retries with other parameters when this gcd equals n, and heuristically needs about √p steps for the smallest prime factor p.

Problem 241 source code

C++

#include <algorithm>
#include <atomic>
#include <chrono>
#include <cstdint>
#include <functional>
#include <future>
#include <iostream>
#include <map>
#include <string>
#include <thread>
#include <utility>
#include <vector>

using u64 = std::uint64_t;
using u128 = __uint128_t;

namespace {

constexpr u64 kDefaultLimit = 1000000000000000000ULL;  // 1e18

u64 gcd_u64(u64 a, u64 b) {
    while (b) {
        u64 t = a % b;
        a = b;
        b = t;
    }
    return a;
}

u128 gcd_u128(u128 a, u128 b) {
    while (b) {
        u128 t = a % b;
        a = b;
        b = t;
    }
    return a;
}

u64 mul_mod_u64(u64 a, u64 b, u64 mod) {
    return static_cast<u64>((static_cast<u128>(a) * b) % mod);
}

u64 pow_mod_u64(u64 a, u64 d, u64 mod) {
    u64 r = 1;
    while (d) {
        if (d & 1) r = mul_mod_u64(r, a, mod);
        a = mul_mod_u64(a, a, mod);
        d >>= 1;
    }
    return r;
}

bool is_prime_u64(u64 n) {
    if (n < 2) return false;
    static u64 small_primes[] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37};
    for (u64 p : small_primes) {
        if (n == p) return true;
        if (n % p == 0) return false;
    }

    u64 d = n - 1, s = 0;
    while ((d & 1) == 0) {
        d >>= 1;
        ++s;
    }

    auto witness = [&](u64 a) -> bool {
        if (a % n == 0) return false;
        u64 x = pow_mod_u64(a, d, n);
        if (x == 1 || x == n - 1) return false;
        for (u64 i = 1; i < s; ++i) {
            x = mul_mod_u64(x, x, n);
            if (x == n - 1) return false;
        }
        return true;
    };

    static u64 bases[] = {2ULL, 325ULL, 9375ULL, 28178ULL, 450775ULL, 9780504ULL, 1795265022ULL};
    for (u64 a : bases) {
        if (witness(a)) return false;
    }
    return true;
}

u64 splitmix64(u64& x) {
    u64 z = (x += 0x9e3779b97f4a7c15ULL);
    z = (z ^ (z >> 30)) * 0xbf58476d1ce4e5b9ULL;
    z = (z ^ (z >> 27)) * 0x94d049bb133111ebULL;
    return z ^ (z >> 31);
}

thread_local u64 rng_state =
    static_cast<u64>(std::chrono::high_resolution_clock::now().time_since_epoch().count()) ^
    0x9e3779b97f4a7c15ULL;

u64 rand_u64(u64 lo, u64 hi) {
    u64 r = splitmix64(rng_state);
    return lo + (hi > lo ? (r % (hi - lo + 1)) : 0);
}

u64 pollard_rho(u64 n) {
    if ((n & 1ULL) == 0) return 2;
    if (n % 3ULL == 0) return 3;

    while (true) {
        u64 c = rand_u64(1, n - 1);
        u64 x = rand_u64(0, n - 1);
        u64 y = x;
        u64 d = 1;

        auto f = [&](u64 v) { return (mul_mod_u64(v, v, n) + c) % n; };

        while (d == 1) {
            x = f(x);
            y = f(f(y));
            u64 diff = (x > y) ? (x - y) : (y - x);
            d = gcd_u64(diff, n);
        }
        if (d != n) return d;
    }
}

u64 min_factor_u64(u64 n);

u64 min_factor_rec(u64 n) {
    if (n == 1) return 1;
    if (is_prime_u64(n)) return n;
    u64 d = pollard_rho(n);
    u64 a = min_factor_rec(d);
    u64 b = min_factor_rec(n / d);
    return std::min(a, b);
}

u64 min_factor_u64(u64 n) {
    if ((n & 1ULL) == 0) return 2;
    return min_factor_rec(n);
}

struct Solver {
    u64 limit;
    u64 target_num;  // A (odd), target ratio is A/2.
    u128 sum = 0;
    std::map<u64, u64> min_factor_cache;

    explicit Solver(u64 limit_, u64 target_num_) : limit(limit_), target_num(target_num_) {
    }

    u64 get_min_factor(u64 n) {
        auto it = min_factor_cache.find(n);
        if (it != min_factor_cache.end()) return it->second;
        u64 res = min_factor_u64(n);
        min_factor_cache.emplace(n, res);
        return res;
    }

    void dfs(u64 n, u64 rn, u64 rs) {
        if (static_cast<u128>(n) * rs > limit) return;
        if (rn == rs) {
            sum += n;
            return;
        }
        if (rn < rs) return;  // overshot target
        if (rs == 1) return;  // no forced prime left; with n even this can't finish within limit

        u64 p = get_min_factor(rs);
        if (n % p == 0) return;  // p already used; would have been chosen with larger exponent earlier

        u64 tmp = rs;
        int e = 0;
        while (tmp % p == 0) {
            tmp /= p;
            ++e;
        }

        u128 p_pow = 1;
        u128 sigma = 1;
        for (int i = 1; ; ++i) {
            if (p_pow > limit / p) break;
            p_pow *= p;
            sigma += p_pow;

            if (i < e) continue;
            u128 n2 = static_cast<u128>(n) * p_pow;
            if (n2 > limit) break;

            u128 rn2 = static_cast<u128>(rn) * p_pow;
            u128 rs2 = static_cast<u128>(rs) * sigma;
            u128 g = gcd_u128(rn2, rs2);
            rn2 /= g;
            rs2 /= g;

            if (rn2 < rs2) break;

            if (n2 * rs2 > limit) continue;
            dfs(static_cast<u64>(n2), static_cast<u64>(rn2), static_cast<u64>(rs2));
        }
    }

    u128 run() {
        dfs(1ULL, target_num, 2ULL);
        return sum;
    }
};

std::string to_string_u128(u128 x) {
    if (x == 0) return "0";
    std::string s;
    while (x > 0) {
        int digit = static_cast<int>(x % 10);
        s.push_back(static_cast<char>('0' + digit));
        x /= 10;
    }
    std::reverse(s.begin(), s.end());
    return s;
}

u128 solve(u64 limit, int threads) {
    std::vector<u64> targets;
    for (u64 a = 3; a <= 13; a += 2) targets.push_back(a);

    if (threads <= 1 || targets.size() == 1) {
        u128 total = 0;
        for (u64 a : targets) {
            Solver solver(limit, a);
            total += solver.run();
        }
        return total;
    }

    int worker_count = std::min<int>(threads, static_cast<int>(targets.size()));
    std::vector<u128> partial(targets.size(), 0);
    std::atomic<std::size_t> idx{0};

    auto worker = [&]() {
        while (true) {
            std::size_t i = idx.fetch_add(1);
            if (i >= targets.size()) return;
            Solver solver(limit, targets[i]);
            partial[i] = solver.run();
        }
    };

    std::vector<std::thread> pool;
    pool.reserve(worker_count);
    for (int i = 0; i < worker_count; ++i) {
        pool.emplace_back(worker);
    }
    for (auto& t : pool) t.join();

    u128 total = 0;
    for (u128 v : partial) total += v;
    return total;
}

void usage(const char* argv0) {
    std::cerr << "Usage: " << argv0 << " [-n LIMIT] [-t THREADS]\n";
}

}  // namespace

int main(int argc, char** argv) {
    std::ios::sync_with_stdio(false);
    std::cin.tie(nullptr);

    u64 limit = kDefaultLimit;
    int threads = static_cast<int>(std::thread::hardware_concurrency());
    if (threads <= 0) threads = 1;

    for (int i = 1; i < argc; ++i) {
        std::string arg = argv[i];
        if (arg == "-n" && i + 1 < argc) {
            limit = std::stoull(argv[++i]);
        } else if (arg == "-t" && i + 1 < argc) {
            threads = std::max(1, std::stoi(argv[++i]));
        } else {
            usage(argv[0]);
            return 1;
        }
    }

    u128 ans = solve(limit, threads);
    std::cout << to_string_u128(ans) << "\n";
    return 0;
}

Python

import math
import random
import multiprocessing

def gcd(a, b):
    while b:
        a, b = b, a % b
    return a

def is_prime(n):
    if n < 2: return False
    if n in (2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37): return True
    if any(n % p == 0 for p in (2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37)):
        return False

    d = n - 1
    s = 0
    while d % 2 == 0:
        d >>= 1
        s += 1

    def witness(a):
        if a % n == 0: return False
        x = pow(a, d, n)
        if x == 1 or x == n - 1: return False
        for _ in range(1, s):
            x = pow(x, 2, n)
            if x == n - 1: return False
        return True

    for a in (2, 325, 9375, 28178, 450775, 9780504, 1795265022):
        if witness(a): return False
    return True

def pollard_rho(n):
    if n % 2 == 0: return 2
    if n % 3 == 0: return 3
    while True:
        c = random.randint(1, n - 1)
        x = random.randint(0, n - 1)
        y = x
        d = 1
        
        def f(v):
            return (pow(v, 2, n) + c) % n

        while d == 1:
            x = f(x)
            y = f(f(y))
            d = gcd(abs(x - y), n)
            
        if d != n:
            return d

def min_factor_rec(n):
    if n == 1: return 1
    if is_prime(n): return n
    d = pollard_rho(n)
    return min(min_factor_rec(d), min_factor_rec(n // d))

def min_factor(n):
    if n % 2 == 0: return 2
    return min_factor_rec(n)

class Solver:
    def __init__(self, limit, target_num):
        self.limit = limit
        self.target_num = target_num
        self.sum = 0
        self.min_factor_cache = {}

    def get_min_factor(self, n):
        if n in self.min_factor_cache:
            return self.min_factor_cache[n]
        res = min_factor(n)
        self.min_factor_cache[n] = res
        return res

    def dfs(self, n, rn, rs):
        if n * rs > self.limit: return
        if rn == rs:
            self.sum += n
            return
        if rn < rs: return
        if rs == 1: return
        
        p = self.get_min_factor(rs)
        if n % p == 0: return
        
        tmp = rs
        e = 0
        while tmp % p == 0:
            tmp //= p
            e += 1
            
        p_pow = 1
        sigma = 1
        i = 1
        while True:
            if p_pow > self.limit // p: break
            p_pow *= p
            sigma += p_pow
            
            if i >= e:
                n2 = n * p_pow
                if n2 > self.limit: break
                
                rn2 = rn * p_pow
                rs2 = rs * sigma
                g = gcd(rn2, rs2)
                rn2 //= g
                rs2 //= g
                
                if rn2 < rs2: break
                if n2 * rs2 <= self.limit:
                    self.dfs(n2, rn2, rs2)
                    
            i += 1

    def run(self):
        self.dfs(1, self.target_num, 2)
        return self.sum

def worker(args):
    limit, target = args
    solver = Solver(limit, target)
    return solver.run()

def solve(limit=10**18):
    targets = list(range(3, 14, 2))
    pool_args = [(limit, t) for t in targets]
    
    threads = min(multiprocessing.cpu_count(), len(targets))
    if threads > 1:
        with multiprocessing.Pool(threads) as pool:
            results = pool.map(worker, pool_args)
    else:
        results = map(worker, pool_args)
        
    return str(sum(results))

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

Java

import java.math.BigInteger;
import java.util.*;
import java.util.concurrent.*;

public class Euler241 {
    static final long LIMIT = 1000000000000000000L;

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

    static BigInteger gcd(BigInteger a, BigInteger b) {
        return a.gcd(b);
    }

    static long mulMod(long a, long b, long mod) {
        return BigInteger.valueOf(a).multiply(BigInteger.valueOf(b)).mod(BigInteger.valueOf(mod)).longValue();
    }

    static long powMod(long a, long d, long mod) {
        long r = 1;
        while (d > 0) {
            if ((d & 1) != 0)
                r = mulMod(r, a, mod);
            a = mulMod(a, a, mod);
            d >>= 1;
        }
        return r;
    }

    static boolean isPrime(long n) {
        if (n < 2)
            return false;
        long[] smallPrimes = { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37 };
        for (long p : smallPrimes) {
            if (n == p)
                return true;
            if (n % p == 0)
                return false;
        }

        long d = n - 1;
        int s = 0;
        while ((d & 1) == 0) {
            d >>= 1;
            s++;
        }

        long[] bases = { 2L, 325L, 9375L, 28178L, 450775L, 9780504L, 1795265022L };
        for (long a : bases) {
            if (a % n == 0)
                continue;
            long x = powMod(a, d, n);
            if (x == 1 || x == n - 1)
                continue;
            boolean composite = true;
            for (int i = 1; i < s; i++) {
                x = mulMod(x, x, n);
                if (x == n - 1) {
                    composite = false;
                    break;
                }
            }
            if (composite)
                return false;
        }
        return true;
    }

    static long splitmix64(long[] x) {
        long z = (x[0] += 0x9e3779b97f4a7c15L);
        z = (z ^ (z >>> 30)) * 0xbf58476d1ce4e5b9L;
        z = (z ^ (z >>> 27)) * 0x94d049bb133111ebL;
        return z ^ (z >>> 31);
    }

    static final ThreadLocal<long[]> rngState = ThreadLocal
            .withInitial(() -> new long[] { System.nanoTime() ^ 0x9e3779b97f4a7c15L });

    static long rand(long lo, long hi) {
        long r = splitmix64(rngState.get());
        if (r < 0)
            r = ~r;
        return lo + (hi > lo ? (r % (hi - lo + 1)) : 0);
    }

    static long pollardRho(long n) {
        if ((n & 1) == 0)
            return 2;
        if (n % 3 == 0)
            return 3;

        while (true) {
            long c = rand(1, n - 1);
            long x = rand(0, n - 1);
            long y = x;
            long d = 1;

            while (d == 1) {
                x = (mulMod(x, x, n) + c) % n;
                y = (mulMod(y, y, n) + c) % n;
                y = (mulMod(y, y, n) + c) % n;
                long diff = Math.abs(x - y);
                d = gcd(diff, n);
            }
            if (d != n)
                return d;
        }
    }

    static long minFactorRec(long n) {
        if (n == 1)
            return 1;
        if (isPrime(n))
            return n;
        long d = pollardRho(n);
        return Math.min(minFactorRec(d), minFactorRec(n / d));
    }

    static long minFactor(long n) {
        if ((n & 1) == 0)
            return 2;
        return minFactorRec(n);
    }

    static class Solver {
        long limit, targetNum;
        BigInteger sum = BigInteger.ZERO;
        Map<Long, Long> minFactorCache = new HashMap<>();

        Solver(long limit, long targetNum) {
            this.limit = limit;
            this.targetNum = targetNum;
        }

        long getMinFactor(long n) {
            return minFactorCache.computeIfAbsent(n, Euler241::minFactor);
        }

        void dfs(long n, long rn, long rs) {
            // Check n * rs > limit handling overflow with BigInteger
            if (BigInteger.valueOf(n).multiply(BigInteger.valueOf(rs)).compareTo(BigInteger.valueOf(limit)) > 0)
                return;
            if (rn == rs) {
                sum = sum.add(BigInteger.valueOf(n));
                return;
            }
            if (rn < rs)
                return;
            if (rs == 1)
                return;

            long p = getMinFactor(rs);
            if (n % p == 0)
                return;

            long tmp = rs;
            int e = 0;
            while (tmp % p == 0) {
                tmp /= p;
                e++;
            }

            BigInteger pPow = BigInteger.ONE;
            BigInteger sigma = BigInteger.ONE;
            BigInteger P = BigInteger.valueOf(p);
            BigInteger Limit = BigInteger.valueOf(limit);

            for (int i = 1;; ++i) {
                if (pPow.compareTo(Limit.divide(P)) > 0)
                    break;
                pPow = pPow.multiply(P);
                sigma = sigma.add(pPow);

                if (i < e)
                    continue;

                BigInteger n2 = BigInteger.valueOf(n).multiply(pPow);
                if (n2.compareTo(Limit) > 0)
                    break;

                BigInteger rn2 = BigInteger.valueOf(rn).multiply(pPow);
                BigInteger rs2 = BigInteger.valueOf(rs).multiply(sigma);
                BigInteger g = gcd(rn2, rs2);
                rn2 = rn2.divide(g);
                rs2 = rs2.divide(g);

                if (rn2.compareTo(rs2) < 0)
                    break;

                if (n2.multiply(rs2).compareTo(Limit) > 0)
                    continue;
                dfs(n2.longValue(), rn2.longValue(), rs2.longValue());
            }
        }

        BigInteger run() {
            dfs(1, targetNum, 2);
            return sum;
        }
    }

    public static String solve() {
        List<Long> targets = new ArrayList<>();
        for (long a = 3; a <= 13; a += 2)
            targets.add(a);

        int threads = Math.min(targets.size(), Math.max(1, Runtime.getRuntime().availableProcessors()));
        ExecutorService executor = Executors.newFixedThreadPool(threads);
        List<Future<BigInteger>> futures = new ArrayList<>();

        for (long target : targets) {
            futures.add(executor.submit(() -> {
                Solver solver = new Solver(LIMIT, target);
                return solver.run();
            }));
        }

        BigInteger total = BigInteger.ZERO;
        for (Future<BigInteger> f : futures) {
            try {
                total = total.add(f.get());
            } catch (Exception e) {
            }
        }
        executor.shutdown();

        return total.toString();
    }

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