「素数 (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 * 3、12 は 22 * 3、30 は 2 * 3 * ５ と素因数分解することができます。素数の並べ方を除けば、素因数分解はただの一通りしかありません。これを「素因数分解の一意性」といいます。大きな数の素因数分解はとても難しい問題で、現代の暗号 (RSA 暗号など) に用いられています。
今回は素数にまつわるプログラムを Python3 で作ることにします。その前に、前提となる数学 (整数論) について簡単に説明しましょう。
数学の世界では、2 つの整数 a, b を 2 以上の整数 n で割ったとき、その剰余が等しいことを「合同」といい、以下の式で表します。
これを「合同式」といい、"a 合同 b モッド n" または "n を法として a と b は合同" と読みます。
合同式の +, -, *, べき乗では、次の規則が成立します。
ただし、合同式の除算には条件があります。
「互いに素」とは、a, n の共通の約数が 1 しかないこと、最大公約数 \(\gcd(a, n)\) の値が 1 であることと同じです。また、\(f(x\)) を整数係数多項式とすると、\(a \equiv b \pmod n\) のとき、\(f(a) \equiv f(b)\) が成立します。
Python の場合、剰余は演算子 % で求めることができます。簡単な例を示しましょう。
>>> 12345678 % 256 78 >>> 123456789 % 256 21 >>> 1234567890 % 256 210 >>> 12345678900 % 256 52 >>> a = 12345678 >>> b = a % 256 >>> b 78 >>> c = 123456789 >>> d = c % 256 >>> d 21 >>> (a + c) % 256 99 >>> (b + d) % 256 99 >>> (a - c) % 256 57 >>> (b - d) % 256 57 >>> (a * c) % 256 102 >>> (b * d) % 256 102 >>> (a ** 5) % 256 224 >>> (b ** 5) % 256 224 >>> 15240 % 256 136 >>> 34440 % 256 136 >>> 15240 % 3 0 >>> 34440 % 3 0 >>> (15240 // 3) % 256 216 >>> (34440 // 3) % 256 216 >>> 15240 % 4 0 >>> 34440 % 4 0 >>> (15240 // 4) % 256 226 >>> (34440 // 4) % 256 162
もう一つ、有名な定理を紹介します。
本当にそうなるか実際に試してみましょう。
>>> a = 10 >>> b = 7 >>> for n in range(1, b): .... print(a * n % b) .... 3 6 2 5 1 4
証明は「背理法」を使うと簡単です。\(ia\) と \(ja\) (\(i \gt j\)) を \(n\) で割った余りが同じ値 \(x\) だと仮定しましょう。すると、次式が成立するはずです。
\((i-j)a\) を \(n\) で割った余りは \(0\) になるので、\((i-j)a\) は \(n\) の倍数でなければいけません。\(i, j\) の前提条件 \(1 \leq i \lt n, \ 1 \leq j \lt n\) より、\(1 \leq i - j \lt n\) になるので、\(i - j\) は \(n\) で割り切れることはありません。そして、\(a\) と \(n\) は互いに素なので、\(a\) も \(n\) で割り切れることはありません。つまり、\((i-j)a\) は \(n\) の倍数ではないので、最初の仮定「\(ia\) と \(ja\) (\(i \gt j\)) を \(n\) で割った余りが同じ値 \(x\) になる」が矛盾していて、余りはすべて異なる値になることが導かれました。
n が素数で a と n が互いに素のとき、a の n - 1 乗を n で割ると余りが必ず 1 になります。これを「フェルマーの小定理」といい、合同式で書くと次のようになります。
本当にそうなるのか、実際に試してみましょう。
>>>> 10 ** (3 - 1) % 3 1 >>>> 100 ** (3 - 1) % 3 1 >>>> 100 ** (7 - 1) % 7 1 >>>> 1000 ** (7 - 1) % 7 1 >>>> 1000 ** (11 - 1) % 11 1 >>>> 10000 ** (11 - 1) % 11 1 >>>> 10000 ** (97 - 1) % 97 1
確かに余りが 1 になりますね。
フェルマーの小定理の証明 (一例) を示します。
このほかにも、いろいろな証明方法があります。興味のある方は 参考 URL 7, 9 をお読みください。
フェルマーの小定理を Python の式で表すと、a ** (n - 1) % n == 1 となります。べき乗の剰余を「べき剰余」といいます。Python の場合、べき剰余は x ** y % n で簡単に計算できるように思いますが、x ** y が大きな値になると、剰余の計算に時間がかかるようになります。Python の関数 pow(x, y) は xy を計算しますが、第 3 引数に n を指定すると、べき剰余を求めることができます。
簡単な例を示しましょう。
>>> import time >>> s = time.time(); print(7 ** 654321 % 10); time.time() - s 7 0.09611320495605469 >>> s = time.time(); print(7 ** 6543210 % 10); time.time() - s 9 3.394315481185913 >>> s = time.time(); print(pow(7, 654321, 10)); time.time() - s 7 7.557868957519531e-05 >>> s = time.time(); print(pow(7, 6543210, 10)); time.time() - s 9 6.031990051269531e-05
このように、x ** y の値が大きくなると、x ** y % n よりも pow(x, y, n) の方が速くなります。
べき剰余を求める関数は私達でも簡単に作ることができます。次のリストを見てください。
リスト : べき剰余
def modpow(x, y, n):
r = 1
while y > 0:
if y & 1 > 0:
r = r * x % n
x = x * x % n
y >>= 1
return r
関数 modpow(x, y, n) は x ** y % n を高速に求めます。基本的には、拙作のページ 再帰定義 累乗 (べき乗) の計算で作成した関数 pow2() に剰余の計算を追加しただけです。剰余の計算回数は増えますが、乗算と剰余の計算は桁数が増えると時間がかかるので、結果的にはこちらのほうが速くなります。
それでは実際に試してみましょう。
>>> s = time.time(); print(modpow(7, 654321, 10)); time.time() - s 7 7.486343383789062e-05 >>> s = time.time(); print(modpow(7, 6543210, 10)); time.time() - s 9 6.890296936035156e-05
Python の標準関数 pow() と同じくらい高速にべき剰余を求めることができました。
整数 n 以下の素数を求めるプログラムは、拙作のページ「新・お気楽 Python プログラミング入門」 第 1 回 と 第 2 回 で取り上げました。ここでもう一つ、素数を求める簡単で高速な方法を紹介しましょう。
最初に、2 から n までの整数列を生成します。先頭の 2 は素数なので、この整数列から 2 で割り切れる整数を取り除き除きます。2 で割り切れる整数が取り除かれたので、残った要素の先頭が素数になります。先頭要素は 3 になるので、今度は 3 で割り切れる整数を取り除けばいいのです。このように、素数を見つけたらそれで割り切れる整数を取り除いていくアルゴリズムを「エラトステネスの篩 (ふるい)」といいます。
これをそのままプログラムすると、次のようになります。
リスト : エラトステネスの篩
def sieve(n):
xs = bytearray(n + 1)
for x in range(2, math.isqrt(n) + 1):
if xs[x] == 0:
for y in range(2 * x, len(xs), x):
xs[y] = 1
return [x for x in range(2, n + 1) if xs[x] == 0]
bytearray() はバイト列を生成する関数です。詳しい説明は拙作のページ お気楽 Python3 プログラミング超入門 bytes と bytearray をお読みください。
プログラムは簡単です。bytearray() で大きさ n + 1 のバイト列を生成し、変数 xs にセットします。xs の添字が数を表します。要素が 0 の数を素数とし、素数でない場合は 1 に書き換えます。バイト列は 0 に初期化されるので、最初はすべての数が素数ということになります。
最初の for ループで、x を 2 から √n まで 1 ずつ増やして、x が素数かチェックします。math.isqrt(n) は引数 n の平方根を整数で返す関数です。xs[x] が 0 ならば x は素数です。2 番目の for ループで、xs から x の倍数を削除します。最後に、内包表記で xs から素数を取り出してリストに格納して返します。
それでは実際に試してみましょう。
>>> import time >>> sieve(100) [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] >>> sieve(500) [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, 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, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499] >>> s = time.time(); print(len(sieve(100000))); time.time() - s 9592 0.014856815338134766 >>> s = time.time(); print(len(sieve(1000000))); time.time() - s 78498 0.14157962799072266 >>> s = time.time(); print(len(sieve(10000000))); time.time() - s 664579 1.4505505561828613
2 以外の偶数は素数ではないので、奇数だけ篩にかけるようにすると、実行速度はもっと速くなります。次のリストを見てください。
リスト : エラトステネスの篩 (2)
def sieve1(n):
xs = bytearray((n - 3) // 2 + 1)
for x in range(3, math.isqrt(n) + 1, 2):
y = (x - 3) // 2
if xs[y] == 0:
for i in range(y + x, len(xs), x):
xs[i] = 1
return [2] + [x * 2 + 3 for x in range(len(xs)) if xs[x] == 0]
バイト列 xs で奇数列 (3, 5, 7, ... ) を表します。奇数を変数 x とし、それに対応する xs の添字を変数 y とすると、変数 x は 3, 5, 7, 9, ... に、それに対応する変数 y は 0, 1, 2, 3, ... になります。この場合、x の倍数に対応する y の値は y + x, y + x * 2, y + x * 3, ... になります。たとえば、3, 5, 7 の倍数は次のようになります。
x | 3 5 7 9 11 13 15 17 19 21 23 25 y | 0 1 2 3 4 5 6 7 8 9 10 11 --+------------------------------------- 3 | O 0 O 0 5 | 0 0 0 7 | 0 0
プログラムは簡単です。最初の while ループで、x を √n まで +2 ずつ増やして素数かチェックします。xs の添字 y は (x - 3) / 2 で求めることができます。xs[y] が 0 ならば x は素数です。x の倍数を xs から削除します。最後に、内包表記で xs から素数を取り出し、リストの先頭に [2] を追加して返します。
実行結果を示します。
>>> s = time.time(); print(len(sieve1(100000))); time.time() - s 9592 0.0068874359130859375 >>> s = time.time(); print(len(sieve1(1000000))); time.time() - s 78498 0.07216382026672363 >>> s = time.time(); print(len(sieve1(10000000))); time.time() - s 664579 0.6974353790283203
sieve1() のほうが約 2 倍速くなりました。
自然数 low 以上 high 以下の区間にある素数を求めることを考えます。この場合、sqrt(high) 以下の素数表があれば、「エラトステネスの篩」を使って高速に素数を求めることができます。これを「区間篩」と呼びます。プログラムは次のようになります。
リスト : 区間篩
def segmented_sieve(low, high):
if low <= 2: return sieve1(high)
xs = bytearray(high - low + 1)
for p in sieve1(math.isqrt(high)):
s = math.ceil(low / p) * p - low
if low <= p: s += p
for i in range(s, len(xs), p):
xs[i] = 1
return [i + low for i in range(len(xs)) if xs[i] == 0]
low が 2 以下であれば区間篩をする必要が無いので sieve1() を呼び出します。変数 xs に low 以上 high 以下の数を表すバイト列をセットします。次の for ループで、√high 以下の素数 p を求め、xs から p の倍数を削除します。変数 s はバイト列 xs で p の倍数の最初の位置 (添字) を表します。
math.ceil(n) は、n の小数点数を切り上げた整数を返します。n が整数であれば、n をそのまま返します。簡単な実行例を示しましょう。
>>> math.ceil(2) 2 >>> math.ceil(2.4) 3 >>> math.ceil(2.6) 3 >>> math.ceil(1000 / 3) * 3 1002 >>> math.ceil(1000 / 7) * 7 1001
たとえば、区間が [1000, 2000] の場合、最初に現れる 3 の倍数は 1000 / 3 の小数点を切り上げて 334 * 3 = 1002 になるので、xs の開始位置は 2 になります。同様に、7 の場合は 1000 / 7 を切り上げて 143 * 7 = 1001 になるので、xs の開始位置は 1 になります。
それから、√high 以下の素数 p と区間 [low, high] が重なる場合があることに注意してください。p >= low の場合、最初に求めた s 番目の位置は素数になるので、xs[s] に 1 をセットしてはいけません。s に p を加算して、次の倍数に移動します。あとは、for ループで倍数を削除し、内包表記で xs にある素数を取り出してリストに格納して返します。
簡単な実行例を示します。
>>> segmented_sieve(0, 100) [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] >>> segmented_sieve(100, 200) [101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199] >>> segmented_sieve(1000, 1100) [1009, 1013, 1019, 1021, 1031, 1033, 1039, 1049, 1051, 1061, 1063, 1069, 1087, 1091, 1093, 1097] >>> segmented_sieve(100000, 100100) [100003, 100019, 100043, 100049, 100057, 100069] >>> segmented_sieve(1000000, 1000100) [1000003, 1000033, 1000037, 1000039, 1000081, 1000099] >>> segmented_sieve(0x80000000, 0x80000100) [2147483659, 2147483693, 2147483713, 2147483743, 2147483777, 2147483783, 2147483813, 2147483857, 2147483867, 2147483869, 2147483887, 2147483893]
正常に動作しているようです。sieve1() のように偶数を省くと、もう少し速くなるかもしれません。興味のある方はプログラムを改造してみてください。
自然数 n が素数であるか判定する一番簡単な方法は、n を √n までの素数 p で割り算してみることです。すべての p で割り切ることができなければ、n は素数と判定することができます。プログラムは次のようになります。
リスト : 素数の判定
def is_prime(n):
for p in sieve1(math.isqrt(n)):
if n % p == 0: return False
return True
sieve1() で √n 以下の素数 p を求め、n % p が 0 ならば n は p で割り切れるので素数ではありません。return で False を返します。for ループが終了したら、すべての素数で割り切れなかったので True を返します。
簡単な実行例を示します。
>>> is_prime(97) True >>> is_prime(111) False >>> is_prime(65537) True >>> 2 ** 31 - 1 2147483647 >>> 2 ** 32 - 1 4294967295 >>> is_prime(2 ** 31 - 1) True >>> is_prime(2 ** 32 - 1) False
正常に動作していますね。ただし、大きな数になると is_prime() で素数を判定するのは困難になるでしょう。このため、確率的に素数判定する方法が考案されています。
フェルマーの小定理の対偶を用いて、確率的な素数判定を行うことができます。これを「フェルマーテスト (Fermat primality test)」といいます。以下にアルゴリズムを示します。
これを複数回繰り返して \(n\) が「確率的素数」と判定されれば、\(n\) は素数である確率が高いといえます。
プログラムは次のようになります。
リスト : フェルマーテスト
def fermat_test(n, k = 20):
if n < 2: return False
elif n == 2: return True
elif n % 2 == 0: return False
else:
for _ in range(k):
a = random.randrange(2, n)
if math.gcd(n, a) != 1 or pow(a, n - 1, n) != 1:
return False
return True
モジュール random の関数 randrange(s, e) は、range(s, e) の範囲から要素を一つランダムに選択します。最大公約数を求める関数 gcd() はモジュール math に定義されています。引数 k はテストの回数を表します。あとは、アルゴリズムをそのままプログラムしただけなので、難しいところはないと思います。
それでは実際に試してみましょう。
>>> fermat_test(2) True >>> fermat_test(97) True >>> fermat_test(111) False >>> fermat_test(2 ** 31 - 1) True >>> fermat_test(2 ** 32 - 1) False >>> for n in range(6): ... x = 2 ** (2 ** n) + 1 ... print(x, end=" ") ... print(is_prime(x), end= " ") ... print(fermat_test(x)) ... 3 True True 5 True True 17 True True 257 True True 65537 True True 4294967297 False False >>> 2 ** (2 ** 6) + 1 18446744073709551617 >>> fermat_test(2 ** (2 ** 6) + 1) False >>> 2 ** (2 ** 7) + 1 340282366920938463463374607431768211457 >>> fermat_test(2 ** (2 ** 7) + 1) False >>> fermat_test(111111111111111111) False >>> fermat_test(1111111111111111111) True >>> fermat_test(11111111111111111111) False
正常に動作していますね。\(F(n) = 2^{2^n} + 1\) を「フェルマー数」といい、素数であるフェルマー数を「フェルマー素数」といいます。現在、判明しているフェルマー素数は F(0) から F(4) までの 5 つだけです。F(5) 以降のフェルマー数に素数があるかどうかはわかっていません。
n が合成数で a と n が互いに素のとき、\(a^{n-1} \equiv 1 \pmod n\) が成立する数を「偽素数」といいます。たとえば、a = 2 のときの偽素数を求めると、次のようになります。
>>> for n in range(2, 5000): ... if not is_prime(n) and pow(2, n - 1, n) == 1: ... print(n) ... 341 561 645 1105 1387 1729 1905 2047 2465 2701 2821 3277 4033 4369 4371 4681
そして、n と互いに素となる a において、すべての値で \(a^{n-1} \equiv 1 \pmod n\) が成り立つ数を「カーマイケル数」といいます。a が素数であれば、条件を満たすのは当然ですね。残りの合成数がすべて条件を満たすとカーマイケル数になります。341 と 561 を調べてみましょう。
>>> for a in range(1, 341): ... if math.gcd(a, 341) == 1 and not is_prime(a) and pow(a, 341 - 1, 341) != 1: ... print(a) ... 6 9 10 ... 略 ... 335 336 338 >>> for a in range(1, 561): ... if math.gcd(a, 561) == 1 and not is_prime(a) and pow(a, 561 - 1, 561) != 1: ... print(a) ... >>>
341 の場合、341 と互いに素で、素数ではなく、\(a^{340} \equiv 1 \pmod {341}\) にはならない a が多数存在するので、カーマイケル数ではありません。561 の場合、561 と互いに素で、素数ではないすべての a が \(a^{560} \equiv 1 \pmod {561}\) を満たすので、561 はカーマイケル数になります。
関数 fermat_test() は条件式 math.gcd(n, a) != 1 or pow(a, n - 1, n) != 1 を満たすものを合成数と判定しています。カーマイケル数の場合、or の右辺式は必ず成立するので、or の左辺式だけで素数判定を行うことになります。つまり、たまたま gcd(a, n) != 1 となる a を選択すれば、n を合成数と判定することができますが、そうでなければ n を「確率的素数」と判定してしまうのです。
カーマイケル数が有限個しかなければ取り除くこともできるでしょうが、残念ながらカーマイケル数は無数に存在することが知られています。
カーマイケル数 561, 1105, 1729, 2465, 2821, 6601, 8911, ...
このため、確率的な素数判定にはフェルマーテストを改良した「ミラーラビン素数判定法」を用いるのが一般的です。
ミラーラビン素数判定法 (以下ミラーラビンテストと記述) は、フェルマーテストと同様に与えられた数が素数か確率的に判定するアルゴリズムです。整数 n を奇素数、a を整数とすると、フェルマーの小定理により次式が成立することを利用します。
n が奇素数でなければ上式は成り立たないので、n は合成数と判定することができます。以下に判定方法の概要を示します。
mod n の世界では、-1 は n - 1 と同じです。具体的には、1 以上 n 未満の範囲で a をランダムに選んで、n が確率的素数か判定します。これを何回か繰り返して確率的素数と判定されれば、n は素数である可能性が高いといえます。ミラーラビンテストの場合、テストを k 回繰り返すと、誤判定する確率は最大で \(4^{-k}\) になることが知られています。たとえば、k = 20 では 4-20 = 9.094947e-13 になります。
また、n の範囲を限定すると、ミラーラビンテストを決定的な素数判定法に変更することができます。たとえば、n が 4759123141 未満の場合、a の値は 2, 7, 61 の 3 つを調べるだけで、素数か否か判定することができます。n が 264 以下であれば、次の値を調べるだけで判定することができます。
2, 325, 9375, 28178, 450775, 9780504, 1795265022
プログラムは次のようになります。
リスト : ミラーラビン素数判定法
# n を (2 ** c) * d に分解する
def factor2(n):
c = 0
while n % 2 == 0:
n //= 2
c += 1
return c, n
def miller_rabin_sub(n, s, d, ks):
for a in ks:
if a >= n: break
if pow(a, d, n) != 1:
r = 0
while r < s:
if pow(a, (2**r) * d, n) == n - 1: break
r += 1
else:
return False
return True
def miller_rabin_test(n, k = 20):
if n < 2: return False
elif n == 2: return True
elif n % 2 == 0: return False
else:
s, d = factor2(n - 1)
if n < 4759123141:
ks = [2, 7, 61]
elif n <= 0x10000000000000000:
ks = [2, 325, 9375, 28178, 450775, 9780504, 1795265022]
else:
ks = [random.randrange(1, n) for _ in range(k)]
return miller_rabin_sub(n, s, d, ks)
今まで説明したことをそのままプログラムしただけなので、特に難しいところはないと思います。それでは実行してみましょう。
>>> miller_rabin_test(2) True >>> miller_rabin_test(97) True >>> miller_rabin_test(111) False >>> miller_rabin_test(2 ** (2 ** 0) + 1) True >>> miller_rabin_test(2 ** (2 ** 1) + 1) True >>> miller_rabin_test(2 ** (2 ** 2) + 1) True >>> miller_rabin_test(2 ** (2 ** 3) + 1) True >>> miller_rabin_test(2 ** (2 ** 4) + 1) True >>> miller_rabin_test(2 ** (2 ** 5) + 1) False >>> miller_rabin_test(2 ** 31 - 1) True >>> miller_rabin_test(2 ** 32 - 1) False >>> miller_rabin_test(111111111111111111) False >>> miller_rabin_test(1111111111111111111) True >>> miller_rabin_test(11111111111111111111) False
正常に動作しているようです。興味のある方はいろいろ試してみてください。
\(2^n - 1\) (n は自然数) の形の自然数を「メルセンヌ数 (Mersenne number)」といい、素数であるメルセンヌ数を「メルセンヌ素数」といいます。n が 32 以下のメルセンヌ素数を求めると、次のようになります。
>>> for n in range(2, 33): ... m = 2 ** n - 1 ... print(m, end=" ") ... print(fermat_test(m), end=" ") ... print(miller_rabin_test(m)) ... 3 True True 7 True True 15 False False 31 True True 63 False False 127 True True 255 False False 511 False False 1023 False False 2047 False False 4095 False False 8191 True True 16383 False False 32767 False False 65535 False False 131071 True True 262143 False False 524287 True True 1048575 False False 2097151 False False 4194303 False False 8388607 False False 16777215 False False 33554431 False False 67108863 False False 134217727 False False 268435455 False False 536870911 False False 1073741823 False False 2147483647 True True 4294967295 False False
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 4 より引用します。
『\(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 である素数の組を「双子素数 (twin prime)」といいます。整数 n 以下の双子素数 (p, p + 2), p <= n は、n が大きくなければエラトステネスの篩で簡単に求めることができます。
リスト : 双子素数
def twin_primes(n):
xs = sieve1(n + 2)
return [(xs[i], xs[i + 1]) for i in range(len(xs) - 1) if xs[i] + 2 == xs[i + 1]]
n + 2 以下の素数を sieve1() で求めて変数 xs にセットします。あとは、xs[i] + 2 == xs[i + 1] を満たす組を内包表記で求めるだけです。それでは実行してみましょう。
>>> twin_primes(100) [(3, 5), (5, 7), (11, 13), (17, 19), (29, 31), (41, 43), (59, 61), (71, 73)] >>> twin_primes(500) [(3, 5), (5, 7), (11, 13), (17, 19), (29, 31), (41, 43), (59, 61), (71, 73), (101, 103), (107, 109), (137, 139), (149, 151), (179, 181), (191, 193), (197, 199), (227, 229), (239, 241), (269, 271), (281, 283), (311, 313), (347, 349), (419, 421), (431, 433), (461, 463)] >>> a = twin_primes(3000000) >>> len(a) 20932 >>> a[-1] (2999831, 2999833)
3,000,000 以下の双子素数は全部で 20932 個で、その最大の組は (2999831, 2999833) になります。なお、参考 URL 3 によると、『双子素数は無数に存在するかという問題、いわゆる「双子素数の予想」や「双子素数の問題」は、いまだに数学上の未解決問題である。無数に存在するだろう、とは、多くの数論学者が予想している。』 とのことです。
#
# prime.py : 素数
#
# Copyright (C) 2023 Makoto Hiroi
#
import math, random
# べき剰余 x ** y % n を求める
def modpow(x, y, n):
r = 1
while y > 0:
if y & 1 > 0:
r = r * x % n
x = x * x % n
y >>= 1
return r
# エラトステネスの篩 (単純版)
def sieve(n):
xs = bytearray(n + 1)
for x in range(2, math.isqrt(n) + 1):
if xs[x] == 0:
for y in range(2 * x, len(xs), x):
xs[y] = 1
return [x for x in range(2, n + 1) if xs[x] == 0]
def sieve1(n):
xs = bytearray((n - 3) // 2 + 1)
for x in range(3, math.isqrt(n) + 1, 2):
y = (x - 3) // 2
if xs[y] == 0:
for i in range(y + x, len(xs), x):
xs[i] = 1
return [2] + [x * 2 + 3 for x in range(len(xs)) if xs[x] == 0]
# 区間篩 (low 以上 high 以下)
def segmented_sieve(low, high):
if low <= 2: return sieve1(high)
xs = bytearray(high - low + 1)
for p in sieve1(math.isqrt(high)):
s = math.ceil(low / p) * p - low
if low <= p: s += p
for i in range(s, len(xs), p):
xs[i] = 1
return [i + low for i in range(len(xs)) if xs[i] == 0]
# 素数の判定
def is_prime(n):
for p in sieve1(math.isqrt(n)):
if n % p == 0: return False
return True
# フェルマーテスト
def fermat_test(n, k = 20):
if n < 2: return False
elif n == 2: return True
elif n % 2 == 0: return False
else:
for _ in range(k):
a = random.randrange(2, n)
if math.gcd(n, a) != 1 or pow(a, n - 1, n) != 1:
return False
return True
# ミラーラビンテスト
# n を (2 ** c) * d に分解する
def factor2(n):
c = 0
while n % 2 == 0:
n //= 2
c += 1
return c, n
def miller_rabin_sub(n, s, d, ks):
for a in ks:
if a >= n: break
if pow(a, d, n) != 1:
r = 0
while r < s:
if pow(a, (2**r) * d, n) == n - 1: break
r += 1
else:
return False
return True
def miller_rabin_test(n, k = 20):
if n < 2: return False
elif n == 2: return True
elif n % 2 == 0: return False
else:
s, d = factor2(n - 1)
if n < 4759123141:
ks = [2, 7, 61]
elif n <= 0x10000000000000000:
ks = [2, 325, 9375, 28178, 450775, 9780504, 1795265022]
else:
ks = [random.randrange(1, n) for _ in range(k)]
return miller_rabin_sub(n, s, d, ks)
# リュカ-レーマー・テスト
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
# 双子素数
def twin_primes(n):
xs = sieve1(n + 2)
return [(xs[i], xs[i + 1]) for i in range(len(xs) - 1) if xs[i] + 2 == xs[i + 1]]