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

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

コード再掲

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$ と表記していきます。¹

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 と実際の数は idxToNum と numToIdx で変換していくことや、 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$ の倍数はすべて消せます。

次に num と num の隣りの数の積を求めると、 $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)$ ということですね。 ↩