C#:エラトステネスの篩を使って素数を求める



エラトステネスの篩 (エラトステネスのふるい)を使った素数を求めるC#のプログラムです。

エラトステネスのふるいについては、Wikipediaの [こちら](http://ja.wikipedia.org/wiki/%E3%82%A8%E3%83%A9%E3%83%88%E3%82%B9%E3%83%86%E3%83%8D%E3%82%B9%E3%81%AE%E7%AF%A9)をみてください。

##ノーマルバージョン
Primesが素数を求めるメソッドです。2から引数で与えた整数(maxnum)の間の素数を列挙しています。
この引数の値で、内部でふるい用の配列(sieve)の要素数を決定します。戻り値は、```IEnumerable<T>```です。

```csharp:SieveOfEratosthenes.cs  
using System;
using System.Collections.Generic;
using System.Linq;

namespace Gushwell.Etude {
    
    class Program {
        static void Main(string[] args) {
            foreach(var n in Primes(200))
                Console.Write("{0,3} ",n);
            Console.WriteLine();
            Console.ReadLine();
        }

        static IEnumerable<int> Primes(int maxnum) {
            int[] sieve = Enumerable.Range(0, maxnum + 1).ToArray();
            sieve[1] = 0;  // 0 : 素数ではない
            int squareroot = (int)Math.Sqrt(maxnum);
            for (int i = 2; i <= squareroot; i++) {
                if (sieve[i] <= 0)
                    continue;
                for (int n = i * 2; n <= maxnum; n += i)
                    sieve[n] = 0;
            }
            return sieve.Where(n => n > 0);
        }
    }
}
```

##省メモリバージョン

もう少し、メモリ効率を良くしたバージョンも載せておきます。intの配列ではなくて、boolの配列にしています。

```csharp:PrimesSavingMemory.cs
 static IEnumerable<int> Primes(int maxnum) {
     // var sieve = new BitArray(maxnum + 1, true);
     bool[] sieve = Enumerable.Repeat(true, maxnum + 1).ToArray();
     int squareroot = (int)Math.Sqrt(maxnum);
     for (int i = 2; i <= squareroot; i++) {
         if (sieve[i] == false)
             continue;
         for (int n = i * 2; n <= maxnum; n += i)
             sieve[n] = false;
     }
     for (int i = 2; i <= maxnum; i++)
         if (sieve[i] == true)
             yield return i;
 }
```

コメントのようにBitArray使えば、もっとメモリ効率は良くなります。

[GitHub](https://github.com/gushwell/CSharpAlgorithmCollection/tree/master/SieveOfEratosthenesPart1)でこのソースコードを公開しています。

##列挙開始を早くする
ところで、ここまで書いて思ったんですが、このメソッドは、列挙が始まるのが遅すぎるのが難点ですね。全部の素数が求まらないと列挙が始まりません。
最初のfor文の中で、```yield return``` を書ければ、列挙を開始するのがもっと早くなるはずです。
いくつまでの素数が必要かわからない場合は、列挙開始が速くなれば、トータルの処理時間も短くなる可能性もあります。

で、書いたのが次のコードです。

```csharp:PrimesStartDash.cs
 static IEnumerable<int> Primes(int maxnum) {
     var sieve = new BitArray(maxnum + 1, true);
     int squareroot = (int)Math.Sqrt(maxnum);
     for (int i = 2; i <= squareroot; i++) {
         if (sieve[i] == false)
             continue;
         for (int n = i * 2; n <= maxnum; n += i)
             sieve[n] = false;
         if (sieve[i])
             yield return i;  // ここで列挙してしまう
     }
     for (int i = squareroot + 1; i <= maxnum; i++)
         if (sieve[i] == true)
             yield return i;
 }
```

うーーん、でも、```Math.Sqrt(maxnum)``` の値までしか列挙開始を早められないから、ほとんどメリット無いかな。
一瞬、良いアイデアと思ったんですが...
もっと突き詰めて考えれば、もう少し実用的なコードになるような気もするのですが、とりあえずこれで止めておきます。

##追記 ("列挙開始を早くする"の改良バージョン)

もう少し、突き詰めて考えてみました。
良く考えてみると、

- 整数2で篩を掛けた後の配列の状態を見ると、2,3,5,7 までの素数が求まっています。
- その後3で篩をかけると、11,13の素数が求まります。
- さらに4で篩をかけると、17, 19, 23の素数が求まります。（実際には、4で篩を掛けても配列の状態は変わらないのだが、説明の都合上、こう書いておく)
- さらに5で篩をかけると、29, 31の素数が求まっています。

つまり、整数nで篩を掛け終わった時には、(n+1) * (n+1) - 1 以下の素数が求まっているということです。実質的に、n+1で篩を掛けることができるのは、(n+1)*(n+1)からであり、これより小さい値は、n以下の値ですでに篩を掛けてあるからです。

そこで書き換えたのが以下のコードとなります。このコードでは、sieveのtrueとfalseの意味を反対にしています。


```csharp:PrimesStartDash2.cs
 static IEnumerable<int> Primes(int maxnum) {
    yield return 2;
    yield return 3;
    var sieve = new BitArray(maxnum + 1);
    int squareroot = (int)Math.Sqrt(maxnum);
    for (int i = 2; i <= squareroot; i++) {
        if (sieve[i] == false) {
            for (int n = i * 2; n <= maxnum; n += i)
                sieve[n] = true;
        }
        for (int n = i * i + 1; n <= maxnum && n < (i + 1) * (i + 1); n++) {
            if (!sieve[n])
                yield return n;
        }
    }
 }
        
```

[GitHub](https://github.com/gushwell/CSharpAlgorithmCollection/tree/master/SieveOfEratosthenesPart2)でこのコードを公開しています。


 
2つ目のfor文がなくなり、```yield return```は、最初のfor文の中に全て移動できた関係で、最初のバージョンから比べると、スムーズに列挙が進むようになりました。

それと、「列挙開始を早くする」で示したコードよりも、コードもそれほど複雑にならずに実現できました。


とは言え、引数に最大値を指定しないといけないという欠点は残ります。「1000までの素数を求めたい」という要求には答えられますが、「小さい順に1000個の素数を求めたい」という要求には、うまく対応できません。
yone64さんがコメントしてくれたような工夫が必要になりそうです。これについては、時間ができた時に、再度考えてみたいと思います。



-----------------------------
この記事は、Gushwell's C# Programming Pageで公開したものを加筆・修正したものです。