Grand Central Dispatchが面白そうだったのでやってみた

Google Code Jam 2009の予選は、なぜか土曜日だと勘違いしていて、実際はもろに仕事中だったので参加できなかった、という言い訳。

さて、我が家でも発売日からSnow Leopardが導入された。
Snow LeopardといえばOpenCLが気になるところだけど、このMacBookにはGeforce積んでないのでまだOpenCLは試せない。
実はIONなATOMマシンを自作したので、いずれLinuxで試したい。

並列実行でもうひとつ実装されたのが、Grand Central Dispatch。
ASCII.jpの記事がとても分かりやすい例を挙げてくれている。

マルチコア時代の新機軸！ Snow LeopardのGCD
http://ascii.jp/elem/000/000/455/455786/

自分でも試したい、ということで円周率を計算してみた。
アルゴリズムはWikipediaに載っていたインドで発明された？方法を採用。

円周率 - Wikipedia
http://ja.wikipedia.org/wiki/%E5%86%86%E5%91%A8%E7%8E%87
14世紀インドの数学者・天文学者であるサンガマグラマのマドハヴァは次のような π の無限級数表示を見いだしている……

まずは比較するために、並列化していないコード。

include <stdio.h>
#include <math.h>
#include <sys/time.h>

#define MAX 500*1000*1000

double getTimeSecond()
{
        struct timeval tv;
        gettimeofday(&tv, NULL);
        return tv.tv_sec + (double)tv.tv_usec * 1e-6;
}

int main() {
        struct timeval startTime, endTime;
        double pi = 0.0;
        double t1, t2;
        int n;

        t1 = getTimeSecond();
        for (n = 0; n < MAX; ++n) {
                pi += pow(-1, n) * (1.0 / (2 * n + 1));
        }
        pi *= 4;
        t2 = getTimeSecond();

        printf("%f\t%f\n", pi, t2 - t1);

        return 0;
}

これを三回ぐらい実行してみた。

π	秒
3.141593	24.381985
3.141593	24.190699
3.141593	24.135658

だいたい24秒ぐらい。
次に、GCDを使ったコード。
計算する範囲は同じだが、それを1000分割して並列化している。

include <stdio.h>
#include <math.h>
#include <sys/time.h>
#include <dispatch/dispatch.h>

#define MAX 500*1000*1000
#define DIV 1000

double getTimeSecond()
{
        struct timeval tv;
        gettimeofday(&tv, NULL);
        return tv.tv_sec + (double)tv.tv_usec * 1e-6;
}

int main() {
        dispatch_group_t g = dispatch_group_create();
        dispatch_queue_t q = dispatch_get_global_queue(DISPATCH_QUEUE_PRIORITY_DEFAULT, 0);
        double result[MAX / DIV];

        struct timeval startTime, endTime;
        double t1, t2;

        t1 = getTimeSecond();
        int i;
        for (i = 0; i < MAX / DIV; ++i) {
                dispatch_group_async(g, q,
                        ^{
                                int j = i;
                                double p = 0.0;
                                int n;
                                for (n = j * DIV; n < j + DIV; ++n) {
                                        p += pow(-1, n) * (1.0 / (2 * n + 1));
                                }
                                result[j] = p;
                        }
                );
        }

        dispatch_group_wait(g, DISPATCH_TIME_FOREVER);

        double pi = 0.0;
        for (i = 0; i < MAX / DIV; ++i) {
                pi += result[i];
        }
        pi *= 4;
        t2 = getTimeSecond();

        printf("%f\t%f\n", pi, t2 - t1);

        return 0;
}

同じく、三回ぐらい実行してみる。

π	秒
3.142592	0.313952
3.142592	0.317246
3.142592	0.314269

むちゃくちゃ速い！
けど、なぜか精度が悪くなっている……。
足し算の順番が違うので、結果が異なることはもちろんありえるけど、ちょっと大きすぎる？
DIVを100にすると、SEGVで落ちたりするので、どこかバグってるのかも。
あとで確認しよう。

GCD簡単で面白いなぁ。