素数無限リスト

先日のHaskellのコードにて。二つのソースで何がそんなに速度を変える原因になっているのだろうか、と思ったのだけれども、僕の書いたほうのソースには欠陥があることが分かった。もう一回そのソースをさらしてみる。

is_prime :: Integer -> [Integer] -> Bool
is_prime x  = True
is_prime x (hd:tl) =
                if (mod x hd == 0) then False else is_prime x tl

nextprime :: Integer -> [Integer] -> Integer
nextprime num  = 2
nextprime num list = if is_prime num list then num 
                     else nextprime (num + 1) list

primelist :: Integer -> [Integer] -> [Integer]
primelist _  = 2:(primelist 3 [2])
primelist num list =
        (np:(primelist (np + 1) (np:list))) where
                np = nextprime num list

primes :: [Integer]
primes = primelist 1 

「すでに素数と分かっている数のリスト」を用いて、そのリストの要素でだけ割ってチェックすれば良い、という方針にてコーディングしてある。ところが、この「素数と分かっている数のリスト」というのの更新を行っているのだけれども、この部分で新しい素数を追加するときに、新しい素数をnpとしてnp:listという形にしてある。つまり、大きい素数ほどリストの最初の方においてあることになる。
基本的に小さい素数で割り切れる確率の方が、大きい素数で割り切れる確率よりも大きい。したがって、最初に小さい素数から割っていくというのが定石だと考えられる。そのため、以下のようにprimelistを書き直してみた。

primelist num list = 
	(np:(primelist (np + 1) (list ++ [np]))) where
		np = nextprime num list

リストの最後尾に追加するのはコストが高いのだけれども、そのあたりは遅延評価とかで何とかなるんじゃないかなとやってみたところ、通常のエラトステネスのふるいと同じ程度の速度が出るようになった。
で、さらにis_primeをちょっと書き直してみた。

is_prime x (hd:tl) =
	if (x < hd * hd) then True
	else
		if (mod x hd == 0) then False else is_prime x tl

これによって、不要な素数で割る作業をやめたところ、今度はこのソースのほうが圧倒的に早くなった。