FastSieve

超高速なエラトステネスの篩の実装です。

アルゴリズム詳細

このプログラムは、エラトステネスの篩に複数の高度な最適化手法を組み合わせることで、極めて高いパフォーマンスを実現しています。以下に、主要な最適化技術を詳述します。

1. ホイール法 (Wheel Factorization)

• 原理: 素数2, 3, 5の倍数をあらかじめ篩の対象から除外する最適化手法です。この3つの素数の最小公倍数は30であるため、「wheel30」と呼ばれます。
• 詳細: 30までの自然数のうち、2, 3, 5のいずれの倍数でもない数は {1, 7, 11, 13, 17, 19, 23, 29} の8個しか存在しません。これは、任意の整数 30k にこの8つの数を加えた数のみが素数候補となることを意味します（1は素数ではないため、実質的には7以上の素数候補がこのパターンに従います）。
• 効果: この手法により、篩うべき数の母集団が元の数の約 8/30 (≈ 26.7%) にまで劇的に減少し、計算量が大幅に削減されます。また、候補が30ごとに8個という周期性を持つため、1バイト（8ビット）で30の範囲をマッピングするビットパッキングと極めて相性が良くなります。

2. ビットパッキング (Bit Packing)

• 原理: 篩の状態を管理する配列において、1つの素数候補を1ビットに対応させることで、メモリ使用量を極限まで削減する手法です。
• 詳細: std::vector<unsigned char> (C++では char は1バイト=8ビット) を使い、各ビットをフラグとして利用します。ホイール法によって選別された30ごとの8個の候補が、配列の1バイトに以下のようにマッピングされます。
  • sieve[k] の0ビット目: 30k + 1
  • sieve[k] の1ビット目: 30k + 7
  • sieve[k] の2ビット目: 30k + 11
  • ...
  • sieve[k] の7ビット目: 30k + 29 ある数が素数でないと判明した場合、対応するビットを 1 から 0 に反転させます。
• 効果: メモリ使用量を単純な bool 配列と比較して1/8に圧縮します。これにより、CPUのキャッシュ（特にL1/L2キャッシュ）に格納できる範囲が8倍に広がり、キャッシュヒット率が大幅に向上します。メモリアクセスは計算処理において主要なボトルネックの一つであるため、このキャッシュ効率の改善は実行速度に絶大な効果をもたらします。

3. 区間篩 (Segmented Sieve)

• 原理: 篩にかける巨大な数値範囲を、CPUのキャッシュサイズに収まるような小さな「セグメント」に分割し、セグメントごとに篩処理を完結させる手法です。
• 詳細:
  1. ベース素数の計算: まず、篩いたい最大値 N の平方根 √N までの素数（ベース素数）を、比較的小さな篩（ベース篩）で全て求めます。なぜなら、N 以下の全ての合成数 C は、必ず √N 以下の素因数 p を持つからです。（もし全ての素因数が >√N なら、その積は N を超えてしまうため）。
  2. セグメント処理: √N から N までの範囲を、SEGMENT_SIZE で定義された固定長の区間（セグメント）に分割します。そして、各セグメントを表す小さな篩配列をメモリ上に確保し、ベース素数を使ってそのセグメント内の合成数を篩い落とします。この処理を最後のセグメントまで繰り返します。
• 効果: セグメントのサイズをCPUのL2キャッシュ（数MB程度）に完全に収まるように設計することで、セグメント内の篩処理中に発生するメモリアクセスがほぼ全てキャッシュヒットするようになります。これにより、メインメモリへのアクセス頻度が劇的に減少し、メモリアクセスのレイテンシが隠蔽されるため、処理速度が飛躍的に向上します。この実装では、SEGMENT_SIZE は 256KB * 30 バイトに設定されており、多くの現代的CPUのL2/L3キャッシュサイズを意識した値となっています。

4. 事前計算テーブルによるループ最適化

• 原理: 篩処理の最深部ループは、プログラム全体の実行時間の大半を占めます。このループ内の計算を極限まで単純化するため、複雑な計算パターンを事前に解析し、その結果を定数テーブルとしてコードに埋め込む手法です。
• 詳細: ホイール法を適用した篩では、ある素数 p = 30k+r1 の倍数 p * q (ここで q = 30m+r2 も素数候補) を篩い落とす必要があります。この積 p*q がビット配列上のどのバイトのどのビットに対応するかを計算するのは、 (p*q)/30 のような高コストな演算が必要となり、非常に非効率です。しかし、この計算結果のパターンを分析すると、p の系列 (r1 の値) ごとに、q の8つの系列 (r2) との積のオフセットが、k (=p/30) に関する一次式で表せることが分かります。 この実装では、その一次式の定数項と係数を offsetIdxConst と offsetIdxCoeff というテーブルに事前計算しています。
  • offsetMod8: p*q の積の剰余 (%30) が8つの候補のうちどれに対応するかのパターンを保持します。これにより、どのビットを 0 にすべきかが瞬時に分かります。
  • offsetIdxConst, offsetIdxCoeff: ビットが含まれるバイトのインデックスを計算するためのオフセットです。p の系列と k (=p/30) が分かれば、8つのパターンのオフセットを オフセット = offsetIdxConst + offsetIdxCoeff * k という極めて単純な計算で求めることができます。
• 効果: 最も実行回数の多い篩い落としのループから、高コストな剰余演算や除算を完全に排除します。ループ内はテーブル参照、加算、乗算、ビット演算といったCPUが得意とする単純な命令のみで構成されることになります。これにより、CPUの命令パイプラインが効率的に動作し、スーパースカラー実行も促進されるため、驚異的なパフォーマンス向上を実現しています。これは、この実装における最も重要なコア最適化です。

使い方

ビルド

g++ -o FastSieve -O3 -march=native src/sieve.cpp

ベンチマーク版をビルドする場合:

g++ -o SieveBench -O3 -march=native src/sieve_bench.cpp

-O3 と -march=native をつけることで、お使いのCPUに合わせた最適なコードが生成され、最高のパフォーマンスが期待できます。

実行方法

2からNまでの素数を列挙する

./FastSieve <N>

例:

./FastSieve 1000000000

実行後、見つかった素数の個数と処理時間が表示されます。

MからNまでの素数を列挙する

./FastSieve <M> <N>

例:

./FastSieve 1000000000 2000000000

実行後、範囲内で見つかった素数の個数と処理時間が表示されます。

注: ベンチマーク版 (SieveBench) も同様の引数で実行できます。こちらは素数リストをメモリに保持せず、個数のみをカウントするため、より高速に動作しメモリ使用量も少ないです。