28
Help us understand the problem. What are the problem?

More than 3 years have passed since last update.

posted at

updated at

細かすぎて伝わらないエラトステネスの篩の高速化

5000兆回考えられたネタ。

エラトステネスの篩とは

指定された数以下の素数の一覧を求める単純なアルゴリズム。
計算量は $O(n\log\log n)$ らしい。

細かいことはググって。

ここでの条件

  • C言語で実装
  • エラトステネスの篩を使う
  • 素数の数を求める(100 なら 25)
  • コンパイラ最適化とかパイプラインとかには適度に期待する
  • 最大値はそれなりに小さいとする

計算時間は手元の Mac で測っている(コンパイラは clang)。コンパイラオプションは -O3 -march=native-march=native なのは popcnt と比較するため。

基本のコード

#include <stdio.h>
#define MAXVAL 200000000

int main() {
    int i, j, count = 0;
    static int flags[MAXVAL + 1] = {};
    for (i = 2; i <= MAXVAL; i++) {
        if (flags[i] == 0) {
            ++count;
            for (j = i * 2; j <= MAXVAL; j += i) {
                flags[j] = 1;
            }
        }
    }
    printf("%d\n", count);
    return 0;
}

5.54秒くらいかかる。おっせ〜。

高速化手法

外周の高速化

2 の倍数をとばす

あるあるの手法。2 の倍数を飛ばす場合、2 の倍数のフラグを立てなくて済むが、実はそこまで性能差はでない。

- int i, j, count = 0;
+ int i, j, count = 1;

