素数を出しておくときの高速化

Question

質問内容

素数である数を後々使うので，先に素数である数をピックアップしようと思って書いたのですが，この部分の計算が遅いことが原因で正解になりませんでした．どのようにすれば高速化できますか？

素数のとき，配列の中身を0にしようと思っています．

コード

int num,num2;
    int box[246913]={1,1};
    int i,j;
for(i=2;i<246913;i++){
        if(box[i]==0){
            for(j=i+1;j<246913;j++){
                if(j%i==0&&box[j]==0) box[j]=1;
            }
        }
    }

Comments

[e-watt]

実行時間を減らせるネタとして、ひとまずこんなのがあります。一通り実装してみてはいかがでしょうか。
・偶数な素数は2だけなので、大外のループは奇数だけで（i=3;i<QQQ;i+=2）回す
　2については、それだけ独立のループで処理する。
（偶数分の表は作らず、あとのコードで「（2以外の）偶数は素数じゃない」というコードで対処する、
　という手もあるにはあるが、「その方が良い」という場合はあまりないかも）
・Nまでの値をふるいにかけるなら、調べる数の最大値は[√N]でよい（[]はガウス記号）
・剰余の計算をしなくても、2i=i+1; (j=i*3;j<RRR;j+=i2) （iの奇数倍）がiの倍数なのだからそいつに印をつければよい(偶数倍は偶数なので前記「2の倍数のループ」ですでにマーク済み)
・（内側のループ内で）box[j]が0かを調べるまでもなく、いきなり1を代入してよい

Answer 1

素数を求めること自体が目的でないなら、固定値でいいんじゃないでしょうか？

Comments

[grape_ll]

確かに固定値でも問題ないように思えますね．ですが，高速化について考えれば他のことにも応用が利くと思いますので今回は高速化をする方法を学ぼうと思います．ご指摘ありがとうございました．

[anndonut]

素数よりも擬素数を保持しておくとメリットが大きいと思います。素数はたくさんありますが擬素数は少ないためです。巨大な素数を計算する、このサイトが参考になると思います。
https://tools.m-bsys.com/calculators/prime_number.php

[grape_ll]

ご助言ありがとうございます．後学の参考にしたいと思います．

Answer 2

結論

hnk_llさんのソースコードが以下の通りですね．

int main()
{
    int num, num2;
    int box[246913] = {1, 1};
    int i, j;
    for (i = 2; i < 246913; i++)
    {
        if (box[i] == 0)
        {
            for (j = i + 1; j < 246913; j++)
            {
                if (j % i == 0 && box[j] == 0)
                    box[j] = 1;
            }
        }
    }
    return 0;
}

これを少しいじって，

int main()
{
    int num, num2;
    int box[246913] = {1, 1};
    int i, j;
    for (i = 2; i < 246913; i++)
    {
        if (box[i] == 0)
        {
            for (j = i + i; j < 246913; j += i)
                box[j] = 1;
        }
    }
    return 0;
}

にすると，大分速くなると思います．

解説1：なぜ従来の手法が遅いのか

従来の手法ですと，

for (j = i + 1; j < 246913; j++)
{
    if (j % i == 0 && box[j] == 0)
        box[j] = 1;
}

の部分で，j を1ずつ増やして，j % i == 0 && box[j] == 0が真のとき（つまり素数でないとき）にbox[j]を0にしていますね．
しかし，これだと，無駄が多いです．
というのも，下の画像を見てみてください．

この画像は従来の手法でした場合のfor文の中で回っている部分です．
緑の部分がfor文の中で回っている部分です．
質問者のhnk_llさんはおそらく，
「if文で条件を絞って，エラトステネスの篩を実装しよう」
と思ったのかもしれません．
しかし，このif文の条件文が計算量の多い計算を含んでいるので，すごい時間がかかっているのです．
ikadzuchiさんのおっしゃるとおり，除算や剰余はたいへん遅い処理なのです…

解説2：どのように改善すればよいか

さて，どのように改善すれば速くなるかということですが，
このif文を使わずに実装することです．

具体的には，
i = 2のとき，j = 4 6 8 10 12 ...のときにbox[j]を0にし，
i = 3のとき，j = 6 9 12 15 18 ...のときにbox[j]を0にし…
というふうにiの倍数（ただし1倍を除く）のときだけfor文を回るようにします．
そうすればif文は不要になります．

画像だとこうなります．（緑の部分がfor文の中で回っている部分です．）
緑の部分が減りましたね！

それが

for (j = i + 1; j < 246913; j++)
{
    box[j] = 1;
}

というふうにすれば速くなる理由です．

Comments

[grape_ll]

確かにこのようにすれば剰余などを使わずにできますし早いですね．仕組みも理解出来た気がします．ご回答ありがとうございました．

Answer 3

一番の問題はここですね。

j%i==0

除算や剰余はたいへん遅い処理ですので可能な限り避けるべきです。
エラトステネスの篩を実装するのが正攻法だと思います。

あと、

box[j]==0

この条件判定は無駄なので無くしたほうが速くなりそうな気がしたのですが、試してみると遅くなったりほぼ変わらなかったりですね…。

Comments

[grape_ll]

