一部の熱狂的な C++ プログラマーは、自分のプログラムは 8 秒しかかからないと主張しています。
それは実時間ですか、それとも CPU 時間ですか?
壁時計で、タスクが 100 個の CPU に分割されている場合、たとえば、あまり印象的ではありません (まともです)。1000 個に分割されている場合は、哀れです。
CPU 時間の場合:
実際に 10 11までふるいにかけても時間に達していないことは確かです。
それまでは 4×10 9素数より数個多く、通常の 2 ～ 3 GHz CPU を想定すると、素数ごとに 4 ～ 6 サイクルになります。
エラトステネスのふるいでも、アトキンのふるいでも、それを達成することはできません. 各素数を検査して数え、各複合体をそのようにマークして検査する必要があります。これにより、配列の初期化、ループ境界のチェック、ループ変数の更新、冗長なマーキングなどを考慮せずに、ふるいの数値ごとに 2 サイクルの理論的な下限が得られます。その理論的な限界に近づくことはありません。
いくつかのデータ ポイント:
Daniel Bernstein の primegen (Sieve of Atkin) は、私の 32KB L1 キャッシュを最大限に活用するようにふるい分けブロックを調整したもので、素数を 10 11にふるい分けてカウントするのに 90 秒かかります (デフォルトのふるいブロック サイズ 8K では 234 秒)。言葉）私のCore i5 2410M（2.3GHz）。(これは 2 32までの範囲で最適化されていますが、それを超えると著しく遅くなります。制限 10 9では、時間は 0.49 対 0.64 秒です。)
私のセグメント化された Eratosthenes のふるい は、公開されていない内部要素を使用してリストの作成を回避し、340 秒で 10 11までふるい分けてカウントします(嗅ぎます :-/ でも、10 9の場合は 2.92 秒かかりました - 近づいており、10 の間のどこかです) 12と 10 13はprimegenを追い越します:) 公開されたインターフェイスを使用して素数のリストを作成すると、32 ビット GHC でコンパイルする場合と同様に、所要時間が約 2 倍になります。
したがって、報告された 8 秒の時間 (CPU 時間の場合) が正しければ、実際に全体をふるいにかけることなく素数の数をカウントするアルゴリズムであることに賭けたいと思います。applicative's answerで示されているように、それははるかに高速に実行できます。
dafis@schwartz:~/Haskell/Repos/arithmoi> time tests/primeCount 100000000000
4118054813
real 0m0.145s
user 0m0.139s
sys 0m0.006s
10^11 は 32 ビット演算には大きすぎることに注意してください。また、10^11 ビット = 12.5 GB であり、Haskell の 32 ビット アドレス空間に収まるにはデータが多すぎます。したがって、一度にビットマップ全体をメモリに格納することはできません。
その範囲をふるいにかけるには、セグメント化されたふるいを使用する必要があります。32 ビットのアドレス空間に制限されていない場合でも、このような大きな配列を使用すると、キャッシュ ミスが頻繁に発生するため、パフォーマンスが低下します。プログラムは、ほとんどの時間をメイン メモリからデータが転送されるのを待つことに費やします。L2キャッシュに収まるチャンクでふるいにかけます（ふるいをL1に合わせて高速化しようとしてもうまくいきませんでした。GHCランタイムのオーバーヘッドが大きすぎて機能しないと思います）。
また、ふるいからいくつかの小さな素数の倍数を削除します。これにより、必要な作業が減り、さらにふるいを小さくすることでパフォーマンスが向上します。偶数の消去は簡単で、3 の倍数は簡単、5 の倍数はそれほど難しくありません。
最後に、10^11 未満の素数の数は 2^32 未満の数にすぎないため、実際の整数を一度にすべて格納することはできません。
ふるいをビット配列のリストとして保存し、2、3、および 5 の倍数を削除すると、チャンクを保存するために約 3.3GB が必要になるため、実際に最大 4GB を使用できる場合は、それが収まります。しかし、不要になったチャンクはすぐにガベージ コレクションする必要があります。
(重要な場合、私の現在の実装では IOUArray Integer Bool と複数のスレッドを使用して、検索空間の独立した部分範囲を処理しています。現在、10MB の配列チャンクから 2 の倍数をすべて削除するには数秒かかります...)
それは問題である。
IntインデックスとunsafeRead/に使用してunsafeWrite、配列の読み取りと変更を行います。Integer計算は計算よりもはるかに遅く、/Intで得られる境界チェックは本当に痛い.readArraywriteArray - 10MB のチャンクは大きすぎるため、キャッシュの局所性が失われます。せいぜい数百 KB のチャンクを使用します (L2 キャッシュからその他の必要なスペースを差し引いたもの)。
- それでも、
Integerインデックス、境界チェック、および 10MB のチャンクを使用しても、2 の倍数を削除するのに数秒かかることはありません。あなたのコードを見ることができますか？
休暇後の更新:
10 11までの素数をふるいにかけるのに8分もあれば、深い魔法がなくても可能です。ここにはキャッシュ効果がないはずなので、コアを 1 つから 4 つにすると 8 倍のスピードアップが得られるかどうかはわかりませんが、それが何であれ、コードを見なければ調査できません。
それでは、コードを見てみましょう。
まず、誤り：
vs <-
mapM
(\ start -> do
let block = (start, start + block_size)
v <- newEmptyMVar
writeChan queue $ do
say $ "New chunk " ++ show block
chunk <- chunk_new block
sieve_top base chunk
c <- chunk_count chunk
putMVar v c
return v
)
(takeWhile (< target) $ iterate (+block_size) base_max)
数値base_max + k*block_sizeはそれぞれ 2 つのチャンクで表示されます。それらのいずれかが素数である場合、その素数は 2 回カウントされます。また、上限を で制限する必要がありますtarget。
次に、パフォーマンスの側面に進みます。
飛び出すことの1つは、それが本当におしゃべりであり、キャッシュに調整すると測定できるほどおしゃべりであるということです（私は512KBのL2キャッシュに256KBのブロックを取りました）-その後、スレッドはメッセージのためにblock_size戦うことで遅くなります.stdoutif prime < 100 then say $ "Sieve " ++ show prime else return ()
(サイレンシングされた) ふるい分けループを見てみましょう。
chunk_sieve :: Chunk -> Integer -> IO ()
chunk_sieve array prime = do
(min, max) <- getBounds array
let n0 = min `mod` prime
let n1 = if n0 == 0 then min else min - n0 + prime
mapM_
(\ n -> writeArray array n (n == prime))
(takeWhile (<= max) $ iterate (+prime) n1)
時間がかかることの 1 つは、各インデックスが倍数がマークされている素数と比較されることです。1 回の比較は安価ですが (1 回の比較よりもかなり高価ですがInt)、膨大な数の比較を行うと、そのうちの 1 つだけが結果を出す可能性がありTrueます。無条件に書き込みFalse、必要に応じTrueてループの後に素数のインデックスに書き込むと、かなりのスピードアップが得られます。
タイミングのために、ターゲットを 10 9に減らし、2 つのコアで実行しました。元のコードは 155 秒 (経過、ユーザーは 292 秒)block_sizeかかり、148 秒が短縮され、143 秒が沈黙しました。比較は省略しますが、
mapM_
(\ n -> writeArray array n False)
(takeWhile (<= max) $ iterate (+prime) n1)
when (min <= prime && prime <= max) $ writeArray array prime True
131秒で実行されます。
今度は、さらに大きなスピードアップを行うときです。境界チェックには多くの時間がかかることは既に述べましたか? ループ条件により、範囲外のアクセスが試行されないことが保証されるため (素数は十分に小さいため、Intオーバーフローが発生することはありません)、実際には未チェックのアクセスを使用する必要があります。
chunk_sieve :: Chunk -> Integer -> IO ()
chunk_sieve array prime = do
(min, max) <- getBounds array
let n0 = min `mod` prime
n1 = if n0 == 0 then min else min - n0 + prime
n2 = fromInteger (n1 - min)
mx = fromInteger (max - min)
pr = fromInteger prime
mapM_
(\ n -> unsafeWrite array n False)
(takeWhile (<= mx) $ iterate (+pr) n2)
when (min <= prime && prime <= max) $ writeArray array prime True
これにより、実行時間が 96 秒に短縮されます。はるかに優れていますが、それでもひどいです。犯人は
takeWhile (<= mx) $ iterate (+pr) n2
GHC はその構成をうまく融合できず、Intトラバースされるボックス化された s のリストを取得します。それを算術シーケンスに置き換えると、 GHC はボックス化されていないs を使用して 37 秒の[n2, n2+pr .. mx]ループを喜んで作成します。Int#
はるかに優れていますが、それでも悪いです。現在、最大の時間を消費しているのは
chunk_count :: Chunk -> IO Integer
chunk_count array = do
(min, max) <- getBounds array
work min max 0
where
work i max count = do
b <- readArray array i
let count' = count + if b then 1 else 0
evaluate count'
let i' = i+1
if i' > max
then return count'
else work i' max count'
繰り返しますが、境界チェックには多くの時間がかかります。と
chunk_count :: Chunk -> IO Integer
chunk_count array = do
(min, max) <- getBounds array
work 0 (fromInteger (max-min)) 0
where
work i max count = do
b <- unsafeRead array i
let count' = count + if b then 1 else 0
evaluate count'
let i' = i+1
if i' > max
then return count'
else work i' max count'
15 秒まで短縮されています。さて、でstrictevaluate count'を作成するには、やや高価な方法です。の代わりに最後の行で使用すると、実行時間が 13 秒に短縮されます。より適切な (少なくとも GHC には) 方法で定義すると、workcountelse work i' max $! count'evaluatework
chunk_count :: Chunk -> IO Integer
chunk_count array = do
(min, max) <- getBounds array
let mx = fromInteger (max-min)
work i !ct
| mx < i = return ct
| otherwise = do
b <- unsafeRead array i
work (i+1) (if b then ct+1 else ct)
work 0 0
時間を 6.55 秒に短縮します。現在、測定可能な違いが生じる状況にありsay $ "New chunk " ++ show block、無効にすると 6.18 秒に短縮されます。
ただし、配列からバイトを読み取り、不要なビットをマスクして、個々のビットごとに 0 と比較することにより、設定されたビットをカウントすることは、最も効率的な方法ではありません。Word配列から ( castIOUArrayを介して) s全体を読み取り、popCountを使用する方が高速です。素数の二乗がチャンクの上限よりも大きくなったときにマーキングを停止する
sieve_top :: Chunk -> Chunk -> IO ()
sieve_top base chunk = work 2
where
work prime = do
chunk_sieve chunk prime
mp <- chunk_next_prime base prime
case mp of
Nothing -> return ()
Just p' -> do
(_,mx) <- getBounds chunk
when (p'*p' <= mx) $ work p'
3.9秒まで。まだ壮観ではありませんが、出発点を考えると悪くありません。他の悪い動作が減った後のキャッシュの局所性の重要性を説明するだけです。元の 10MB ブロック サイズの同じコードは 8.5 秒かかります。
コードのもう 1 つの小さな問題は、すべてのスレッドがふるい分けに小さな素数の同じ可変配列を使用することです。これは変更可能であるため、アクセスを同期する必要があり、オーバーヘッドが少し追加されます。スレッドが 2 つしかない場合、オーバーヘッドはそれほど大きくありません。不変のコピーを使用してふるい分けを行っても、ここでは時間が 3.75 秒に短縮されるだけですが、スレッドが多いほど効果が大きくなると予想されます。(ハイパースレッディングを使用する物理コアは 2 つしかないため、2 つ以上のスレッドを使用して同じ種類の作業を行うと、そこから導き出された結論が無効になる可能性があるスローダウンが発生しますが、4 つのスレッドを使用すると、5.3 秒に対して不変配列で 4.55 秒が得られますこれは、増加する同期オーバーヘッドを裏付けているようです。)
GHC のオプティマイザー (より多くのワーカー/ラッパーの変換、いくつかの静的引数の変換) のコードを書くことで、より多くの計算を排除することで得られるものはまだ少しありますがInteger、それほど多くはなく、おそらく 10-15% です。
次の大きな改善は、ふるいから偶数を取り除くことです。これにより、作業、割り当て、および実行時間が半分以下に削減されます。素数ふるいは偶数を考慮すべきではありません。実際、それは無意味な無駄な作業です。