素数列生成。
前回のエントリで、ProjectEular 10があまりに遅かったので素数列生成について色々考えてみた。
結構速くなったのだけど、多分それでもまだまだ遅いカモ。
調べれば色々アルゴリズムあるのだろうけど、うちの弱い頭だと論文とか読めないので独力で頑張ってみた。
そんなこんな。
書いている間になんどもやる気を失って、2〜3週間かけて書いたのですが、当然よい文章になるわけもなく、超駄文です。
一応、色々頭悩ませて書いたのでナニカの足しになればいいな、と思って駄文のままPOST。最速(?)な素数列生成はエントリの一番最後に。
あと、小ネタ、というか完全ネタですが、dcスクリプトによる最大最小除去平均を求めてます。
ナニカの参考になれば。
まず。
素数生成について。
基本的にはエラトステネスのふるい使います。というかそれ以外方法知りません。
で、個人的な欲望としては、1) n以下の素数列、2) n個の素数、二種をちゃんと出せるものが欲しい。
たとえばC言語なんかで(筆者が)実装すると、前者はn個のbool配列かなにか用意してふるいをそのまま使うけど、後者はn個のintなりlongな配列を用意して、素数をぶちこんでいってふるいをかける、みたいに方法ががらっと変わる。
それはいやだなぁ、ということで、アルゴリズムは同じで、以下に用意した関数を定義する。
takePrimes :: (Integral a) => a -> [a] underPrimes :: (Integral a) => a -> [a]
の二つ。*1。
で、これで速さを競います。
速度の比較はtakeの方は100個、200個、...、10000個の、underの方は1000以下、2000以下、...、100000以下の素数をそれぞれ生成させてグラフで比較します。
あ、ちげ。underは「未満」だ。いけない、いけない。全然意識しないで書いてた。
求めた素数列は適当にsumして評価を走らせます。
値の方は誤差を減らすため複数回走らせて、そのうち最大と最小を抜いた値の平均を使います。
測定用スクリプトはまいどおなじみシェル+dcスクリプトという変態仕様。
dcスクリプト書きつつ、RubyやHaskell...せめてシェルスクリプトで書いた方がどんだけ楽か、思ったけど、まぁ、気にしない。
測定用シェルスクリプト(test.sh)
#!/bin/bash prog=$1 init=$2 finite=$3 step=$4 dat=$(basename ${prog}).dc rm -f ${dat} echo "10k" >> ${dat} echo "[dsm]sM [dsn]sN [dlm<M dln>N Ss li1-si li0<A]sA [LsLs+Ss li1-si li0<B]sB" >> ${dat} for ((i = ${init}; i < ${finite}; i += ${step})); do # without dc # echo -n "${i} " >> ${dat} # /usr/bin/time -p ./timeout.rb 10 ./${prog} ${i} 2>&1 1>/dev/null | sed "s/^real *//p;d" >> ${dat} # with dc echo "[${i} ]n" >> ${dat} for ((j=0;j<7;j++)); do /usr/bin/time -p ./timeout.rb 10 ./${prog} ${i} 2>&1 1>/dev/null | sed "s/^real *//p;d" >> ${dat} done echo "zsz zsi dsm dsn lAx lmd0r-Ss lnd0r-Ss lz2+1-si lBx Lslz2-/p" >> ${dat} done
問題となるdcスクリプトは上記の
[dsm]sM [dsn]sN [dlm<M dln>N Ss li1-si li0<A]sA [LsLs+Ss li1-si li0<B]sB zsz zsi dsm dsn lAx lmd0r-Ss lnd0r-Ss lz2+1-si lBx Lslz2-/p
の二つ。
いらないと思うけど一応説明いれると、使うレジスタがマクロ用にM,N,A,B、値保持用にm,n,s,i,zを使用。
M,Nはそれぞれ最大/最小を記録するためのマクロ。大小比較でのマクロ実行に即値が使えないのでレジスタに。
m,nがそのM,Nで使う最大/最小を保持するレジスタ。
Aはメインスタックからsスタックへ値を移しつつも最大/最小を見つけるループ。
Bはそのsスタックの平均を求めるループ。
iはA,Bで使うiteratorで、zは要素数。メインスタックはzコマンドで要素数が拾えるのだけどレジスタスタックは拾えないので保持。
一行目が共通部分で、マクロの定義。
M: [dsm] d ... dup sm ... store to register m N: [dsn] d ... dup sn ... store to register n A: [dlm<M dln>N Ss li1-si li0<A] // 最大値の更新 d ... dup lm ... load from register m <M ... if (lm) < (dup) then call M // 最小値の更新 d ... dup ln ... load from register n >N ... if (lm) > (dup) then call N // レジスタsへ待避 Ss ... push to register s // レジスタiのデクリメント li ... load from register i 1- ... decrement si ... store to register i // next li ... load from register i 0 ... push 0 <A ... if (0) < (li) then recursive call A B: [LsLs+Ss li1-si li0<B] // 加算 Ls ... pop from register s Ls ... pop from register s + ... add Ss ... push to register s // レジスタiのデクリメント li ... load from register i 1- ... decrement si ... store to register i // next li ... load from register i 0 ... push 0 <B ... if (0) < (li) then recursive call B
二行目のが実際平均を出すスクリプト
// 要素数の記憶 z ... push stack depth sz ... store to register z // レジスタiの初期化 z ... push stack depth si ... store to register i // レジスタmの初期化 d ... dup sm ... store to register m // レジスタnの初期化 d ... dup sn ... store to register n // 最大/最小をもとめつつレジスタsに待避 lAx ... call A // 最大値の除去(符号反転させてpush) lm ... load from register m d ... dup 0 ... push 0 r ... swap - ... subtract Ss ... push to register s //最小値の除去 ln ... load from register n d ... dup 0r- ... invert sign Ss ... push to register s // レジスタiの初期化(最大/最小の除去分で+2、加算が二稿演算なので-1) lz ... load from register z 2+ ... add 2 1- ... subtract 1 si ... store to register i // レジスタsを総和 lBx ... call B // 総和から平均を求める(最大/最小を抜いたので-2) Ls ... pop from register s lz ... load from register z 2- ... subtract 2 / ... divide p ... print
うわぁ...
自分で書いたんだけどね...。
上にも書いたけどメインスタックの要素数はzコマンドで取れるのだけれど、レジスタスタックは深さが分からない。
ということでレジスタzを使ったのだけれど、今になってよくよく考えるとsレジスタすら使わなくてよかったことに気がついた。
加算しながら最大/最小みつけて、最後総和から引くだけじゃん...orz
まぁ、いいや、いまから書き直すのもあれだし。
あと、シェルスクリプトにtimeout.rbというのがあるけど、これはご想像の通りRubyスクリプト。
久々にRuby書いたのでつっこみは無しの方向で。
#!/usr/bin/ruby timeout = ARGV.shift prog = ARGV.shift args = ARGV p1 = fork { exec [prog,prog], *args } p2 = fork { sleep timeout.to_f Process.kill( 9, p1 ) } Process.wait p1 Process.kill( 9, p2 )
全盛期でもそんな変わらない、というオチがあったり...。
ということで、シェル+dc+Rubyという超変態スクリプトではかっていきます。
あと、テストプログラムのひな形( Take.hs )
import System import {-# SOURCE #-} Primes main :: IO () main = do args <- getArgs case args of (x:_) -> print $ sum $ takePrimes $ (read :: String -> Integer) $ x [] -> do progName <- getProgName putStrLn $ "usage: " ++ progName ++ " n"
こちらはunder用(Under.hs)
import System import {-# SOURCE #-} Primes main :: IO () main = do args <- getArgs case args of (x:_) -> print $ sum $ underPrimes $ (read :: String -> Integer) $ x [] -> do progName <- getProgName putStrLn $ "usage: " ++ progName ++ " n"
アルゴリズム部分だけ変えられるようにということでPrimes.hs-boot
module Primes (primes, sieve, takePrimes, underPrimes ) where primes :: (Integral a) => [a] sieve :: (Integral a) => [a] -> [a] takePrimes :: (Integral a) => a -> [a] underPrimes :: (Integral a) => a -> [a]
あとMakefile
HC=ghc HFLAGS=-Wall TYPE=FilterSieve FilterSieve2 FilterSieve3 DiffSieve DiffSieve2 DiffSieve3 DiffSieve4 DiffSieve5 Strides LazySieve Sqrt SRCS=Take.hs Under.hs $(patsubst %,Primes.%.hs,$(TYPE)) OBJS=$(SRCS:.hs=.o) INTS=$(SRCS:.hs=.hi) BINS=$(patsubst %,t%,$(TYPE)) $(patsubst %,u%,$(TYPE)) DCS=$(BINS:=.dc) DATS=$(BINS:=.dat) BOOTSRCS=Primes.hs-boot BOOTINTS=$(BOOTSRCS:.hs-boot=.hi-boot) PHONY= PHONY+=all all: $(BINS) PHONY+=clean clean: $(RM) $(OBJS) $(INTS) $(BINS) PHONY+=dc dc: $(DCS) PHONY+=dat dat: $(DATS) result.pdf: $(DATS) plot.gpi gnuplot plot.gpi %.dat: %.dc dc $^ > $@ t%.dc: t% test.sh ./test.sh ./t$* 100 10000 100 u%.dc: u% test.sh ./test.sh ./u$* 1000 100000 1000 t%: Primes.%.o Take.o $(BOOTINTS) $(HC) $(HFLAGS) -o $@ $^ u%: Primes.%.o Under.o $(BOOTINTS) $(HC) $(HFLAGS) -o $@ $^ %.hi: %.o @: %.o: %.hs $(HC) $(HFLAGS) -c -o $@ $^ %.o: %.lhs $(HC) $(HFLAGS) -c -o $@ $^ %.hi-boot: %.o-boot @: %.o-boot: %.hs-boot $(HC) $(HFLAGS) -c -o $@ $^ %.o-boot: %.lhs-boot $(HC) $(HFLAGS) -c -o $@ $^ .PHONY: $(PHONY) .PRECIOUS: %.dc
ふぅ。やっと下準備。
ここまで前置き。
と、いうわけで、まず比較対象となる一番遅い、と思われるものから。
Primes.FilterSieve.hs
module Primes (primes, sieve, takePrimes, underPrimes ) where import List primes :: (Integral a) => [a] primes = sieve [2..] sieve :: (Integral a) => [a] -> [a] sieve (x:xs) = x:sieve (filter ((/=0) . (`mod`x)) xs) sieve [] = [] takePrimes :: (Integral a) => a -> [a] takePrimes n = genericTake n $ primes underPrimes :: (Integral a) => a -> [a] underPrimes n = takeWhile (< n) $ primes
入門書とかでよく見る形。
filterを使ってsieve関数を実装しているのでFilterSieveと勝手に名付けました。
ネーミングセンスが(ry
で、前回エントリでも使ったdiffを使った定義。
Primes.DiffSieve.hs
module Primes (primes, sieve, takePrimes, underPrimes ) where import List diff :: (Ord a) => [a] -> [a] -> [a] diff xxs@(x:xs) yys@(y:ys) | x > y = diff xxs ys | x < y = x:diff xs yys | otherwise = diff xs ys diff [] _ = [] diff xxs [] = xxs primes :: (Integral a) => [a] primes = sieve [2..] sieve :: (Integral a) => [a] -> [a] sieve (x:xs) = x:sieve (xs `diff` (iterate (+x) x)) sieve [] = [] takePrimes :: (Integral a) => a -> [a] takePrimes n = genericTake n $ primes underPrimes :: (Integral a) => a -> [a] underPrimes n = takeWhile (< n) $ primes
ということで、結果がこちら
だいぶいいかんじ。
2の倍数、3の倍数、nの倍数の特殊化。
さて、diffを使うでもなく、なにか高速化する方法はないか？と考えてまず出てくるのが2の倍数を特殊扱いする方法。
C言語なんかでも使う方法だけど、どれほど効果があがるのか。
考え方としては
sieve [2..]
よりは
2:sieve [3,5..]
の方が単純2倍速いはず。...だよね？
あと、さらに3の倍数を特殊化することもできる。
2の倍数と違って[3,5..]のような単純な形で表せないのでけど、2と3の最小公倍数が6であることを考えれば6周期ごとに同じ差で表せるはず。
| 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | ... | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 2で割った余り | 1 | 0 | 1 | 0 | 1 | 0 | 1 | 0 | 1 | 0 | 1 | 0 | 1 | ... |
| 3で割った余り | 1 | 2 | 0 | 1 | 2 | 0 | 1 | 2 | 0 | 1 | 2 | 0 | 1 | ... |
| 2でも3でも割り切れない数 | 1 | 5 | 7 | 11 | 13 | ... |
ということで、4つ飛ばして2つ飛ばして、をワンセットで考えるとループすることがわかる。
3の次の素数5から考えると2,4,2,4,2,4,...と飛ばしていけばいいことが分かる。
よって、
primes = 2:3:sieve (scanl (+) 5 $ cycle [2,4])
となります。具体的には
(scanl (+) 5 $ cycle [2,4]) -- [5,7,11,13,17,19,23,25...
と、2の倍数と3の倍数が取り除かれてるのがわかる。
よりつっこんで見てみると32 = 9以下に関しては素数と言い切ることができる。
同様に
primes :: (Integral a) => [a] primes = 2:3:5:sieve (scanl (+) 7 $ cycle [4,2,4,2,4,6,2,6])
2,3,5の倍数をあらかじめ除去。
(scanl (+) 7 $ cycle [4,2,4,2,4,6,2,6]) -- [7,11,13,17,19,23,29,31,37,41,43,47,49,53,59,61,67,71,73,77,...
ちなみにこれは、52=25までは素数。
さらにさらに、
primes :: (Integral a) => [a] primes = 2:3:5:7:sieve (scanl (+) 11 $ cycle [2,4,2,4,6,2,6,4,2,4,6,6,2,6,4,2,6,4,6,8,4,2,4,2,4,8,6,4,6,2,4,6,2,6,6,4,2,4,6,2,6,4,2,4,2,10,2,10])
やたら爆発的に伸びますね。
まぁ、さきのと比べると周期が7倍になりますので...。
(scanl (+) 11 $ cycle [2,4,2,4,6,2,6,4,2,4,6,6,2,6,4,2,6,4,6,8,4,2,4,2,4,8,6,4,6,2,4,6,2,6,6,4,2,4,6,2,6,4,2,4,2,10,2,10]) -- [11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,...
これは72=49まで素数のはず。
ということで、ここまでやれば速くなるだろう、という予測の下、実際の結果がこちら。
無印がそのままsieveを使い、2は2の倍数、3は2,3の倍数をそれぞれ除外した物。
あ れ？
全然変わらない。というよりはむしろ遅くなってる。
当然DiffSieveで同様のことをやっても全然速くならない。
無印がそのままsieveを使い、2は2の倍数を、3は2,3の倍数を、4は2,3,5の倍数、5は2,3,5,7の倍数をそれぞれ除外した物。
むぅ。もっと大がかりな変更が必要な様子。
ということで、新しい方法を考える。
上のscanl...の件で「52=25までは素数」「72=49までは素数」、という表現をちょこちょこ書いていたのですが、これがナニカ使えるかも。
と思っていじくり回した結果が↓のもの。
module Primes (primes, sieve, takePrimes, underPrimes ) where import Data.List import Control.Arrow diff :: (Ord a) => [a] -> [a] -> [a] diff xxs@(x:xs) yys@(y:ys) | x > y = diff xxs ys | x < y = x:diff xs yys | otherwise = diff xs ys diff [] _ = [] diff xxs [] = xxs primes :: (Integral a) => [a] primes = sieve [2..] sieve :: (Integral a) => [a] -> [a] sieve (x:xs) = x:sieve (xs `diff` (iterate (+x) x)) sieve [] = [] multiples :: (Num a) => a -> [a] multiples x = iterate (+x) x strides :: (Num a, Ord a, Enum a) => [a] -> a -> [a] strides primeList nextPrime = diffList $ nonMultiples primeList where lccm = foldl1' (*) primeList nonMultiples = foldl' (\acc -> (acc `diff`) . multiples) [nextPrime..nextPrime+lccm] diffList xxs@(x:xs) = zipWith subtract xxs xs diffList [] = [] growPrimes :: (Integral a) => ([a],a) -> ([a],a) growPrimes (primeList,nextPrime) = second head $ span (<n) $ scanl (+) nextPrime $ cycle $ strides primeList nextPrime where n = (nextPrime - 1) ^ 2 underPrimes :: (Integral a) => a -> [a] underPrimes n = takeWhile (<n) $ fst $ underPrimes' n where underPrimes' n' | n' < 500 = second head $ span (<n') $ primes | otherwise = let (primeList,nextPrime) = underPrimes' $ isqrt n' in first (primeList++) $ growPrimes $ (primeList,nextPrime) takePrimes :: (Integral a) => a -> [a] takePrimes n = takePrimes' (second head $ span (<500) $ primes) n where takePrimes' :: (Integral a) => ([a],a) -> a -> [a] takePrimes' (primeList,nextPrime) n' = let len = genericLength primeList in if n' < len then genericTake n' primeList else (primeList ++) $ takePrimes' (growPrimes (primeList, nextPrime)) $ n' - len isqrt :: (Integral a, Integral b) => a -> b isqrt = ceiling . sqrt . fromIntegral
なんか込み入ってしまったので、軽く説明。
ポイントとなるのがstridesで、これはさっきの[2,4]だとか[4,2,4,2,4,6,2,6]を素数列から生成する関数。(引数は素数列+次の素数)
このstridesの返値をscanlすると、基にした素数列の最後の数の二乗までは正しい素数となるはず。
と、いうことは、例えば[2,3] 5から始めたとすると
strides [2,3] 5 -- [2,4] -- scanlを使って、 [2,3,(5),(5+2),(5+2+4),(5+2+4+2),(5+2+4+2+4),...]を生成。 -- つまり、 [2,3,5,7,11,13,17,..]。これは3^2=9までは正しい素数。ということで、 strides [2,3,5,7] 11. -- [2,4,2,4,6,2,6,4,2,4,6,6,2,6,4,2,6,4,6,8,4,2,4...] -- 同様に、[2,3,5,7,(11),(11+2),(11+2+4),(11+2+4+2),(11+2+4+2+4),(11+2+4+2+4+6),..] -- つまり、[2,3,5,7,11,13,17,19,23,29,...]。これは7^2=49までは正しい素数。
の用に、素数列を「成長」させていく、という考え方。やたら七面倒なことをやってますが、紆余曲折の末です。
で、この方法の最大の特徴は、最初は3まで、次は7まで、その次は47まで、という風に「素数である」と言える範囲が爆発的に増えることと、その上限までは一定時間で求まる(?)ということ。
なので、実行時間のグラフは階段形になるかな、とか予想しつつ実行してみたのですが、実際はこんな感じ。
みごとに真っ平ら。
などと喜んでいたのですが、ある一定の値を超すと途端に落ちる...orz。
初期値を適当に決めてるのでそのあたりをどうにかすればいいのだろうけど、いい方法が思いつかない。
あといちいちテストする気力もやる気もなかった。
とどめにもっと速い方法見つけた、という諸事情につきこいつは放棄！！！
まぁ、書いていたときから無駄が多いなぁ、とか、書いている途中に分け分からなくなってきた、と色々思っていたし、むしろここまで数字が出たことに驚き、という感じ。
他の奴と値が離れ過ぎていてどれくらい速いか分かりづらいのですが、
underPrimes :: (Integral a) => a -> [a] underPrimes n = foldl' (\acc -> (acc `diff`) . tail . multiples) [2..n] [2..isqrt n]
↑のような非常な純粋な実装と比較しても十分勝てるぐらいの速さ。
SqrtはtakePrimesを書けないないのでupper limitなケースだけ
あとは一定値で落ちることをなくせれば...。
これが最速？
ということで、まずグラフがこんな感じ。Stridesと比較しています。
upper limitの方に関してはもはやほぼ誤差。
takeの方が随分ひらいた。むぅ。これはStrideの初期値の取り方が悪いということか。
とはいえどうすればいいか分からない。
で、肝心のソースコードが↓のような感じ。
module Primes (primes, sieve, takePrimes, underPrimes ) where import Data.List diff :: (Ord a) => [a] -> [a] -> [a] diff xxs@(x:xs) yys@(y:ys) | x > y = diff xxs ys | x < y = x:diff xs yys | otherwise = diff xs ys diff [] _ = [] diff xxs [] = xxs merge :: (Ord a) => [a] -> [a] -> [a] merge xxs@(x:xs) yys@(y:ys) | x < y = x:merge xs yys | x > y = y:merge xxs ys | otherwise = x:merge xs ys merge [] yys = yys merge xxs [] = xxs primes :: (Integral a) => [a] primes = sieve [2..] multiples :: (Num a) => a -> [a] multiples x = iterate (+x) x sieve :: (Num a, Ord a) => [a] -> [a] sieve = sieve' [2] sieve' :: (Num a, Ord a) => [a] -> [a] -> [a] sieve' [] (x:xs) = sieve' [x] xs sieve' (p:ps) xxs = let (p1,p2) = span (< p^2) ((xxs `diff`) $ tail $ multiples p) in (p1 ++) $ sieve' (ps ++ p1) p2 underPrimes :: (Integral a) => a -> [a] underPrimes n = takeWhile (< n) $ sieve [2..] takePrimes :: (Integral a) => a -> [a] takePrimes n = genericTake n $ sieve [2..]
戦略としてはsieveにもう一つ引数を増やして、その時点ですでに素数と言える数列、つまり篩いに使った数の二乗までの素数を引数として保持しつつもそれは「素数だ！」と言って返してしまう。
で、必要になれば引数として取っておいた(?)素数を使って篩いをかけ直す。
必要になるまで篩いにかけない、という気分で名前はLazySieve。ネーミングセン(ry
とりあえず、前回のエントリの10番目の問題は3秒ちょいで終了。
Enterキーたたいてから少し間を置くけど、十分でしょう。
これ以上いい実装が思い浮かばないので、このあたりでFA。
多分最新の研究とかそういうの使うともっと速くなるのだろうけど、個人的には満足。
最後に全部のグラフを重ねてお別れです。ノシ
*1:ネーミングセンスがひどい、というのは仕様ということで...
