ここでは、d(n)をn自身を除いたnの約数の総和とします。
完全数
d(n) = n
を満たすとき、nを完全数（Perfect number）と呼びます。
p = 2m - 1
で、pが素数のとき、pをメルセンヌ素数と呼びます。メルセンヌ素数には特別な高速の素数判定法があり、現在見つかっている大きな素数はすべてメルセンヌ素数です。
さて、p = 2m - 1がメルセンヌ素数のとき、
n = 2m-1(2m - 1)
とすると、nは完全数です。なぜなら、約数のところで紹介した計算法より、
d(n) = (2m - 1)(1 + 2m - 1) - n = n
を確認できます。
逆に、偶数の完全数はこの形であることが証明されています。
奇数の完全数についてはあまりよくわかっていません。
社交数
n2 = d(n1)
n3 = d(n2)
…
n1 = d(nm)
m > 2
となるとき、n1, ..., nmを社交数（Social number）と言います。
完全数・友愛数とあわせて社交数を求めるプログラムを示しましょう。エラトステネスのふるい的に素因数分解すればよいのですが、素因数分解してから約数の和を求めるのではなく、直接約数の和を求めます。すなわち、最初に配列をすべて1にしておいて、例えば、32で割り切れるなら、配列の要素に1 + 3 + 32 = 1 + 3 * (1 + 3) = 13を掛けるといった処理を施していけばよいです。
ここでは全て1000万以内の社交数を求めています。
[12496, 14288, 15472, 14536, 14264]
[14316, 19116, 31704, 47616, 83328, 177792, 295488, 629072, 589786, 294896, 358336, 418904, 366556, 274924, 275444, 243760, 376736, 381028, 285778, 152990, 122410, 97946, 48976, 45946, 22976, 22744, 19916, 17716]
[1264460, 1547860, 1727636, 1305184]
[2115324, 3317740, 3649556, 2797612]
[2784580, 3265940, 3707572, 3370604]
[4938136, 5753864, 5504056, 5423384]
[7169104, 7538660, 8292568, 7520432]
def make_primes(max_p):
def is_prime(n):
for p in primes:
if p * p > n:
return True
elif n % p == 0:
return True
return True
primes = [ 2 ]
for n in range(3, max_p + 1, 2):
if is_prime(n):
primes.append(n)
return primesdef calc_p_factor(n, p):
m = 1
while n % p == 0:
m = 1 + m * p
n /= p
return m, ndef calc_sum_divs(max_n):
a = range(max_n + 1)
d = [ 1 ] * (max_n + 1)
for p in primes:
for k in xrange(p, max_n + 1, p):
m, a[k] = calc_p_factor(a[k], p)
d[k] *= m
for k in xrange(2, max_n + 1):
if a[k] != 1:
d[k] *= (a[k] + 1)
d[k] -= k
return ddef print_social(n):
a = [ n ]
m = d[n]
while m != n:
a.append(m)
m = d[m]
print adef find_social_number(n):
s = set()
m = n
while m not in s:
s.add(m)
m = d[m]
if m > N:
return 0
if m == n:
a = list(s)
a.sort()
if a[0] == n:
print_social(n)N = 10 ** 7
M = int(N ** 0.5 + 0.5)
primes = make_primes(N)
d = calc_sum_divs(N)
for n in xrange(2, N + 1):
find_social_number(n)
