0. はじめに
「笑わない数学者」のビリヤードの問題は魔円陣（Magic Circle）とも呼ばれており，古くから知られていた数学上の問題のようだ。今回は過去の文献を調査した結果をレポートしたい。
参考文献[1]によると，魔円陣（Magic Circle）とはフィリピンの Pablo Manalastas 博士が提唱した問題と述べており，近現代に提案された問題と思いきや，文献[3]によると歴史はもっと古く，初出は1857年まで遡るとのこと。まあ，これはその後の研究で情報が刷新されたのだろう。また既に当時から岩堀氏（おそらく岩堀長慶氏）が以下の予想を挙げていた。
- 魔円陣が作れる必要条件は，魔円陣の個数 $n$ が素数 $p$ のべき乗に１を加えた値であること。すなわち $n = 1 + p^m$ であること。
- 魔円陣が存在するような $n$ について，相異なる魔円陣の個数は $\phi(n^2 - n + 1)/6m$ であること。ここで $\phi(x)$ はオイラーのφ関数である。
ただし，当時から40年以上経過した現在でも予想は予想のままのようで，あまり進展していないようだ。
文献[2]は，当時中学教諭の下林山稔氏が魔円陣を作成する BASIC プログラムを公開し，$n = 14$ までの結果を掲載している。当時最新機種だった 8bit パソコン，日立製 S1/40 を駆使し，都合100時間を費やして求めたそうで，いろいろ泣けてくる。
さて，このプログラムを今の PC で実行させたらどのくらいの速さになるのだろうか？衝撃的な結末に期待して頂きたい。まあタイトルでネタバレしているが。
1. プログラム解説
公開された BASIC プログラムはコメントを除くと僅か40行足らずのコンパクトなコードではあるが，GOTO文でアチコチに飛んでいる＊いわゆる＊スパゲティコードである。同時に掲載されたフローチャートを読んでもなかなか構造を理解するのが難しいと思われる。また，残念ながら掲載プログラムをそのまま実行できる環境を今でも保有している方も少ないであろう。とはいえ，ソースコードを眺めるだけではこのプログラムの凄さが伝わらないだろう。逆に動かしてみると凄さが分かる。このため最新の BASIC 環境である Visual Basic 向けに移植作業を行った。
Visual Basic を採択した理由は，昨今の Windows は Visual Basic コンパイラを内蔵しており，新たに開発ツールを購入する必要なく，フリーソフトすらインストールする必要もなく Visual Basic コンパイラを使用することができるからだ。もちろんオリジナルの BASIC プログラムに敬意を示す意味もある。ただし，使用するにはちょっとした手順が必要なので次章で説明する。
移植に際して，まず再帰関数を使えば非常にすっきりとした構造に収められる。あとは Visual Basic の文法に沿った変更を行い，コンソール入出力関数を DOTNET 環境向けに書き換えれば良い。加えてバッチ処理向けに魔円陣の個数をコマンドラインで指定できるようにした。
なお，オリジナルプログラムとは完全に同じ結果が得られることを確認している。
OPTION EXPLICIT ON
'*******************************************************************
' 魔円陣プログラム
'
' ORIGINAL PROGRAM by M.S 1985/12/20
' PORTED TO VB.NET by TETSURO 2024/03/11
'*******************************************************************
MODULE MAENJIN
DIM A() AS INTEGER ' 魔円陣を作る配列
DIM B() AS BOOLEAN ' 表される数をチェックする配列
DIM N AS INTEGER ' 魔円陣の個数
DIM H AS INTEGER ' 表される最大数
'***************************************************************
' メイン関数
'***************************************************************
SUB MAIN(ARGS() AS STRING)
'***********************************************************
' 準備
'***********************************************************
IF ARGS.LENGTH <> 0 THEN
N = CINT(ARGS(0))
ELSE
CONSOLE.WRITE("何個の魔円陣ですか? ")
N = CINT(CONSOLE.READLINE)
END IF
H = N * (N - 1) + 1
REDIM A(N)
REDIM B(H)
'***********************************************************
' １の数は先頭に固定する
' ２の数は二番目から中央まで移動する
'***********************************************************
A(1) = 1
RCHECK(2, 1 + (N >> 1))
CONSOLE.WRITELINE("END")
END SUB
'***************************************************************
' 再帰的チェック関数，非表示の最小数 M，検索範囲の上限 E とする
'***************************************************************
SUB RCHECK(M AS INTEGER, E AS INTEGER)
DIM K, S AS INTEGER
'***********************************************************
' 非表示の最小数 M を二番目から E まで移動する
'***********************************************************
FOR K = 2 TO E
IF A(K) <> 0 THEN CONTINUE FOR
A(K) = M
'*******************************************************
' 加算チェック
'*******************************************************
IF ACHECK THEN
'***************************************************
' 次の非表示の最小数を探す
'***************************************************
FOR S = 2 TO H
IF NOT B(S) THEN EXIT FOR
NEXT
IF S >= H THEN
'***********************************************
' 魔円陣の表示
'***********************************************
OUTPUT
ELSE
'***********************************************
' 再帰的に次の数をチェックする
'***********************************************
RCHECK(S, N)
END IF
END IF
A(K) = 0
NEXT
END SUB
'***************************************************************
' 加算チェック関数
'***************************************************************
FUNCTION ACHECK
DIM I, J, L, T AS INTEGER
FOR I = 1 TO H: B(I) = FALSE: NEXT
FOR I = 1 TO N
IF A(I) = 0 THEN CONTINUE FOR
L = 0
T = I
FOR J = 1 TO N - 1
IF A(T) = 0 THEN EXIT FOR
L += A(T)
IF L > H ORELSE B(L) THEN RETURN FALSE
B(L) = TRUE
T += 1
IF T > N THEN T -= N
NEXT
NEXT
RETURN TRUE
END FUNCTION
'***************************************************************
' 魔円陣の表示関数
'***************************************************************
SUB OUTPUT
DIM I AS INTEGER
CONSOLE.WRITE("[")
FOR I = 1 TO N
CONSOLE.WRITE(" " & A(I))
NEXT
CONSOLE.WRITELINE(" ]")
END SUB
END MODULE
当初，文献[2]に掲載されたオリジナルプログラムをそのまま転載しようと思ったが，著作権者の許可を得ようにも連絡先が分からないこと，引用であれば許可を得る必要はないが，全50行余りの小さなプログラムとはいえプログラム全文を引用することが必要最小限と認められるのか法的判断が難しいため，オリジナルプログラムをそのまま掲載することは控えた。
Visual Basic に移植したプログラムは移植とは名ばかりでゼロから新規で書き起こしており，オリジナルのフローチャートに沿ってはいるが，近代的な構造化プログラミングの形で再構築され，元のプログラムの原型を留めていない。すなわち移植したプログラムは新たな創作物であり，公開しても著作権侵害には当たらないと判断した。
2. ビルド手順
Windows 内蔵の Visual Basic コンパイラを使うためのバッチファイルを以下に示す。
@echo off
for /F "usebackq tokens=1,2,*" %%I in ( `reg query "HKLM\SOFTWARE\Microsoft\NET Framework Setup\NDP\v4\Full" /v InstallPath` ) do (
if "%%I"=="InstallPath" set "DOTNETPATH=%%K"
)
if "%DOTNETPATH:~-1%"=="\" set "DOTNETPATH=%DOTNETPATH:~0,-1%"
set "PATH=%PATH%;%DOTNETPATH%"
echo PATH に %DOTNETPATH% を追加しました。
コマンドプロンプト上で上記のバッチファイルを実行して Visual Basic コンパイラにパスを通した後，以下のコマンドを実行してコンパイルする。MAENJIN.EXE という実行ファイルが作成されていたらビルド成功である。
vbc /optimize MAENJIN.VB
3. 実行例
実行例を以下に示す。引数なしで実行すると魔円陣の個数を尋ねてくる。以下は 6 個を入力した場合である。
maenjin
何個の魔円陣ですか? 6
[ 1 2 5 4 6 13 ]
[ 1 2 7 4 12 5 ]
[ 1 3 2 7 8 10 ]
[ 1 3 6 2 5 14 ]
[ 1 7 3 2 4 14 ]
[ 1 14 4 2 3 7 ]
[ 1 14 5 2 6 3 ]
END
魔円陣の個数を引数として与えることもできる。以下は 8 個の例である。
maenjin 8
[ 1 2 10 19 4 7 9 5 ]
[ 1 8 2 3 21 4 12 6 ]
[ 1 12 2 6 3 7 22 4 ]
[ 1 4 2 10 18 3 11 8 ]
[ 1 3 8 2 16 7 15 5 ]
[ 1 3 5 11 2 12 17 6 ]
[ 1 6 17 12 2 11 5 3 ]
END
4. 実行結果
以下のバッチファイルを用いて 2～17 個の範囲で魔円陣を探索した。
@echo off
for /L %%I in (2,1,17) do (
echo n = %%I のとき
maenjin %%I
echo.
)
実行結果を以下に示す。$n = 4$ の場合を見ると一目瞭然だが，[1 3 2 7] と [ 1 7 2 3] の左右対称解が両方含まれており，左右対称解を排除し切れていない。ただし，全ての解を漏れなく得られていることは確認できた。
n = 2 のとき
[ 1 2 ]
END
n = 3 のとき
[ 1 2 4 ]
END
n = 4 のとき
[ 1 2 6 4 ]
[ 1 3 2 7 ]
[ 1 7 2 3 ]
END
n = 5 のとき
[ 1 5 2 10 3 ]
END
n = 6 のとき
[ 1 2 5 4 6 13 ]
[ 1 2 7 4 12 5 ]
[ 1 3 2 7 8 10 ]
[ 1 3 6 2 5 14 ]
[ 1 7 3 2 4 14 ]
[ 1 14 4 2 3 7 ]
[ 1 14 5 2 6 3 ]
END
n = 7 のとき
END
n = 8 のとき
[ 1 2 10 19 4 7 9 5 ]
[ 1 8 2 3 21 4 12 6 ]
[ 1 12 2 6 3 7 22 4 ]
[ 1 4 2 10 18 3 11 8 ]
[ 1 3 8 2 16 7 15 5 ]
[ 1 3 5 11 2 12 17 6 ]
[ 1 6 17 12 2 11 5 3 ]
END
n = 9 のとき
[ 1 2 4 8 16 5 18 9 10 ]
[ 1 16 22 2 3 4 6 8 11 ]
[ 1 8 12 2 3 13 24 4 6 ]
[ 1 14 8 2 28 3 6 7 4 ]
END
n = 10 のとき
[ 1 2 6 18 22 7 5 16 4 10 ]
[ 1 36 2 12 7 8 3 13 4 5 ]
[ 1 4 2 20 8 9 23 10 3 11 ]
[ 1 18 3 2 8 4 29 11 9 6 ]
[ 1 4 3 10 2 9 14 16 6 26 ]
[ 1 18 28 5 2 8 6 11 9 3 ]
END
n = 11 のとき
END
n = 12 のとき
[ 1 2 14 4 37 7 8 27 5 6 13 9 ]
[ 1 2 9 8 14 4 43 7 6 10 5 24 ]
[ 1 2 12 31 25 4 9 10 7 11 16 5 ]
[ 1 2 14 12 32 19 6 5 4 18 13 7 ]
[ 1 7 2 22 4 16 3 11 29 21 12 5 ]
[ 1 11 2 8 7 25 6 3 27 20 19 4 ]
[ 1 6 2 33 17 14 12 10 15 3 16 4 ]
[ 1 14 3 2 4 7 21 8 25 10 12 26 ]
[ 1 31 10 2 4 3 8 13 5 20 14 22 ]
[ 1 9 48 2 4 18 5 11 3 12 13 7 ]
[ 1 8 5 2 18 11 10 22 6 24 23 3 ]
[ 1 6 28 2 13 16 23 19 5 9 8 3 ]
[ 1 3 12 34 21 2 8 9 5 6 7 25 ]
[ 1 4 7 3 16 2 6 17 20 9 13 35 ]
[ 1 35 26 3 11 2 4 21 7 5 10 8 ]
[ 1 15 5 3 25 2 7 4 6 12 14 39 ]
[ 1 41 8 4 17 2 3 6 7 20 10 14 ]
[ 1 7 6 16 15 2 9 10 40 3 20 4 ]
END
n = 13 のとき
END
n = 14 のとき
[ 1 2 21 17 11 5 9 4 26 6 47 15 12 7 ]
[ 1 2 13 7 5 14 34 6 4 33 18 17 21 8 ]
[ 1 2 28 14 5 6 9 12 48 18 4 13 16 7 ]
[ 1 10 2 16 32 20 3 21 33 5 4 22 8 6 ]
[ 1 5 2 24 15 29 14 21 13 4 33 3 9 10 ]
[ 1 15 34 2 3 17 7 6 8 4 19 9 48 10 ]
[ 1 4 20 2 12 3 6 7 33 11 8 10 35 31 ]
[ 1 12 48 6 2 38 3 22 7 10 11 5 4 14 ]
[ 1 33 14 36 2 16 6 5 8 4 3 25 21 9 ]
[ 1 12 11 6 2 18 16 35 21 4 3 40 5 9 ]
[ 1 9 8 28 2 24 5 6 44 14 20 7 12 3 ]
[ 1 25 23 4 10 2 3 17 11 7 6 45 21 8 ]
[ 1 9 14 26 4 2 11 5 3 12 27 34 7 28 ]
[ 1 19 16 5 8 2 7 55 11 14 12 6 24 3 ]
[ 1 12 22 7 17 2 18 10 23 32 25 6 5 3 ]
[ 1 4 6 31 3 13 2 7 14 12 17 46 8 19 ]
[ 1 17 8 5 24 14 2 4 3 12 27 34 22 10 ]
[ 1 8 16 10 12 19 2 11 4 3 42 27 23 5 ]
[ 1 18 7 17 15 14 2 4 56 5 23 10 3 8 ]
[ 1 4 8 52 3 25 18 2 9 24 6 10 7 14 ]
[ 1 14 7 10 6 24 9 2 18 25 3 52 8 4 ]
END
n = 15 のとき
END
n = 16 のとき
END
n = 17 のとき
[ 1 2 4 8 16 32 27 26 11 9 45 13 10 29 5 17 18 ]
[ 1 7 31 2 11 3 9 36 17 4 22 6 18 72 5 10 19 ]
[ 1 23 11 3 2 13 33 22 6 4 17 9 41 25 44 12 7 ]
[ 1 27 8 7 2 3 26 25 16 14 4 6 13 39 50 11 21 ]
[ 1 3 12 10 31 7 27 2 6 5 19 20 62 14 9 28 17 ]
[ 1 7 3 15 33 5 24 68 2 14 6 17 4 9 19 12 34 ]
END
以上，バッチファイルが終了するまで筆者の古いノート PC で約７時間要した。
5. 計算時間の比較
あれれ？このプログラム，めっちゃ速いぞ！というか
我 が プ ロ グ ラ ム は 圧 倒 的 に 負 け て い る
拙作の C プログラム billiard との比較結果を以下に示す。maenjin は VB.NET のプログラムのため最初に JIT コンパイラによるコンパイルを一定時間要するが，ボールの個数（魔円陣の個数）が増えてコンパイル時間を無視できるようになると，ざっと40～50倍速い結果となった。計算量は $O(8.5^n)$ 程度に見えるが，下林山氏の予測 $O(n!/2^n)$ も正しいように見える。
以下計算時間の比較を示す。単位は秒である。
| n | billiard | maenjin |
|---|---|---|
| 4 | 0.01 | 0.03 |
| 5 | 0.01 | 0.03 |
| 6 | 0.01 | 0.03 |
| 7 | 0.01 | 0.03 |
| 8 | 0.01 | 0.03 |
| 9 | 0.03 | 0.03 |
| 10 | 0.15 | 0.04 |
| 11 | 1.17 | 0.11 |
| 12 | 11.39 | 0.65 |
| 13 | 112.67 | 4.22 |
| 14 | 1,124.96 | 36.15 |
| 15 | 11,178.60 | 280.58 |
| 16 | 112,038.24 | 2,562.51 |
| 17 | 1,122,292.52 | 21,700.87 |
グラフにすると綺麗に対数スケールで計算時間が増加していることが分かる。
6. アルゴリズム考察
6.1 探索範囲の比較
なぜ，こんなに速いのか？というと探索範囲が桁違いに狭いからである。以下，ボールの個数（魔円陣の個数）に対する再帰関数の呼び出し回数を比較するが，maenjin のほうが圧倒的に少ない。もちろん再帰関数内での処理負荷はまったく異なる（maenjin のほうが重い）ので公正な比較ではないが，一つの目安になろう。
| $n$ | billiard | maenjin |
|---|---|---|
| 4 | 3 | 6 |
| 5 | 24 | 17 |
| 6 | 165 | 72 |
| 7 | 1,139 | 242 |
| 8 | 8,428 | 1,114 |
| 9 | 67,794 | 4,642 |
| 10 | 553,483 | 25,846 |
| 11 | 4,774,154 | 130,133 |
| 12 | 41,279,613 | 828,252 |
| 13 | 366,012,625 | 4,903,408 |
| 14 | 3,304,581,480 | 35,414,975 |
| 15 | 29,801,016,280 | 238,493,464 |
| 16 | 272,351,149,456 | 1,906,617,239 |
| 17 | 2,511,320,203,782 | 14,350,975,350 |
グラフにすると一目瞭然である。
6.2 maenjin の探索アルゴリズム
maenjin のアルゴリズムは人間が手計算で行う手続きに似ており分かり易い。必ず１と２の数値は必要なのでまず１の数値を先頭に固定し，２の数値は左右対称解を避けて二番目から「中央」まで移動させる。「中央」とは魔円陣の個数 $n$ が奇数のときは $(n + 1)/2$ とし，偶数のときは $n/2 + 1$ とする。これらをまとめると $\lfloor n / 2 \rfloor + 1$ となる。
次に空いたスペースに「表されていない最小数」を順に置いていく。「表されていない最小数」とは，たとえば２の数値が１の隣にある場合は４となり，隣にない場合は３となる。
２の数値の置く位置は二番目から中央過ぎ $\lfloor n / 2 \rfloor + 1$ までだから「場合の数」は $\lfloor n / 2 \rfloor$ 通り，以降，「表されていない最小数」を空いているスペースに置く「場合の数」は $n - 2$，$n - 3$，$n - 4$ 通りと一つ置くたびに減っていき，これらの積が $n$ 個の魔円陣における組み合わせの数 $c_n$ と考えると
c_n = \left\lfloor \frac{1}{2} n \right\rfloor \times (n - 2)!
となるが，実際の再帰関数の呼び出し回数はこれよりもずっと少なくなっている。「表されていない最小数」を順に置いていく途中で既にある数と重複してしまうケースが次第に増えてくるからだ。
6.3 billiard の探索アルゴリズム
一方，billiard では１の数値を先頭に固定する，すなわち $a(1) = 1$ とする点は同じだが，以降の数値は二番目，三番目と先頭から順に埋めていく点が異なる。ボールの個数を $n$ とおくと二番目に置く数 $a(2)$ の範囲は左右対称解を避けて $2 \le a(2) \le \lfloor (n + 2)(n - 1)/4 \rfloor$ であり，「場合の数」は $\lfloor (n + 2)(n - 1)/4 \rfloor - 1$ 通りであるが，まずこの時点で maenjin の $n / 2$ 通りよりも圧倒的に多いことに気づく。
次に三番目に置く数 $a(3)$ の範囲は二番目の数値 $a(2)$ によって変わり，
- $a(2) \le n - 2$ のとき，$2 \le a(3) \le n(n - 1)/2 + 1$
- $a(2) \ge n - 1$ のとき，$2 \le a(3) \le n(n + 3)/2 - 3 - 2a(2)$
となる。$a(3)$ の取り得る「場合の数」は $a(2)$ の値を除いて
- $a(2) \le n - 2$ のとき，$n(n - 1)/2 - 1$ 通り
- $a(2) \ge n - 1$ のとき，$n(n + 3)/2 - 5 - 2a(2)$ 通り
となる。これらの平均値が実質的な $a(3)$ の取り得る「場合の数」となるが，概算でも $n^2 / 4$ 以上あり，maenjin の $n - 2$ 通りと比べるとやはり圧倒的に多い。
これ以降に置く数の「場合の数」についても同様に圧倒的な差があると考えられる。これらの差の積み重ねより，billiard の探索範囲は maenjin に比べて桁違いに多くなっているものと考える。
なお billiard の探索アルゴリズムについては過去の記事で詳しく述べている。
7. まとめ
おそらく下林山氏の作成した魔円陣プログラムは人間が行う手計算を忠実にトレースしたものと思われる。加算チェックのたびに配列をクリアして計算し直す点が非効率的だと思うが，その点を差し置いても全体としては計算量の少ない非常に優れたアルゴリズムであることに気づかされる。
当初，温故知新というタイトルで昔の文献を紹介しようと思ったのだが，40年前のBASICプログラムをリスペクトする記事になるとは思わなかった。だがしかし，40年前のBASICプログラムに負けてばかりでは遺憾なので次回リベンジしたい。
8. 参考文献
参考文献[1]は今読んでも面白い。$3n+1$ 問題，いわゆるコラッツ予想を取り上げている。
[1] 島内剛一，λ倍取りゲーム，数学セミナー1982年1月号，10～14頁，日本評論社
[2] 下林山稔，パソコンで魔円陣を，数学セミナー1987年7月号，55～59頁，日本評論社
[3] 秋山茂樹，魔円陣と有限幾何，数学セミナー2013年7月号，25～31頁，日本評論社 
[4] 小林吹代，ガロアの数学「体」入門～魔円陣とオイラー方陣を例に～，2018年5月26日，技術評論社
[5] J. Singer, A Theorem in Finite Projective Geometry and Some Applications to Number Theory,Trans. American Mathematical Society, 43-3 (1938), pp.377-385.