スパコンで約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で試す – パラレルに恋して

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

  1. あまり詳しくはない者ですが、元のニュースを読んだ時に「高校生なのに凄いな~」と印象に残っていたので興味深く読ませて頂きました。
    本ブログの内容ですと、
    ①5×5の全解を求めるのは(工夫は必要だけど)それほど難しいことではない
    ②スパコンを使わなくても普通のPCでできる、それも短時間に
    と思えるので、元記事のインパクトがないような気がしたのですが、そうなのでしょうか? 高校生がスパコンを使いこなして頑張っているよ!というレベルなのでしょうか。それとも、やはり記事に取り上げられるな何かインパクトやオリジナリティーがあるのでしょうか。結果だけ見るとYABUKI Taroさんの方が全然凄いような気がするのですが。

    • ①は時間をあまり気にしなくていいならそのとおりです。

      ②は「短時間」をどう考えるかによります。
      この記事のような単純な方法でも10分でできるので、
      アルゴリズムを考えるのが得意な方にとっては、
      1分とか10秒とかの話になるのではないでしょうか。

      「高校生が・・・」というのは問題の本質とは無関係ですが、
      注目を集める効果はあるでしょう。

      「スパコン・・・」というこの記事の書き方も、
      同じ効果を狙ったもので、褒められるものではありません。

  2. 「解の、すべての数xを(26-x)に置き換えたもの」は回転・対称解の方に入っているのではないでしょうか

  3. ピンバック: 2014/03/16(Sun) | gattya.run

  4. 自分でもプログラムを書いてみて、i7-4.4GHzで、TDE-GCC481で、35分かかりました。コア数を考慮しても自分のプログラムは、2倍以上遅いですね。(キャッシュ量の差が効くのかは不明ですが。)

    以下の部分ですが、
    >(3)すべての場合を調べていること、の確認が必要です。
    ですが、意味がよくわかりませんでした。条件を満たさない物の枝狩りは問題ないと思うのですが、多重ループを廻していて、正しい枝狩りであれば、すべての場合を調べている事にはならないのでしょうか?

    自分の場合、異なるアルゴリズムを使っており、対象解を除く判定に時間がかかっています。もう一度挑戦してみます。

    これからコードを読ましていただくのですが、中央の数字に注目して、たとえば中央が1の物と25の物は同数あり、どちらか一方が求まれば、他方は生成できますが、このような性質は使っていないと思うのですが、如何でしょうか。

    では、失礼します。

    • masaさん

      >>(3)すべての場合を調べていること、の確認が必要です。
      >ですが、意味がよくわかりませんでした。

      枝刈りの正しさに疑問を持たれたときに、それをどう証明するかという意味でした。

      >中央の数字に注目して、たとえば中央が1の物と25の物は同数あり、どちらか一方が求まれば、他方は生成できますが、このような性質は使っていないと思うのですが、如何でしょうか。

      その知識で省けるものは、対称性の考慮によって省けているので、あまり効果はないと思います(本文中に書いたように、一応試してはみました)。

  5. 早速のご返信有難うございます。
    1個目を中央に置く場合、1個目は、1~25ではなく、1~13、または、13~25の場合を調べれば個数を求められると思うのですが、如何でしょうか。
    中央の数字を i として その場合の個数をn(i) とすると、n(i)= n(26-i) ですので、総数=n(1)+n(2)+…+n(25)= ( n(1)+n(2)+…n(12) )x2 + n(13) となります。
    個数のみを数えるのであれば、計算時間は、1/2強になりそうです。

    ソースコード拝見しました。速いです。PIIx6 3.8GHz 6コアで、13分強で終わりました。
    自分のプログラムより、2倍以上速いです。自分も、もう少しがんばって見ます。

    • masaさん

      その件については、本文中、大森清美『魔方陣の世界』を紹介したあたりにも書きましたが、xを(26-x)に置き換えたものは、対称性を考慮して省いたたものと重複するのです。

    • 正しいと思います。
      中央のマスを1から12までに制限して数えた解数と、14から25に制限して数えた解数と同じになるので、重複して数える必要はないと言えます。
      なお、中央が1から12までの解はどのように回転と対象の変換を行っても、中央値が14から25のいずれかになることはありません。

    • masaさん・酒井さん

      ご指摘の点を考慮したら、10分だった計算時間が6分になりました。この結果を踏まえて、本文を更新しました。

      • こんばんは。
        inquisitorさんの対称解の除去の判断が効率的である事に気が付き、採用させていただきました。
        (以前は、自分のアルゴリズムには適さないと思っていました。)

        数字を置く順番を含め変更し、その結果 436秒から、242秒 と大幅に短縮できました。

        HDが遅いので、1個目(中心)と2個目(左上)に置く数値とその場合の解の個数(312個)のみ記録しています。1個目は、1から13に限定しています。
        CPUは、i7-ivy 4.3GHz です。コンパイラは、TDM-GCC 4.9.2 を使いました。

        以前は、対称解の除去に、特定の3つの数字の位置関係を使っていたのですが、コードが複雑で、速度低下の原因になっていたようです。

        自己満足にすぎませんが満足できる結果となりました。

        では、では。

        • masaさん

          コンパイラによる違いはいろいろあるみたいですね。

          解を出力することについては、元ネタにならっただけです。
          その部分にも工夫の余地はあるかと思います。

  6. 酒井様、コードをダウンロードいたしました。参考にさせていただきます。

    自分は、HP等持っていないので コードをご紹介できないのですが、i7-has 4.5GHz で 8 threads で、196 秒でした。i7-ivy と i7-has で、同一クロックでも差がありました。

    現在は、アルゴリズムの改良を中止し、クラスタシステムで実行するためのインターフェースのコードを書き始めました。(物量で勝負ですね。) 4PC,32スレッドとかで動かすことが目標です。

    そのあと、時間があれば、CUDA に移植しようと思っています。

    • PC 3台でクラスタ構成で実行する事に成功しました。
      PC1: Fx9370 4.4GHz master x1, slave x8
      PC2: i7-ivy 4.3GHz slave x 8
      PC3: i7-has 4.45GHz slave x 8
      以上の構成で、90.2秒で終了しました。

      Z:\5×5-801-Claster>mt702-n-13-01-clstr.exe 24
      *** n1= 24
      Pass= z:\temp00.grd
      Pass= z:\temp01.grd
      (中略)
      Pass= z:\temp22.grd
      Pass= z:\temp23.grd
      13/13
      12/13
      11/13
      (中略)
      3/13
      2/13
      1/13
      Ans= 275305224
      Time= 90.199 sec
      Z:\5×5-801-Claster>
      Windowsの共有ファイル設定の勉強になりました。
      今回は、slaveを24個手動で起動するのが大変でしたが、OpenMPを使えば、PC毎に1コマンドで起動できそうな気がします。

      では、では。   masa

      • PC毎にSlaveを必要な数だけ一括起動できる様になったので、5 PCでの動作を試してみました。Slave番号の割り当てを変えると実行時間が変わります。以下は、3回通り試し最も時間の短かった場合の結果です。
        i7-4770 4.45 GHz Slave 0 to 7
        i7-3770 4.3 GHz Slave 8 to 15
        Fx-9370 4.4 GHz Slave16 to 23, Master
        Fx-8350 4.0-4.1 GHz Slave 24 to 31
        i7-4702MQ 2.9 GHz Slave 32 to 40 ( Note PC )
        Ans= 275305224
        Time= 63.241 sec
        Start time= Wed Feb 25 00:57:43 2015
        End time= Wed Feb 25 00:58:46 2015
        もう一台PCを増やせば、60秒を切れそうです。(Fx-8150がパーツとして有ります。。。)
        masa

コメントは停止中です。