書記が数学やるだけ#894 素数判定アルゴリズム
素数判定アルゴリズムについていくつか例を示す。
問題
説明
エラトステネスの篩 (Sieve of Eratosthenes) は，古代ギリシアの科学者エラトステネスが考案したとされる，指定された整数以下の全ての素数を発見するための単純なアルゴリズムである。素数2に対し2の倍数を消し，次の素数3に対し3の倍数を消し，を繰り返す構成となっている。
確率的素数判定については以前に扱った：
別の確率的素数判定として，ルジャンドル記号の性質を利用したソロベイ–シュトラッセン素数判定法がある。
AKS素数判定法は，与えられた自然数が素数であるかどうかを決定的多項式時間で判定できる世界初のアルゴリズムであり，AKS素数判定法は2002年8月6日に "PRIMES is in P" と題された論文で発表された。
解答
素数を2から順に愚直に調べる方法では，√k以下の素数で割り算をする必要があり，その合計の計算量はO(n^(3/2))で見積もられる。
エラトステネスの篩では，1ステップごとに約n/p個の数を消すことができることから，計算量はO(nloglogn)まで減少する。
ソロベイ–シュトラッセン素数判定法を用いて21の素数判定を実施する。まずa=2の時の手計算について等式が成り立たないことからEuler witnessesであり，合成数であることがわかる。
他のaについても計算した結果，12個のaのうち9個がEuler witnessesであった。一般に，任意の奇数の合成数nに対して全ての底のうち少なくとも1/2の底がEuler witnessであることが知られている。
import math
def jacobi_symbol(a, n):
"""
ヤコビ記号(a/n)を計算。
nは正の奇数。
"""
if n <= 0 or n % 2 == 0:
raise ValueError("n must be a positive odd integer")
if math.gcd(a, n) != 1:
return 0
a %= n
res = 1
while a != 0:
# 2を因数分解
while a % 2 == 0:
a //= 2
if n % 8 == 3 or n % 8 == 5:
res = -res
# 平方剰余の相互法則を適用
a, n = n, a
if a % 4 == 3 and n % 4 == 3:
res = -res
a %= n
if n == 1:
return res
else:
return 0
# 検証する合成数 n
n = 21
half_n_minus1 = (n - 1) // 2
# ヘッダーの表示
print(f"{'a':>2} {'gcd(a,n)':>8} {'a^10 mod 21':>12} {'Jacobi(a/21)':>14} {'Euler test':>12}")
print("-" * 60)
# 1からn-1までのaについて検証
for a in range(1, n):
g = math.gcd(a, n)
mod_exp = pow(a, half_n_minus1, n)
if g > 1:
jacobi = 0
euler_result = "gcd>1"
else:
jacobi = jacobi_symbol(a, n)
# オイラーの規準 a^((n-1)/2) === (a/n) (mod n) を確認
if mod_exp == jacobi % n:
euler_result = "liar"
else:
euler_result = "witness"
# 結果の表示
print(f"{a:>2} {g:>8} {mod_exp:>12} {jacobi:>14} {euler_result:>12}") a gcd(a,n) a^10 mod 21 Jacobi(a/21) Euler test
------------------------------------------------------------
1 1 1 1 liar
2 1 16 -1 witness
3 3 18 0 gcd>1
4 1 4 1 witness
5 1 16 1 witness
6 3 15 0 gcd>1
7 7 7 0 gcd>1
8 1 1 -1 witness
9 3 9 0 gcd>1
10 1 4 -1 witness
11 1 4 -1 witness
12 3 9 0 gcd>1
13 1 1 -1 witness
14 7 7 0 gcd>1
15 3 15 0 gcd>1
16 1 16 1 witness
17 1 4 1 witness
18 3 18 0 gcd>1
19 1 16 -1 witness
20 1 1 1 liarAKS素数判定法の鍵となるのは次の合同式である。
$$(X+a)^n \not\equiv X^n+a\mod n$$
この合同式はフェルマーの小定理の対偶を改良したものであり、厳密に定式化することにより、素数であることの必要十分条件を与えている。この合同式を適切に評価すれば、素数性を判定する決定性アルゴリズムを構築できる。しかし、実際には最悪の場合、n 個の係数を評価する必要があり、指数関数的な計算時間がかかる。そこで、ある小さな r を選び、多項式を X^r−1 で割った剰余で評価することで、評価すべき係数の数を大幅に削減することができる。
本記事のもくじはこちら：
