「約数 (divisor)」は整数 n を割り切る整数 d のことです。これを d | n と表すことがあります。たとえば、24 の正の約数は 1, 2, 3, 4, 6, 8, 12, 24 の 8 個あります。ここでは正の約数を考えることにしましょう。整数 n を素因数分解すると、約数の個数や合計値、約数の集合を簡単に求めることができます。
約数の個数を求めるのは簡単です。たとえば、n の素因数分解が \(p^a \times q^b \times r^c\) だとしましょう。n の約数を d とすると、d は次の式で表すことができます。
\(p^i\) が表す数は \(p^a, \ p^{a-1}, \ \cdots, \ p, 1\) の a + 1 通りあります。同様に、\(q^j\) は b + 1 通り、\(r^k\) は c + 1 通りあるので、約数の個数は次式のようになります。
たとえば、12 は \(2^2 \times 3^1\) になるので、約数の個数は 3 * 2 = 6 になります。実際、12 の約数は 1, 2, 3, 4, 6, 12 の 6 個です。24 は \(2^3 \times 3^1\) になるので、約数は (3 + 1) * (1 + 1) = 8 個になります。
拙作のページ「素因数分解」で作成した関数 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 が素数の場合、pa の約数は次のように簡単に求めることができます。
n の素因数分解が \(p^a \times q^b\) だったとすると、その約数は次のようになります。
たとえば、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):
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]
数学の整数論では、約数の 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
参考 URL『完全数 - Wikipedia』によると、『完全数（かんぜんすう，perfect number）とは、その数自身を除く約数の和が、その数自身と等しい自然数のことである。』 とのことです。約数関数を使って定義すると、\(\sigma(n) = 2n\) を満たす整数 \(n\) のことになります。完全数を求める perfect_number() は簡単です。次のリストを見てください。
リスト : 完全数
def perfect_number(n):
for x in range(2, n + 1):
if divisor_total(x) - x == x:
print(x, end=' ')
print()
x の約数の合計値を divisor_total() で求め、その値から x を引いた値が x と等しければ完全数です。print() で x を表示します。
実行結果を示します。
>>> perfect_number(10000) 6 28 496 8128
ところで、参考 URL 『完全数 - Wikipedia』によると、メルセンヌ素数を \(M_n\) とすると、偶数の完全数は \(2^{n-1} \times M_n\) で表すことができるそうです。この式を使うと偶数の完全数は次のようになります。
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 『過剰数 - Wikipedia』と『不足数 - Wikipedia』によると、『その数自身を除く約数の総和が元の数より大きい数』を「過剰数 (abundant number)」といい、『その数自身を除く約数の総和が元の数より小さい数』を「不足数 (deficient number)」というそうです。約数関数を使って定義すると、\(\sigma(n) \gt 2n\) を満たす整数 n が過剰数で、\(\sigma(n) \lt 2n\) を満たす整数 n が不足数になります。
過剰数と不足数の個数を求めるプログラムも簡単に作ることができます。次のリストを見てください。
リスト : 過剰数と不足数
# 過剰数
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『過剰数 - Wikipedia』によると、『自然数のうち過剰数が占める割合は0.2474から0.2480の間であると証明されている。』 とのことで、1000000 以下の過剰数の個数は確かにこの範囲内に入っています。
参考 URL『友愛数 - Wikipedia』によると、『友愛数（ゆうあいすう）とは、異なる2つの自然数の組で、自分自身を除いた約数の和が、互いに他方と等しくなるような数をいう。』 とのことです。友愛数を求める関数 amicable_number() も簡単です。次のリストを見てください。
リスト : 友愛数
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()
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 + p^1 + \cdots + p^q\) を求め、それを 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.2 秒 (実行環境 : Python 3.10.12, Ubuntu 22.04 LTS (WSL2), Intel Core i5-6200U 2.30GHz) かかりましたが、amicable_number1() では約数の和を求める処理を含めて 0.48 秒ですみました。興味のある方はいろいろ試してみてください。
数学の整数論では、整数 n に対して、n 以下の自然数で n と互いに素となる個数を求める関数を、「オイラーのトーシェント関数 (Euler's totient function)」といいます。慣例的に \(\varphi(n)\) と表記されることから、オイラーのファイ関数 (phi function) と呼ばれることもあります。
たとえば、1, 2, 3, 4 の中で 4 と互いに素となる数は 1 と 3 の 2 つあるので、totient(4) は 2 になります。totient(5) は、5 が素数なので、1, 2, 3, 4 と 5 は互いに素となるため、totient(5) は 4 になります。p を素数とすると、totient(p) = p - 1 になります。
totient(n) を gcd() を使って力任せにプログラムすると次のようになります。
リスト : オイラーのトーシェント関数 (単純版)
import math
def totient0(n):
c = 1
for x in range(2, n):
if math.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)\) は以下の式で求めることができます。
この他に以下の式が成り立ちます。
ファイ関数の和は、n の約数を d とすると、\(\varphi(d)\) の合計値を求めることです。たとえば、8 の約数は 1, 2, 4, 8 の 4 通りありますが、\(\varphi(1) + \varphi(2) + \varphi(4) + \varphi(8)\) を計算すると、1 + 1 + 2 + 4 = 8 になります。オイラーの定理で n が素数の場合が「フェルマーの小定理」になります。詳しい説明は参考 URL『オイラーのファイ関数のイメージ』, 『約数関数 - Wikipedia』をお読みください。
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)\) を求めることができました。
メビウス関数 (Möbius function) は、正の整数 n を与えると -1, 0, 1 のいずれかの値を返します。
メビウス関数 \(\mu(n)\) は引数 n が 1 の場合は 1 を返します。1 より大きい場合、n を素因数分解した結果で値が決まります。同じ素数で 2 回割り切れる (平方因子を持つ) 場合は 0 を返します。それ以外の場合、n は異なる素数の積で表されるので、その個数 k が奇数ならば -1 を、偶数ならば 1 を返します。
たとえば、24 は \(2^3 \times 3\) と素因数分解できますが、2 で 2 回割り切れるので \(\mu(24)\) は 0 になります。つまり、素因数分解して 2 以上の指数がある場合は 0 を返します。\(\mu(30)\) は \(30 = 2 \times 3 \times 5\) なので、\((-1)^3 = -1\) を返します。n が素数の場合も -1 になります。\(\mu(35)\) は \(35 = 5 \times 7\) なので、\((-1)^2 = 1 \) を返します。
関数 factorization() を使うと、プログラムは次のようになります。
リスト : メビウス関数
def mobius(n):
if n == 1: return 1
ps = factorization(n)
for _, q in ps:
if q > 1: return 0
return (-1) ** len(ps)
定義をそのままプログラムしただけなので、難しいところはないと思います。試し割りをしながら値を求めることも可能です。興味のある方はプログラムを作ってみてください。
それでは実際に試してみましょう。
>>> mobius(1) 1 >>> mobius(2) -1 >>> mobius(4) 0 >>> mobius(6) 1 >>> mobius(30) -1 >>> for n in range(1, 30): print(mobius(n), end=" ") ... 1 -1 -1 0 -1 1 -1 0 0 1 -1 0 -1 1 1 0 -1 0 -1 0 1 1 -1 0 0 1 0 0 -1
n の約数を d とすると、n != 1 のとき \(\mu(d)\) の合計値は 0 になります。
それでは実際に試してみましょう。
>>> divisors(24) [1, 2, 3, 4, 6, 8, 12, 24] >>> [mobius(x) for x in divisors(24)] [1, -1, -1, 0, 1, 0, 0, 0] >>> sum(mobius(x) for x in divisors(24)) 0 >>> sum(mobius(x) for x in divisors(12345678)) 0 >>> sum(mobius(x) for x in divisors(123456789)) 0 >>> sum(mobius(x) for x in divisors(1234567890)) 0 >>> sum(mobius(x) for x in divisors(1111111111)) 0
自然数 n を引数に受け取り、数を返す関数 f(n), g(n) を考えます。
式 1 を満たす場合、式 2 が成り立ちます。これを「メビウスの反転公式」といいます。式 2 の場合、d は n の約数なので、n / d も n の約数です。d が決まれば n / d は一意に定まるので、式 2 の \(\mu(n/d) g(d)\) の和と \(\mu(d) g(n/d)\) の和は等しくなります。たとえば、n = 6 とすると、約数は 1, 2, 3, 6 の 4 通りあります。式 2 を計算すると、以下のように総和は等しくなることがわかります。
詳しい説明は参考 URL『メビウスの反転公式の証明と応用』, 『メビウス関数 - Wikipedia』をお読みくださいませ。
簡単な例として、ファイ関数の和の式にメビウスの反転公式を適用してみましょう。
実際に \(\varphi(n)\) を計算することができるのか、Python3 で試してみましょう。
>>> def phi(n): ... return sum(mobius(n // d) * d for d in divisors(n)) >>> phi(1000) 400 >>> totient(1000) 400 >>> phi(12345) 6576 >>> totient(12345) 6576
phi() と totient() は同じ値を返します。メビウスの反転公式は正しいことがわかります。
#
# divisor.py : 約数とオイラーのトーシェント関数
#
# Copyright (C) 2023 Makoto Hiroi
#
import math
# 素因数分解
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 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
# 約数
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 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 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 totient0(n):
c = 1
for x in range(2, n):
if math.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
# メビウス関数
def mobius(n):
if n == 1: return 1
ps = factorization(n)
for _, q in ps:
if q > 1: return 0
return (-1) ** len(ps)