k + sigma(k) is perfect

k_plus_sigma_k_is_perfect.py

"""
    k_plus_sigma_k_is_perfect.py
    edited and modified by MFH 
    forked from D.Radcliffe

Problem:
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?

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

Theorem: (Euclid & Euler)
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.


## Idea of our approach:

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

We search for solutions of the form k = q*r, where q is a prime not dividing r.
We will scan "all" small values of r (e.g., all even r less than 10^8). 

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

Known examples:  k = q*r = 5*2 ; s = 3 ; r+s = 5
                 P = 28 = 2^2 * (2^3 - 1)  <=>  p = 3
                 P - s = 28 - 3 = 25 ; q = (P-s)/(r+s) = 5 
    (Probably) next larger:
        p = 31 <=> P =  2^30*(2^31-1) = 2305843008139952128
        r = 187268 = 2^2 * 46817, s = 327726 = 2 * 3^4 * 7 * 17^2 ;
        r+s = 514994 = 2 * 257497 ;  P-s = (r+s) * q ; q = 4477417228433

    No other with p < 90 000 and r < 10^7

Radcliff's 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).

"""

import numpy as np
from sympy import isprime, divisor_sigma as sigma
from datetime 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]

def perfect_number(p):
    "Return the perfect number 2^(p-1) * (2^p - 1) associated with the Mersenne exponent p."
    return ((1<<p) - 1) << (p-1)
    #P = 2 ** (p - 1)
    #return P * (2 * P - 1)

