Szarny.io

There should be one-- and preferably only one --obvious way to do it.

Pythonで実装する素数判定アルゴリズムと実行効率の比較

はじめに

Pythonで5種の素数判定アルゴリズムを実装します.

後で効率を比較するために,各関数は  [2, limit] の範囲において,各数が素数であるかを計算します.そして,その範囲における素数の総数を戻り値として設定します.

ごり押し判定法

説明

ある数 k素数かどうか判定するために, [2, k-1]の範囲において約数が存在するかどうか調べます.
もし,その範囲内に約数が存在しないなら,素数であると判断します.

ソースコード

def normal(limit):
    n = 0
    for k in range(2, limit+1):
        factor = 0
        
        for divisor in range(2, k):
            # 約数が存在するか?
            if k % divisor == 0:
                factor += 1
                
        # 約数が存在しないなら素数
        if factor == 0:
            n += 1
            
    return n

ちょっと改善した判定法

説明

ある数 k について, [\lceil\frac{k+1}{2}\ \rceil, k-1]の範囲に約数が存在しないことは明らかです.
そのため,k素数かどうか判定するために, [2, \lfloor\frac{k}{2}\ \rfloor]の範囲において約数が存在するかどうか調べます.
もし,その範囲内に約数が存在しないなら,素数であると判断します.

加えて,2以外の偶数は確実に素数ではないため,スキップします.

ソースコード

def improve(limit):
    n = 0
    for k in range(2, limit+1):
        factor = 0
        
        # 2以外の偶数は素数ではないので無視する
        if k % 2 == 0 and k != 2:
            continue
        
        # 繰り返しの上限を半分にする
        for divisor in range(2, k//2):
            if k % divisor == 0:
                factor += 1
                
        if factor == 0:
            n += 1
            
    return n

素数の性質を利用した判定法

説明

合成数 xp <= \sqrt{x} を満たす素因子 p をもつ

素数判定 | アルゴリズムとデータ構造 | Aizu Online Judge

という性質を利用します.

つまり,ある数 k素数かどうか判定するために, [2, \sqrt{k}]の範囲において約数が存在するかどうか調べれば十分であるということです.
もし,その範囲内に約数が存在しないなら,素数であると判断します.

前項で実装した,偶数をスキップするアルゴリズムも採用します.

ソースコード

def using_sqrt(limit):
    n = 0
    for k in range(2, limit+1):
        factor = 0
        
        # 2以外の偶数は素数ではないので無視する
        if k % 2 == 0 and k != 2:
            continue
        
        # 繰り返しの上限を対象の平方根にする
        for divisor in range(2, math.floor(math.sqrt(k))+1):
            if k % divisor == 0:
                factor += 1
                
        if factor == 0:
            n += 1
            
    return n

フェルマーの小定理を利用した判定法

説明

ap が互いに素な自然数のとき
a^{p-1} \equiv 1 \mod p

フェルマーの小定理の証明と例題 | 高校数学の美しい物語

というフェルマーの小定理を利用して,素数判定を行います.

非常に高速ですが,確率的素数判定法であるため,素数でない数(擬素数)を素数であると判断することがあります.

ソースコード

def fermat(limit):
    n = 0
    for k in range(2, limit+1):
        # 2以外の偶数は素数ではないので無視する
        if k % 2 == 0 and k != 2:
            continue
            
        if pow(2, k-1, k) == 1:
            n += 1
            
    return n

エラトステネスの篩

説明

  1.  [2, limit-1]の範囲に含まれる整数を列挙した集合Aを用意する
  2. Aの中で最小の数 p素数とし,素数リストPに追加する
  3. Aの中で,pの倍数である要素を削除する
  4. p > \sqrt{limit}となった時点で終了し,残ったAの全要素をPに追加する

以上のようにして,素数を篩い分けるアルゴリズムです.

ソースコード

def eratosthenes(limit):
    A = [i for i in range(2, limit+1)]
    P = []
    time = 0
    
    while True:
        prime = min(A)
        
        if prime > math.sqrt(limit):
            break
            
        P.append(prime)
            
        i = 0
        while i < len(A):
            if A[i] % prime == 0:
                A.pop(i)
                continue
            i += 1
            
    for a in A:
        P.append(a)
            
    return len(P)

実行時間の比較

説明

以上,5つのアルゴリズムを,上限を変えながらそれぞれ実行し,実行時間の比較を行います.
今回は上限を 1000, 5000, 10000, 100000 と設定しました.

実験用のソースコード

# 上限
limits = [1000, 5000, 10000, 100000]

# 適用する関数名
labels = ["Eratosthenes", "Fermat", "SQRT", "Improve", "Normal"]

# 関数本体
functions = [normal, improve, using_sqrt, fermat, eratosthenes]

for answer, limit in zip(answers, limits):
    records = []
    
    for t_f, f in enumerate(functions):
        start = time.time()
        result = f(limit)
        end = time.time()
        records.append(round(end - start, 5))
    
    records = list(reversed(records))
    
    plt.barh(range(5), records)
    plt.yticks(range(5), (labels))
    plt.title("limit is " + str(limit))
    plt.ylabel("Algorithm")
    plt.xlabel("Time(sec)")
    plt.axis(xmax=max(records)*1.2)
    
    for x,y in zip(records, [0,1,2,3,4]):
        plt.text(x*1.01, y, records[y])
    
    plt.show()

実行結果

結果はそれぞれ以下のようになりました.
上限が小さいときはそれほど差はありませんが,大きくなるにつれて差がどんどん広がっていきました.

グラフからわかる通り,実行時間が最も短いのはフェルマーの小定理を利用したアルゴリズムですが,やはり結果は不正確です.
正確かつ高速な判定が行いたいときは,平方根を上限として判定を行うアルゴリズムが最適かと思います.

f:id:Szarny:20170921231446p:plain
f:id:Szarny:20170921231457p:plain
f:id:Szarny:20170921231504p:plain
f:id:Szarny:20170921231940p:plain