Odd Perfect Numbers

Introduction

Euler proved that an odd perfect number \(N\) must be of the form \(N = q^{\alpha}p_1^{2e_1}\ldotp \ldotp \ldotp p_k^{2e_k}\), where \(q\), \(p_1\), ..., \(p_k\) are distinct odd primes and \(q\)≡\(α\)≡1(mod 4). Very comprehensive searches have not found any numbers which satisfy these conditions: they either do not exist or are very rare. We have searched instead for spoof odd perfect numbers which satisfy a weaker condition:
\(N = qp_1^2\ldotp \ldotp \ldotp p_k^2\), where \(p_1\), ..., \(p_k\) are distinct odd primes.

For a spoof odd perfect number \(N\), \(\sigma (N) = 2N\), where \(\sigma\) is the sum-of-divisors function. Let \(a = \prod_{i = 1}^k\sigma (p_i^2)\) and \(b = \prod_{i = 1}^kp_i^2\). Then \(a(q + 1) = 2bq\) and \(q = a/ (2b - a)\). We can check each of the \(C(n,k)\) combinations of k primes from a list of the first \(n\) primes, checking whether \(a \equiv 0 mod(2b - a)\). With k=4 and n >= 5 this finds \(q=22021\) with primes [3, 7, 11, 13], which is the famous Descartes number. No other such numbers are currently known.

Search Combinations of Primes

We can relax the conditions even further, and check whether the remainder \(r\) of \(a/ (2b - a)\) is close to 0 or \(2b - a -1\). Let \(d = 2b - a\), then find the combinations where \(r / d\) have minimum and maximum values. Increasing the number of odd primes searched can find better results, at the cost of increasing the number of combinations and the execution time. In the results below, the number of odd primes has been set as large as possible to give no more than 1,000, 000, 000 combinations.

Choose smallest primes which keep \(\sigma (N) \leq 2N\)

There is another approach to finding these numbers which does not depend on searching combinations of a rather arbitrary set of primes. Start with a list of one or more primes and a value of \(q\) for which the fraction \(\varphi = \frac{a(q - 1)}{2bq} < 1\). Values of \(\frac{\sigma (p^2)}{p^2}\) are always greater than 1, getting closer to 1 as \(p\) increases. So there is a prime \(p\) for which \(\varphi < \varphi\frac{\sigma (p^2)}{p^2} \leq 1\). This can be found by solving \[\frac{\varphi(p^3 - 1)}{(p - 1)p^2} \leq 1,\] \[p^2(\frac{1}{\varphi} - 1) - p - 1 \geq 0.\] Taking the positive root of this quadratic equation gives \[p = \frac{1 + \sqrt{1 + 4(\frac{1}{\varphi} - 1)}}{2(\frac{1}{\varphi} - 1)}.\] It is convenient to calculate in Python with \(\varphi\) as a gmpy2 rational type gmpy2.mpq(). This precludes the use of a square root, however, so we need to make a substitution in this formula. if \(\varepsilon = \frac{1}{\varphi } - 1\), \[\left(1 + 2\varepsilon \right)^2 = 1 + 4\varepsilon + 4\varepsilon^2\] \[1 + 2\varepsilon > \sqrt{1 + 4\varepsilon }\] \[p < \frac{1 + 1 + 2\varepsilon }{2\varepsilon }\] \[p < \frac{1}{1 - \varphi }\] As \(p\) will almost certainly not be an integer, we can set \(p\) to the next prime after \(\left\lfloor p \right\rfloor\). (We also need to try setting \(p\) to the next prime after \(\left\lfloor p-1 \right\rfloor\) in case the substitution has caused \(p\) to be slightly overestimated.) The fraction \(\varphi\) can repeatedly be increased towards 1 by this procedure. It may never ever reach 1, but in practice the value of \(p\) soon becomes so huge that finding the next prime (with gmpy2.next_prime(p), which actually find the next probable prime) takes a very long time. Starting with primes [5, 13] and q = 3, this procedure added primes as follows:-
[5, 13, 11, 67, 853, 7330607, 2519305090553, 14426102101776053284860379, 1927211959374477020846289210499909161066541397893, 41187444261697319336775171368660819898740194508980459349407138663903744584302715436780243183729]
Subsequent primes had 188, 372, 742, 1479, and 2954 digits. The value of \(1 - \varphi\) was then less than \(3 \times 10^{-5905}\), so \(\varphi\) was extremely close to 1. Here are the Python functions which perform these calculations:-

def σ(p, e):
    return (p ** (e + 1) - 1) // (p - 1)

def odd_perfect_number(primes, q):
    import gmpy2
    φ = gmpy2.mpq(q + 1, 2 * q)
    for p in primes: φ = φ * σ(p, 2) / p ** 2
    assert φ < 1
    while True:
        p = gmpy2.mpz(1 / (1 - φ))
        for p in [p - 1, p]:
            p = gmpy2.next_prime(p)
            if p != 2 and φ * σ(p, 2) / p ** 2 < 1:
                break
        else:
            print('φ >= 1')
            break
        digits = p.num_digits()
        print('next prime', p if digits <= 100 else '' , digits, 'digits')
        φ = φ * σ(p, 2) / p ** 2
        primes.append(p)
    return primes

Use primes which give cancellation

The fraction \(\varphi\) will never reach 1 unless all the prime factors in the numerator and denominator cancel out, which won't happen with the above approach. Here is better approach, where any prime factor \(p\) found in the numerator is immediately applied to \(\varphi\) by \(\varphi = \varphi \frac{\sigma (p^2)}{p^2}\), causing some cancellation. Initially, \(\varphi = \frac{(q + 1)}{2q}\). This often results in \(\varphi > 1\), which means that the value \(q\) cannot result in a perfect number or spoof, as \(\varphi\) can only get bigger. There is a complication if a value of \(p\) has already been used, as \(\sigma\) needs to be calculated differently for exponents other than 2. So we keep track of all primes and their exponents previously used, and set \(\varphi = \varphi \frac{\sigma (p^{e + 2})}{\sigma (p^e)p^2}\), where \(e\) is the exponent previously used. This approach can quite rapidly check a range of odd values of q increasing by 2. In the following table, all odd values of \(q\) between those in the rows have been checked and rejected. It is encouraging that the table includes Euler's even perfect numbers, spoof even perfect numbers when \(q\) is not prime, and Descartes number, when these have not been explicitly coded for.

Christopher B. Jones 2025-02-13