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 おすすめ

でした。