「素数 (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
今回は素数に関連する問題を出題します。
素数を求めるいちばん簡単な方法は、奇数 3, 5, 7, 9, ... をそれまでに見つかった素数で割ってみることです。見つかった素数はリストに格納しておけばいいでしょう。プログラムは次のようになります。使用するプログラミング言語は Python3 (PyPy3, ver 7.3.1) です。
リスト : 素数を求める (prime.py)
import time
def solver0(n):
primes = [2]
for x in range(3, n + 1, 2):
for p in primes:
if p * p > x:
primes.append(x)
break
if x % p == 0: break
return primes
def test1(func, n):
s = time.time()
ps = func(n)
print(len(ps), ps[-1])
print('{:.3f}'.format(time.time() - s))
変数 primes は素数のリストで [2] に初期化します。奇数の生成は range を使うと簡単です。最初の for 文で奇数を取り出して変数 x にセットします。次の for 文で primes から素数を取り出して変数 p にセットします。
x が素数か判別する場合、√x より小さい素数を調べるだけで OK です。p * p が x よりも大きくなったら、x は素数であることがわかります。append() で p を primes に追加して break で for ループを脱出します。x % p が 0 であれば、x は素数 p で割り切れたので x は素数ではありません。break で for ループを脱出します。最後に return で primes を返します。
それでは実行してみましょう。
>>>> from prime import * >>>> test1(solver0, 3000000) 216816 2999999 1.257 >>>> test1(solver0, 3000000) 216816 2999999 1.236
素数の個数は 216816 個、最大の素数は 2999999 で、実行時間は約 1.25 秒 (実行環境 : Ubuntu 20.04 LTS (WSL1), Intel Core i5-6200U 2.30GHz) でした。
ところで、素数 p は 2, 3 を除外すると p % 6 の値は 1 または 5 になります。
: 5 7 11 13 17 19 23 29 31 37 41 43 47 ...
------+--------------------------------------------
p % 6 : 5 1 5 1 5 1 5 5 1 1 5 1 5 ...
つまり、6 * i ± 1 (i = 1, 2, 3, ...) の数値だけを調べればよいわけです。これをプログラムすると次のようになります。
リスト : 高速化
def solver1(n):
primes = [2, 3]
x = 5
while x <= n:
for p in primes:
if p * p > x:
primes.append(x)
break
if x % p == 0: break
if x % 6 == 5:
x += 2
else:
x += 4
return primes
x % 6 == 5 ならば x に 2 を加算し、そうでなければ 4 を加算します。これで 6 * i ± 1 (i = 1, 2, 3, ...) の数値だけ調べることができます。
それでは実際に試してみましょう。
>>>> test1(solver1, 3000000) 216816 2999999 0.835 >>>> test1(solver1, 3000000) 216816 2999999 0.824
実行時間は 0.824 秒になりました。残念ながら、劇的には速くなりませんでした。
もう一つ、素数を求める簡単な方法を紹介しましょう。
最初に、2 から N までの整数列を生成します。先頭の 2 は素数なので、この整数列から 2 で割り切れる整数を取り除き除きます。2 で割り切れる整数が取り除かれたので、残った要素の先頭が素数になります。先頭要素は 3 になるので、今度は 3 で割り切れる整数を取り除けばいいのです。このように、素数を見つけたらそれで割り切れる整数を取り除いていくアルゴリズムを「エラトステネスの篩 (ふるい)」といいます。
プログラムは次のようになります。
リスト : エラトステネスの篩
def sieve(n):
primes = [True] * (n // 2)
ps = [2]
x = 3
while x * x <= n:
y = (x - 3) // 2
if primes[y]:
ps.append(x)
y += x
while y < len(primes):
primes[y] = False
y += x
x += 2
while x <= n :
if primes[(x - 3) // 2]: ps.append(x)
x += 2
return ps
リスト primes で奇数列 (3, 5, 7, ... ) を表します。True で素数を表し、素数でない場合は False に書き換えます。primes は True で初期化されるので、最初はすべての数が素数ということになります。
奇数を変数 x とし、それに対応する primes の添字を変数 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 ずつ増やして素数かチェックします。primes の添字 y は (x - 3) / 2 で求めることができます。primes[y] が True ならば x は素数です。x を ps に追加して、x の倍数を primes から削除します。次の while ループで、√n よりも大きい素数を求めます。最後に return で ps を返します。
それでは実際に試してみましょう。
>>>> test1(sieve, 3000000) 216816 2999999 0.149 >>>> test1(sieve, 3000000) 216816 2999999 0.106 >>>> test1(sieve, 3000000) 216816 2999999 0.082 >>>> test1(sieve, 3000000) 216816 2999999 0.082
PyPy は JIT コンパイラを使用しているので、初回の実行には時間がかかることがあります。2 回目で約 0.1 秒、3 回目以降は 0.1 秒未満で実行することができました。
ところで、エラトステネスの篩は次のようにプログラムすることもできます。
リスト : エラトステネスの篩 (2)
def sieve2(n):
ps = [2]
xs = list(range(3, n + 1, 2))
while xs[0] * xs[0] <= n:
m = xs[0]
ps.append(m)
xs = [x for x in xs if x % m != 0]
return ps + xs
素数の倍数を削除するとき、内包表記で新しいリストを作っています。プログラムは簡単ですが、新しいリストを作る分だけ時間がかかります。実際試してみると、実行時間は約 2.5 秒になり、一番遅くなってしまいました。他のプログラミング言語では、この方法でも高速に実行できるかもしれません。興味のある方は試してみてください。
素因数分解はエラトステネスの篩と同じ考え方で行うことができます。プログラムは次のようになります。
リスト : 素因数分解
def factor_sub(n, m):
c = 0
while n % m == 0:
c += 1
n //= m
return c, n
def factorization(n):
buff = []
c, m = factor_sub(n, 2)
if c > 0: buff.append((2, c))
c, m = factor_sub(m, 3)
if c > 0: buff.append((3, c))
x = 5
while m >= x * x:
c, m = factor_sub(m, x)
if c > 0: buff.append((x, c))
if x % 6 == 5:
x += 2
else:
x += 4
if m > 1: buff.append((m, 1))
return buff
最初に 2 と 3 で割り算します。それから、６＊i ± 1 (i = 1, 2, 3, ...) で割り算していきます。割り算するときは、その数で割り切れるあいだは割り算を続けることに注意してください。たとえば、27 を素因数分解すると 3 * 3 * 3 になりますが、3 を一回だけしか割り算しないと、結果は 3 * 9 のように素数ではない数が含まれてしまいます。この処理を関数 factor_sub() で行っています。
あとは、factor_sub() の返り値をチェックして、割り算した回数 c が 0 よりも大きければ、素数と c をタプルに格納して、それを buff に追加します。これを変数 m が x * x 以下のあいだ繰り返します。最後に m が 1 よりも大きければ、(m, 1) を buff に追加して、return で buff を返します。
それでは実行してみましょう。
>>>> factorization(3000001) [(853, 1), (3517, 1)]
3000001 は 853 * 3517 と素因数分解できるので素数ではありません。なお、これはとても単純なアルゴリズムなので、大きな整数の素因数分解には適していません。巨大な合成数の素因数分解はとても難しい問題です。興味のある方は素因数分解について調べてみてください。
大きな整数ではないので、素数のチェックは factorization() を使いましょう。プログラムは次のようになります。
リスト : フィボナッチ素数
# フィボナッチ数
def fibo(n):
a1, a2 = 0, 1;
for _ in range(n):
a1, a2 = a2, a1 + a2
return a1
# フィボナッチ素数
def prime_fibo(n):
cnt = 0
i = 2
while True:
x = fibo(i)
if x > n: break
xs = factorization(x)
if len(xs) == 1 and xs[0][1] == 1:
cnt += 1
print(i, ":", x)
i += 1
print(cnt)
関数 fibo() でフィボナッチ数を求め、それを factorization() で素因数分解します。返り値のリスト xs の長さが 1 で因子の個数が 1 であれば素数です。print() でフィボナッチ数を表示します。実行結果は次のようになりました。
>>>> prime_fibo(3000000) 3 : 2 4 : 3 5 : 5 7 : 13 11 : 89 13 : 233 17 : 1597 23 : 28657 29 : 514229 9
フィボナッチ素数の個数はそれほど多くはないようです。なお、参考 URL 2 によると、『フィボナッチ素数が無限にあるかどうかは分かっていない。』 とのことです。
3,000,000 以下の双子素数であれば、単純な方法でも簡単に求めることができます。次のプログラムを見てください。
リスト : 双子素数
# 素数ジェネレータ
def make_primes(n):
primes = [True] * (n // 2)
yield 2
x = 3
while x * x <= n:
y = (x - 3) // 2
if primes[y]:
yield x
y += x
while y < len(primes):
primes[y] = False
y += x
x += 2
while x <= n :
if primes[(x - 3) // 2]: yield x
x += 2
# 双子素数
def twin_prime(n):
p = 1
q = 0
count = 0
for x in make_primes(n):
if x - p == 2:
count += 1
q = p
p = x
print(count)
print(q, q + 2)
関数 make_primes() は素数を生成するジェネレータです。これはエラトステネスの篩で作成した関数 sieve() をジェネレータに改造しただけです。あとは、ジェネレータから素数を順番に取り出して変数 x にセットし、一つ前の素数 p との差が 2 であれば双子素数であることがわかります。p を q にセットして、次の双子素数を探します。最後に、双子素数の個数 count と双子素数 q, q + 2 を出力します。
実行結果は次のようになりました。
>>>> twin_prime(3000000) 20932 2999831 2999833
3,000,000 以下の双子素数は全部で 20932 個で、その最大の組は (2999831, 2999833) になりました。なお、参考 URL 3 によると、『双子素数は無数に存在するかという問題、いわゆる「双子素数の予想」や「双子素数の問題」は、いまだに数学上の未解決問題である。無数に存在するだろう、とは、多くの数論学者が予想している。』 とのことです。
n が 32 以下のメルセンヌ素数は大きな値ではないので、単純な factorization() でも素数判定が可能です。実際、factorization() でメルセンヌ数を素因数分解すると次のようになります。
2 [(3, 1)] 3 [(7, 1)] 4 [(3, 1), (5, 1)] 5 [(31, 1)] 6 [(3, 2), (7, 1)] 7 [(127, 1)] 8 [(3, 1), (5, 1), (17, 1)] 9 [(7, 1), (73, 1)] 10 [(3, 1), (11, 1), (31, 1)] 11 [(23, 1), (89, 1)] 12 [(3, 2), (5, 1), (7, 1), (13, 1)] 13 [(8191, 1)] 14 [(3, 1), (43, 1), (127, 1)] 15 [(7, 1), (31, 1), (151, 1)] 16 [(3, 1), (5, 1), (17, 1), (257, 1)] 17 [(131071, 1)] 18 [(3, 3), (7, 1), (19, 1), (73, 1)] 19 [(524287, 1)] 20 [(3, 1), (5, 2), (11, 1), (31, 1), (41, 1)] 21 [(7, 2), (127, 1), (337, 1)] 22 [(3, 1), (23, 1), (89, 1), (683, 1)] 23 [(47, 1), (178481, 1)] 24 [(3, 2), (5, 1), (7, 1), (13, 1), (17, 1), (241, 1)] 25 [(31, 1), (601, 1), (1801, 1)] 26 [(3, 1), (2731, 1), (8191, 1)] 27 [(7, 1), (73, 1), (262657, 1)] 28 [(3, 1), (5, 1), (29, 1), (43, 1), (113, 1), (127, 1)] 29 [(233, 1), (1103, 1), (2089, 1)] 30 [(3, 2), (7, 1), (11, 1), (31, 1), (151, 1), (331, 1)] 31 [(2147483647L, 1)] 32 [(3, 1), (5, 1), (17, 1), (257, 1), (65537L, 1)]
n が 32 以下のメルセンヌ素数 Mn は n = 2, 3, 5, 7, 13, 17, 19, 31 の 8 個あります。メルセンヌ素数の場合、n = 2 を除いて n の値は素数 (奇素数) になります。ただし、n が奇素数だからといって、メルセンヌ数 Mn が素数になるとは限りません。実際、n = 11 は素数ですが、M11 はメルセンヌ素数ではありません。
ただし、これ以上大きなメルセンヌ数を単純な素因数分解で素数判定するのは時間がかかって現実的ではありません。幸いなことに、メルセンヌ素数は「リュカ-レーマー・テスト (Lucas-Lehmer primality test)」という高速な素数判定法が確立しています。参考 URL 4 より引用します。
『p が奇素数のとき、Mp が素数となるための必要十分条件は、S0 = 4, Sn = Sn−12 − 2 (n ≧ 1) と定義したときに Sp−2 が Mp で割り切れることである』
簡単な例を示しましょう。
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): .... if lucas_lehmer_test(n): print(n) .... 3 5 7 13 17 19 31 61 89 107 127 521 607
半素数は 2 つの素数の積なので、n / 2 以下の素数を求めて、それを掛け合わせていくことで求めることができます。また、エラトステネスの篩のアルゴリズムを改造して、素数を求めながら半素数を求めることもできます。この場合、新しく求めた素数で割り算した回数を記憶しておいて、それが 2 ならば半素数であることがわかります。どちらが速いか実際にプログラムを作って試してみましょう。
プログラムは次のようになります。
リスト : 半素数
# 半素数を求める
def semiprime(n):
ps = [0] * (n + 1)
semi = []
x = 2
while x <= n:
if ps[x] == 0:
# 素数
for y in range(x + x, n + 1, x):
z = y
while ps[y] < 3 and z % x == 0:
ps[y] += 1
z /= x
elif ps[x] == 2:
# 半素数
semi.append(x)
x += 1
return semi
# 半素数を求める
def semiprime1(n):
ps = sieve(n // 2)
semi = []
i = 0
while ps[i] * ps[i] <= n:
for j in range(i, len(ps)):
x = ps[i] * ps[j]
if x > n: break
semi.append(x)
i += 1
return sorted(semi)
def test6(func, n):
s = time.time()
print(func(n))
print('{:.3f}'.format(time.time() - s))
関数 semiprime() はエラトステネスの篩を改造したものです。配列 ps を 0 に初期化します。ここに素数で割った回数を記憶します。これがふるいの役目を果たしていて、ps[x] の値が 0 ならば x は素数であることがわかります。そうであれば、x の倍数を x で割り算していきます。このとき、割り算の回数が 2 を超えたならば、その数は半素数ではないことがわかるので、ここで割り算を終了します。ps[x] が 2 ならば半素数なので、x を append() で配列 semi に追加します。
関数 semiprime1() は sieve() で n / 2 以下の素数を配列 ps に求めます。あとは、ps の先頭から順番に素数を取り出し、それ以降の要素と掛け算して、その値 x を配列 semi に追加していくだけです。x が n より大きくなったら、次の素数を取り出して掛け算します。この場合、semi の要素は昇順に並ばないので、最後に sorted() で semi をソートします。
実行結果は次のようになりました。
>>>> semiprime(100) [4, 6, 9, 10, 14, 15, 21, 22, 25, 26, 33, 34, 35, 38, 39, 46, 49, 51, 55, 57, 58, 62, 65, 69, 74, 77, 82, 85, 86, 87, 91, 93, 94, 95] >>>> semiprime1(100) [4, 6, 9, 10, 14, 15, 21, 22, 25, 26, 33, 34, 35, 38, 39, 46, 49, 51, 55, 57, 58, 62, 65, 69, 74, 77, 82, 85, 86, 87, 91, 93, 94, 95] >>>> test6(semiprime, 10000000) 1904324 2.861 >>>> test6(semiprime1, 10000000) 1904324 0.434
semiprime1() のほうが速くなりました。なお、ソートしないでよければ、semiprime1() はもう少し速くなります。また、半素数の個数をカウントするだけでよければ、ps[i] * ps[j] が n 以下になる j の位置を ps の後ろから求めると、さらに速くすることができます。興味のある方はいろいろ試してみてください。
素数の数が無限にあることは古くから知られていて、ユークリッドが「背理法」を用いて証明しています。背理法はある命題が真であることを証明するために、それが偽であることを仮定して、そこから矛盾を導くことで命題が真であることを証明する方法です。ここでは、素数の数が無限にあることを証明するわけですから、それを否定する命題、つまり「素数の数は有限個しかない」と仮定して、それが矛盾することを示せばいいわけです。
まず最初に「素数の数が有限個 (n 個) しかない」と仮定し、その最大値を p とします。2 から p までの素数を乗算して 1 を足した整数を a としましょう。a は 2 から p までの素数で割ると、割り切れずに 1 あまります。最初の仮定により a は素数ではないので素因数分解できるはずですが、2 から p までの素数では割り切れないので、これ以外の素数が存在することになり、最初の仮定と矛盾します。つまり、最初の仮定「素数の数は有限個しかない」が間違いだったわけで、素数は無数にあることが証明されました。
ところで、ユークリッドの証明でよくある誤解に「2 から p までの素数を乗算した値に 1 を足すことで素数を生成できる」というものがあります。これは間違いで、a = 2 * 3 * ... * p + 1 とすると、p から √a までの間に素数が存在する場合、a を割り切る素数があるかもしれません。実際に試してみましょう。
2 + 1 = 3 2 * 3 + 1 = 7 2 * 3 * 5 + 1 = 31 2 * 3 * 5 * 7 + 1 = 211 2 * 3 * 5 * 7 * 11 + 1 = 2311 2 * 3 * 5 * 7 * 11 * 13 + 1 = 30031 = 59 * 509 2 * 3 * 5 * 7 * 11 * 13 * 17 + 1 = 510511 = 19 * 97 * 277
30031 や 510511 は素因数分解できるので素数ではありません。このように、2 から p までの素数を乗算して 1 を足しても、それが素数になるとは限らないわけです。
#
# prime.py : 素数
#
# Copyright (C) 2015-2022 Makoto Hiroi
#
import time
# n 以下の素数を求める
def solver0(n):
primes = [2]
for x in range(3, n + 1, 2):
for p in primes:
if p * p > x:
primes.append(x)
break
if x % p == 0: break
return primes
def solver1(n):
primes = [2, 3]
x = 5
while x <= n:
for p in primes:
if p * p > x:
primes.append(x)
break
if x % p == 0: break
if x % 6 == 5:
x += 2
else:
x += 4
return primes
# エラトステネスの篩
def sieve(n):
primes = [True] * (n // 2)
ps = [2]
x = 3
while x * x <= n:
y = (x - 3) // 2
if primes[y]:
ps.append(x)
y += x
while y < len(primes):
primes[y] = False
y += x
x += 2
while x <= n :
if primes[(x - 3) // 2]: ps.append(x)
x += 2
return ps
def sieve2(n):
ps = [2]
xs = list(range(3, n + 1, 2))
while xs[0] * xs[0] <= n:
m = xs[0]
ps.append(m)
xs = [x for x in xs if x % m != 0]
return ps + xs
def test1(func, n):
s = time.time()
ps = func(n)
print(len(ps), ps[-1])
print('{:.3f}'.format(time.time() - s))
# 素因数分解
def factor_sub(n, m):
c = 0
while n % m == 0:
c += 1
n //= m
return c, n
def factorization(n):
buff = []
c, m = factor_sub(n, 2)
if c > 0: buff.append((2, c))
c, m = factor_sub(m, 3)
if c > 0: buff.append((3, c))
x = 5
while m >= x * x:
c, m = factor_sub(m, x)
if c > 0: buff.append((x, c))
if x % 6 == 5:
x += 2
else:
x += 4
if m > 1: buff.append((m, 1))
return buff
# フィボナッチ数
def fibo(n):
a1, a2 = 0, 1;
for _ in range(n):
a1, a2 = a2, a1 + a2
return a1
# フィボナッチ素数
def prime_fibo(n):
cnt = 0
i = 2
while True:
x = fibo(i)
if x > n: break
xs = factorization(x)
if len(xs) == 1 and xs[0][1] == 1:
cnt += 1
print(i, ":", x)
i += 1
print(cnt)
# 素数ジェネレータ
def make_primes(n):
primes = [True] * (n // 2)
yield 2
x = 3
while x * x <= n:
y = (x - 3) // 2
if primes[y]:
yield x
y += x
while y < len(primes):
primes[y] = False
y += x
x += 2
while x <= n :
if primes[(x - 3) // 2]: yield x
x += 2
# 双子素数
def twin_prime(n):
p = 1
q = 0
count = 0
for x in make_primes(n):
if x - p == 2:
count += 1
q = p
p = x
print(count)
print(q, q + 2)
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 semiprime(n):
ps = [0] * (n + 1)
semi = []
x = 2
while x <= n:
if ps[x] == 0:
# 素数
for y in range(x + x, n + 1, x):
z = y
while ps[y] < 3 and z % x == 0:
ps[y] += 1
z /= x
elif ps[x] == 2:
# 半素数
semi.append(x)
x += 1
return semi
# 半素数を求める
def semiprime1(n):
ps = sieve(n // 2)
semi = []
i = 0
while ps[i] * ps[i] <= n:
for j in range(i, len(ps)):
x = ps[i] * ps[j]
if x > n: break
semi.append(x)
i += 1
return sorted(semi)
def test6(func, n):
s = time.time()
print(len(func(n)))
print('{:.3f}'.format(time.time() - s))