約数 (divisor) は整数 n を割り切る整数のことです。たとえば、24 の正の約数は 1, 2, 3, 4, 6, 8, 12, 24 の 8 個あります。ここでは正の約数を考えることにします。
それでは問題です。
本ページでは Python3 (ver 3.10.12) を使ってプログラムを作ることにします。
最大公約数は「ユークリッドの互除法」で求めることができます。
詳しい説明は拙作のページ 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 gcd1(a, b):
while b > 0:
a, b = b, a % b
return a
引数 a, b の値を書き換えることで最大公約数を求めています。
簡単な実行例を示します。
>>> from divisor import * >>> gcd(12345678, 123456789) 9 >>> gcd1(12345678, 123456789) 9 >>> gcd(1234321, 12345654321) 121 >>> gcd1(1234321, 12345654321) 121
複数の整数の最大公約数を Python で求める場合は標準ライブラリ functools の関数 reduce() を使うと簡単です。
>>> import functools >>> functools.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 で求める場合は functools.reduce() を使うと簡単です。
>>> functools.reduce(lcm, range(2, 21)) 232792560 >>> functools.reduce(lcm, range(2, 31)) 2329089562800
n の素因数分解ができると、約数の個数を求めるのは簡単です。\(n = p^a \times q^b \times r^c\) とすると、約数の個数は \((a + 1) \times (b + 1) \times (c + 1)\) になります。たとえば、12 は \(2^2 * 3^1\) になるので、約数の個数は 3 * 2 = 6 になります。実際、12 の約数は 1, 2, 3, 4, 6, 12 の 6 個です。
拙作のページ 素数 で作成した関数 factorization() を使うと、プログラムは次のようになります。
リスト : 約数の個数
def divisor_count(n):
a = 1
for _, x in factorization(n):
a *= x + 1
return a
関数 divisor_count() は for ループでリストの要素 (タプル) を順番に取り出し、x + 1 を a に掛け算していくだけです。
簡単な実行例を示します。
>>> divisor_count(24) 8 >>> divisor_count(12345678) 24 >>> divisor_count(123456789) 12 >>> divisor_count(1234567890) 48 >>> divisor_count(1111111111) 16
n の素因数分解ができると、約数の合計値を求めるのは簡単です。素因数分解した結果が \(p^a\) の場合、その約数は \(1, p, p^2, \cdots, p^a\) となり、初項が 1 で公比が p の等比数列になります。約数の合計値を \(s(p, a)\) で表すことにすると、等比数列の和の公式から \(s(p, a)\) は以下の式になります。
たとえば、8 の素因数分解は \(2^3\) になり、素数の合計値は 8 + 4 + 2 + 1 = 15 になります。整数 n の素因数分解が \(p^a \times q^b \times r^c\) だとすると、\(s(p, a) \times s(q, b) \times s(r, c)\) は約数の合計値になります。実際に式を展開すると次のようになります。
このように、各項は約数と一対一に対応しているので、その和は約数の合計値になるわけです。たとえば、12 は \(2^2 \times 3\) に素因数分解できますが、その合計値は (4 + 2 + 1) * (3 + 1) = 28 となります。12 の約数は 1, 2, 3, 4, 6, 12 なので、その合計値は確かに 28 になります。
拙作のページ 素数 で作成した関数 factorization() を使うと、プログラムは次のようになります。
リスト : 約数の合計値
def divisor_total(n):
a = 1
for p, q in factorization(n):
a *= (p ** (q + 1) - 1) // (p - 1)
return a
for ループの中で、等比数列の和の公式を使って s(p, q) を計算します。あとはその値を累積変数 a に掛け算していくだけです。
簡単な実行例を示します。
>>> divisor_total(24) 60 >>> divisor_total(12345678) 27319968 >>> divisor_total(123456789) 178422816 >>> divisor_total(1234567890) 3211610688 >>> divisor_total(1111111111) 1246404096
p が素数の場合、\(p^a\) の約数は次のように簡単に求めることができます。
n の素因数分解が \(p^a \times q^b\) だったとすると、その約数は次のようになります。
たとえば、12 の約数は \(2^2\) = (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):
return [p ** i for i in range(q + 1)]
def divisors(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 のリストの要素と掛け合わせていくだけです。
簡単な実行例を示します。
>>> divisors(24) [1, 2, 3, 4, 6, 8, 12, 24] >>> divisors(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] >>> divisors(123456789) [1, 3, 9, 3607, 3803, 10821, 11409, 32463, 34227, 13717421, 41152263, 123456789] >>> divisors(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] >>> divisors(1111111111) [1, 11, 41, 271, 451, 2981, 9091, 11111, 100001, 122221, 372731, 2463661, 4100041, 27100271, 101010101, 1111111111]
リスト : 完全数
def perfect_number(n):
for x in range(2, n + 1):
if divisor_total(x) - x == x:
print(x, end=' ')
print()
完全数を求める perfect_number() は簡単です。x の約数の合計値を divisor_total() で求め、その値から x を引いた値が x と等しければ完全数です。print() で x を表示します。
実行結果を示します。
>>> perfect_number(10000) 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 range(1, n + 1):
y = divisor_total(x) - x
if y > x: count += 1
return count
# 不足数
def deficient_number(n):
count = 0
for x in range(1, n + 1):
y = divisor_total(x) - x
if y < x: count += 1
return count
>>> abundant_number(1000000) 247545 >>> deficient_number(1000000) 752451
1000000 以下の過剰数は 247545 個、不足数は 752451 個、完全数は 4 個なので、合計で 1000000 になります。参考 URL 3 によると、『自然数のうち過剰数が占める割合は0.2474から0.2480の間であると証明されている。』 とのことで、1000000 以下の過剰数の個数は確かにこの範囲内に入っています。
リスト : 友愛数
def amicable_number(n):
for x in range(2, n + 1):
m = divisor_total(x) - x
if m < x and x == divisor_total(m) - m:
print(m, x)
print()
友愛数を求める関数 amicable_number() も簡単です。divisor_total() で x の約数の合計値を求め、その値から x を引いた値を変数 m にセットします。m の約数の合計値から m を引いた値が x と等しければ、x と m は友愛数です。print() で x と m を表示します。同じ組を表示しないようにするため、m < x を条件に入れています。
実行結果を示します。
>>> amicable_number(100000) 220 284 1184 1210 2620 2924 5020 5564 6232 6368 10744 10856 12285 14595 17296 18416 66928 66992 67095 71145 63020 76084 69615 87633 79750 88730
なお、友愛数が無数に存在するか否かは、未解決な問題だそうです。
perfect_number() と amicable_number() は、divisor_total() を呼び出して約数の和を求めていますが、あらかじめ約数の和を計算してリスト (配列) に格納しておく方法もあります。この場合、約数の和の計算にちょっと時間がかかりますが、完全数と友愛数を求める処理は高速になります。
基本的な考え方は簡単です。約数の和を格納する配列 sum_table を用意します。sum_table の要素は 1 に初期化します。2 から順番に素数 p で割り算して 1 + p1 + ... + pq を求め、それを sum_table の値に掛け算します。素数は「エラトステネスの篩」と同じ方法で求めることができます。素数の倍数であれば、sum_table の値は 1 よりも大きくなっているはずです。つまり、sum_table が 1 であれば、その整数は素数であることがわかります。
プログラムは次のようになります。
リスト : n 以下の整数の約数の和をリストに格納して返す
def make_divisor_total(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 range(2, n + 1):
if sum_table[i] != 1: continue
for j in range(i, n + 1, i):
sum_table[j] *= factor_sub(j, i)
return sum_table
局所関数 factor_sub() は n を 素数 p で割り算して、1 + p1 + ... + pq を求めます。make_divisor_total() は sum_table を初期化してから、for ループで素数を探します。sum_table[i] が 1 でない場合、i は素数ではありません。contiune で次の数をチェックします。素数の場合は、その倍数に対応する sum_table の値を更新します。i の倍数を 変数 j にセットし、sum_table[j] に factor_sub(j, i) の値を乗算するだけです。最後に sum_table を返します。
make_divisor_total() を使うと amicable_number() は次のようになります。
リスト : 友愛数
def amicable_number1(n):
table = make_divisor_total(n)
for x in range(2, n + 1):
y = table[x] - x
if y < x and table[y] - y == x:
print(y, x)
実際に 200,000 以下の友愛数を求めてみたところ、20 個の友愛数を出力するのに amicable_number() では 3.22 秒 (実行環境 : Ubuntu 22.04 LTS (WSL2), Intel Core i5-6200U 2.30GHz) かかりましたが、amicable_number1() では約数の和を求める処理を含めて 0.52 秒ですみました。興味のある方はいろいろ試してみてください。
数学の整数論では、約数の k 乗和を求める関数を「約数関数 (divisor function)」といい、\(\sigma_k(n)\) と表記します。\(\sigma_0(n)\) は約数の個数を表し、\(d(n)\) や \(\tau(n)\) と表記されることもあります。\(\sigma_1(n)\) は約数の合計値を表し、1 を省略して \(\sigma(n)\) と表記されることもあります。
n を素因数分解すると、約数関数の値を求めるのは簡単です。素因数分解した結果が \(p^a\) の場合、その約数の k 乗は \(1, \ p^k, \ p^{k*2}, \ \cdots, \ p^{k*a}\) となり、初項が 1 で公比が \(p^k\) の等比数列になります。約数の k 乗の合計値を s(p, a, k) で表すことにすると、等比数列の和の公式から s(p, a, k) は以下の式になります。
たとえば、8 の素因数分解は 23 になり、約数の 2 乗の合計値は 64 + 16 + 4 + 1 = 85 になります。整数 n の素因数分解が \(p^a \times q^b \times r^c\) だとすると、約数関数の値は \(s(p, a, k) \times s(q, b, k) \times s(r, c, k)\) になります。
プログラムは次のようになります。
リスト : 約数関数
def divisor_sigma(n, k):
if k == 0: return divisor_count(n)
a = 1
for p, q in factorization(n):
a *= (p ** (k * (q + 1)) - 1) // (p ** k - 1)
return a
k = 0 の場合、公式の分母 p ** k - 1 が 0 になってしまい、エラーが送出されます。そこで、k = 0 のときは divisor_count() を呼び出すようにしています。あとは公式をそのままプログラムしただけなので、難しいところはないと思います。それでは実際に試してみましょう。
>>> divisor_sigma(24, 0) 8 >>> divisor_sigma(12345678, 0) 24 >>> divisor_sigma(123456789, 0) 12 >>> divisor_sigma(1234567890, 0) 48 >>> divisor_sigma(1111111111, 0) 16 >>> divisor_sigma(24, 1) 60 >>> divisor_sigma(12345678, 1) 27319968 >>> divisor_sigma(123456789, 1) 178422816 >>> divisor_sigma(1234567890, 1) 3211610688 >>> divisor_sigma(1111111111, 1) 1246404096 >>> divisor_sigma(24, 2) 850 >>> divisor_sigma(12345678, 2) 214137553857500 >>> divisor_sigma(123456789, 2) 17123257639169500 >>> divisor_sigma(1234567890, 2) 2226023493092035000 >>> divisor_sigma(1111111111, 2) 1245528410223519376
数学の整数論では、整数 n に対して、n 以下の自然数で n と互いに素となる個数を求める関数を、「オイラーのトーシェント関数 (Euler's totient function)」といいます。慣例的に \(\varphi(n)\) と表記されることから、オイラーのファイ関数 (phi function) と呼ばれることもあります。gcd() を使って力任せにプログラムすると次のようになります。
リスト : オイラーのトーシェント関数 (単純版)
def totient0(n):
c = 1
for x in range(2, n):
if gcd(n, x) == 1: c += 1
return c
gcd(n, 1) は 1 になるので、2 から n - 1 までの数をチェックします。そして、gcd() の返り値が 1 となる数の個数をカウントします。とても簡単ですね。
それでは実際に試してみましょう。
>>> for n in range(1, 31): ... print(totient0(n), end=" ") ... 1 1 2 2 4 2 6 4 6 4 10 4 12 6 8 8 16 6 18 8 12 10 22 8 20 12 18 12 28 8
とても単純なプログラムですが、n が大きくなると時間がかかるようになります。
>>> s = time.time(); print(totient0(100000)); time.time() - s 40000 0.1185765266418457 >>> s = time.time(); print(totient0(1000000)); time.time() - s 400000 1.345728874206543 >>> s = time.time(); print(totient0(10000000)); time.time() - s 4000000 15.34862232208252
幸いなことに、\(\varphi(n)\) を求める公式があるので、それを使うと高速に計算することができます。\(n\) の素因数分解が \({p_i}^a \times {p_j}^b \times \cdots \times {p_k}^c\) とすると、\(\varphi(n)\) は以下の式で求めることができます。
詳しい説明は 参考 URL 6, 8 をお読みください。
factorization() を使うと次のようになります
リスト : オイラーのトーシェント関数 (公式)
def totient(n):
d = n
for p, _ in factorization(n):
d = d * (p - 1) // p
return d
1 - 1/p を (p - 1) / p に変形して計算していることに注意してください。あとは特に難しいところはないと思います。
それでは実際に試してみましょう。
>>> for n in range(1, 31): ... print(totient(n), end=" ") ... 1 1 2 2 4 2 6 4 6 4 10 4 12 6 8 8 16 6 18 8 12 10 22 8 20 12 18 12 28 8 >>> s = time.time(); print(totient(1000000)); time.time() - s 400000 7.224082946777344e-05 >>> s = time.time(); print(totient(10000000)); time.time() - s 4000000 8.869171142578125e-05 >>> s = time.time(); print(totient(100000000)); time.time() - s 40000000 0.00010657310485839844
公式を使うことで高速に \(\varphi(n)\) を求めることができました。
#
# divisor.py : 約数
#
# Copyright (C) 2015-2023 Makoto Hiroi
#
# 素因数分解
def factor_sub(n, m):
c = 0
while n % m == 0:
c += 1
n //= m
return c, n
# wheel factorization
def factorization(n):
buff = []
wheel = [1, 2, 2, 4, 2, 4, 2, 4, 6, 2, 6]
s = 3
i, x, m = 0, 2, n
while m >= x * x:
c, m = factor_sub(m, x)
if c > 0: buff.append((x, c))
x += wheel[i]
i += 1
if i >= len(wheel): i = s
if m > 1: buff.append((m, 1))
return buff
# 最大公約数
def gcd(a, b):
if b == 0: return a
return gcd(b, a % b)
def gcd1(a, b):
while b > 0:
a, b = b, a % b
return a
# 最小公倍数
def lcm(a, b):
return a * b // gcd(a, b)
# 約数の個数
def divisor_count(n):
a = 1
for _, x in factorization(n):
a *= x + 1
return a
# 約数の和
def divisor_total(n):
a = 1
for p, q in factorization(n):
a *= (p ** (q + 1) - 1) // (p - 1)
return a
# 約数
# pq の約数を求める
def divisor_sub(p, q):
return [p ** i for i in range(q + 1)]
def divisors(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)
# 完全数
def perfect_number(n):
for x in range(2, n + 1):
if divisor_total(x) - x == x:
print(x, end=' ')
print()
# 過剰数
def abundant_number(n):
count = 0
for x in range(1, n + 1):
y = divisor_total(x) - x
if y > x: count += 1
return count
# 不足数
def deficient_number(n):
count = 0
for x in range(1, n + 1):
y = divisor_total(x) - x
if y < x: count += 1
return count
# 友愛数
def amicable_number(n):
for x in range(2, n + 1):
m = divisor_total(x) - x
if m < x and x == divisor_total(m) - m:
print(m, x)
print()
# n 以下の整数の約数の和をリストに格納して返す
def make_divisor_total(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 range(2, n + 1):
if sum_table[i] != 1: continue
for j in range(i, n + 1, i):
sum_table[j] *= factor_sub(j, i)
return sum_table
# 高速化
def amicable_number1(n):
table = make_divisor_total(n)
for x in range(2, n + 1):
y = table[x] - x
if y < x and table[y] - y == x:
print(y, x)
# 約数関数
def divisor_sigma(n, k):
if k == 0: return divisor_count(n)
a = 1
for p, q in factorization(n):
a *= (p ** (k * (q + 1)) - 1) // (p ** k - 1)
return a
# オイラーのトーシェント関数 (単純版)
def totient0(n):
c = 1
for x in range(2, n):
if gcd(n, x) == 1: c += 1
return c
# 公式を使う
def totient(n):
d = n
for p, _ in factorization(n):
d = d * (p - 1) // p
return d