「約数 (divisor)」は整数 n を割り切る整数のことです。たとえば、24 の正の約数は 1, 2, 3, 4, 6, 8, 12, 24 の 8 個あります。ここでは正の約数を考えることにします。
それでは問題です。
最大公約数は「ユークリッドの互除法」で求めることができます。
詳しい説明は拙作のページ Python 入門 第 3 回 再帰定義と高階関数 : ユークリッドの互除法 をお読みください。
プログラムは次のようになります。
リスト : 最大公約数 (末尾再帰)
def gcd(a, b):
if b == 0: return a
return gcd(b, a % b)
関数 gcd() は引数 a と b の最大公約数を求めます。b が 0 の場合は a を返します。これが再帰呼び出しの停止条件になります。そうでなければ、gcd() を再帰呼び出しして、b と a % b の最大公約数を求めます。
残念ながら、Python は末尾再帰最適化をサポートしていません。繰り返しに変換すると次のようになります。
リスト : 最大公約数 (繰り返し)
def gcd(a, b):
while b > 0:
a, b = b, a % b
return a
引数 a, b の値を書き換えることで最大公約数を求めています。
簡単な実行例を示します。
gcd(12345678, 123456789) => 9 gcd(1234321, 12345654321) =>121
複数の整数の最大公約数を Python で求める場合は reduce() を使うと簡単です。
reduce(gcd, [123, 12345, 12345678]) => 3
最小公倍数は最大公約数を使って簡単に求めることができます。プログラムは次のようになります。
リスト : 最大公倍数
def lcm(a, b):
return a * b / gcd(a, b)
整数 a と b の最小公倍数は a * b / gcd(a, b) で求めることができます。
簡単な実行例を示します。
lcm(5, 7) => 35 lcm(14, 35) => 70
複数の整数の最小公倍数を Python で求める場合は reduce() を使うと簡単です。
reduce(lcm, range(2, 21)) => 232792560 reduce(lcm, range(2, 31)) => 2329089562800
n の素因数分解ができると、約数の個数を求めるのは簡単です。n = pa * qb * rc とすると、約数の個数は (a + 1) * (b + 1) * (c + 1) になります。たとえば、12 は 22 * 31 になるので、約数の個数は 3 * 2 = 6 になります。実際、12 の約数は 1, 2, 3, 4, 6, 12 の 6 個です。
拙作のページ 素数 で作成した関数 factorization() を使うと、プログラムは次のようになります。
リスト : 約数の個数
def divisor_num(n):
a = 1
for _, x in factorization(n):
a *= x + 1
return a
関数 divisor_num() は for ループでリストの要素 (タプル) を順番に取り出し、x + 1 を a に掛け算していくだけです。
簡単な実行例を示します。
divisor_num(24) => 8 divisor_num(12345678) => 24 divisor_num(123456789) => 12 divisor_num(1234567890) => 48 divisor_num(1111111111) => 16
n の素因数分解ができると、約数の合計値を求めるのは簡単です。n の素因数分解が pa だった場合、その約数の合計値は次の式で求めることができます。
σ(p, a) = pa + pa-1 + ... + p2 + p + 1
たとえば、8 の素因数分解は 23 になり、素数の合計値は 8 + 4 + 2 + 1 = 15 になります。
pa の約数の合計値を σ(p, a) で表すことにします。n = pa * qb * rc の場合、n の約数の合計値は σ(p, a) * σ(q, b) * σ(r, c) になります。たとえば、12 は 22 * 3 に素因数分解できますが、その合計値は (4 + 2 + 1) * (3 + 1) = 28 となります。12 の約数は 1, 2, 3, 4, 6, 12 なので、その合計値は確かに 28 になります。
拙作のページ 素数 で作成した関数 factorization() を使うと、プログラムは次のようになります。
リスト : 約数の合計値
# σ(p, n) の計算
def div_sum_sub(p, n):
a = 0
while n > 0:
a += p ** n
n -= 1
return a + 1
def divisor_sum(n):
a = 1
for p, q in factorization(n):
a *= div_sum_sub(p, q)
return a
関数 div_sum_sub() は σ(p, n) を計算します。あとは for ループで div_sum_sub() の返り値を累積変数 a に掛け算していくだけです。
簡単な実行例を示します。
divisor_sum(24) => 60 divisor_sum(12345678) => 27319968 divisor_sum(123456789) => 178422816 divisor_sum(1234567890) => 3211610688 divisor_sum(1111111111) => 1246404096
p が素数の場合、pa の約数は次のように簡単に求めることができます。
pa, pa-1, ... p2, p, 1
n の素因数分解が pa * qb だったとすると、その約数は次のようになります。
(pa, pa-1, ... p2, p, 1) * qb,
(pa, pa-1, ... p2, p, 1) * qb-1,
.....
(pa, pa-1, ... p2, p, 1) * q2,
(pa, pa-1, ... p2, p, 1) * q,
(pa, pa-1, ... p2, p, 1) * 1
たとえば、12 の約数は 24 = (1, 2, 4) と 3 = (1, 3) から、(1, 2, 4) * 1 と (1, 2, 4) * 3 のすべての要素 (1, 2, 4, 3, 6, 12) になります。
プログラムは次のようになります。
リスト : 約数をすべて求める
# pq の約数を求める
def divisor_sub(p, q):
a = []
for i in xrange(0, q + 1):
a.append(p ** i)
return a
def divisor(n):
xs = factorization(n)
ys = divisor_sub(xs[0][0], xs[0][1])
for p, q in xs[1:]:
ys = [x * y for x in divisor_sub(p, q) for y in ys]
return sorted(ys)
関数 divisor_sub() は pq の約数をリストに格納して返します。引数 n を factorization() で素因数分解して変数 xs にセットします。xs の先頭要素を divisor_sub() に渡してリストに変換して変数 ys にセットします。あとは for ループで xs の 1 番目から要素を順番に取り出し、pq を divisor_sub() でリストに変換して、それを内包表記で累積変数 ys のリストの要素と掛け合わせていくだけです。
簡単な実行例を示します。
divisor(24) => [1, 2, 3, 4, 6, 8, 12, 24] divisor(12345678) => [1, 2, 3, 6, 9, 18, 47, 94, 141, 282, 423, 846, 14593, 29186, 43779, 87558, 131337, 262674, 685871, 1371742, 2057613, 4115226, 6172839, 12345678] divisor(123456789) => [1, 3, 9, 3607, 3803, 10821, 11409, 32463, 34227, 13717421, 41152263, 123456789] divisor(1234567890) => [1, 2, 3, ,5 ,6 9, 10, 15, 18, 30, 45, 90, 3607, 3803, 7214, 7606, 10821, 11409, 18035, 19015, 21642, 22818, 32463, 34227, 36070, 38030, 54105, 57045, 64926, 68454, 108210, 114090, 162315, 171135, 324630, 342270, 13717421, 27434842, 41152263, 68587105, 82304526, 123456789, 137174210, 205761315, 246913578, 411522630, 617283945, 1234567890] divisor(1111111111) => [1, 11, 41, 271, 451, 2981, 9091, 11111, 100001, 122221, 372731, 2463661, 4100041, 27100271, 101010101, 1111111111]
リスト : 完全数
def perfect_number(n):
for x in xrange(2, n + 1):
if divisor_sum(x) - x == x:
print x,
print
perfect_number(10000)
完全数を求める perfect_number() は簡単です。x の約数の合計値を divisor_sum() で求め、その値から x を引いた値が x と等しければ完全数です。print で x を表示します。
実行結果を示します。
6 28 496 8128
ところで、参考 URL 1 によると、メルセンヌ素数を Mn とすると、偶数の完全数は 2n-1 * Mn で表すことができるそうです。この式を使うと偶数の完全数は次のようになります。
n : メルセンヌ素数 : 完全数 ---+----------------+---------------------- 2 : 3 : 6 3 : 7 : 28 5 : 31 : 496 7 : 127 : 8128 13 : 8191 : 33550336 17 : 131071 : 8589869056 19 : 524287 : 137438691328 31 : 2147483647 : 2305843008139952128
なお、奇数の完全数はまだ発見されておらず、偶数の完全数 (つまりメルセンヌ素数) が無数に存在するか否かも未解決な問題だそうです。
ところで、参考 ULR 3, 4 によると、『その数自身を除く約数の総和が元の数より大きい数』を「過剰数 (abundant number)」といい、『その数自身を除く約数の総和が元の数より小さい数』を「不足数 (deficient number)」というそうです。過剰数と不足数の個数を求めるプログラムも簡単に作ることができます。
リスト : 過剰数と不足数
# 過剰数
def abundant_number(n):
count = 0
for x in xrange(1, n + 1):
y = divisor_sum(x) - x
if y > x: count += 1
return count
# 不足数
def deficient_number(n):
count = 0
for x in xrange(1, n + 1):
y = divisor_sum(x) - x
if y < x: count += 1
return count
print abundant_number(1000000)
print deficient_number(1000000)
247545 752451
1000000 以下の過剰数は 247545 個、不足数は 752451 個、完全数は 4 個なので、合計で 1000000 になります。参考 URL 3 によると、『自然数のうち過剰数が占める割合は0.2474から0.2480の間であると証明されている。』 とのことで、1000000 以下の過剰数の個数は確かにこの範囲内に入っています。
リスト : 友愛数
def yuuai_number(n):
for x in xrange(2, n + 1):
m = divisor_sum(x) - x
if m < x and x == divisor_sum(m) - m:
print m, x
print
友愛数を求める関数 yuuai_number() も簡単です。divisor_sum() で x の約数の合計値を求め、その値から x を引いた値を変数 m にセットします。m の約数の合計値から m を引いた値が x と等しければ、x と m は友愛数です。print で x と m を表示します。同じ組を表示しないようにするため、m < x を条件に入れています。
実行結果を示します。
220 284 1184 1210 2620 2924 5020 5564 6232 6368 10744 10856 12285 14595 17296 18416 63020 76084 66928 66992 67095 71145 69615 87633 79750 88730
なお、友愛数が無数に存在するか否かは、未解決な問題だそうです。
perfect_number() と yuuai_number() は、divisor_sum() を呼び出して約数の和を求めていますが、あらかじめ約数の和を計算してリスト (配列) に格納しておく方法もあります。この場合、約数の和の計算にちょっと時間がかかりますが、完全数と友愛数を求める処理は高速になります。
基本的な考え方は簡単です。約数の和を格納する配列 sum_table を用意します。sum_table の要素は 1 に初期化します。2 から順番に素数 p で割り算して 1 + p1 + ... + pq を求め、それを sum_table の値に掛け算します。素数は「エラトステネスの篩」と同じ方法で求めることができます。素数の倍数であれば、sum_table の値は 1 よりも大きくなっているはずです。つまり、sum_table が 1 であれば、その整数は素数であることがわかります。
プログラムは次のようになります。
リスト : n 以下の整数の約数の和をリストに格納して返す
def make_divisor_sum(n):
def factor_sub(n, p):
a = 1
q = 1
while n % p == 0:
a += p ** q
q += 1
n /= p
return a
#
sum_table = [1] * (n + 1)
for i in xrange(2, n + 1):
if sum_table[i] != 1: continue
for j in xrange(i, n + 1, i):
sum_table[j] *= factor_sub(j, i)
return sum_table
局所関数 factor_sub() は n を 素数 p で割り算して、1 + p1 + ... + pq を求めます。make_divisor_sum() は sum_table を初期化してから、for ループで素数を探します。sum_table[i] が 1 でない場合、i は素数ではありません。contiune で次の数をチェックします。素数の場合は、その倍数に対応する sum_table の値を更新します。i の倍数を 変数 j にセットし、sum_table[j] に factor_sub(j, i) の値を乗算するだけです。最後に sum_table を返します。
make_divisor_sum() を使うと yuuai_number() は次のようになります。
リスト : 友愛数
def yuuai_number(n):
table = make_divisor_sum(n)
for x in xrange(2, n + 1):
y = table[x] - x
if y < x and table[y] - y == x:
print y, x
実際に 1,000,000 以下の友愛数を求めてみたところ、40 個の友愛数を出力するのに最初のプログラムでは 2.9 秒 (実行環境 : Windows 7, Core i7-2670QM 2.20GHz, PyPy 4.0.1) かかりましたが、make_divisor_sum() を使ったプログラムは約数の和を求める処理を含めて 0.28 秒ですみました。約 10 倍高速化することができました。