四則演算の計算速度についてはあまり気にしていなかったのですが，今度からはそういった視点でみ見てみるようにしたいと思います．ご指摘ありがとうございました．

Answer 4

エラトステネスの篩では高速化出来ませんか？
手持ちの機械では↓のように成ります。
usr ~/Project/NewProject % soe 246913
Time:0.000495

*** prime count: 21787 ***

≒5/10000秒です。
GitHubに上げている『エラトステネスの篩』です。
最初に見つかった非０を残してその倍数を（０で）消去しています。
最後には穴だらけの配列に成りますが・・・

/*
 * エラトステネスの篩本体 O(n log log n)
 */
void SieveOfEratosthenes(long arry[], long siz)
{
    long lim = (long)(ceil(sqrt(siz+2)));  // 最大値の平方根+αまで調べればいい
    for (long i = 2; i < lim; i++) {
        if (arry[i]) {
            long tmp = arry[i];  // 素数を取り出す
            for (long j = i + i; j < siz; j += tmp) {
                arry[j] = 0;  // 素数の倍数をクリア
            }
        }
    }
}

最後に別の配列を用意して、非０だけを抽出すれば、素数の表が作れます。

Comments

[cateye]

参考：https://github.com/neko2313/soe/blob/master/soe.c

[grape_ll]

wikipediaの動きの動画のようなものを見たのですが，知識がない私からでは私が書いたコードのようなことをしていると思ってしまったのですが，具体的にはどのように違いますか？もしよろしければ教えていただけますと幸いです．

[e-watt]

私が質問につけたコメントの（大部分が実装されているという）相違があります。

[grape_ll]

すいません，気づきませんでした．
確認させていただいて，上から二つ目のポイントまでは理解できたのですが，そこから下がよく分かりませんでした．ひとまず2の倍数は飛ばすという風に考えて書いてみたのですが，まだ足りていないと思いますのでご指摘いただければ幸いです．よろしくお願いいたします。
int box[246913]={1,1};
int i,j;
for(i=2;i<=123456;i++) box[2*i]=1;
for(i=2;i<=123456;i++){
if(box[2*i-1]==0){
for(j=i+1;j<=123456;j++){
if((2*j-1)%(2*i-1)==0&&box[2*j-1]==0) box[2*j-1]++;
}
}
}

[cateye]

ceil(sqrt(siz+2))などという面倒くさいことをしているのは、平方根の値が素数だった時に取りこぼしをしないためです。

[grape_ll]

難しいことではないのかもしれませんが，なにぶん経験がないものでして，しっかり考えてみたいと思います．ご指摘ありがとうございました．

Answer 5

素数の計算を行う場合、巨大な配列を扱うことになるかと思いますが、スタック・ヒープ・キャッシュの仕組みについて理解しておく必要があるかと思います。

スタックというのはCのローカル変数などを保持しておくための領域で、スタック領域のサイズは小さめにとってあります。なのでグローバル変数・ローカル変数で巨大な固定長の配列を用意するとコンパイラがエラーを出します。スタック領域のサイズは大きくすることもできますが、下述のキャッシュの仕組みにより推奨されません。

ヒープはスタックよりも多くのデータを扱うことができます。ただしヒープは動的にメモリを確保する必要があります。最近のPCはRAM(主記憶メモリ)を4GB以上搭載していることが多いのでヒープを大量に使ってもメモリ不足になることはあまりないかと思います。ヒープがRAMを食いつぶしてもOSにメモリスワップの仕組みがあるのでRAMの容量以上のメモリをプログラムが使用してもプログラムが停止することはないかと思います。ただしメモリスワップは急激なパフォーマンスの低下を招きます。

PCに搭載されているプロセッサ(CPU)にはキャッシュ(キャッシュメモリ)という仕組みがあります。キャッシュはデータの局在性を用いた高速化の仕組みです。スワップ領域のサイズが小さい理由の一つとしては、わざと近くのメモリ領域だけで計算をさせるという意図があります。ヒープを使うとデータの局在性が失われ、計算速度が低下します。

こんなことを言われてもわからないと思いますが、一つわかりやすい例をあげます。複雑なアルゴリズムのほうが簡単なアルゴリズムよりも高速に計算しそうに思われます。理論上はオーダー…つまりO(n)やO(n^2)などのこと、をさげるなど高速な計算に寄与しそうに思われますが、複雑なアルゴリズムはプログラムサイズと使用するデータ量を増加させやすく、結果高速化にはつながらないということがありえます。

素数の計算には様々な手法がありますが、理論的な計算速度の他にもデータの局在性による影響も絡んでいるということを頭の隅にいれておくとよいかと思います。

Comments

[cateye]

たしかにそうですね。私のエラトステネスの篩は（素数をそのまま格納するために）10億まででlong(8byte)×10億≒８Gをヒープから引っ張ってきます。
・・・事前に『封筒の裏の計算』はちゃんとしましょう。
最近のAMDのCPUのキャッシュラインは、1ラインあたりのバイト数はどれくらいなんだろう？

[grape_ll]

まだ知識が浅いので完全に理解することは出来ませんが，データの局在性というものも絡んでくるかもしれないということは覚えておこうと思います．ご回答ありがとうございました．