魔方陣の解の数を一瞬で求めるプログラム

以前、「スパコンで約2時間36分かかったという、5×5の魔方陣の全解列挙を、パソコンで試す(C++)」という記事を書きましたが、解を列挙するのではなく、解の数を数えるだけなら一瞬で終わらせるプログラムがあります。

4×4の魔方陣の場合

clang++ -O3 -std=c++11 magicsquare4.cc」などとしてコンパイルしてから実行します。実行時間は一瞬です。

このプログラムの欠点は、コンパイルにかなりの時間がかかることです(TANSTAAFL)。手元の環境では、コンパイルに4分くらいかかりました(説明するのは野暮ですが、これはネタです)。

C++のconstexprを使っています。複数のコンパイラがconstexprをサポートしていますが、このコードをコンパイルできるのは、私が試したところではClang 3.3のみです(Ubuntuでは「sudo apt-get install clang-3.3」などとしてインストールできます)。ステップ数に上限が導入されたClang 3.4以降や、途中経過を記憶してメモリを圧迫するgccではコンパイルできません。

Clang 3.3でサポートされるconstexprの関数には、大ざっぱに言えばreturn文しか書けないので、マスに数が入るかどうかをチェックする関数を作り、「数がすでに使われていなければ次のマスを試し、使われていれば1だけ大きい数を再帰的に試す」ということを、return文の条件演算子で実現しています。

このコードは私が直接書いたものではなく、私が書いたプログラムが生成したものです(手で直接書くのは大変です)。同じプログラムで5×5の場合も生成しましたが、試すのはやめた方がいいでしょう(constexprを消してコンパイル・実行すれば、正しいことは確認できます)。

スパコンで約2時間36分かかったという、5×5の魔方陣の全解列挙を、パソコンで試す(C++)

追記:対称性の活用法をコメントで教えていただきました。それを採用すると本稿の結果よりもかなり高速になります(10分→2分30秒。解を数えるだけなら約90秒)。(コード。本文にはまだ反映していません。)

筑波大学のスパコン「T2K-Tsukuba」で約2時間36分かけて5×5の魔方陣(≠魔法陣)の全解を求めたというニュースがありました。

「スパコン」ではなく「パソコン」だったらどうなのだろう、と思って試してみたら、10分でした(Core i7 4930K)。

解の「数」がわかればいいだけ、つまり全解をディスクに書き出さなくてももいいならもう少し早くなります。

コードはこちら(最初のバージョン)。頭の悪そうなコードですが、スクリプトで生成した枠組みを手直しして作れば、そんなに大変ではないかと。理想は全自動生成ですが。Visual Studio 2013 ExpressとICC 14.0.2、GCC 4.6.3、Clang 3.4で動作を確認しています。ClangはOpenMPに対応させるのが面倒なうえに、ちょっと試したところでは、GCCより遅かったです。(実際に試す場合は、OpenMPを有効にしてください。計算中のCPU使用率は100%になるはずです。)

多重for文(笑)ですべての場合を探索するというのが基本方針ですが、それではさすがに終わらないので、次のような工夫をします(工夫というほどのものではありませんが)。

  • 各行・各列・対角線の和は65。(1から25までの和 / 5)
  • (カンで)効率が良くなりそうなところから順番に埋めていく(下図の番号順)。
  • 回転と裏返しで同一になるものを重複して数えないように、下図の緑部分に、「左上<右上<左下」かつ「左上<右下」のような制約を入れる。左上は22以下としてよい。
  • a+b+c+d+eのうち、a,b,c,dが決まったら、e=65-a-b-c-dになる(下図の紫部分)。
  • a+b+c+d+eのうち、a,b,cが決まったら、dはmax(1,65-a-b-c-25)からmin(25,65-a-b-c-1)の範囲で調べればよい。max(un,65-a-b-c-ux)からmin(ux,65-a-b-c-un)にするとさらに速くなる(unは未使用の最小数,uxは未使用の最大数。これを調べるのに時間をかけると効果はなくなる)。

この方針の場合、ループを関数で表現したり、盤面をコンテナで表現したりすると遅くなります。

