Problem 241: Perfection Quotients

Mathematical Development

Definition 1 (Abundancy index). For a positive integer $n$, the abundancy index is $\sigma(n)/n$. We say $n$ is half-integer perfect when $\sigma(n)/n \in \mathbb{Z} + \tfrac{1}{2}$.

Theorem 1 (Characterisation). A positive integer $n > 1$ is half-integer perfect if and only if the rational number $2\sigma(n)/n$ is an odd integer $\ge 3$.

Proof. Suppose $\sigma(n)/n = k + \tfrac{1}{2}$ for some non-negative integer $k$. Then $2\sigma(n)/n = 2k + 1$, which is an odd integer. Since $\sigma(n) \ge 1 + n > n$ for $n > 1$ (as $n$ has at least two divisors), we have $\sigma(n)/n > 1$, hence $2k + 1 \ge 3$.

Conversely, if $2\sigma(n)/n = 2k + 1$ for some integer $k \ge 1$, then $\sigma(n)/n = k + \tfrac{1}{2}$. $\blacksquare$

Theorem 2 (Multiplicative decomposition). Write $n = 2^a \cdot m$ with $m$ odd and $a \ge 0$. Since $\sigma$ is multiplicative and $\gcd(2^a, m) = 1$:

The half-integer perfect condition $2\sigma(n)/n \in \{3, 5, 7, \ldots\}$ becomes:

Proof. The identity $\sigma(2^a) = 2^{a+1} - 1$ follows from the geometric series $\sum_{j=0}^{a} 2^j$. Multiplicativity of $\sigma$ for coprime arguments is standard. The algebraic manipulation is direct. $\blacksquare$

Lemma 1 (Divisibility constraint on $\sigma(m)$). For $a \ge 1$, the expression $(2^{a+1} - 1)\,\sigma(m) / (2^{a-1}\,m)$ can be an integer only if $2^{a-1} \mid \sigma(m)$.

Proof. Observe that $2^{a+1} - 1$ is odd and $m$ is odd, so $\gcd(2^{a-1},\, (2^{a+1} - 1) \cdot m) = 1$. For the fraction to be integral, the denominator $2^{a-1} \cdot m$ must divide the numerator $(2^{a+1} - 1)\,\sigma(m)$. Since $\gcd(2^{a-1}, (2^{a+1}-1)m) = 1$, the factor $2^{a-1}$ must divide $\sigma(m)$. $\blacksquare$

Lemma 2 (Abundancy upper bound). For $n = p_1^{e_1} \cdots p_k^{e_k}$,

Proof. Each factor $\sigma(p^e)/p^e = 1 + p^{-1} + \cdots + p^{-e}$ is a partial sum of the geometric series $\sum_{j=0}^{\infty} p^{-j} = p/(p-1)$, hence strictly less than $p/(p-1)$. $\blacksquare$

Corollary 1 (Pruning criterion). During a depth-first search that constructs $n$ by choosing odd prime factors in increasing order, if the product $\prod_{p \mid m,\, p \ge p_0} p/(p-1)$ (taken over all primes $\ge p_0$ dividing $m$ and all larger potential primes) is insufficient to bring $2\sigma(n)/n$ to any odd integer $\ge 3$, the branch can be pruned.

Editorial

The search also tries “denominator-targeted” primes: for each divisor $d$ of the current denominator, the values $d \pm 1$ and $2d \pm 1$ are tested for primality, since such primes are likely to cancel denominator factors. We build n by multiplying by prime powers p^e in increasing order of p. We then update fraction: new_num/new_den = (num * sig_pe) / (den * pe). Finally, reduce by gcd.

Pseudocode

Begin search: n = 1, ratio 2*sigma(1)/1 = 2/1
Build n by multiplying by prime powers p^e in increasing order of p
num/den = 2*sigma(n)/n in lowest terms
Update fraction: new_num/new_den = (num * sig_pe) / (den * pe)
Reduce by gcd

Complexity Analysis

• Time. The search tree explores factorisations of odd parts $m$ with $2^a m \le 10^{18}$. The abundancy bound (Lemma 2) and divisibility constraint (Lemma 1) provide strong pruning: empirically, $O(10^6)$ nodes are visited across all branches.
• Space. $O(\log N)$ for the recursion stack. The maximum depth equals the number of distinct prime factors of $n$, which is at most $\lfloor \log_3(10^{18}) \rfloor = 38$.

Answer

/*
 * Project Euler Problem 241: Perfection Quotients
 *
 * Find all n <= 10^18 with sigma(n)/n a half-integer, i.e., 2*sigma(n)/n
 * is an odd integer >= 3.  We build n by DFS over its prime factorisation,
 * tracking the ratio 2*sigma(n)/n as a reduced fraction num/den.
 */

#include <bits/stdc++.h>
using namespace std;

typedef long long ll;
typedef __int128 lll;

static const ll LIMIT = 1000000000000000000LL; // 10^18

ll my_gcd(ll a, ll b) {
    while (b) { a %= b; swap(a, b); }
    return a;
}

vector<ll> results;
vector<int> small_primes;

