C#:エラトステネスの篩を使った素数の計算

  • 3
    Like
  • 2
    Comment
More than 1 year has passed since last update.

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

エラトステネスのふるいについては、Wikipediaの こちらをみてください。

ノーマルバージョン

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

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の配列にしています。

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使えば、もっとメモリ効率は良くなります。

列挙開始を早くする

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

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

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の意味を反対にしています。

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;
        }
    }
 }

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

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

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


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