数を埋める順番
6 10 11 12 8
20 1 18 2 21
22 13 5 16 23
24 3 19 4 25
9 14 17 15 7

魔方陣の解の列挙は並列化しやすそうな問題ですが、ここでの方針では、探索効率を上げるためには条件分岐が不可欠なので、(「数」を求めるだけだとしても)GPGPUでうまくやる方法がわかりません。そこで、CPUに載っているコアのみで並列化します(Xeon Phiなら簡単なのでしょうか→追記参照)。

一番外側の、0から(1<<25)-1まで変化する変数iのループをOpenMPで並列化します(schedule(guided)では遅くなります。schedule(auto)はVisual C++でサポートされたら試します)。変数iは上の図の緑の部分(カンで5個にしました)を各数5ビットで表現し、つなげたものです。マスに入りうる数は1から25までなので、5ビットというのはちょっと冗長ですが、とりあえずはよしとしましょう。

出力はバイナリ形式で、1つの解に25バイト使います(1つのマスに入る数を1バイトで表現する)。これもちょっと冗長ですが、テキスト形式よりは小さくて速いはずなので、とりあえずはよしとしましょう。

解は全部で2億7530万5224個、1つの解を25バイトで表現するので、(スレッド数と同数できる)出力ファイルのサイズの合計は275305224 x 25 = 6882630600バイト(約6.4GB)になります。

最初のバージョンの実行時間は以下の通りです。最終的には、先述の通り、もう少し速くなります。

