Project Euler - Problem 433 を並列化して強引に解こう

問題文はこちら。

• http://projecteuler.net/problem=433
• http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%20433

E(x, y) = E(nx, ny) である事を利用し、互いに素な x, y のみを考えるとして、Stern-Brocot tree を巡回する方法で、Gauche を使ってゴリ押しするならばこのようになる。

N = 1e4 で 1m23.882s も掛かっているので、とてもではないがこれでは解けそうになく、E(x, y) の計算コストを削減しなければならない。そこで、E(x, y) を求める際の手順を逆にして考え、E(x, y) + 1 = E(nx+y, x) である事を利用する。また、x と y が互いに素なら nx+y, x も同様に互いに素であるので、逆算する要領で x, y を辿っていけば、互いに素なペアを全て網羅しつつ、E(x, y) の計算コストを殆ど無視出来る。これを C++ で書くとこのようになる。

このアルゴリズムを使った時の計算量は O(N^2) となっており、N = 1e5 の時 14.350s なので、N = 5e6 ならば約 36000 秒つまり約 10 時間程で処理を終えられそうである。

OpenMP で並列化してみた

さて、という事である程度の目処が付いたので、次は上記のプログラムを OpenMP を利用し並列化させてみよう。

6 コアの CPU を使っているので 6 スレッド使ったところ、シングルスレッドでは 14.350s 掛かっていたものが 3.890s で処理出来るようになった。が、しかし理想的に並列化出来たならば 14.350s / 6 = 2.3916666...s にもっと近似しても良さそうである。

main() から直接呼ばれる solve(n,1,5,0) それぞれについて計測を行うと、一番初めに呼ばれる solve(2,1,5,0) の実行時間が全体の約半分を占める事が分かる。逆算を行う際の初期値が小さければ小さい程、再帰呼び出しを行う回数が増えるからだ。そこで今度は、main() から直接行う最初の幾つかの solve() の呼び出しを、個別に並列化してみた。

このコードは solve(2,1,5,0), solve(3,1,5,0), solve(4,1,5,0) の 3 つをそれぞれ個別に並列処理するようにし、またループ内の処理を終えたスレッドは nowait で次のループへ取り掛かるように変更したもの。これにより N = 1e5 の時 1.941s で処理出来るようになり、3.890s より約 2 倍に高速化出来た。N = 5e6 でもぶん回せば 1 時間半以内に答えが出る。

という事で

• Core i7-4930K おすすめ
• OpenMP の nowait おすすめ

でした。

Project Euler - Problem 229 を強引に素早く解こう

問題文はこちら。

• http://projecteuler.net/problem=229
• http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%20229

内容は、2*10^9 までの整数のうち、以下の 4 つの式を成立させられる n はいくつあるかというもの。

• n = a₁² +   b₁²
• n = a₂² + 2 b₂²
• n = a₃² + 3 b₃²
• n = a₇² + 7 b₇²

この問題はとても簡単に力任せで解く事が可能で、2*10^9 個の要素を持つフラグテーブルを作って、a と b を網羅し対応する n のフラグを立て、全てのフラグが立っているものをカウントするだけで済んでしまう。

$ g++ -std=c++11 -O4 229.cpp; time ./a.out
********

real    1m0.754s
user    1m0.140s
sys     0m0.516s

ちなみに上記の計測に使った CPU は Core i7-4930K で、以前 Core2Duo E8500 で計測した時でも 4 分で処理出来ていた。

高速化してみた

先程のプログラムは、2*10^9 バイトのおよそ 2GB という広い領域を持つフラグテーブルを、a と b によって塗り潰すものとなっており、参照の局所性というものを全く考えていない。そこでそのテーブルを、キャッシュメモリに乗るだけの区間に分けて処理し、ついでに OpenMP を使いマルチスレッドで動作するように書き直した。

すると OpenMP を使わずシングルスレッドで動かした場合でも

$ g++ -std=c++11 -Wall -O4 229-segmented.cpp; time ./a.out 
229-segmented.cpp:44:0: 警告: #pragma omp parallel を無視します [-Wunknown-pragmas]
********

real    0m13.339s
user    0m13.328s
sys     0m0.000s

と、如何に DRAM が遅くキャッシュメモリが速いのかが分かる結果となった。 そして更に OpenMP を用いてマルチスレッドで処理させると

$ g++ -std=c++11 -Wall -O4 -fopenmp 229-segmented.cpp; time ./a.out
********

real    0m1.648s
user    0m19.632s
sys     0m0.004s

このようになる。というわけで

• 参照の局所性を考えるのおすすめ
• OpenMP を使ったお手軽な並列処理おすすめ

でした。