素数列挙プログラムのスピードアップ

以下のようなプログラムがあった。

このプログラムだと，10**8 以下の素数を探索するのに 18 日くらいかかる。

import time

primes_list = []
upper_lim   = 10**3
start_time  = time.time()

for integer in range(2, upper_lim + 1, 1):
    if len(primes_list) < 1 :
        primes_list.append(integer)
    else:
        is_divisible = False
        for prime in primes_list:
            if integer % prime == 0:
                is_divisible = True
                break
        if not is_divisible:
            primes_list.append(integer)

#print(primes_list)
print(len(primes_list), "prime numbers foud.")

elapsed_time = time.time() - start_time # 経過時間 = 現在時刻 - 開始時刻
print ("run time: {0}".format(elapsed_time) + " seconds")

偶数の素数は 2 だけ。それ以外は奇数。range(2, upper_lim + 1, 1) でぶん回すのは少なくとも2倍の時間が掛かる。

for ループの中に len(primes_list) < 1 を毎回判定する if ブロックは無駄。

is_divisible という変数を使っているが，後で not is_divisible として使っているのなら，is_prime という変数を作ったほうがわかりやすい。

ある数が，primes_list の中にある素数で割りきれるか primes_list の全体を走査しているが，その数の平方根以下の要素まで調べればよい。そこまでで割り切れなければ素数。

Python の癖であるが，関数にしたほうが速い。

＋＋＋＋＋＋＋＋＋＋

他の人からアドバイスがあって，書き直されたプログラムが3とおりあった。

関数化したこと
2 を特別に扱い，range(3, upper_limit + 1, 2) で探索すること。
all(integer % prime != 0 for prime in primes_list) を使うこと。
しかし，これはあまりよい提案ではないことがわかる。
素数かどうかの走査を早めに終えるためにはこの方法はよくない。

def primes_up_to_(upper_limit):
    if upper_limit < 2:
        return []
    primes_list = [2]
    for integer in range(3, upper_limit + 1, 2):
        if all(integer % prime != 0 for prime in primes_list):
            primes_list.append(integer)
    return primes_list

＋＋＋＋＋＋＋＋＋＋

もう一つのバージョン

for integer in range(2, upper_lim + 1, 1):
    if len(primes_list) < 1 :
        primes_list.append(integer)
    else:
        is_divisible = False
        for prime in primes_list:
            # if prime >= math.sqrt(integer): # このやり方は実によい！！
            if prime > math.sqrt(integer): # 但し，> でないと間違い
                break
            elif integer % prime == 0:
                is_divisible = True
                break
        if not is_divisible:
            primes_list.append(integer)

＋＋＋＋＋＋＋＋＋＋

もう一つのバージョンだが，何をやっているか理解しておらず，処理の重複のため長時間かかることになってしまった。

追加したものは if all() と重複する余分なもの。百害あって一利なし。それに，前のバージョンの prime > math.sqrt(integer) がなくなってしまった。

def primes_up_to_3(upper_limit): # 関数化
    if upper_limit < 2:
        return []
    primes_list = [2]
    for integer in range(3, upper_limit + 1, 2):
        temp_primes_list = []                                       # 追加
        for prime in primes_list:                                   # 追加
            if prime <= math.sqrt(integer):                         # 追加
                temp_primes_list.append(prime)                      # 追加
        if all(integer % prime != 0 for prime in temp_primes_list):
            primes_list.append(integer)
    return primes_list

＋＋＋＋＋＋＋＋＋＋

私が考えていたバージョン。

def gen_primes2(upper_lim = 10000):
    primes_list = [2]
    for integer in range(3, upper_lim + 1, 2):
        is_prime = True
        upper_lim_sqrt = int(integer**0.5)
        for prime in primes_list:
            if prime > upper_lim_sqrt: # 時短にかなり有効
                break
            elif integer % prime == 0:
                is_prime = False
                break
        if is_prime:
            primes_list.append(integer)
    return primes_list

これだと，最初のバージョンでは 10**8 までの検索に 18 日かかっていたであろうものが，10 分で終わる。

＋＋＋＋＋＋＋＋＋＋

しかし，以下の Julia で書いたプログラムだと，さらにその 10 倍速い。

function gen_primes2(upper_lim = 10000)
    primes_list = [2]
    for integer = 3:2:upper_lim
        upper_lim_sqrt = ceil(sqrt(integer))
        is_prime = true
        for prime in primes_list
            if prime >  upper_lim_sqrt
                break
            elseif integer % prime == 0
                is_prime = false
                break
            end
        end
        if is_prime
            append!(primes_list, integer)
        end
    end
    println("$(length(primes_list)) prime numbers foud.")
    # println(primes_list)
end

@time gen_primes2(10^8)
# $ julia gen_primes2.jl
# 5761455 prime numbers foud.
# 54.732656 seconds (35 allocations: 65.001 MiB, 0.00% gc time)

10**8 までの検索が 1分で終わる。

つまり，元のプログラムの2万6千倍速いのだ。