「素数 (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.9) です。
リスト : 素数を求める (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.240 >>>> test1(solver0, 3000000) 216816 2999999 1.228
素数の個数は 216816 個、最大の素数は 2999999 で、実行時間は約 1.24 秒 (実行環境 : Ubuntu 22.04 LTS (WSL2), 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.812 >>>> test1(solver1, 3000000) 216816 2999999 0.797
実行時間は約 0.8 秒になりました。残念ながら、劇的には速くなりませんでした。
もう一つ、素数を求める簡単な方法を紹介しましょう。
最初に、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.162 >>>> test1(sieve, 3000000) 216816 2999999 0.121 >>>> test1(sieve, 3000000) 216816 2999999 0.085 >>>> test1(sieve, 3000000) 216816 2999999 0.086
PyPy は JIT コンパイラを使用しているので、初回の実行には時間がかかることがあります。2 回目で約 0.12 秒、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 と素因数分解できるので素数ではありません。なお、これはとても単純なアルゴリズムなので、大きな整数の素因数分解には適していません。巨大な合成数の素因数分解はとても難しい問題です。興味のある方は素因数分解について調べてみてください。
ところで、これまでのプログラムは最初に 2 と 3 の倍数を篩にかけていましたが、2, 3, 5 の倍数を篩にかけることを考えてみましょう。次の図を見てください。
:8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 -+---------------------------------------------------------------------------------------- 2:X X X X X X X X X X X X X X X 3: X X X X X X X X X X 5: X X X X X X
8 から 37 までの区間で調べる数は 11, 13, 17, 19, 23, 29, 31, 37 になります。素数 7 を先頭として差分をとると、[4, 2, 4, 2, 4, 6, 2, 6] になります。つまり、7 のあとはリストをグルグル回して要素を加算していけば、2, 3, 5 の倍数を省くことができます。これを wheel (車輪) factorization と呼びます。プログラムは次のようになります。
リスト : wheel factorization
def wheel_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
変数 wheel が車輪を表すリストです。実際に循環するのは wheel[3:] の部分になります。先頭位置 3 を変数 s にセットします。wheel の位置は変数 i で表していて、i が wheel の末尾に到達したら、i の値を s に戻します。これで車輪 (wheel) が機能します。
それでは実際に試してみましょう。
>>>> 2**64 + 1 18446744073709551617 >>>> s = time.time(); factorization(2**64 + 1); print(time.time() - s) [(274177, 1), (67280421310721, 1)] 0.076416015625 >>>> s = time.time(); wheel_factorization(2**64 + 1); print(time.time() - s) [(274177, 1), (67280421310721, 1)] 0.0694725513458252 >>>> 2**67 - 1 147573952589676412927 >>>> s = time.time(); factorization(2**67 - 1); print(time.time() - s) [(193707721, 1), (761838257287, 1)] 3.716808319091797 >>>> s = time.time(); wheel_factorization(2**67 - 1); print(time.time() - s) [(193707721, 1), (761838257287, 1)] 2.9785001277923584
wheel_factorization の方がちょっとだけ速くなりました。エラトステネスの篩に車輪を使うと、もう少し速くなるかもしれません。興味のある方は試してみてください。
大きな整数ではないので、素数のチェックは 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 「フィボナッチ素数 - Wikipedia」によると、『フィボナッチ素数が無限にあるかどうかは分かっていない。』 とのことです。
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 「双子素数 - Wikipedia」によると、『双子素数は無数に存在するかという問題、いわゆる「双子素数の予想」や「双子素数の問題」は、いまだに数学上の未解決問題である。無数に存在するだろう、とは、多くの数論学者が予想している。』 とのことです。
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 「メルセンヌ数 - Wikipedia」より引用します。
『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 の後ろから求めると、さらに速くすることができます。興味のある方はいろいろ試してみてください。
陳素数を求めるプログラムは、sieve() と factorization() を使うと簡単です。次のリストを見てください。
リスト : 陳素数
def is_chenprime(n):
xs = factorization(n)
if len(xs) == 1:
return xs[0][1] <= 2
else:
return len(xs) == 2 and xs[0][1] == 1 and xs[1][1] == 1
def chen_prime(n):
cs = []
for p in sieve(n):
if is_chenprime(p + 2):
cs.append((p, p + 2))
return cs
chen_prime() は sieve() で n 以下の素数を生成し、それを順番に取り出して変数 p にセットします。そして、p + 2 が陳素数の条件を満たしているか、is_chenprime() でチェックします。そうであれば、(p, p + 2) を cs にセットします。最後に cs を返します。
is_chenprime() では、引数 n を factorization() で素因数分解します。xs の長さが 1 であれば、xs[0][1] が 1 または 2 であれば陳素数の条件を満たします。そうでなれば、xs の長さが 2 で、xs[0][1] と xs[1][1] が 1 であれば、半素数になるので True を返します。
それでは実行してみましょう。
>>>> chen_prime(100) [(2, 4), (3, 5), (5, 7), (7, 9), (11, 13), (13, 15), (17, 19), (19, 21), (23, 25), (29, 31), (31, 33), (37, 39), (41, 43), (47, 49), (53, 55), (59, 61), (67, 69), (71, 73), (83, 85), (89, 91)] >>>> chen_prime(500) [(2, 4), (3, 5), (5, 7), (7, 9), (11, 13), (13, 15), (17, 19), (19, 21), (23, 25), (29, 31), (31, 33), (37, 39), (41, 43), (47, 49), (53, 55), (59, 61), (67, 69), (71, 73), (83, 85), (89, 91), (101, 103), (107, 109), (109, 111), (113, 115), (127, 129), (131, 133), (137, 139), (139, 141), (149, 151), (157, 159), (167, 169), (179, 181), (181, 183), (191, 193), (197, 199), (199, 201), (211, 213), (227, 229), (233, 235), (239, 241), (251, 253), (257, 259), (263, 265), (269, 271), (281, 283), (293, 295), (307, 309), (311, 313), (317, 319), (337, 339), (347, 349), (353, 355), (359, 361), (379, 381), (389, 391), (401, 403), (409, 411), (419, 421), (431, 433), (443, 445), (449, 451), (461, 463), (467, 469), (479, 481), (487, 489), (491, 493), (499, 501)]
陳素数は双子素数と似ていますが、その個数は双子素数よりも多いようです。「双子素数は無数に存在するか」という問題はまだ未解決ですが、陳素数は無数に存在することが陳景潤氏によって証明されているそうです。
素数の数が無限にあることは古くから知られていて、ユークリッドが「背理法」を用いて証明しています。背理法はある命題が真であることを証明するために、それが偽であることを仮定して、そこから矛盾を導くことで命題が真であることを証明する方法です。ここでは、素数の数が無限にあることを証明するわけですから、それを否定する命題、つまり「素数の数は有限個しかない」と仮定して、それが矛盾することを示せばいいわけです。
まず最初に「素数の数が有限個 (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-2023 Makoto Hiroi
#
import time, math, random
# 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 wheel_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 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))
# 陳素数
def is_chenprime(n):
xs = factorization(n)
if len(xs) == 1:
return xs[0][1] <= 2
else:
return len(xs) == 2 and xs[0][1] == 1 and xs[1][1] == 1
def chen_prime(n):
cs = []
for p in sieve(n):
if is_chenprime(p + 2):
cs.append((p, p + 2))
return cs