素数の列を返すプログラム、素数かどうか判定するプログラム

素数の列を返すプログラムを書いてみる。

def prime(n): # nまでの素数のリストを返す
    if n == 2:
        return [2]
    elif n > 2:
        array = [2] # 素数を集める
        for i in range(3, n + 1): # 3～nが素数かどうか
            for p in array:
                if i % p == 0: # 素数で割り切れたとき
                    break
            else:
                array.append(i)
        return array

2以上の整数nが素数であることの定義は、nが1とn以外のいずれの整数でも割り切れないことなので、
3以上の整数については、2～n-1の範囲で割れるか調べればよい。
のだが、それでは効率が悪い、実際は2～n-1の中の素数で割れるか調べるだけで十分。
（素数以外の整数（合成数という）は素数の積に分解できるので、素数で割れないならその素数を素因数にもつ合成数でも割れない）

なので、素数を集めておくリストを作っておき、素数判定のたびにリストに追加していき、そのリストを用いて次の数の素数判定を行っている。
素数の列を返して欲しいプログラムなのでどうせこのリストは必要になり、この方法が効率がよい。

だが、単にある数nが素数かどうか判定したいだけのときに、まずn-1までの数が素数かどうか調べるこの方法がいいのかどうかは知らん。
というか、ちょっと考えただけでもめっちゃ効率悪そう。
「9999までの素数をすべて列挙せよ」ならともかく、「9999が素数かどうか判定せよ」なんていう本来一瞬で解き終わる問題が、この方法だと途方もない時間がかかりそうだ。

アルゴリズムの発想はあまり変えず、微妙に書き方を変えてみる。

def ismyprime(n): # nが素数かどうか判定する
    if n == 2:
        return True
    elif n > 2:
        # 「nが素数である」とは、「nがn未満の素数で割り切れない」
        array = filter(ismyprime, range(2, n))
        for p in array:
            if n % p == 0:
                return False
        else:
            return True
    else:
        return False

def myprime(n): # nまでの素数のリストを返す
    return list( filter(ismyprime, range(2, n + 1)) )

先ほどのプログラムでは2重ループでarrayにどんどん追加していったが、
ここでは再帰を用いて宣言的に書き換えた。

このプログラムでも動くが、代償として実行速度がお察し状態になった。
myprime(70)を試すと一瞬で動作が終わるprime(70)と比べて遅く感じられる。
myprime(100)は待ち切れないので強制終了させた。

インタプリタだから遅いとか元のアルゴリズムが非効率とかいうより、
酷い再帰の書き方で実行速度がお察しになったのだろう。

同じアルゴリズムでHaskell版。

isprime :: Int -> Bool
isprime n
    | n <= 1     = False
    | n == 2     = True
    | otherwise  = foldr (\p acc -> if n `mod` p == 0 then False else acc) True (filter isprime [2..n - 1])

prime :: Int -> [Int]
prime n = filter isprime [2..n]

畳み込み関数を使うとコンパクトにまとまるが、慣れていないと読みにくい。しかも無駄な計算してる気がする。
原始的な再帰で書くとこうなる。

not_divisible :: Int -> [Int] -> Bool
not_divisible n [] = True
not_divisible n (x:xs)
    | n `mod` x == 0    =   False
    | otherwise         =   not_divisible n xs

isprime :: Int -> Bool
isprime n
    | n <= 1     = False
    | n == 2     = True
    | otherwise  = not_divisible n (filter isprime [2..n - 1])

部分適用とelem関数を上手く使うともっとわかりやすくなる。遅延評価の言語だから途中で無駄な計算は発生しないはず？

isprime2 :: Int -> Bool
isprime2 n
    | n <= 1 = False
    | n == 2 = True
    | otherwise  = not $ 0 `elem` (map (n `mod`) (filter isprime [2..n - 1]))

prime2 :: Int -> [Int]
prime2 n = filter isprime2 [2..n]

いずれのプログラムもprime 70くらいから遅い。これならfilter isprime [2..n - 1]とか使わずに[2..n - 1]をそのまま使っとけって感じだ。

というわけで、エラトステネスの篩ってどう実装すればいいんだっけとググる。

正しくはこうらしい。

prime :: [Int] -> [Int]
prime [] = []
prime (x:xs) = x : prime [y | y <- xs, mod y x /= 0]

なるほど、頭いいな。これなら先の酷い再帰アルゴリズムと違って、一瞬で動作が終わる。
Pythonならこう。こちらもさっきの酷い再帰アルゴリズムと違って、一瞬で動作が終わる。

def myprime2(n):
    array = list(range(2, n + 1))

    def primes(array):
        if not array:
            return []
        else:
            return [array[0]] + primes([x for x in array if x % array[0] != 0])

    return primes(array)

再帰を使わずに書くとこうか。

def myprime3(n):
    array = list(range(2, n + 1))
    answer = []

    while array:
        answer.append(array[0])
        array = [x for x in array if x % array[0] != 0]

    return answer

効率のよいアルゴリズムってむずかしい

＜補足＞
2番目の再帰でなぜこんなに遅くなったか？

多分だが、isprime(n)を定義するのに、isprime(2)～isprime(n-1)を使っていて、・・・（ア）
isprime(n-1)を求めるにはisprime(2)～isprime(n-2)を使っていて、
isprime(n-2)を求めるにはisprime(2)～isprime(n-3)を使っていて・・・となっているのだが、
求めたisprime(i)を共用できていないからではないか？

アの各isprimeの下にisprimeのツリーが広がって、その下にまたツリーが広がって、と
おぞましいことになっている気がする。
計算結果を共有できてればいいのだが、できてないからとんでもない時間がかかっている、ということだろうか。