Skip to content

Instantly share code, notes, and snippets.

@minoki

minoki/result.txt

Created April 21, 2019 10:45
Show Gist options
  • Save minoki/333fb9c2d82e2279f5bfda54b51b9279 to your computer and use it in GitHub Desktop.
Save minoki/333fb9c2d82e2279f5bfda54b51b9279 to your computer and use it in GitHub Desktop.
Download ZIP
Haskellでもエラトステネスの篩がしたい！
Raw
result.txt
benchmarked trial division (2)/1000
time 337.3 μs (333.0 μs .. 342.5 μs)
0.996 R² (0.987 R² .. 0.999 R²)
mean 344.1 μs (339.5 μs .. 353.6 μs)
std dev 22.82 μs (7.337 μs .. 41.88 μs)
variance introduced by outliers: 42% (moderately inflated)
benchmarked trial division (2)/10000
time 22.05 ms (21.57 ms .. 22.77 ms)
0.996 R² (0.989 R² .. 0.999 R²)
mean 21.55 ms (21.31 ms .. 21.92 ms)
std dev 757.0 μs (481.0 μs .. 1.056 ms)
benchmarked trial division (6)/1000
time 327.0 μs (322.3 μs .. 332.3 μs)
0.995 R² (0.987 R² .. 0.999 R²)
mean 329.0 μs (325.8 μs .. 337.3 μs)
std dev 15.77 μs (7.806 μs .. 29.09 μs)
variance introduced by outliers: 27% (moderately inflated)
benchmarked trial division (6)/10000
time 21.65 ms (20.90 ms .. 22.34 ms)
0.996 R² (0.993 R² .. 0.999 R²)
mean 21.40 ms (21.08 ms .. 22.28 ms)
std dev 1.150 ms (526.8 μs .. 2.231 ms)
variance introduced by outliers: 18% (moderately inflated)
benchmarked eratos sieve (2)/1000
time 4.049 μs (3.965 μs .. 4.148 μs)
0.989 R² (0.976 R² .. 0.997 R²)
mean 4.388 μs (4.317 μs .. 4.493 μs)
std dev 295.4 ns (222.7 ns .. 399.9 ns)
variance introduced by outliers: 41% (moderately inflated)
benchmarked eratos sieve (2)/10000
time 44.48 μs (43.14 μs .. 45.50 μs)
0.993 R² (0.984 R² .. 0.998 R²)
mean 48.08 μs (46.94 μs .. 50.15 μs)
std dev 5.083 μs (3.207 μs .. 7.514 μs)
variance introduced by outliers: 66% (severely inflated)
benchmarked eratos sieve (2)/10^5
time 447.0 μs (440.7 μs .. 453.0 μs)
0.998 R² (0.997 R² .. 0.999 R²)
mean 455.4 μs (451.2 μs .. 465.0 μs)
std dev 19.39 μs (10.55 μs .. 31.05 μs)
variance introduced by outliers: 22% (moderately inflated)
benchmarked eratos sieve (2)/10^6
time 4.733 ms (4.683 ms .. 4.776 ms)
0.999 R² (0.997 R² .. 1.000 R²)
mean 4.779 ms (4.744 ms .. 4.851 ms)
std dev 147.6 μs (78.98 μs .. 247.0 μs)
variance introduced by outliers: 13% (moderately inflated)
benchmarked eratos sieve (2, bit)/1000
time 5.299 μs (5.105 μs .. 5.583 μs)
0.993 R² (0.987 R² .. 0.999 R²)
mean 5.378 μs (5.299 μs .. 5.545 μs)
std dev 378.8 ns (237.9 ns .. 623.4 ns)
variance introduced by outliers: 46% (moderately inflated)
benchmarked eratos sieve (2, bit)/10000
time 56.15 μs (53.29 μs .. 59.02 μs)
0.990 R² (0.985 R² .. 0.998 R²)
mean 54.92 μs (54.19 μs .. 55.87 μs)
std dev 2.879 μs (2.161 μs .. 3.596 μs)
variance introduced by outliers: 31% (moderately inflated)
benchmarked eratos sieve (2, bit)/10^5
time 535.5 μs (503.0 μs .. 570.4 μs)
0.988 R² (0.980 R² .. 0.999 R²)
mean 512.6 μs (508.7 μs .. 521.9 μs)
std dev 19.15 μs (11.41 μs .. 35.28 μs)
variance introduced by outliers: 18% (moderately inflated)
benchmarked eratos sieve (2, bit)/10^6
time 5.051 ms (4.907 ms .. 5.201 ms)
0.996 R² (0.994 R² .. 0.998 R²)
mean 5.159 ms (5.083 ms .. 5.372 ms)
std dev 370.3 μs (187.6 μs .. 697.6 μs)
variance introduced by outliers: 43% (moderately inflated)
benchmarked eratos sieve (6)/1000
time 2.843 μs (2.752 μs .. 2.958 μs)
0.993 R² (0.989 R² .. 0.998 R²)
mean 2.772 μs (2.744 μs .. 2.815 μs)
std dev 114.3 ns (82.90 ns .. 153.1 ns)
variance introduced by outliers: 22% (moderately inflated)
benchmarked eratos sieve (6)/10000
time 31.84 μs (30.64 μs .. 33.09 μs)
0.995 R² (0.992 R² .. 0.999 R²)
mean 31.11 μs (30.81 μs .. 31.46 μs)
std dev 1.077 μs (864.2 ns .. 1.419 μs)
variance introduced by outliers: 18% (moderately inflated)
benchmarked eratos sieve (6)/10^5
time 331.1 μs (319.5 μs .. 343.5 μs)
0.995 R² (0.993 R² .. 0.999 R²)
mean 325.2 μs (323.1 μs .. 328.7 μs)
std dev 8.970 μs (6.684 μs .. 12.23 μs)
variance introduced by outliers: 11% (moderately inflated)
benchmarked eratos sieve (6)/10^6
time 3.276 ms (3.086 ms .. 3.603 ms)
0.973 R² (0.950 R² .. 1.000 R²)
mean 3.109 ms (3.078 ms .. 3.193 ms)
std dev 170.3 μs (49.93 μs .. 318.6 μs)
variance introduced by outliers: 33% (moderately inflated)
benchmarked eratos sieve (30)/1000
time 3.637 μs (3.515 μs .. 3.838 μs)
0.989 R² (0.977 R² .. 0.999 R²)
mean 3.571 μs (3.539 μs .. 3.650 μs)
std dev 154.7 ns (93.19 ns .. 265.3 ns)
variance introduced by outliers: 24% (moderately inflated)
benchmarked eratos sieve (30)/10000
time 26.29 μs (25.29 μs .. 27.77 μs)
0.990 R² (0.983 R² .. 0.999 R²)
mean 25.83 μs (25.57 μs .. 26.25 μs)
std dev 1.113 μs (766.9 ns .. 1.529 μs)
variance introduced by outliers: 23% (moderately inflated)
benchmarked eratos sieve (30)/10^5
time 305.1 μs (296.7 μs .. 318.2 μs)
0.993 R² (0.983 R² .. 0.999 R²)
mean 299.8 μs (297.6 μs .. 304.3 μs)
std dev 10.15 μs (5.584 μs .. 18.25 μs)
variance introduced by outliers: 15% (moderately inflated)
benchmarked eratos sieve (30)/10^6
time 2.965 ms (2.885 ms .. 3.058 ms)
0.995 R² (0.989 R² .. 0.999 R²)
mean 2.907 ms (2.884 ms .. 2.941 ms)
std dev 94.54 μs (67.93 μs .. 149.3 μs)
variance introduced by outliers: 14% (moderately inflated)
Raw
sieve.hs
{-# LANGUAGE BangPatterns #-}
import Gauge
import Data.Word
import Data.Bits
import Control.Monad
import qualified Data.Vector.Unboxed as U
import qualified Data.Vector.Unboxed.Mutable as UM
infixr 5 !:
(!:) :: a -> [a] -> [a]
(!x) !: xs = x : xs
-- 試し割りを使って、与えられた自然数以下の素数のリストを返す
-- 2の倍数はあらかじめ除外しておく
trialDivision2 :: Int -> [Int]
trialDivision2 !max = 2 : sieve [3,5..max]
where
sieve [] = []
sieve (p:xs) = p : sieve (filter (\x -> x `rem` p /= 0) xs)
-- 試し割りを使って、与えられた自然数以下の素数のリストを返す
-- 2または3の倍数はあらかじめ除外しておく
trialDivision6 :: Int -> [Int]
trialDivision6 !max = 2 : 3 : sieve list
where
-- 6n + 1, 6n + 5 の形の自然数だけを考える
list = 5 : listX 1
where listX !n | 6 * n + 1 > max = []
| otherwise = 6*n+1 !: 6*n+5 !: listX (n+1)
sieve [] = []
sieve (p:xs) = p : sieve (filter (\x -> x `rem` p /= 0) xs)
-- エラトステネスの篩で、与えられた自然数以下の素数のリストを返す
-- 2の倍数はあらかじめ除外しておく
eratosSieve2 :: Int -> [Int]
eratosSieve2 !max = 2 : U.ifoldr (\i isPrime xs -> if isPrime then (2 * i + 1) !: xs else xs) [] vec
where
vec = U.create $ do
-- vectorの最大のindex n に関して、 2 * n + 1 <= max < 2 * n + 3 となるように選ぶ
-- 2 * len - 1 <= max < 2 * len + 1
vec <- UM.replicate ((max - 1) `quot` 2 + 1) True
UM.write vec 0 False -- 1 is not a prime
let clear !p = forM_ [3*p,5*p..max] $ \n -> UM.write vec (n `quot` 2) False
factorBound = floor (sqrt (fromIntegral max) :: Double) -- max <= 2^53 を仮定する
loop !i | 2 * i + 1 > factorBound = return ()
| otherwise = do b <- UM.read vec i
when b $ clear (2 * i + 1)
loop (i + 1)
loop 1
-- vec ! i is True if (2 * i + 1) is prime
return vec
-- エラトステネスの篩で、与えられた自然数以下の素数のリストを返す
-- 2の倍数はあらかじめ除外しておく
eratosSieve2bit :: Int -> [Int]
eratosSieve2bit !max = 2 : U.ifoldr (\i bb xs -> let add !j xs | j >= 32 = xs
| testBit bb j = (64 * i + 2 * j + 1) !: add (j+1) xs
| otherwise = add (j+1) xs
in add 0 xs
) [] vec
where
vec = U.create $ do
-- vectorの最大のindex n に関して、 2 * n + 1 <= max < 2 * n + 3 となるように選ぶ
-- 2 * len - 1 <= max < 2 * len + 1
let maxIndex = (max - 1) `quot` 64
max' = 64 * maxIndex + 2 * 31 + 1
vec <- UM.replicate (maxIndex + 1) (0xffffffff :: Word32)
UM.modify vec (\x -> clearBit x 0) 0 -- 1 is not a prime
let clear !p = forM_ [3*p,5*p..max'] $ \n -> do
-- let (!q,!r) = n `quotRem` 64
let !q = n `shiftR` 6
!r = (n .&. 63) `shiftR` 1
UM.modify vec (\x -> clearBit x r) q
factorBound = floor (sqrt (fromIntegral max) :: Double) -- max <= 2^53 を仮定する
loop !i !n
-- n = 64 * i + 1
| n > factorBound = return ()
| otherwise = do
b <- UM.read vec i
forM_ [0..31] $ \j ->
when (testBit b j) $ clear (n + 2 * j)
loop (i + 1) (n + 64)
loop 0 1
-- testBit (vec ! i) j is True <=> (64 * i + 2 * j + 1) is prime
return vec
-- エラトステネスの篩で、与えられた自然数以下の素数のリストを返す
-- 2または3の倍数はあらかじめ除外しておく
eratosSieve6 :: Int -> [Int]
eratosSieve6 !max = 2 : 3 : U.ifoldr (\i bb xs -> case bb of
0x22 -> (6 * i + 1) !: (6 * i + 5) !: xs
0x02 -> (6 * i + 1) !: xs
0x20 -> (6 * i + 5) !: xs
_ -> xs
) [] vec
where
vec = U.create $ do
-- vectorの最大のindex n に関して、 6 * n + 1 <= max < 6 * n + 7 となればOK
let !maxIndex = (max - 1) `quot` 6
max' = 6 * maxIndex + 5
vec <- UM.replicate (maxIndex + 1) (0x22 :: Word8) -- 0b0100010
UM.write vec 0 0x20 -- 1 is not a prime, 5 is a prime
-- (6n+1, 6n+5)
let -- clear p : p の倍数を消す。ただし、 p の2倍、3倍はすでに消えているので、 p*(6n+1), p*(6n+5) の形の数を消す。
clear !p = let !(!p5q,!p5r) = (5*p) `quotRem` 6
!(!p7q,!p7r) = (7*p) `quotRem` 6
clearLoop !n !nq !nq'
-- n == 6*nq + p5r
-- n + 2 * p == 6*nq' + p7r
| n+2*p > max' = if n > max'
then return ()
else UM.modify vec (\x -> clearBit x p5r) nq
| otherwise = do
-- n = 6*nq+p5r
UM.modify vec (\x -> clearBit x p5r) nq
-- n+2*p = 6*nq'+p7r
UM.modify vec (\x -> clearBit x p7r) nq'
clearLoop (n + 6 * p) (nq + p) (nq' + p)
in clearLoop (5 * p) p5q p7q
factorBound = floor (sqrt (fromIntegral max) :: Double) -- max <= 2^53 を仮定する
-- loop i n : n == 6 * i + 1
loop !i !n | n > factorBound = return ()
| otherwise = do bb <- UM.read vec i
when (testBit bb 1) $ clear n -- 6 * i + 1
when (testBit bb 5) $ clear (n + 4) -- 6 * i + 5
loop (i + 1) (n + 6)
clear 5
loop 1 7
return vec
-- エラトステネスの篩で、与えられた自然数以下の素数のリストを返す
-- 2または3または5の倍数はあらかじめ除外しておく
eratosSieve30 :: Int -> [Int]
eratosSieve30 !max = 2 : 3 : 5 : U.ifoldr (\i bb xs -> let add j xs | testBit bb j = (30 * i + j) !: xs
| otherwise = xs
in add 1 $ add 7 $ add 11 $ add 13 $ add 17 $ add 19 $ add 23 $ add 29 xs
) [] vec
where
vec = U.create $ do
-- vectorの最大のindex n に関して、 30 * n + 1 <= max < 30 * n + 31 となればOK
let !maxIndex = (max - 1) `quot` 30
max' = 30 * maxIndex + 29
vec <- UM.replicate (maxIndex + 1) (0x208a2882 :: Word32) -- 0b10_0000_1000_1010_0010_1000_1000_0010
UM.write vec 0 0x208a2880 -- 0b10_0000_1000_1010_0010_1000_1000_0000, 7, 11
-- 30n+1, 30n+7, 30n+11, 30n+13, 30n+17, 30n+19, 30n+23, 30n+29
let clear !p = let !(!p7q,!p7r) = (7*p) `quotRem` 30
!(!p11q,!p11r) = (11*p) `quotRem` 30
!(!p13q,!p13r) = (13*p) `quotRem` 30
!(!p17q,!p17r) = (17*p) `quotRem` 30
!(!p19q,!p19r) = (19*p) `quotRem` 30
!(!p23q,!p23r) = (23*p) `quotRem` 30
!(!p29q,!p29r) = (29*p) `quotRem` 30
!(!p31q,!p31r) = (31*p) `quotRem` 30
!max'' = max - 24*p
clearLoop !n !n7q !n11q !n13q !n17q !n19q !n23q !n29q !n31q
-- n == (30 * k + 7) * p for some k
-- n == 30 * n7q + p7r
-- n + 4 * p == 6 * n11q + p11r
-- ...
-- n + 24 * p == 6 * n31q + p31r
| n > max'' = if n+12*p <= max'
then do UM.modify vec (\x -> clearBit x p7r) n7q
UM.modify vec (\x -> clearBit x p11r) n11q
UM.modify vec (\x -> clearBit x p13r) n13q
UM.modify vec (\x -> clearBit x p17r) n17q
when (n+12*p <= max') $ do
UM.modify vec (\x -> clearBit x p19r) n19q
when (n+16*p <= max') $ do
UM.modify vec (\x -> clearBit x p23r) n23q
when (n+22*p <= max') $ do
UM.modify vec (\x -> clearBit x p29r) n29q
else do when (n <= max') $ do
UM.modify vec (\x -> clearBit x p7r) n7q
when (n+4*p <= max') $ do
UM.modify vec (\x -> clearBit x p11r) n11q
when (n+6*p <= max') $ do
UM.modify vec (\x -> clearBit x p13r) n13q
when (n+10*p <= max') $ do
UM.modify vec (\x -> clearBit x p17r) n17q
| otherwise = do
-- n = 30*n7q+p7r
UM.modify vec (\x -> clearBit x p7r) n7q
-- n+4*p = 30*n11q+p11r
UM.modify vec (\x -> clearBit x p11r) n11q
UM.modify vec (\x -> clearBit x p13r) n13q
UM.modify vec (\x -> clearBit x p17r) n17q
UM.modify vec (\x -> clearBit x p19r) n19q
UM.modify vec (\x -> clearBit x p23r) n23q
UM.modify vec (\x -> clearBit x p29r) n29q
UM.modify vec (\x -> clearBit x p31r) n31q
clearLoop (n + 30 * p) (n7q + p) (n11q + p) (n13q + p) (n17q + p) (n19q + p) (n23q + p) (n29q + p) (n31q + p)
in clearLoop (7 * p) p7q p11q p13q p17q p19q p23q p29q p31q
factorBound = floor (sqrt (fromIntegral max) :: Double) -- max <= 2^53 を仮定する
loop !i !n | n > factorBound = return ()
| otherwise = do
-- n == 30 * i + 1
bb <- UM.read vec i
when (testBit bb 1) $ clear n -- 30 * i + 1
when (testBit bb 7) $ clear (n + 6) -- 30 * i + 7
when (testBit bb 11) $ clear (n + 10) -- 30 * i + 11
when (testBit bb 13) $ clear (n + 12) -- 30 * i + 13
when (testBit bb 17) $ clear (n + 16) -- 30 * i + 17
when (testBit bb 19) $ clear (n + 18) -- 30 * i + 19
when (testBit bb 23) $ clear (n + 22) -- 30 * i + 23
when (testBit bb 29) $ clear (n + 28) -- 30 * i + 29
loop (i + 1) (n + 30)
loop 0 1
return vec
main :: IO ()
main = defaultMain
[ bgroup "trial division (2)"
[ bench "1000" $ nf trialDivision2 1000
, bench "10000" $ nf trialDivision2 10000
-- , bench "10^5" $ nf trialDivision2 (10^5)
]
, bgroup "trial division (6)"
[ bench "1000" $ nf trialDivision6 1000
, bench "10000" $ nf trialDivision6 10000
-- , bench "10^5" $ nf trialDivision6 (10^5)
]
, bgroup "eratos sieve (2)"
[ bench "1000" $ nf eratosSieve2 1000
, bench "10000" $ nf eratosSieve2 10000
, bench "10^5" $ nf eratosSieve2 (10^5)
, bench "10^6" $ nf eratosSieve2 (10^6)
]
, bgroup "eratos sieve (2, bit)"
[ bench "1000" $ nf eratosSieve2bit 1000
, bench "10000" $ nf eratosSieve2bit 10000
, bench "10^5" $ nf eratosSieve2bit (10^5)
, bench "10^6" $ nf eratosSieve2bit (10^6)
]
, bgroup "eratos sieve (6)"
[ bench "1000" $ nf eratosSieve6 1000
, bench "10000" $ nf eratosSieve6 10000
, bench "10^5" $ nf eratosSieve6 (10^5)
, bench "10^6" $ nf eratosSieve6 (10^6)
]
, bgroup "eratos sieve (30)"
[ bench "1000" $ nf eratosSieve30 1000
, bench "10000" $ nf eratosSieve30 10000
, bench "10^5" $ nf eratosSieve30 (10^5)
, bench "10^6" $ nf eratosSieve30 (10^6)
]
]
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment