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

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

The residual quotient is the real search invariant

Instead of testing $\sigma(n)/n$ from scratch at every node, the search keeps the reduced residual quotient

$$Q(n)=T\frac{n}{\sigma(n)}=\frac{u}{v},\qquad \gcd(u,v)=1.$$

At the root, $n=1$, so $Q(1)=T$. If we extend the current partial factorization by a new prime power $p^e$, then

$$Q(np^e)=Q(n)\cdot \frac{p^e}{\sigma(p^e)}.$$

Since $\sigma(p^e)/p^e>1$, each update multiplies $Q$ by a factor strictly smaller than 1. The target condition becomes

$$Q(n)=1 \iff T\frac{n}{\sigma(n)}=1 \iff \frac{\sigma(n)}{n}=T.$$

So the entire search is a controlled descent from $Q=T$ down to $Q=1$.

Why the denominator forces the next prime

Assume the current state is $Q(n)=u/v$ in lowest terms, and suppose that some unused cofactor $m$ finishes the construction. Then

$$1=Q(n)\frac{m}{\sigma(m)}=\frac{u}{v}\frac{m}{\sigma(m)},$$

which is equivalent to

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

Because $\gcd(u,v)=1$, this forces

$$v \mid m.$$

That divisibility is the key invariant. Every prime that appears in the current denominator must still appear later in the unfinished part of the number. Let $p$ be the smallest prime divisor of $v$. Then any completion must contain some power $p^e$, so the recursion branches only on that prime $p$.

There is also a lower bound on the exponent. If $p^a \mid v$, then the new numerator must contribute at least $p^a$ in order to cancel that denominator power after reduction. But in the update factor

$$\frac{p^e}{\sigma(p^e)},$$

the only source of $p$-adic valuation in the numerator is $p^e$ itself, because $\sigma(p^e)\equiv 1 \pmod p$. Therefore the branch must start at $e \ge a$.

Canonical prime order avoids duplicate factorizations

Once a prime has been inserted into the current partial factorization, the search never comes back to that prime later. This is not an arbitrary restriction. If a valid solution contains $p^e$, then that full exponent must be chosen precisely when the denominator first forces the prime $p$. Revisiting $p$ later would only recreate the same final number in a different order. The recursion therefore explores each admissible exponent exactly once, in a canonical order dictated by the changing denominator of $Q$.

Two strong pruning rules come directly from the invariant

First, if the reduced residual quotient satisfies $u<v$, then $Q(n)<1$. From that point onward every extra prime power multiplies by another factor $p^e/\sigma(p^e)<1$, so the branch can only move farther below 1. Such a branch is already lost.

Second, the divisibility $v\mid m$ implies $m\ge v$. Hence every completion $N=nm$ satisfies

$$N \ge n\,v.$$

So whenever $n\,v>10^{18}$, no completion can stay inside the problem limit, and the branch can be discarded immediately.

There is also monotonicity inside a fixed prime loop. For one prime $p$, the factor

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

decreases as $e$ grows. Therefore, once a chosen exponent pushes the reduced residual below 1, larger exponents cannot repair the branch, so the loop may stop at once.

Worked example: reaching the target $5/2$

Take $T=\frac52$. Start from

$$Q(1)=\frac52.$$

The denominator is 2, so the next forced prime is $2$. Choose the exponent $e=3$. Then

$$Q(2^3)=\frac52\cdot\frac{8}{1+2+4+8}=\frac52\cdot\frac{8}{15}=\frac43.$$

Now the denominator is 3, so the next forced prime is $3$. Taking exponent $e=1$ gives

$$Q(2^3\cdot 3)=\frac43\cdot\frac{3}{1+3}=\frac43\cdot\frac34=1.$$

Thus

$$\frac{\sigma(24)}{24}=\frac52,$$

so $24$ is a valid term. This example shows the exact rhythm of the search: read the denominator, force its smallest prime, try admissible exponents, reduce the residual quotient, and repeat until either $Q=1$ or a pruning rule fires.

How the Code Works

The C++, Python, and Java implementations follow the same number-theoretic plan. They solve the six target quotients independently and add the six partial sums.

Residual states and denominator factorization

Each recursive state stores the current partial number together with the reduced numerator and denominator of $Q$. To keep the denominator-driven search fast, the implementations repeatedly ask for the smallest prime factor of the current denominator. Those denominators can still be large 64-bit integers, so primality is tested with deterministic Miller-Rabin bases, and composite denominators are split with Pollard-Rho. A small cache avoids refactoring the same denominator many times.

Depth-first search over forced prime powers

At each node, the implementation extracts the smallest prime dividing the current denominator and counts its multiplicity. It then walks through exponents starting at that multiplicity, while updating $p^e$, $\sigma(p^e)$, and the reduced residual quotient. The branch stops as soon as one of the mathematical impossibility tests is triggered: the partial number exceeds the limit, the residual quotient drops below 1, or the lower bound $n\,v\le 10^{18}$ fails.

The implementation also rejects any attempt to reuse a prime that is already present in the partial factorization. That enforces the canonical prime order and removes duplicates automatically.

Independent target searches and parallel execution

The six half-integer targets are independent tasks, so they can be processed separately. The compiled implementations distribute them across worker threads, and the Python implementation uses multiple processes when possible. The arithmetic logic is the same in all three languages; only the surrounding execution model differs.

Complexity Analysis

There is no simple closed-form running time, because the algorithm is a heavily pruned search tree rather than a regular table-filling recurrence. The practical cost is determined by how many residual states survive the pruning rules and by how expensive it is to factor the denominators that appear along those states.

Memory use is modest. The recursion depth is the number of distinct primes chosen on the current branch, and the auxiliary storage is mainly a cache of previously factored denominators plus a few larger integer temporaries for overflow-safe arithmetic. The main gain is structural: instead of iterating through all $10^{18}$ candidates, the search visits only states that are compatible with the exact divisor-sum constraints.

Footnotes and References

Project Euler problem page: https://projecteuler.net/problem=241
Divisor function: Wikipedia - Divisor function
Multiplicative function: Wikipedia - Multiplicative function
Pollard's rho algorithm: Wikipedia - Pollard's rho algorithm
Miller-Rabin primality test: Wikipedia - Miller-Rabin primality test

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