k + sigma(k) is perfect

k_plus_sigma_k_is_perfect.py

"""
The number 10 has a special property. When 10 is added to the sum 
of its divisors, the result is 28, which is a perfect number.

    10 + sigma(10) = 10 + 1 + 2 + 5 + 10 = 28

Are there other perfect numbers that can be expressed as the sum 
of a number and the sum of its divisors?

Recall that an integer n is called perfect if sigma(n) = 2n, 
where sigma(n) is the sum of the divisors of n.

An even number n is perfect if and only if it can be represented as 
(2^p - 1) * 2^(p-1) for some prime p, where 2^p - 1 is also prime. 
No odd perfect numbers are known.

Let P = 2^(p-1) * (2^p - 1) be an even perfect number. 
We wish to find k such that sigma(k) + k = P.

The value of k may be very large, so we will search for solutions 
of the form k = q*r, where r is small (e.g. less than 10^8) and 
q is a prime not dividing r.

Let s = sigma(r). Then P = q*r + (q+1)*s = q*(r+s) + s, 
so q = (P - s) / (r + s).

Our search algorithm is as follows:

1. Precompute the values of sigma(r) for r = 1 to max_r.
2. For each Mersenne exponent p, compute the perfect number 
   P = 2^(p-1) * (2^p - 1).
3. For each r between 1 and max_r, check whether 
   P  - sigma(r) is divisible by r + sigma(r).
4. If so, compute q = (P - sigma(r)) / (r + sigma(r)).
5. If q is prime and q does not divide r, then we have found a solution. 
   Report the values of p, r, and sigma(r).

The value of k is not displayed, since it may be very large. 
But it can be computed easily from p, r, and sigma(r).
Specifically, k = q*r, where q = (P - sigma(r)) / (r + sigma(r)) 
and P = 2^(p-1) * (2^p - 1).

Thanks to Maximillian Hasler for suggestions and corrections to the original code.
"""

import numpy as np
from sympy import isprime
import datetime


# Mersenne exponents: Primes p such that 2^p - 1 is prime.
# The list of values below is complete (as of 2025-03-07) up to 57885161.
mersenne_exponents = [
    2, 3, 5, 7, 13, 17, 19, 31, 61, 89, 107, 127, 521, 607, 1279, 
    2203, 2281, 3217, 4253, 4423, 9689, 9941, 11213, 19937, 21701, 
    23209, 44497, 86243, 110503, 132049, 216091, 756839, 859433, 
    1257787, 1398269, 2976221, 3021377, 6972593, 13466917, 
    20996011, 24036583, 25964951, 30402457, 32582657, 37156667, 
    42643801, 43112609, 57885161, 74207281, 77232917, 82589933, 136279841]

max_r = 1_000_000  # Check all values of r less than max_r.
max_p = None # Check all exponents up to max_p. If None, check all exponents.

print(f'Run started {datetime.datetime.now()}')

def perfect_number(p, modulus=None):
    """Compute the perfect number 2^(p-1) * (2^p - 1) 
    associated with the Mersenne exponent p."""
    return  pow(2, 2*p - 1) - pow(2, p - 1)


def perfect_number_mod(p, modulus):
    """Compute the perfect number 2^(p-1) * (2^p - 1) 
    associated with the Mersenne exponent p, modulo modulus."""
    modulus = int(modulus)
    return (pow(2, 2*p - 1, modulus) - pow(2, p - 1, modulus)) % modulus

print(f'Precomputing sigma(r) for r < {max_r:_}')
sigma = np.ones(max_r, dtype=np.int64)
for r in range(2, max_r):
    sigma[r: max_r: r] += r


for i, p in enumerate(mersenne_exponents):
    if max_p is not None and p > max_p:
        break
    print(f'Testing {p=:_}')
    perfect = perfect_number(p)
    for r in range(2, max_r, 2):
        s = sigma[r]
        perfect_mod = perfect_number_mod(p, modulus=r+s)
        if perfect_mod == s:
            q = (perfect - s) // (r + s)
            if q % 2 and r % q and isprime(q):
                print(f'\tSolution Found: {p=:_}, {r=:_}, {s=:_}')


print(f'Run ended {datetime.datetime.now()}')