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

それでは [Haskell で「エラトステネスの篩」その１](https://qiita.com/little_Haskeller/items/e73d54835b5772a55e81) の続きです。

コード再掲
==========

``` {.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$ と表記していきます。[^1]

``` {.haskell}
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` と同じです。

``` {.haskell}
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)$ ということですね。