- for (i = 2; i <= MAXVAL; i++) {
+ for (i = 3; i <= MAXVAL; i += 2) {

3 の倍数もとばす

2、4、2、4、… と飛ばしていけばよい。

- int i, j, count = 1;
+ int i, j, f, count = 2;

- for (i = 3; i <= MAXVAL; i += 2) {
+ for (i = 5, f = 4; i <= MAXVAL; i += (f = 6 - f)) {

最大値のルートまでしか回さない

100 までの場合、11 とか 13 とかについてはフラグを操作することがないので、これらを無視する。これで随分はやくなる。

#include <stdio.h>
#include <math.h>
#define MAXVAL 200000000

int main() {
    int i, j, f, count = 2;
    int max = (int)sqrt(MAXVAL) + 1;
    static int flags[MAXVAL + 1] = {};
    for (i = 5, f = 4; i <= max; i += (f = 6 - f)) {
        if (!flags[i]) {
            ++count;
            for (j = i * 2; j <= MAXVAL; j += i) {
                flags[j] = 1;
            }
        }
    }
    for (; i <= MAXVAL; i += (f = 6 - f))
        if (!flags[i])
            ++count;
    printf("%d\n", count);
    return 0;
}

内周の高速化

開始インデックスと 2 の倍数

冷静に考えると、i * i から計算すれば十分である。なぜなら、たとえば 7 の倍数を消すとき、14 とか 35 とかは、それぞれ 2 と 5 を見ているときに消しているはずだからだ。

あと、2、4、…の倍数はフラグを起こさなくていいので、これも飛ばす。

これでもかなり速くなる。

-             for (j = i * 2; j <= MAXVAL; j += i) {
+             for (j = i * (i | 1); j <= MAXVAL; j += i * 2) {

後ろから見る

知っている合成数 × 今見ている素数 のパターンは、既にフラグが立っているはずなので、これはフラグを起こさなくてよい。
ただ、下から見てしまうと、このフラグを勘違いしてしまうことがある。
たとえば、11 を見ているとき、121 の合成数フラグを立ててしまうと、1331 のフラグを立てなくていいと勘違いしてしまう。
これを避けるために、後ろから for 文を回す。

- int i, j, f, count = 2;
+ int i, j, f, s, count = 2;

-             for (j = i * 2; j <= MAXVAL; j += i) {
-                 flags[j] = 1;
+             s = MAXVAL / i;
+             for (j = s - !(s & 1); j >= i; j -= 2) {
+                 if (!flags[j]) flags[j * i] = 1;

ちなみにこの高速化は、フラグをビット演算すると無駄になる(キャッシュの方が勝るので、どんどん書いてしまう方が速い)。

フラグをビット演算化する

フラグをビット演算化する

フラグは 0 と 1 しかないんだから、ビット演算にしろ。かなり高速になる。

#include <stdio.h>
#include <math.h>
#define MAXVAL 200000000
#define BITS (sizeof(int) * 8)
#define FLAGS_NUM (MAXVAL / BITS + 1)

static int flags[FLAGS_NUM] = {};

inline static void upflag(int i) {
    flags[i / BITS] |= 1 << (i % BITS);
}

inline static int getflag(int i) {
    return (flags[i / BITS] >> (i % BITS)) & 1;
}

int main() {
    int i, j, f, s, count = 2;
    int max = (int)sqrt(MAXVAL) + 1;
    for (i = 5, f = 4; i <= max; i += (f = 6 - f)) {
        if (!getflag(i)) {
            ++count;
            s = MAXVAL / i;
            for (j = s - !(s & 1); j >= i; j -= 2) {
                if (!getflag(i))
                    upflag(j * i);
            }
        }
    }
    for (; i <= MAXVAL; i += (f = 6 - f))
        if (!getflag(i))
            ++count;
    printf("%d\n", count);
    return 0;
}

popcnt を使う

1 になっているビットの数を数える関数。2 と 3 の倍数のフラグを立てる必要があるので、これはテーブルを作ってしまう。

平方根までは count をインクリメントしてもいいかな、と思ったけどどっちが速いんだろう。

count = -2 は、0・2・3 を合成数扱いしてしまっているのと、1 を素数扱いしてしまっている分の埋め合わせ。


#include <stdio.h>
#include <math.h>
#include <x86intrin.h>
#define MAXVAL 200000000
#define BITS (sizeof(int) * 8)
#define FLAGS_NUM (MAXVAL / BITS + 1)

static int defaults[] = { 0x5D75D75D, 0xD75D75D7, 0x75D75D75 };
static int flags[FLAGS_NUM];

inline static void upflag(int i) {
    flags[i / BITS] |= 1 << (i % BITS);
}

inline static int getflag(int i) {
    return (flags[i / BITS] >> (i % BITS)) & 1;
}

int main() {
    int i, j, f, count = -2;
    int max = (int)sqrt(MAXVAL) + 1;

    for (i = 0; i < FLAGS_NUM; i++)
        flags[i] = defaults[i % 3];
    flags[FLAGS_NUM - 1] &= (1 << (MAXVAL % BITS + 1)) - 1;

    for (i = 5, f = 4; i <= max; i += (f = 6 - f)) {
        if (!getflag(i)) {
            for (j = i * (i | 1); j <= MAXVAL; j += i * 2)
                upflag(j);
        }
    }

    for (i = 0; i < FLAGS_NUM; i++)
        count += _popcnt32(flags[i]);

    printf("%d\n", MAXVAL - count);
    return 0;
}

手元の環境では、これで 0.41 秒までになった。

そのほか

あまり意味のなかった手法

  • 数値型を変える
    • unsigned にしたり 64bit にしたりしたけど、ほとんど変わらなかった
    • 組み合わせによっては若干速くなることがある
    • コンパイラ最適化の問題?
  • flags の初期化
    • 実は {} はいらない
    • とるのは未初期化っぽい雰囲気が出るのであまり好きじゃない(わがまま)
    • ちなみにあってもなくてもほとんど変わらない(ないほうがほんの少し速い)

ダメだった手法

  • while 文にする
    • for 文より while 文の方が速いことがあるが、スクリプト言語じゃないしまあ無意味
    • むしろ遅い
  • 平方根の改良
    • math.h が意外と速い
    • 開平法を実装しても無意味だった

まだやってない手法

  • 5 の倍数も飛ばす
    • 計算が複雑化して遅くなりそう(?)
    • 飛ばし方が 4, 2, 4, 2, 4, 6, 2, 6 となる(煩雑!)
    • このパターンを高速に作れれば速くなるかも
  • 2 の倍数のフラグを作らない
    • 直接インデックスでアクセスできなくなって遅くなりそう(?)

まとめ

環境によっては選び方によって違うかも。なんか他にアイデアがあったらください。

Special Thanks: @dsk_saito@D_Plius@KRiver1@lpha_z@eukaryo

Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
Sign upLogin
28
Help us understand the problem. What are the problem?