受験の役に立たない（と思う）数学シリーズ
 

５×５の魔方陣の総数を求めるプログラム

アルゴリズムとプログラム自体は大森清美氏の「魔方陣」という書籍にも載っているのですが、ここでは、私の考えたより高速でかつ能率のよいアルゴリズムを紹介します。工学系なので、こういうのは得意分野です。

基本的に魔方陣の総数を勘定するには、しらみつぶしに頼るしかないようです。 最初に説明を簡単にするために用語を定義します。

定義

魔方陣の各マスに次のようにラベルを付けます。

00-11-22-33-44のラインを対角要素と呼ぶことにします。逆の04-13-22-31-40のラインを副対角要素と呼ぶことにします。 5次の魔方陣の各ラインの和は60です。0〜24の数字を5個選び足して60になる組み合わせは 

0 1 12 23 24
0 1 13 22 24
0 1 14 21 24
....
10 11 12 13 14
の1394種類です。これを「組み」と呼ぶことにし、 
[0 1 12 23 24]
のように表現します。「組み」の段階ではその並びは重要であはありません。「どの数字が含まれているか」だけが重要です。そこで組みとして表現するときは小さい順番に数字を並べて書くことにします。

次にこれらの組みに順番をつけます。順番はどうでもよいのですが、ここではプログラム的楽になるよう、次のようなルールで順番をつけることにします。

ルール

  1. 組みの数字を小さい順番に並べたとき、一番目の数字が小さいほうが順番的に先となる
  2. 一番目の数字が等しい場合は二番目の数字で比較する。二番目が等しい場合は三番目の数字で比較する。
以上で全て順番が定義されます。 順番が定義できたので、組みの大小比較を 
[0 2 13 22 23] < [1 5 14 19 21]
のように不等号で表現することにします。

アルゴリズムの流れ

5次の魔方陣の総数を求めるアルゴリズムは何種類かあるようですが、ここでは 
1) 中央の数字が0の5次の魔方陣の総数
2) 中央の数字が1の5次の魔方陣の総数
....
25) 中央の数字が24の5次の魔方陣の総数
を求めその合計を求めるというアルゴリズムを採用します。このアルゴリズムはＰＣが複数台あれば並列処理が可能であり、高速化が期待できるという利点があります。ところで、中央の数字が24の魔方陣の総数は中央の数字が0の総数に等しいという事実があります。というのも 
0<->24
1<->23
2<->22
3<->21
という置換を行えば、その結果は魔方陣に必ずなるからです。そこで実際には 
1) 中央の数字が0の5次の魔方陣の総数
2) 中央の数字が1の5次の魔方陣の総数
....
13) 中央の数字が12の5次の魔方陣の総数
を求めれば十分で，よって，このアルゴリズムは最高13台のPCで並列処理可能ということになります。

アルゴリズムの流れ２

それでは、中央(22)の数字が Nであるような魔方陣を探すアルゴリズムを説明します。 まず、最初に互いに共通の数字(独立と呼ぶ)のない 5つのユニークな組み、C0〜C4を選びます。その際C0は Nを含むものに限定します。例えば N=3ならば 
C0 [3 10 14 15 18]
C1 [0 6 11 20 23]
C2 [1 8 12 17 22]
C3 [2 5 13 19 21]
C4 [4 7 9 16 24]
などがその例です。 次にC0〜C4のそれそれと共通要素を1つしか持たない組み(交差と呼ぶ)でかつ、互いに独立であるような組みL0〜L4を選びます。L0はだけはNを含むものとします。上のC0〜C4に対してこの条件を満足するものとしては例えば以下のようになります。 
L0 [3 11 13 16 17]
L1 [0 6 11 20 23]
L2 [1 8 12 17 22]
L3 [2 5 13 19 21]
L4 [5 8 9 18 20]
次にC0〜C4、L0〜L4の全てと交差する組みD0を探します。D0はN を含んでいなければなりません。上のC0〜C4、L0〜L4に対して満足するものとしては例えば以下のようになります。 
D0 [3 5 6 22 24]
D0を適当に並び換えて魔方陣の対角要素とし、C0〜C4を列、L0〜L4を行とする と方陣の数字が全てユニークに決定します。例えば、対角要素として D0を 6-5-3-22-24と並べるとします。 
 6  *  *  *  * 
 *  5  *  *  * 
 *  *  3  *  * 
 *  *  * 22  * 
 *  *  *  * 24
00を含む行と列を考えます。6は C1と L1に含まれているので、候補はC1とL1です。01の列は5を含まなければならないので、C3が候補です。C3とL1の共通の数字は21ですから、01は21になります。このようにして埋めていくと 
      C0
      V
 6 21 14 12  7 
20  5 18  8  9 
11 13  3 17 16 <-L0
 0 19 15 22  4 
23  2 10  1 24
となります。この方陣は全ての行、列および対角の和が60になっています。もしこのとき副対角の和が60ならば方陣は間違いなく魔方陣です。例えば上の例では副対角の和は60ですから、間違いなく魔方陣です。このときに副対角の組みをD1で表現します。上の場合は 
D1 [3 7 8 19 23]
です。

なお、L0、C0、D0はN(例では3)を含んでいるので、真中が Nとなるように並び変えをする場合、中央の行および列の要素はC0と L0に必ずなることに注意してください。

D0の並び方としては5!=120通りあるのですが、真中が Nとなるものは4!= 24通りしかありません． つまり，C0〜C4、L0〜L4、D0を全ての組み(1394種類)とD0の並び方(24種)に対してしらみつぶしに魔法陣になるかチェックを行い個数をカウントすれば全ての魔方陣の総数が求まることになります。

最適化

その1

