Haskell
素数
エラトステネスの篩

Haskell で「エラトステネスの篩」 その2

それでは Haskell で「エラトステネスの篩」その1 の続きです。


コード再掲

import Control.Monad

import qualified Data.Vector.Unboxed as U
import qualified Data.Vector.Unboxed.Mutable as UM

--
-- エラトステネスの篩 ( 2 の倍数をあらかじめ除去 )
--
primeList :: Int -> [Int]
primeList n = 2 : U.ifoldr f [] sieve
where
f i v xs = if v then idxToNum i : xs else xs

idxToNum i = 2 * i + 3
numToIdx v = div (v - 3) 2
lastIdx = numToIdx n

sieve = U.create $ do
mVec <- UM.replicate (lastIdx + 1) True
forM_ [0 .. numToIdx (floor $ sqrt $ fromIntegral n)] $ \i -> do
v <- UM.read mVec i
when v $ do
let num = idxToNum i
s = numToIdx (num * num)
forM_ [s, s + num .. lastIdx] $ \j -> do
UM.write mVec j False
return mVec

--
-- エラトステネスの篩 ( 2 または 3 の倍数をあらかじめ除去 )
--
primeList2 :: Int -> [Int]
primeList2 2 = [2]
primeList2 n = 2 : 3 : U.ifoldr f [] sieve
where
f i v xs = if v then num : xs else xs
where num = idxToNum i

idxToNum i = if (even i) then 6 * (div i 2) + 5 else 6 * (div i 2) + 7
numToIdx n = if (mod n 6 == 5) then (div n 6) * 2 else (div n 6) * 2 - 1
lastIdx
| mod6 < 1 = div n 6 * 2 - 2
| mod6 < 5 = div n 6 * 2 - 1
| otherwise = div n 6 * 2
where mod6 = mod n 6

sieve = U.create $ do
mVec <- UM.replicate (lastIdx + 1) True
forM_ [0 .. numToIdx (floor $ sqrt $ fromIntegral n)] $ \i -> do
v <- UM.read mVec i
when v $ do
let num = idxToNum i
s1 = numToIdx (num * num)
s2 = numToIdx (num * (idxToNum (i + 1)))
forM_ [s1, s1 + num * 2 .. lastIdx] $ \j -> do
UM.write mVec j False
forM_ [s2, s2 + num * 2 .. lastIdx] $ \j -> do
UM.write mVec j False
return mVec


"2 または 3 の倍数をあらかじめ除去したバージョン" の説明

前回同様、 primeList2 の中心となる sieve について書いていきます。

以下の説明では、「6 で割ると 1 余る数」を $m_1$ , 「6 で割ると 5 余る数」を $m_5$ と表記していきます。1

mVec <- UM.replicate (lastIdx + 1) True

mVec は「2 または 3 の倍数を取り除いた自然数 (5, 7, 11, 13, 17, 19 ...)」を想定しています。

"5, 7, 11, 13, 17, 19 ..." という数列は、よく見ると " ${m_5, m_1, m_5, m_1 ...}$ " というように、 $m_5$ と $m_1$ が交互に並んでい ます。ここが primeList とは違うところで、後の処理に影響してきます。

index と実際の数は idxToNumnumToIdx で変換していくことや、 mVec の中 の値の意味は primeList と同じです。

v <- UM.read mVec i

when v $ do
let num = idxToNum i
s1 = numToIdx (num * num)
s2 = numToIdx (num * (idxToNum (i + 1)))
forM_ [s1, s1 + num * 2 .. lastIdx] $ \j -> do
UM.write mVec j False
forM_ [s2, s2 + num * 2 .. lastIdx] $ \j -> do
UM.write mVec j False

v の値が True ならば num = idxToNum i は素数なので num の倍数を消して いきます。

mVec は $m_5$ と $m_1$ が交互に並んでいるため num の倍数が規則的には並んで いるわけではありません。しかし偶数 index は $m_5$ だけ、奇数 index は $m_1$ だけが並んでいるので、偶数 index と奇数 index を分けて考えれば primeList と同じ方法が使えます。

ここでひとつ確認しておきたいのが、 $m_5 \times m_5 = m_1$ , $m_1 \times m_1 = m_1$ , $m_1 \times m_5 = m_5$ という関係があるということです。

ですから primeList のように num の 2 乗から num * 2 跳びに num の倍数を消していくと $m_1$ の倍数はすべて消せます。

次に numnum の隣りの数の積を求めると、 $m_5 \times m_1$ または $m_1 \times m_5$ となるので、 $m_5$ の列で消すべき倍数の起点を見付けることができます。





  1. すなわち $m_1 ≡ 1 \ (\mathrm{mod} \ 6), m_5 ≡ 5 \ (\mathrm{mod} \ 6)$ ということですね。