ぬあんと

素数10億まで3秒 - JGeek Log
danさんの場合, 1bit でフラグを記憶してるのでメモリが1/8 で済む。そこでメモリアクセスの時間が効いてるんだろう。それならキャッシュに収まる位のブロックに計算を分割しその内側で素数pのループ回せばもっと速くなるかも？と思いやってみた。見事3秒で終わった!

これを上回ることは出来るか?

出来ました。

以下、証拠。

% gcc -W -Wall -O6 -fopenmp primes_ita.c
primes_ita.c: In function ‘main’:
primes_ita.c:69: warning: unused variable ‘p’
% time ./a.out
06.983u 0.084s 0:01.54 458.4%	0+0k 0+5io 0pf+0w
% time ./a.out 9999000000000
10.223u 0.090s 0:02.05 502.9%	0+0k 1+3io 0pf+0w
% time ./a.out 0 10000000000
74.587u 0.906s 0:17.11 441.1%	0+0k 4+13io 0pf+0w

最後のものは、探査範囲を10億ではなく100億に広げた例。これくらいでないとActivity Monitorでバーが8つフルに延びるのを確認できないので。

Utilizationが400%を超えているのは、Core i7のHyperThreadingが効いてる証拠。

しかし今回の主題が、速さではなく早さにあることを示すのが、このdiff。

--- primes_ita.c.orig	2010-08-06 04:59:37.000000000 +0900
+++ primes_ita.c	2010-08-06 05:00:21.000000000 +0900
@@ -3,6 +3,7 @@
 #include <limits.h>
 #include <errno.h>
 #include <math.h>
+#include <omp.h>
 
 #include <stdint.h>
 typedef uint32_t U32;
@@ -68,6 +69,7 @@
     U64 p, i;
     U64 block_size = BLOCKSIZE*2;
 
+#pragma omp parallel shared(b,sieve)
     for(i = 0; i*block_size < size; i++)
     {
 	U64 size1 = block_size;

そう。OpenMPを使っただけ。

恐ろしいまでの手軽さ。

asin:4873114357
並行コンピューティング技法
Clay Breshears / 千住 治郎訳
[原著:The Art of Concurrency]

しかしこのOpenMP、以外に使われてないのは、多分いい資料がなかったから。というわけでこの場を借りて「並行コンピューティング技法」を紹介。以前献本いただいたのだけど、出すタイミングを逸していた。OpenMPって#pragmaつけるだけで既存のコードがマルチスレッドになるという優れものである一方、きちんと使わないと結果がへんてこになっちゃりかえって遅くなったり。たとえばこの例だと、shared(b,sieve)を指定するとコンパイルが通らない。そのあたりのことを知るのにこの本はうってつけ。余裕があれば「インテル スレッディング・ビルディング・ブロック」も目を通しておくといいかも。

OpenMP++

Dan the Casual Parallel Programmer

P.S. BLOCKSIZEはこの場合4MiBより、コア数で割った1MiBにした方がよさげ。L3は共有なので。

% time ./a.out 0 10000000000
42.609u 0.741s 0:12.48 347.2%	0+0k 0+11io 0pf+0w