Core i7 4930K(12スレッド・出力先はRAMディスク)
  • Windows 7 64 bit, Visual Studio 2013 Express, 593秒(約10分)(ターゲットはx64。/openmp /Ox /arch:AVX
  • Fedora 20 64 bit, Intel C++ Compiler 14.0.2, 471秒(約8分)-fopenmp -O3 -xHost
  • Fedora 20 64 bit, GCC 4.8.2, 583秒(約10分)-fopenmp -O3 -march=native
Core i7 4700MQ(8スレッド・出力先はSSD)
  • Windows 7 64 bit, Visual Studio 2013 Express, 1156秒(約20分)/openmp /Ox /arch:AVX

4535786569大森清美『魔方陣の世界』(日本評論社, 2013)(参考文献リストあり・索引あり)を見たら、a+b+c+d+e=65のとき、(26-a)+(26-b)+(26-c)+(26-d)+(26-e)=65なので、解の、すべての数xを(26-x)に置き換えたものもやはり解になるということが載っていました。コメントでご指摘いただいたように、この知識を使って、真ん中のセルに入る数字を1から13までに限定し、そこが13ではない解が見つかったら、すべての数xを(26-x)に置き換えたものも解として出力するようにするとさらに速くなります(同条件でちょうど6分)。

『魔方陣の世界』では、この問題のためのコードも紹介されていますが、並列化されていないこともあって、解の「数」を調べる場合は約5000秒(確認済み)、全解を出力する場合は約3日(未確認)と、あまりふるいません。

この問題は、パソコン・プログラミングに最適の題材である。(p.111)

今日では、5次方陣の総数検索問題は、パソコンに最適な題材となった。(p.124)

パソコンに最適とも、スパコンに不適とも思いませんが、興味深い題材であることは確かです。

追記

  • 正解の数(2億7530万5224個)を知っているから速いのか? 正解の数は1981年にはすでに知られていましたが、ここで試した方法では、その知識は使っていません。ディスクの空き容量は小さくてよいことがわかっているとか、下に書くような理由で、正解が既知だと試すのが気楽だということはあります。
  • 結果の正しさをどうやって確認したのか? 確認していません。解の数は既知なので、それと同じ数が出てくればまあいいかなと。正しさを確認するためには、(1)得られたものが本当に解であること、(2)結果に重複がないこと、(3)すべての場合を調べていること、の確認が必要です。(1)はループの最も内側で解になっているかどうかをチェックすることによって、(2)は得られた解をハッシュテーブル等に格納することによって確認できますが、ある程度のメモリと時間(上のマシンだと約200秒)が余分に必要です(コード)。(3)はここでは難しいですね。
  • アルゴリズム? 穴埋め問題の穴を多重for文(笑)で埋めているだけなので、アルゴリズムというほどのものはありません。(参考:Puzzles for Hackers
  • 10分くらいならやってみようかなという人が、きれいな書き方やGPUで解く方法、効率的な方法を見つけてくれればと思います。(参考:Toy problemsは役立たずか
  • Xeon Phiなら(ほぼ)そのままのコードが動いて、計算は5分で終わるそうです。(スパコンで約2時間36分かかったという、5×5の魔方陣の全解列挙をXeon Phiで試す – パラレルに恋して

Toy problemsは役立たずか(Floydの問題・ネタバレ)

U = {√1, √2, …, √50}を2つのグループに分け、グループ毎の和を求める。和の差がなるべく小さくなるようなグループ分けを、制限時間10秒で見つけよ。

1881526917先日紹介したFloyd問題です。

10秒という制限時間で、2の49乗個の候補の中から最良のものを見つけ出すのは、1977年当時のコンピュータでは難しかったはずです。そこでKnuthは、以下のような工夫をしました。

  1. まず小数部分を近づけることを目指し、その後で整数部分を近づける。(この問題では結果オーライ。√40までの場合はうまくいかない。)
  2. I = {√1, √4, √9, √16, √25, √36, √49}は、和の整数部分にしか関与しないから後で考えることにする。
  3. Uを2つのグループAとBに分け、Aの部分集合とBの部分集合の合計の小数部分を、全体の合計の半分の小数部分に近づける。
  4. Aの部分集合の小数部分をすべて記録しておき、Bの部分集合に対して、上の条件に会うようなAの部分集合がすぐに見つかるようにする。
  5. 部分集合の生成にはグレイコードを使う。
  6. X = {√2, √8, √18, √32, √50}の組み合わせでできる和は16通りしかないことを利用する。同様に、Y = {√3, √12, √27, √48}の組み合わせでできる和は11通りしかないことも利用する。
  7. 集合をA = X ∪ Y ∪ {√5, √6, √7, √10, √11, √13, √14, √15}、B = U – I – Aとする。

今日のコンピュータは当時と比べればはるかに強力なので、これらの工夫のうちの1から4を採用するだけで、あとは富豪的なプログラムでも最良解を見つけるのに5秒とかかりません(コード)。本稿執筆時点のideone.comで約2.5秒(double)あるいは約2.8秒(long double)、私のデスクトップ(Core i7 950 3.07Ghz)のVisual C++で約1.3秒でした(いずれもシングルスレッド)。

4774157155Aの部分集合の記録には、C++11で導入されたunordered_multimapを使いました。キーAのサブセット和の小数部分8桁を、値はサブセットを表現する整数とサブセット和のpairです(この実装では、mapよりunordered_mapのほうがかなり速いです)。sizeAを20にしたり、keyを小数点以下8桁にすると速いのは、動かしてみてわかったことです。

OpenCLに対応するデバイスの列挙(C言語・Mathematica)

GPUをグラフィック処理ではなく汎用計算に利用しようというGPGPUのためには、CUDAかOpenCLを利用するのが一般的です。NVIDIA的には、CUDAはGPGPUの開発環境であり、プログラミング言語としてC for CUDAかOpenCLを選べる、つまりCUDAはOpenCLの上位概念らしいのですが、一般にはCUDAとOpenCLは対立するものとして認識されているような気がします。

CUDAはNVIDIAのGPUにしか対応していないのに対して、OpenCLはNVIDIAのGPUとAMDのGPU、IntelとAMDのCPUにも対応しているので便利です(対応していないGPUやCPUもあります)。

IntelのCPUとHD GraphicsでOpenCLを利用できるようにする、Intel SDK for OpenCL Applicationsの新版(2013)が出たので、ちょっと使ってみます。

OpenCLに対応したデバイスを列挙することで、複数のデバイスに対応しているというOpenCLの特徴を確認してみましょう。規格で定まっているわけではないようですが、複数のOpenCL環境をインストールしているときには、複数のOpenCLプラットフォームがOpenCLのAPIから認識できるようになるようです(後述の参考書p.71)。

実行するためには、以下のいずれかが必要です。

4844331728開発は、WindowsとGNU/Linux、Macのいずれでも可能です。株式会社フィックスターズ『OpenCL入門 1.2対応 マルチコアCPU・GPUのための並列プログラミング』(インプレスジャパン, 改訂新版, 2012)サポートサイト)などを参考にすると、開発環境を比較的簡単に構築できます。(参考:Visual Studio + Intel SDK for OpenCL

適当な開発環境を用意したら、以下のプログラムをビルド・実行します(RSSリーダーでは見られないかもしれません)。

OpenCLの実装は複数のプラットフォームを認識できるものになっているので、上のコードでは、プラットフォームを列挙しつつ、各プラットフォームが持つデバイスを列挙しています。

私のデスクトップPCでの実行結果はこんな感じです。プラットフォーム毎にデバイスが1つあるというわかりやすい構成です。

Platform: NVIDIA CUDA
CL_PLATFORM_VERSION: OpenCL 1.1 CUDA 4.2.1

  Device: GeForce GTX 580
  CL_DEVICE_VERSION: OpenCL 1.1 CUDA
------------------------------------------------------------------------------
Platform: Intel(R) OpenCL
CL_PLATFORM_VERSION: OpenCL 1.2

  Device: Intel(R) Core(TM) i7 CPU         950  @ 3.07GHz
  CL_DEVICE_VERSION: OpenCL 1.2 (Build 63463)
------------------------------------------------------------------------------

私のノートPCでの実行結果はこんな感じです。2番目のプラットフォームにはデバイスが2つあり、そのうち1つは1番目のプラットフォームのデバイスと同じという、わかりにくい構成になっています。

Platform: AMD Accelerated Parallel Processing
CL_PLATFORM_VERSION: OpenCL 1.2 AMD-APP (923.1)

  Device:       Intel(R) Core(TM) i7-3612QM CPU @ 2.10GHz
  CL_DEVICE_VERSION: OpenCL 1.2 AMD-APP (923.1)
------------------------------------------------------------------------------
Platform: Intel(R) OpenCL
CL_PLATFORM_VERSION: OpenCL 1.2

  Device:       Intel(R) Core(TM) i7-3612QM CPU @ 2.10GHz
  CL_DEVICE_VERSION: OpenCL 1.2 (Build 63463)

  Device: Intel(R) HD Graphics 4000
  CL_DEVICE_VERSION: OpenCL 1.1
------------------------------------------------------------------------------

プラットフォームについて得られる情報は、http://www.khronos.org/registry/cl/sdk/1.0/docs/man/xhtml/clGetPlatformInfo.htmlにまとまっています。

デバイスについて得られる情報は、http://www.khronos.org/registry/cl/sdk/1.0/docs/man/xhtml/clGetDeviceInfo.htmlにまとまっています。

Mathematicaなら、すべての情報を以下の2行で得られて便利です(参考:OpenCLInformation)。

Needs["OpenCLLink`"]
OpenCLInformation[]

マンデルブロ集合をOpenCLで描く(Mathematica)

CUDAを使ってマンデルブロ集合の描画を3桁速くする方法を以前紹介したのですが、同じことをOpenCLでやってみます。NVIDIAのGPUを搭載していないノートPCを使うことが多くなった自分用のメモでもあります。書き換え方は、OpenCLLink プログラミングで紹介されています。

書き換えたコードは以下の通りです(RSSリーダーでは表示されないかもしれません)。

最初に作成したコードをCore i7-3612QM 2.1GHz上で実行するのと比べて、ここで作成したコードをIntel HD Graphics 4000上で実行すると、3桁くらい速くなります。

Manipulateを使ってインタラクティブに描くコードは以下の通りです(RSSリーダーでは表示されないかもしれません)。

似たような話がMathematicaのマニュアルにも載っていますが、ここで書いたコードなら、描画の中心もインタラクティブに変わります。

realtime