以上で全ての個数が求まるのですが、これだけでは対称形など無駄な魔方陣も探してしまいます。鏡像や回転はあまりに自明なので、普通はそれらはカウントアップしません。

先程の説明ではD0の並び方は24通りと書いていますが、それを書き出してみると 

1 2 0 3 4←→4 3 0 2 1
1 2 0 4 3←→3 4 0 2 1
1 3 0 2 4←→4 2 0 3 1
1 3 0 4 2←→2 4 0 3 1
1 4 0 2 3←→3 2 0 4 1
1 4 0 3 2←→2 3 0 4 1
2 1 0 3 4←→4 3 0 1 2
2 1 0 4 3←→3 4 0 1 2
3 1 0 2 4←→4 2 0 1 3
3 1 0 4 2←→2 4 0 1 3
4 1 0 2 3←→3 2 0 1 4
4 1 0 3 2←→2 3 0 1 4
のように半分は中央に対して対象です。D0の並びが対象のとき魔方陣が点対象の形になるのは自明ですから、実際は12通りのみカウントアップすればよいことがわかります。 ところで5次の魔方陣には次のような事実があります。

(A) D0を適当に並び換えて魔方陣となっているのなら、00と44を入れ換えても魔方陣となる

例えば 
17 24  2  3 14          21  6 13  4 16  
 8  7 22 18  5           5  7 22 18  8
10 11  0 20 19 →       19 11  0 20 10
 9 12 23 15  1           1 12 23 15  9
16  6 13  4 21          14 24  2  3 17
のように右側も魔方陣となっているのがわかると思います(証明は簡単です)。この他にも

(B) 00<->11および 33<->44を入れ換える
(C) 11<->33を入れ変える
(D) 00<->33および 11<->44を入れ変える

などという自明な入れ換えが存在します。(C)は(A)の点対称形、(D)は(A)と(B)から導かれるので、実際は (A)と(B)のみ考えれば十分です。よってさらにチェックする並び換えが1/4となり、最終的には 

1 2 0 3 4
1 2 0 4 3
1 3 0 4 2
の3通りのみをチェックすれば十分ということになります。

その2

あきらかに C0とL0を交換した魔方陣の総数は交換しても魔方陣の総数はかわりません。 なぜなら，C0 > L0の魔方陣は必ずその対称形が C0 < L0の魔方陣としてカウントアップされているからです．よって， 
C0 < L0
としても問題はありません。

同様の考え方で D0 < D1としてしまっても問題ありません。なぜなら D0 > D1の魔方陣は対称形が D0 < D1の魔方陣としてカウントアップされているからです．よって、カウントする必要もありません。

C1〜C4、L1〜L4はD0の選び方および並び方により魔方陣に収まる位置が一意に定まるため、その順番に意味はありません。よって， 

C1 < C2 < C3 < C4
L1 < L2 < L3 < L4
としてしまっても問題ありません。つまり， 
C0 < L0
D0 < D1
C1 < C2 < C3 < C4
L1 < L2 < L3 < L4
の場合だけを検索すれば十分で，これにより検索時間がかなり節約できます．

プログラム

これです。簡単に説明します。

mk_sqtab()

組みを作成する関数です。実際のプログラムでは演算を高速にするため、 
0  1UL<<24
1  1UL<<23
2  1UL<<22
...
24 1UL<<0
と現わすこととし、各組をこれらの数字のorで現わしたものを「組み」としています。例えば「0 1 12 23 24」なら 
(1UL<<24) | (1UL<<23) | (1UL<<12) | (1UL<<1) | (1UL<<0) =  25169923(0x1801003)
です。言うまでもなくこの数字を2進数で表現すると 1が5つ現われる計算になります。そして、この数字の大きい順番に組みを並べかえると「組み」の順番のルールを満します。

このように選ぶと、C0と C1が互いに独立というのは 

if(C0 & C1 == 0)
と単に論理積を取ることで求められるので、演算を高速にすることができます。 また，順序づけも単純な整数の大小比較で求めることができるので，とても高速になります．

is_cross()

交差のチェックをする関数です。少なくとも25回のシフト演算と比較演算で求めることが可能です。

sq_test()

D0を並びかえ、副対角要素の和が60であるかどうかをチェックする関数です。チェックした結果魔方陣ならば個数をカウントアップします。また100個毎に魔方陣を表示します。100個毎に表示するのは全てを表示しようとすると出力のオーバヘッドが無視できなくなるからです。

コンパイル＆実行

コンパイルをする際に NUMをオプションで #defineしてください。真中が10の魔方陣の総数をチェックしたいのなら、 

% cc -O3 -DNUM=10 msq5.c -o msq5
などとすればよいでしょう。以上のプログラムで NUM=0からNUM=12を計算で求めると以下のようになります。 
0 found 1,091,448 magic squares
1 found 1,366,179 magic squares
2 found 1,914,984 magic squares
3 found 1,958,837 magic squares
4 found 2,431,806 magic squares
5 found 2,600,879 magic squares
6 found 3,016,881 magic squares
7 found 3,112,161 magic squares
8 found 3,472,540 magic squares
9 found 3,344,034 magic squares
10 found 3,933,818 magic squares
11 found 3,784,618 magic squares
12 found 4,769,936 magic squares
13〜24の総数は 0〜11の総数に等しいのでトータルでは 
68,826,306
となります。ただしD0を決定するときに除外した、対称形ではない自明な交換 
(A) 00<->44
(B) 00<->11, 33<->44
があるので、それらを考慮すると、 
68,826,306 * 4 =  275,305,224
となります。この個数は別のアルゴリズムで求めた大森清美氏の個数と一致します。

前へ

次へ

To Menuへ By Issei