Perfect Numbers Introduction
Time for some cool mathematics! I recently was exploring the work and life of [Paul Erdős](https://nl.wikipedia.org/wiki/Paul_Erd%C5%91s), who was one of the most prolific, productive and eccentric mathematicians ever. He seemed to have a special fondness for [perfect numbers](https://en.wikipedia.org/wiki/Perfect_number). This sentiment was shared by Pythagoras, who believed that perfect numbers have a mystical significance. In this TIL I'll explain what perfect numbers are and how to generate them.
What are perfect numbers?
A perfect number is a positive integer $n \in \mathbb{N}$ (1, 2, 3, ...) that is equal to the sum of its [proper divisors](https://en.wikipedia.org/wiki/Divisor). Proper divisor just means that we can divide by that number without having a remainder and exclude the number itself. For example, $6$ is a perfect number because its divisors are $1$, $2$, $3$, and $1 + 2 + 3 = 6$. $28$ is a perfect number because its proper divisors are $1$, $2$, $4$, $7$, $14$, and $1 + 2 + 4 + 7 + 14 = 28$. And so on. These numbers become large pretty quickly. The next perfect numbers are $496$, $8128$, and $33550336$. One fascinating fact about perfect numbers is that no odd perfect number has ever been found and it is unknown if any exist.
Generating perfect numbers
Let's write a simple Python function to generate perfect numbers. One of the most efficient ways to generate perfect numbers is by using the [Euclid-Euler theorem](https://en.wikipedia.org/wiki/Euclid%E2%80%93Euler_theorem). This states that if a number $n$ is a [Mersenne prime](https://en.wikipedia.org/wiki/Mersenne_prime):
$$n = 2^p - 1$$
Then we can find the corresponding perfect number $P$ using $p$ by computing:
$$P = 2^{p-1} \times (2^p - 1)$$
To find perfect numbers we therefore have to loop over $p > 1$ and can compute a perfect number from $p$ each time $2^p -1$ is a prime number. Evaluating if a number is prime can be done efficiently using [sympy's isprime function](https://docs.sympy.org/latest/modules/ntheory.html#sympy.ntheory.primetest.isprime). Here is how we can discover perfect numbers with Python:
from sympy import isprime
def perfect_numbers(n: int) -> list[int]:
'''
Get the first n perfect numbers using the Euclid-Euler theorem.
:param n: How many perfect numbers to return.
:return: List of perfect numbers.
'''
# Start from first non-trivial Mersenne prime, namely 2^2-1=3
p = 2
perfects = []
# Loop over all p until we have n perfect numbers
while len(perfects) < n:
# Compute the Mersenne prime
mersenne = 2**p - 1
# Check if we have a valid p to compute a perfect number for
if isprime(mersenne):
# Compute the perfect number
euclid_euler = 2**(p-1) * mersenne
perfects.append(euclid_euler)
p += 1
return perfects
# Get the first 15 perfect numbers
print(perfect_numbers(15))
# [6, 28, 496, 8128, ...]
Thats it! Feel free to try out the function and see how many perfect numbers you can generate. The 15th perfect number has 770 digits and should look something like this:
5416252628436584741265446537439131614085649053903169578460392081838720699415853485919899992105671992191
9057390080263646159280013827605439746262788903057303445505827028395139475207769044924431494861729435113
1262808379049304627406817179604658673487209925721905694655452996299198234310310926242444635477896354414
8139171981644160558678809214788667732139875666162471455172696430221755428178425481731961195165985555357
3937788923405146222324506715979193757372820860878214322052227584537552897476256179395176624426314480313
4469350852036575847982475360211728804037830486028736212593137899949003366739415037472249669840282408060
4210869007767039525923189466627361521277560353576470795225017385830517102860302123489664785136394992890
4973292145107505979911456221519899345764984291328Thank you for reading and exploring the fascinating world of perfect numbers with me! Until next time!
