約数 (divisor) は整数 n を割り切る整数 d のことです。これを \(d \ | \ n\) と表すことがあります。たとえば、24 の正の約数は 1, 2, 3, 4, 6, 8, 12, 24 の 8 個あります。ここでは正の約数を考えることにしましょう。
素数 (prime number) は正の約数が 1 と自分自身しかない自然数のことです。1 は素数に含めません。1 でない正整数で、素数ではないものを合成数 (composite number) といいます。たとえば、100 以下の素数は 25 個あります。
2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97
合成数は素数の積の形に書き表すことができます。これを素因数分解 (prime factorization) といいます。たとえば、\(6\) は \(2 \times 3\) に、\(12\) は \(2^2 \times 3\) に、\(30\) は \(2 \times 3 \times 5\) に素因数分解することができます。素数の並べ方を除けば、素因数分解はただの一通りしかありません。これを素因数分解の一意性といいます。大きな数の素因数分解はとても難しい問題で、現代の暗号 (RSA 暗号など) に用いられています。
素数の数が無限にあることは古くから知られていて、ユークリッドが「背理法」を用いて証明しています。背理法はある命題が真であることを証明するために、それが偽であることを仮定して、そこから矛盾を導くことで命題が真であることを証明する方法です。ここでは、素数の数が無限にあることを証明するわけですから、それを否定する命題、つまり「素数の数は有限個しかない」と仮定して、それが矛盾することを示せばいいわけです。
まず最初に「素数の数が有限個 (n 個) しかない」と仮定し、その最大値を p とします。2 から p までの素数を乗算して 1 を足した整数を a としましょう。a は 2 から p までの素数で割ると、割り切れずに 1 あまります。最初の仮定により a は素数ではないので素因数分解できるはずですが、2 から p までの素数では割り切れないので、これ以外の素数が存在することになり、最初の仮定と矛盾します。つまり、最初の仮定「素数の数は有限個しかない」が間違いだったわけで、素数は無数にあることが証明されました。
SymPy には約数を求める関数 divisors(n), i 番目の素数を求める関数 prime(i), 素因数分解を行う関数 factorint(n)、素数を判定する関数 isprime(n) など、便利な関数が用意されています。簡単な使用例を示します。SymPy をインポートしてください。
>>> import sympy as sy
>>> for i in range(25): print(sy.prime(i + 1), end=" ")
...
2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97
>>> for i in range(25): print(sy.prime(i + 26), end=" ")
...
101 103 107 109 113 127 131 137 139 149 151 157 163 167 173 179 181 191 193 197 199 211 223 227 229
>>> sy.divisors(24)
[1, 2, 3, 4, 6, 8, 12, 24]
>>> sy.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]
>>> sy.divisors(111111111)
[1, 3, 9, 37, 111, 333, 333667, 1001001, 3003003, 12345679, 37037037, 111111111]
>>> sy.factorint(111111111)
{3: 2, 37: 1, 333667: 1}
>>> sy.factorint(2**32-1)
{3: 1, 5: 1, 17: 1, 257: 1, 65537: 1}
>>> for n in range(100):
... if sy.isprime(n + 1): print(n + 1, end=" ")
...
2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97
\(2^n - 1\) (\(n\) は自然数) の形の自然数を「メルセンヌ数 (Mersenne number)」といい、素数であるメルセンヌ数を「メルセンヌ素数」といいます。\(n\) が 32 以下のメルセンヌ素数を求めると、次のようになります。
>>> for n in range(2,33): ... m = 2**n - 1 ... if sy.isprime(m): print(n, m) ... 2 3 3 7 5 31 7 127 13 8191 17 131071 19 524287 31 2147483647
n が 32 以下のメルセンヌ素数 \(M_n\) は n = 2, 3, 5, 7, 13, 17, 19, 31 の 8 個あります。メルセンヌ素数の場合、n = 2 を除いて n の値は素数 (奇素数) になります。ただし、n が奇素数だからといって、メルセンヌ数 \(M_n\) が素数になるとは限りません。実際、\(n = 11\) は素数ですが、\(M_{11}\) はメルセンヌ素数ではありません。
メルセンヌ素数には「リュカ-レーマー・テスト (Lucas-Lehmer primality test)」という高速な素数判定法があります。参考 URL 2 より引用します。
『\(p\) が奇素数のとき、\(M_p\) が素数となるための必要十分条件は、\(S_0 = 4, S_n = {S_{n−1}}^2 − 2 \ (n \geq 1)\) と定義したときに \(S_{p−2}\) が \(M_p\) で割り切れることである』
簡単な例を示しましょう。
M5 = 25 - 1 = 31 x0 = 4 x1 = (4 * 4 - 2) % 31 = 14 x2 = (14 * 14 - 2) % 31 = 8 x3 = (8 * 8 - 2) % 31 = 0 x3 % 31 = 0 => 素数である M9 = 29 - 1 = 511 x0 = 4 x1 = (4 * 4 - 2) % 511 = 14 x2 = (14 * 14 - 2) % 511 = 194 x3 = (194 * 194 - 2) % 511 = 331 x4 = (331 * 331 - 2) % 511 = 205 x5 = (205 * 205 - 2) % 511 = 121 x6 = (121 * 121 - 2) % 511 = 331 x7 = (331 * 331 - 2) % 511 = 205 x7 % 511 = 205 => 素数ではない
これを Python でプログラムすると次のようになります。
リスト : リュカ-レーマー・テスト
def lucas_lehmer_test(p):
m = 2 ** p - 1
x = 4
for _ in range(p - 2):
x = (x * x - 2) % m
return x % m == 0
この方法を使うと、n が 1000 未満であれば Python3 (PyPy3) でも高速にメルセンヌ素数を求めることができます。実際に n が奇素数で試してみると、メルセンヌ素数 Mn の n の値は次のようになります。
>>> for n in range(3, 1000, 2): ... if lucas_lehmer_test(n): print(n) ... 3 5 7 13 17 19 31 61 89 107 127 521 607 >>> 2 ** 607 - 1 531137992816767098689588206552468627329593117727031923199444138200403559860852242739162502265229285 668889329486246501015346579337652707239409519978766587351943831270835393219031728127 >>> miller_rabin_test(2 ** 607 - 1) True >>> miller_rabin_test(2 ** 608 - 1) False
2 つの整数 \(a, b\) について、両方とも割り切る数 (公約数) の中で最大なものを「最大公約数 (greatest common divisor)」といい、これを \(\gcd(a, b)\) と表します。負でない整数 \(a, b\) の最大公約数は次の手順で求めることができます。
これを「ユークリッドの互除法」といいます。
ユークリッドの互除法は簡単に証明できます。\(a\) と \(b\) の割り算を式 (1) で表します。
ここで、\(a\) と \(b\) の最大公約数を \(m\) とすると、\(a = m \times a', \ b = m \times b'\) となります。すると、式 (1) は式 (2) で表すことができます。
左辺は \(m\) で割り切れるので、右辺も \(m\) で割り切れる必要があります。\(q \times m \times b'\) は \(m\) で割り切れるので、\(r\) も \(m\) で割り切れることになります。つまり、\(m\) は \(b\) と \(r\) の公約数であることがわかります。\(b\) と \(r\) の最大公約数を \(m'\) とすると、式 (3) が成り立ちます。
次に、\(b = m' \times b'', \ r = m' \times r' \)として式 (1) に代入すると、式 (4) が成り立ちます。
右辺は \(m'\) で割り切れるので、\(a\) も \(m'\) で割り切れる必要があります。つまり、\(m'\) は \(a\) と \(b\) の公約数であることがわかります。したがって、式 (5) が成り立ちます。
式 (3) と (5) より \(m = m'\) となり、\(a\) と \(b\) の最大公約数は \(b\) と \(r\) の最大公約数に等しいことが証明されました。
あとは b を a とし、r を b にして同じ計算をすればいいわけです。この計算を繰り返し行うと、a と b はどんどん小さくなっていき、r = 0 になったときの b が最大公約数になります。
Python のモジュール math には最大公約数を求める関数 gcd() が用意されています。簡単な使用例を示します。
>>> import math >>> math.gcd(30, 18) 6 >>> math.gcd(123456789, 12345678) 9 >>> math.gcd(1234567890, 12345678) 18
約数の個数を求めるのは簡単です。たとえば、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 個になります。
SymPy には約数の個数を求める関数 divisor_count(n) が用意されています。簡単な実行例を示します。
>>> sy.divisor_count(24) 8 >>> sy.divisor_count(12345678) 24 >>> sy.divisor_count(123456789) 12 >>> sy.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 になります。
数学の整数論では、約数の 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)\) になります。
SymPy には約数関数 divisor_sigma(n, k = 1) が用意されています。簡単な実行例を示します。
>>> sy.divisor_sigma(24, 0) 8 >>> sy.divisor_sigma(24, 1) 60 >>> sy.divisor_sigma(24, 2) 850 >>> sy.divisor_sigma(12345678, 0) 24 >>> sy.divisor_sigma(12345678) 27319968 >>> sy.divisor_sigma(12345678, 2) 214137553857500 >>> sy.divisor_sigma(111111111, 0) 12 >>> sy.divisor_sigma(111111111) 164831992 >>> sy.divisor_sigma(111111111, 2) 13879968251176300
参考 URL 4 によると、『完全数（かんぜんすう，perfect number）とは、その数自身を除く約数の和が、その数自身と等しい自然数のことである。』 とのことです。約数関数を使って定義すると、\(\sigma(n) = 2n\) を満たす整数 \(n\) のことになります。
SymPy には整数 n が完全数か判定する関数 is_perfect(n) が用意されています。簡単な実行例を示します。
>>> for n in range(1, 10000): ... if sy.is_perfect(n): print(n) ... 6 28 496 8128
ところで、参考 URL 4 によると、メルセンヌ素数を \(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
なお、奇数の完全数はまだ発見されておらず、偶数の完全数 (つまりメルセンヌ素数) が無数に存在するか否かも未解決な問題だそうです。
参考 URL 6, 7 によると、『その数自身を除く約数の総和が元の数より大きい数』を「過剰数 (abundant number)」といい、『その数自身を除く約数の総和が元の数より小さい数』を「不足数 (deficient number)」というそうです。約数関数を使って定義すると、\(\sigma(n) \gt 2n\) を満たす整数 n が過剰数で、\(\sigma(n) \lt 2n\) を満たす整数 n が不足数になります。
SymPy には過剰数を判定する関数 is_abundant() と不足数を判定する関数 is_deficient() が用意されています。簡単な実行例を示します。
>>> sy.is_abundant(20) True >>> sy.is_abundant(15) False >>> sy.is_deficient(20) False >>> sy.is_deficient(15) True >>> a, d, p = 0, 0, 0 >>> for n in range(1, 100001): ... if sy.is_abundant(n): a += 1 ... elif sy.is_deficient(n): d += 1 ... else: p += 1 ... >>> (a, d, p) (24795, 75201, 4)
100000 以下の過剰数は 24795 個、不足数は 75201 個、完全数は 4 個なので、合計で 100000 になります。参考 URL 6 によると、『自然数のうち過剰数が占める割合は0.2474から0.2480の間であると証明されている。』 とのことで、100000 以下の過剰数の個数は確かにこの範囲内に入っています。
参考 URL 5 によると、『友愛数（ゆうあいすう）とは、異なる2つの自然数の組で、自分自身を除いた約数の和が、互いに他方と等しくなるような数をいう。』 とのことです。SymPy には友愛数を判定する関数 is_ amicable(a, b) が用意されています。
簡単な実行例を示します。
>>> sy.is_amicable(220, 284) True >>> sy.is_amicable(220, 285) False >>> for n in range(2, 100001): ... m = sy.divisor_sigma(n) - n ... if m < n and sy.is_amicable(m, n): ... print(m, n) ... 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
なお、友愛数が無数に存在するか否かは、未解決な問題だそうです。