''' This seems memory vorace to me, I'll avoid,
    using the loop over r as outer loop and computing sigma on the fly.

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

max_r = 10_000_000  # Check all values of r less than max_r.
max_p = 2000 # 9941 # Check all exponents up to max_p.
             # Note that P = 2^(p-1) (2^p - 1) ~ 2^2p ~ 10 ^ 0.6p

def run(max_p = 9000, max_r = 1000_000):
    print(f'Run started {datetime.now()}')
    even_odd=[0,0]
    for r in range(2, max_r, 2): # check only even numbers
        if r % 10**5 == 0: print(end=f"[r={r/10**6}], ")
        s = int(sigma(r))
        min_P = 2*r + 3*s

        for p in mersenne_exponents:
            if p > max_p: break # go on with next r
            # print(f'Testing {p=:_}')
            P = perfect_number(p)
            if P >= min_P and (P - s) % (r + s) == 0: # this is the bottleneck
                # the above condition is very rarely satisfied:
                # for p < 5000 and r < 10^6, the following code is evaluated only 49 times
                q = (P - s) // (r + s)
                even_odd[q & 1] += 1    # in 27 cases it's even
                if not ( q & 1 or q == 2): continue
                # the test below was " q % r != 0 " which means q not a multiple of r
                # That's impossible when q is prime.
                # Also, q must be odd but r is even, so q % 1 (> 0) seems the appropriate test.
                # Of course, most q are  >> r, so r % q = 0 seems ridiculous
                # while q % r = 0 seems more probable...
                if r % q and isprime(q):
                    print(f"\tFound: {p=:_}, {r=:_}, {s=:_}")

    print(f'\nRun ended {datetime.now()}')
    print(f"{even_odd = :}" 
"""
(PARI/gp)
run(max_p = 1e6, max_r = 1e7, min_r = 2, check = 0, step=2)={ system("echo Start: %TIME%"); my( candidates = List(), two );
    forstep( r=min_r, max_r, step, /* step = 2 | min_r  <=>  check only even numbers */
        r % 10^5 || printf("[r=%.1fe6]",r/1e6); /* progress indicator */
        s = sigma(r) ; two = Mod(2,2*(r+s)) ;
        foreach ( mersenne_exponents, p, p > max_p && break ;
            two^(p-1)*(two^p-1)+r && next; /* continue only if  2(r+s) | P-s +-(r+s) : odd multiple of (r+s) */
            if ( check, /* then retain only candidates that don't have a small factor */
                my(f = factor(q = ((2^p-1)<<(p-1)) \ (r+s), 0)~); #f > 1 && next );
            listput(candidates, [r,p]); print1([r, p]", ");
            if ( check < 2, next );
            ispseudoprime ( ((2^p-1)<<(p-1)) \ (r+s) ) || next;
            printf("\n\tFound: p = %d, r = %d, s = %d\n",p,r,s);
    )   ); /*foreach ; forstep*/
    printf("\nFound %d cases where r+s | P-s.\n", #candidates); system("echo End: %TIME%"); candidates}
    
 Yields:
 candidates = [[2, 3], [2, 7], [2, 19], [2, 31], [2, 107], [2, 127], [2, 607], [2, 1279], [2, 2203], [2, 4423], [2, 86243], [2, 110503], [2, 216091], [2, 756839],
 [8, 17], [8, 31], [8, 61], [8, 127], [8, 9689], [8, 21701], [8, 756839], [12, 3], [22, 17], [22, 521], [46, 21701], [46, 23209], 
 [70, 86243], [144, 17], [156, 21701], [160, 132049], [254, 17], [320, 13], [328, 2203], [358, 107], 
 [1136, 132049], [1462, 1279], [1748, 61], [1786, 107], [2728, 13], [2728, 4253], [2876, 9689], 
 [3022, 4253], [3024, 132049], [4064, 19937], [5014, 61], [5728, 23209], [5952, 13], [6614, 1279], [8944, 2281], [9934, 44497], 
 [12250, 21701], [15132, 19], [15736, 17], [17270, 859433], [20728, 31], [21952, 19], [37088, 86243], [38748, 9941], [45012, 2203], 
 [50264, 756839], [57692, 19], [58144, 9941], [61336, 11213], [71258, 86243], [83416, 21701], [96844, 44497],
 /*[r=0.1]*/
 [103168, 19937][107456, 9941][153824, 127][176096, 859433][187268, 31][188080, 3217][191822, 44497][r=0.2][211392, 756839][r=0.3][314236, 607][314492, 607][327504, 859433][344786, 31][397544, 23209][r=0.4][426478, 1279][444544, 127][r=0.5][581500, 17][585488, 859433][594032, 19][r=0.6][692008, 127][697516, 107][r=0.7][r=0.8][853024, 44497][r=0.9][r=1.0][1056742, 23209][1072402, 756839][1096412, 607][r=1.1][1129760, 9941][1171308, 132049][1175462, 107][r=1.2][1275356, 4253][r=1.3][r=1.4][r=1.5][r=1.6][r=1.7][r=1.8][r=1.9][r=2.0][r=2.1][r=2.2][2240560, 23209][2270968, 607][r=2.3][r=2.4][r=2.5][r=2.6][r=2.7][r=2.8][r=2.9][r=3.0][r=3.1][r=3.2][r=3.3][r=3.4][r=3.5][r=3.6][r=3.7][r=3.8][r=3.9][r=4.0][r=4.1][r=4.2][r
 [103168, 19937], [107456, 9941], [153824, 127], [176096, 859433], [187268, 31], [188080, 3217], [191822, 44497], 
 [211392, 756839], [314236, 607], [314492, 607], [327504, 859433], [344786, 31], [397544, 23209], [426478, 1279], [444544, 127], 
 [581500, 17], [585488, 859433], [594032, 19], [692008, 127], [697516, 107], [853024, 44497], 
 [1056742, 23209], [1072402, 756839], [1096412, 607], [1129760, 9941], [1171308, 132049], [1175462, 107], [1275356, 4253], 
 [2240560, 23209], [2270968, 607], [4422472, 1279], [5392666, 31], [5576362, 23209], [6340640, 216091], [6579116, 9689], 
 [8306216, 86243], [8930512, 107], [9400130, 1279], [9566544, 61]]

M=Map(); foreach(%, c, mapput(M,c[2],concat(iferr(mapget(M,c[2]),E,[]),c[1]))); M
%223 = Map([
      3, [2, 12]; 7, [2]; 13, [320, 2728, 5952]; 17, [8, 22, 144, 254, 15736, 581500]; 19, [2, 15132, 21952, 57692, 594032];
     31, [2, 8, 20728, 187268, 344786, 5392666]; 61, [8, 1748, 5014, 9566544]; 107, [2, 358, 1786, 697516, 1175462, 8930512];
    127, [2, 8, 153824, 444544, 692008]; 521, [22]; 607, [2, 314236, 314492, 1096412, 2270968]; 
   1279, [2, 1462, 6614, 426478, 4422472, 9400130]; 2203, [2, 328, 45012]; 2281, [8944]; 3217, [188080]; 
   4253, [2728, 3022, 1275356]; 4423, [2]; 9689, [8, 2876, 6579116]; 9941, [38748, 58144, 107456, 1129760]; 
  11213, [61336]; 19937, [4064, 103168]; 21701, [8, 46, 156, 12250, 83416]; 23209, [46, 5728, 397544, 1056742, 2240560, 5576362]; 
  44497, [9934, 96844, 191822, 853024]; 86243, [2, 70, 37088, 71258, 8306216];  110503, [2]; 132049, [160, 1136, 3024, 1171308]; 
  216091, [2, 6340640]; 756839, [2, 8, 50264, 211392, 1072402]; 859433, [17270, 176096, 327504, 585488]])
  
Selected candidates (where q = 2 or odd):
[P(p) = (2^p-1)<<(p-1), rps(r)=r+sigma(r)]
CANDID = [ rp | rp <- candidates, (q = P(rp[2]) \ rps(rp[1]))%2 || q==2 ]
%247 = [[2, 3], [2, 7], [2, 19], [2, 31], [2, 107], [2, 127], [2, 607], [2, 1279], [2, 2203], [2, 4423], [2, 86243], 
    [2, 110503], [2, 216091], [2, 756839], [8, 17], [8, 31], [8, 61], [8, 127], [8, 9689], [8, 21701], [8, 756839], 
    [144, 17], [160, 132049], [320, 13], [328, 2203], 
    [1136, 132049], [5728, 23209], [8944, 2281], [21952, 19], [37088, 86243], [83416, 21701], 
    [103168, 19937], [107456, 9941], [187268, 31], [211392, 756839], [444544, 127], [692008, 127], [4422472, 1279]]

#CANDID \\ = 38 !
for(i=1,#CANDID,print1(i", ");rp=CANDID[i];ispseudoprime(q=P(rp[2])\rps(rp[1]))&& print("/* YES! */ "))
1, /* YES! */
2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 
30,  *** ispseudoprime: user interrupt after 4,751 ms
31,  *** ispseudoprime: user interrupt after 4,767 ms
32, 33, 34, /* YES! */
35,  *** ispseudoprime: user interrupt after 8,970 ms
36, 37, 38,

vecextract(CANDID, [30, 31, 35])
%246 = [[37088, 86243], [83416, 21701], [211392, 756839]]

foreach(%246,rp, printf("%s: len(q)=%d, prime: ",rp,#Str(q=P(rp[2])\rps(rp[1]))); print(if(ispseudoprime(q),YES,NO)))
[37088, 86243]: len(q)=51919, prime: NO (small factors :  1091 · 7487)
[83416, 21701]: len(q)=13060, prime: NO (small factor : 3373)
[211392, 756839]: len(q)=455657, prime: NO (small factor : 143287)

rp=CANDID[30]; if(1 < #F=factor(q=P(rp[2])\rps(rp[1]),0)~, F[1..-2], "no small factor")

"""

PARI/GP changed to possibly produce only list of candidates (r, p), i.e., pairs for which r+s | P-s.

added more investigation in candidates