C#:エラトステネスの篩(ふるい)再び

  • 8
    いいね
  • 2
    コメント

背景

C#:エラトステネスの篩を使った素数の計算という記事を書いたら、yone64さんからコメントがあって、列挙開始も速く、最大値を指定する必要もない「少しずつふるいにかけるメソッド」を提示してもらいました。

でも、コードを読む力がない僕は30秒で挫折...(^^;;
リンク先のページを読んでおぼろげながら何をやろうとしているのかが分かりました。

コード読むより自分で考えた方が面白そう、ということで、「少しずつふるいにかける」をヒントに自力でアルゴリズムを考え、C#のコードにしてみることにしました。

最終版エラトステネスの篩

以下、出来上がったC#のコードを示します。
結果的に、yone64さんのコードと同じことをしているコードになったような気がします。たぶん...(^^;;

求める値は、intの範囲内です。long型まで対応するとなるとさらに工夫が必要になりますし、longよりも大きな素数も求めたいとなると、また違った工夫が必要になってきますが、今の僕にとっては、これが最良の素数を求めるメソッドだと思います。

PrimesBest.cs
 // 下限、上限が指定できるbool型配列
 class BoundedBoolArray {
     private BitArray _array;
     private int _lower;

     public BoundedBoolArray(int lower, int upper) {
         _array = new BitArray(upper - lower + 1);
         _lower = lower;
     }

     public bool this[int index] {
         get {
             return _array[index - _lower];
         }
         set {
             _array[index - _lower] = value;
         }
     }
 }

 static IEnumerable<int> Primes() {
     // 2,3は既知の素数とする
     var primes = new List<int>() { 2, 3 };
     foreach (var p in primes)
         yield return p;

     // 4以上の整数から素数を列挙する。int.MaxValueを超えたときには対処していない
     int ix = 0;
     while (true) {
         int prime1st = primes[ix];
         int prime2nd = primes[++ix];
         // ふるい用の配列の下限、上限を求め、配列を確保する。
         var lower = prime1st * prime1st;
         var upper = prime2nd * prime2nd - 1;
         // ふるいは、[4:8], [9:24], [25:48], [49:120]... と変化する。
         // []内の数値は、配列の下限と上限
         var sieve = new BoundedBoolArray(lower, upper);

         // 求まっている素数を使い、ふるいに掛ける 
         foreach (var prime in primes.Take(ix)) {
             var start = (int)Math.Ceiling((double)lower / prime) * prime;
             for (int index = start; index <= upper; index += prime)
                 sieve[index] = true;
         }

         // ふるいに掛けられて残った値が素数。これを列挙する。
         // 併せて、求まった素数は、primesリストに記憶していく。
         // この素数が次のステップ以降で、ふるいに掛ける際に利用される。
         for (int i = lower; i <= upper; i++) {
             if (sieve[i] == false) {
                 primes.Add(i);
                 yield return i;
             }
         }
     }
 }

簡単な説明

コードにコメント書いておきましたが、備忘のため、簡単にこのコードが何をやっているかを説明しておきます。

まず、既知の素数 2 を使い篩にかけます。この時の篩の下限と上限は、2^2, 3^2-1 です。つまり、4~8までの素数を求めます。この時の 3^2-1 の3は既知の素数として扱っています。これで、素数 5,7が求まります。

次に、篩に掛ける素数を一つ増やし、2,3で篩に掛けます。篩の範囲は、9 (3^2) ~ 24 (5^2-1) です。これで素数 11,13,17,19,23が求まります。

次は、2,3,5で篩に掛け、29,31,37,41,43,47 が素数として求まります。この時の篩の範囲は、25 (5^2) ~ 48 (7^2-1) です。

このステップを繰り返せば、どんどん素数が求まることになります。

なお、下限、上限を指定できる配列もどきのクラスBoundedBoolArrayを作成しています。こうすれば、頭の悪い僕でも、篩にかける部分(trueをセットしているところと、素数を列挙しているところ)が苦労なく書くことができます。

少しずつ篩に掛けるので、「1億以下の素数をすべて求めよ」という問いに対して、トータルの処理時間は通常のエラトステネスの篩のコードよりは、わずかに遅くなりますが、列挙開始までの待ちがほぼゼロ、篩のサイズを指定する必要が無いという大きな利点があります。
「100000番目の素数を求めよ」という問題があった場合、エラトステネスの篩を使った通常のやり方では、篩のサイズをいくつにしたら良いかわからないですから、今回作成したメソッドはとても有用であると思います。

使い方

戻り値を、IEnumerable<int>にすれば、LINQ使って、

    // 1000個の素数を取得
    var primes = Primes().Take(1000);

とか

    // 500以下の素数を取得
    var primes = Primes().TakeWhile(p => p <= 500);

とか

    // 500より大きい最初の素数を取得
    var prime = Primes().First(p => p > 500);

とか書くことができます。

以下のようなコードだと、永久に素数を求め続けようとしますので、このような書き方は基本的にできません。メモリオーバーになるか、オーバーフローするかしてしまいます。

    var primes = Primes();

補足

ちなみに、このプログラムでは、篩のサイズは、[4:8], [9:24], [25:48], [49:120]... と変化するわけですが、理論的には、


  [4:8]
  [9:48]       48 = 7 * 7 - 1   (7は、9未満の最大の素数)
  [49,2208]    2208 = 47 * 47 - 1  (47は、49未満の最大の素数)
  ……

としていくことも可能です。

しかし、そうすると、あっというまに大きなサイズの配列が必要になり、必要のない大きな数の素数まで求めることになってしまいます。これでは、返って速度の遅いプログラムになってしまいすし、列挙を早くするという最初の趣旨にも沿いません。