bool is_prime_check(ll n) {
    if (n < 2) return false;
    if (n < 4) return true;
    if (n % 2 == 0 || n % 3 == 0) return false;
    for (ll i = 5; i * i <= n; i += 6)
        if (n % i == 0 || n % (i + 2) == 0) return false;
    return true;
}

void search(ll n_val, ll num, ll den, ll min_prime) {
    if (den == 1 && num % 2 == 1 && num >= 3 && n_val > 1)
        results.push_back(n_val);

    if (num > 200LL * den) return;

    ll max_p = LIMIT / n_val;
    if (max_p < min_prime) return;

    auto try_prime = [&](ll p) {
        if (p < min_prime || p > max_p) return;
        ll pe = 1;
        ll sig_pe = 1;
        for (int e = 1; pe <= max_p / p; e++) {
            pe *= p;
            sig_pe += pe;
            ll new_n = n_val * pe;
            if (new_n > LIMIT) break;

            ll g1 = my_gcd(num, pe);
            ll g2 = my_gcd(sig_pe, den);
            ll a_num = num / g1, a_pe = pe / g1;
            ll a_sig = sig_pe / g2, a_den = den / g2;

            if (a_num > 0 && a_sig > (ll)4e18 / a_num) break;
            ll new_num = a_num * a_sig;
            if (a_den > 0 && a_pe > (ll)4e18 / a_den) break;
            ll new_den = a_den * a_pe;

            ll g3 = my_gcd(new_num, new_den);
            new_num /= g3;
            new_den /= g3;

            search(new_n, new_num, new_den, p + 1);
        }
    };

    for (int p : small_primes) {
        if (p > max_p) break;
        if (p < min_prime) continue;
        try_prime(p);
    }

    if (den > 1) {
        set<ll> candidates;
        for (ll d = 1; d * d <= den; d++) {
            if (den % d == 0) {
                for (ll c : {d - 1, d + 1, 2*d - 1, 2*d + 1,
                             den/d - 1, den/d + 1, 2*(den/d)-1, 2*(den/d)+1}) {
                    if (c >= min_prime && c <= max_p && c > 2000000)
                        candidates.insert(c);
                }
            }
        }
        for (ll c : candidates)
            if (is_prime_check(c))
                try_prime(c);
    }
}

int main() {
    ios_base::sync_with_stdio(false);

    vector<bool> sieve(2000001, true);
    for (int i = 2; i * i <= 2000000; i++)
        if (sieve[i])
            for (int j = i*i; j <= 2000000; j += i)
                sieve[j] = false;
    for (int i = 2; i <= 2000000; i++)
        if (sieve[i]) small_primes.push_back(i);

    search(1, 2, 1, 2);

    ll ans = 0;
    for (ll x : results) ans += x;
    cout << ans << endl;
    return 0;
}

"""
Project Euler Problem 241: Perfection Quotients

Find the sum of all n <= 10^18 such that sigma(n)/n is a half-integer,
i.e., 2*sigma(n)/n is an odd integer >= 3.

We perform a backtracking search over the prime factorisation of n,
maintaining 2*sigma(n)/n as a reduced fraction num/den.
"""

from math import gcd
from sympy import isprime

LIMIT = 10**18
results = []


def sigma_pe(p, e):
    """Return sigma(p^e) = 1 + p + ... + p^e."""
    return (p**(e + 1) - 1) // (p - 1)


def search(n_val, num, den, min_prime):
    """
    n_val: current partial product
    num/den: 2*sigma(n_val)/n_val in lowest terms
    min_prime: smallest prime allowed as next factor
    """
    if den == 1 and num % 2 == 1 and num >= 3 and n_val > 1:
        results.append(n_val)

    if num > 200 * den:
        return

    max_p = LIMIT // n_val
    if max_p < min_prime:
        return

    candidates = set()

    for p in small_primes:
        if p > max_p:
            break
        if p >= min_prime:
            candidates.add(p)

    if den > 1:
        divisors_of_den = []
        d = 1
        while d * d <= den:
            if den % d == 0:
                divisors_of_den.append(d)
                if d != den // d:
                    divisors_of_den.append(den // d)
            d += 1
        for d in divisors_of_den:
            for c in [d - 1, d + 1, 2 * d - 1, 2 * d + 1]:
                if min_prime <= c <= max_p and c > 1 and isprime(c):
                    candidates.add(c)

    for p in sorted(candidates):
        pe = 1
        sig = 1
        for e in range(1, 200):
            if pe > max_p // p:
                break
            pe *= p
            sig += pe

            new_n = n_val * pe
            if new_n > LIMIT:
                break

            g1 = gcd(num, pe)
            g2 = gcd(sig, den)
            new_num = (num // g1) * (sig // g2)
            new_den = (den // g2) * (pe // g1)
            g3 = gcd(new_num, new_den)
            new_num //= g3
            new_den //= g3

            search(new_n, new_num, new_den, p + 1)


def sieve(limit):
    is_p = [True] * (limit + 1)
    is_p[0] = is_p[1] = False
    for i in range(2, int(limit**0.5) + 1):
        if is_p[i]:
            for j in range(i * i, limit + 1, i):
                is_p[j] = False
    return [i for i in range(2, limit + 1) if is_p[i]]


small_primes = sieve(2_000_000)

search(1, 2, 1, 2)

print(sum(results))