下面是一个(更有效的)筛子版本：# search all numbers in [2..limit] for perfect numbers

# (ones whose proper divisors sum to the number)

limit = int(input("enter upper limit for perfect number search: "))

# initialize - all entries are multiples of 1

# (ignore sieve[0] and sieve[1])

sieve = [1] * (limit + 1)

n = 2

while n <= limit:

# check n

if sieve[n] == n:

print(n, "is a perfect number")

# add n to all k * n where k > 1

kn = 2 * n

while kn <= limit:

sieve[kn] += n

kn += n

n += 1

找到10000个6 is a perfect number

28 is a perfect number

496 is a perfect number

8128 is a perfect number

将这些因素分解可以得到一个有趣的模式：6 3 * 2 ( 4 - 1) * ( 4 / 2)

28 7 * 2 * 2 ( 8 - 1) * ( 8 / 2)

496 31 * 2 * 2 * 2 * 2 ( 32 - 1) * ( 32 / 2)

8128 127 * 2 * 2 * 2 * 2 * 2 * 2 (128 - 1) * (128 / 2)

其中第一个因子(3，7，31，127)是一个素数，它比二的幂小一倍，并且乘以二的相同幂的一半。此外，所涉及的幂是素数(2**2，2**3，2**5，2**7)。

事实上，Euclid证明了当2**p - 1是素数时(2**p - 1) * 2**(p - 1)是一个完美数，而只有当p是素数时才可能(尽管不能保证)。欧拉进一步证明了所有的偶数都必须是这种形式。

这意味着一个更高效的版本-我将继续使用循环，可以随意重写它。首先，我们需要一个素数源和一个is-u素数测试：def primes(known_primes=[7, 11, 13, 17, 19, 23, 29]):

"""

Generate every prime number in ascending order

"""

# 2, 3, 5 wheel

yield from (2, 3, 5)

yield from known_primes

# The first time the generator runs, known_primes

# contains all primes such that 5 < p < 2 * 3 * 5

# After each wheel cycle the list of known primes

# will be added to.

# We need to figure out where to continue from,

# which is the next multiple of 30 higher than

# the last known_prime:

base = 30 * (known_primes[-1] // 30 + 1)

new_primes = []

while True:

# offs is chosen so 30*i + offs cannot be a multiple of 2, 3, or 5

for offs in (1, 7, 11, 13, 17, 19, 23, 29):

k = base + offs # next prime candidate

for p in known_primes:

if not k % p:

# found a factor - not prime

break

elif p*p > k:

# no smaller prime factors - found a new prime

new_primes.append(k)

break

if new_primes:

yield from new_primes

known_primes.extend(new_primes)

new_primes = []

base += 30

def is_prime(n):

for p in primes():

if not n % p:

# found a factor - not prime

return False

elif p * p > n:

# no factors found - is prime

return True

然后搜索看起来像# search all numbers in [2..limit] for perfect numbers

# (ones whose proper divisors sum to the number)

limit = int(input("enter upper limit for perfect number search: "))

for p in primes():

pp = 2**p

perfect = (pp - 1) * (pp // 2)

if perfect > limit:

break

elif is_prime(pp - 1):

print(perfect, "is a perfect number")

发现enter upper limit for perfect number search: 2500000000000000000

6 is a perfect number

28 is a perfect number

496 is a perfect number

8128 is a perfect number

33550336 is a perfect number

8589869056 is a perfect number

137438691328 is a perfect number

2305843008139952128 is a perfect number

不到一秒钟；-)