泥庭

2010年10月17日

素数生成アルゴリズム

Filed under: Project Euler — タグ: — yone64 @ 5:43 AM

Project Eulerに挑戦していると、素数に関する問題によく出会います。(次はフィボナッチ数列かな?)
で、適当に作成した素数生成プログラムもパフォーマンス的に限界が近いので、見直すことに。

ちなみに元のプログラムは↓な感じです。

public static IEnumerable<int> PrimeNumbers()
{
    yield return 2;
    var list = new List<int>();
    var current = 1;

    while (true)
    {
        current += 2;
        if (current < 0) yield break;
        int sqrt = (int)Math.Sqrt(current);
        if (list.TakeWhile(x => x <= sqrt).All(x => current % x != 0))
        {
            list.Add(current);
            yield return current;
        }
    }
}

intがOverflowするまで、素数を求め続けてくれます。また、intをlongに変更するとlongの限界まで演算してくるはず。多分。
でも、数が大きくなってくると結構遅い。CPUががんばってくれてたけど結構限界。
なお計算速度は下表の通り。

N

時間

N番目の素数

100,000

0.8秒

1299709

1,000,000

20秒

15485863

10,000,000

8分30秒

179424673

(Windows7 x64, Core i7 920, 12GB Memory)

1000万個求めてると8分とかかかるので、結構使えない感じです。
そこで、アルゴリズムを見直し。ということでは、速いと噂(?)の、エラストテネスのふるいを採用することに。
アルゴリズムの詳細はWikipediaでも参照してもらうとして、ざっくり書いてみまると↓なかんじです。

public static IEnumerable<int> PrimeNumbers()
{
    int maxnum = 200000000;
    var ba = new bool[maxnum];
             
    int squareroot = (int)Math.Sqrt(maxnum);
    int i;
    for (i = 2; i <= squareroot; i++)
    {
        if (!ba[i])
        {
            yield return i;
            for (int n = i * i; n < maxnum; n += i)
                ba[n] = true;
        }
    }
    for (; i < ba.Length; i++)
    {
        if (!ba[i]) yield return i;
    }
}

エラストテネスのふるいの弱点(?)は、Maxを決めてあげないといけない点ですよね。
とりあえず、2億をMaxとして指定しました。(1000万個目の素数が、2億弱だからですが。。。)

N

時間

N番目の素数

100,000

6.7秒

1299709

1,000,000

6.9秒

15485863

10,000,000

9.0秒

179424673

劇的に速くなりましたね。でも、Nが小さい場合も相応の時間がかかるのが残念な感じ。
というわけで、段階的にふるいをかけるように修正してみました。これにより、とりあえずメモリの限界まで計算してくれるはず。

public static IEnumerable<long> PrimeNumbers()
{
    var list = new List<int>();
    yield return 2;
    list.Add(2);
    yield return 3;
    list.Add(3);
    int ci = 0;
    while (true)
    {
        int current = list[ci];
        int next = list[ci + 1];
        var ba = new bool[(next - current) * (current + next)];
        ba[0] = true;
        var bn = (long)current * (long)current;
        for (int i = 0; i <= ci; i++)
        {
            var item = list[i];

            int a = (int)(bn % item);
            for (int n = item - a; n < ba.Length; n += item)
            {
                ba[n] = true;
            }
        }

        for (int index = 0; index < ba.Length; index++)
        {
            if (!ba[index])
            {
                long result = bn + index;
                yield return result;
                if (result < int.MaxValue)
                {
                    list.Add((int)result);
                }
            }
        }
        ci++;
    }
}

さらに速くなったところで、とりあえず満足。あとは、偶数を除外とかするとさらに速くなるんだろうなぁ。

N

時間

N番目の素数

100,000

0.05秒

1299709

1,000,000

0.5秒

15485863

10,000,000

6.6秒

179424673

コメントする »

まだコメントはありません。

RSS feed for comments on this post. TrackBack URI

コメントを残す

以下に詳細を記入するか、アイコンをクリックしてログインしてください。

WordPress.com ロゴ

WordPress.com アカウントを使ってコメントしています。 ログアウト / 変更 )

Twitter 画像

Twitter アカウントを使ってコメントしています。 ログアウト / 変更 )

Facebook の写真

Facebook アカウントを使ってコメントしています。 ログアウト / 変更 )

Google+ フォト

Google+ アカウントを使ってコメントしています。 ログアウト / 変更 )

%s と連携中

WordPress.com Blog.

%d人のブロガーが「いいね」をつけました。