今さら森博嗣「笑わない数学者」のビリヤードの問題を解く(5)～未解決問題に挑む

0. はじめに

前回の記事で一般化したビリヤードの問題のボールの個数最大（世界最高記録）を目指したが，よくよく考えてみると計算が容易な点（ボールの個数が１＋素数であるとき）のみを計算しており，計算が困難な点（ボールの個数が１＋素数のべき乗であるとき）をスキップしていた。

まるで五合目や八合目を飛び越えて頂上にワープしたようなもので，これでは世界最高峰に登山したとは言えないのではないか？

残念ながら，参考文献のサイトにはボールの個数が１＋素数のべき乗の場合については詳しく書かれていないようだ。本記事は参考文献のサイトのヒントを手掛かりに，手探りでコーディングしていたら，なんとなく１＋素数のべき乗の場合でも上手く出来てしまったというものである。

1. アルゴリズム解説

素数 $p$，自然数 $m$ に対して，ボールの個数 $n = 1 + p^m$ とする。素数のべき乗数（指数）$m$ に応じて $3m$ 次ボナッチ数列 $s_i$ を作る。初項 $s_i$（$0 \le i < 3m$）は

s_i = \left\lbrace \begin{array}{ll}
0 & 0 \le i < 3m - 1 \\
1 & i = 3m - 1
\end{array} \right.

とし，一般項 $s_i$（$i \ge 3m$）は

s_i = \left[ \displaystyle \sum_{k = 1}^{3m} a_{3m - k} \times s_{i - k} \right] \bmod p

とする。ここで数列を決定する重要なパラメータ $0 \le a_k < p$ である。

ここでボールの数値の総和 $\sigma = n(n - 1) + 1$ として，$s_i$ から $s_{i + (m - 1)\sigma}$ まで間隔 $\sigma$ で数列の値がすべてゼロになる $i$ の位置を読み取り，これを数列 $t_j$ とする。

表１　$3m$ 次ボナッチ数列

$i$	$0$	$\cdots$	$t_0$	$\cdots$	$t_1$	$\cdots$	$t_2$	$\cdots$	$2\sigma - 1$
$s_i$	$\cdots$	$\cdots$	$0$	$\cdots$	$0$	$\cdots$	$0$	$\cdots$	$\cdots$
$s_{i + \sigma}$	$\cdots$	$\cdots$	$0$	$\cdots$	$0$	$\cdots$	$0$	$\cdots$	$\cdots$
$s_{i + 2\sigma}$	$\cdots$	$\cdots$	$0$	$\cdots$	$0$	$\cdots$	$0$	$\cdots$	$\cdots$
$\cdots$	$\cdots$	$\cdots$	$0$	$\cdots$	$0$	$\cdots$	$0$	$\cdots$	$\cdots$
$s_{i + (m - 1)\sigma}$	$\cdots$	$\cdots$	$0$	$\cdots$	$0$	$\cdots$	$0$	$\cdots$	$\cdots$

ちなみに数列 $s_i$ をいくつまで計算する必要があるかというと，$t_j = 1$ となる位置を探す（頭出しする）必要から $t_j$ は二周期 $2\sigma$ 分得るとして$s_i$ はプラスして $m - 1$ 周期分 $(m - 1)\sigma$ だけ求める必要がある。すなわち合計 $m + 1$ 周期分 $(m + 1)\sigma$ 必要である。

次に数列 $u_j = t_{j + 1} - t_j$ を求める。$t_j$ は純単調増加の数列のため $u_j$ は必ず自然数の数列となる。

パラメータ $a_k$ を正しく設定すると，$u_j$ の任意の位置から連続する $n$ 個の項を抜き出せばビリヤード問題の解になるが，回転対称の解を除くため $u_j = 1$ となる位置から連続する $n + 1$ 個の項を抜き出し，これを数列 $v_h$（$0 \le h \le n$）とおく。正しくビリヤード問題の解を得られていれば必ず $v_n = 1$ となる。

表２　抜き出した $n + 1$ 個の項

$h$	$0$	$1$	$\cdots$	$n - 1$	$n$
$v_h$	$1$	$v_1$	$\cdots$	$v_{n - 1}$	$1$

左右対称の解を除くため，$v_1 < v_{n - 1}$ の場合は $v_0$ から連続する $n$ 個の項を抜き出し，$v_1 > v_{n - 1}$ の場合は $v_n$ から逆向きに連続する $n$ 個の項を抜き出す。

ちなみにこれで何故ビリヤード問題の解が得られるのかはさっぱり理屈が分からない。素数 $p$ で割った余りの整数列であることから取りうる値は制限されるため，何らかの周期的な数列が作られるところまでは分かる。割る数が素数なのも意味があるのだろう。しかし，これから完全差集合が作られる仕組みのほうはさっぱり分からない。

2. 実装コード

前回と同様 JavaScript（WSH）で作成した。ボールの個数 $n = 1 + p^m$ とおくと，今回，べき乗数（指数）$m$ 別にプログラムを作り分けた。以下，一例としてボールの個数 $n = 1 + p^2$ の場合のプログラムを示す。

OnePlusPrime2.js

var	args = WScript.Arguments.Unnamed;
var	ret = main(args);
try {
	WScript.Quit(ret);
} catch(e) {
	/* 何もしない */
}
//------------------------------------------------------------------------------
// メイン関数
//------------------------------------------------------------------------------
function main(args) {
	//--------------------------------------------------------------------------
	// ヘルプメッセージ
	//--------------------------------------------------------------------------
	if(args.Count == 0) {
		WScript.Echo("Usage: OnePlusPrime2(.js) [prime-number]");
		return false;
	}
	//--------------------------------------------------------------------------
	// コマンドライン解析
	//--------------------------------------------------------------------------
	var	p = parseInt(args(0));
	//--------------------------------------------------------------------------
	// ヘキサナッチ数列のパラメータ a[0]～a[5] を全探索する
	//--------------------------------------------------------------------------
	var	num = p * p + 1;
	var	sum = num * (num - 1) + 1;
	var	a = [0, 0, 0, 0, 0, 0];
	var	d = [];
	function search() {
		//----------------------------------------------------------------------
		// ヘキサナッチ数列より完全差集合を作成する
		//----------------------------------------------------------------------
		var	s = [];
		s.push(0); s.push(0); s.push(0); s.push(0); s.push(0); s.push(1);
		for(var i = 6; i < sum * 3; i++) {
			var	v = a[0] * s[i - 6]
				  + a[1] * s[i - 5]
				  + a[2] * s[i - 4]
				  + a[3] * s[i - 3]
				  + a[4] * s[i - 2]
				  + a[5] * s[i - 1];
			s[i] = v % p;
		}
		var	t = [];
		for(var i = sum; i < sum * 3; i++)
			if(s[i] == 0 && s[i - sum] == 0)
				t.push(i - sum);
		if(t.length != 2 * num) return;
		//----------------------------------------------------------------------
		// 完全差集合をボール配置に変換する
		//----------------------------------------------------------------------
		var	u = [];
		for(var i = 1; i < t.length; i++)
			u.push(t[i] - t[i - 1]);
		//----------------------------------------------------------------------
		// ボール配置を正規化（回転対称，左右対称を除く）
		//----------------------------------------------------------------------
		for(var i = 0; i < u.length; i++)
			if(u[i] == 1) break;
		if(i >= u.length) return;
		var	v = [];
		for(j = 0; j < num + 1; j++)
			v.push(u[i + j]);
		if(v[1] > v[num - 1]) v.reverse();
		v.pop();
		//----------------------------------------------------------------------
		// ボール配置のチェック＆重複するボール配置を取り除く
		//----------------------------------------------------------------------
		if(check(v)) {
			var	w = v.join(" ");
			if(d[w] == null) d[w] = a.join(" ");
		}
	}
	//--------------------------------------------------------------------------
	// ヘキサナッチ数列のパラメータ a[0]～a[5] を全探索する
	//--------------------------------------------------------------------------
	for(a[0] = 0; a[0] < p; a[0]++)
		for(a[1] = 0; a[1] < p; a[1]++)
			for(a[2] = 0; a[2] < p; a[2]++)
				for(a[3] = 0; a[3] < p; a[3]++)
					for(a[4] = 0; a[4] < p; a[4]++)
						for(a[5] = 0; a[5] < p; a[5]++)
							search();
	//--------------------------------------------------------------------------
	// ボール配置をソートして出力する
	//--------------------------------------------------------------------------
	var	e = [];
	for(var val in d)
		e.push({ key:d[val], val:val, arr:val.split(" ") });
	e.sort(function(a, b) {
		for(var i = 0; i < a.arr.length; i++)
			if(a.arr[i] != b.arr[i]) return a.arr[i] - b.arr[i];
		return 0;
	});
	for(var i = 0; i < e.length; i++)
		WScript.Echo([p, e[i].key, e[i].val].join("\t"));
	//--------------------------------------------------------------------------
	// 解の個数を出力する
	//--------------------------------------------------------------------------
	WScript.Echo(e.length);
	return true;
}
//------------------------------------------------------------------------------
// ボール配置のチェック
//------------------------------------------------------------------------------
function check(ball) {
	var	num = ball.length;
	var	sum = num * (num -  1) + 1;
	for(var i = 0, s = 0; i < ball.length; i++)
		s += ball[i];
	if(sum != s) return false;
	var	flag = [];
	for(var i = 0; i <= sum; i++)
		flag.push(false);
	for(var i = 0; i < ball.length - 1; i++) {
		for(var j = 0, s = 0; i + j >= 0; j--) {
			s += ball[i + j];
			if(flag[s] || flag[sum - s]) return false;
			flag[s] = true;
			flag[sum - s] = true;
		}
	}
	return true;
}

以下，べき乗数（指数）$m$ が増えるにつれて，全探索するループのネストが深くなるという非常に見苦しいプログラムになっている。誰か上手いコーディングの仕方があれば教えて欲しい。Ruby あたりだと順列・組み合わせのパターンを自動生成してくれる機能があるらしいが，それ以外の言語だとどう書けば良いのだろうか？

パフォーマンスは落ちるかもしれないが，再帰関数で書けばシンプルに書けるかもしれない。

3. 実行例

ボールの個数 $n = 1 + p^4$ の場合のプログラム OnePlusPrime4.js の実行結果を示す。一例としてボールの個数 $n = 17 = 1 + 2^4$ の場合の実行例を示す。引数には素数 $p = 2$ の値を与える。

最初の $2$ は与えた素数の数，続く１２個の数列はパラメータ $a_k$ であり，その次の数列がボールの配置となる。最後の行の $6$ が解の個数である。

ボールの個数が17個の場合の実行例

OnePlusPrime4 2
2	1 0 0 1 0 0 1 1 1 1 1 1	1 2 4 8 16 32 27 26 11 9 45 13 10 29 5 17 18
2	1 0 0 0 0 0 1 0 1 0 0 1	1 3 12 10 31 7 27 2 6 5 19 20 62 14 9 28 17
2	1 0 0 0 1 0 1 0 0 1 1 1	1 7 3 15 33 5 24 68 2 14 6 17 4 9 19 12 34
2	1 0 0 0 1 0 0 1 1 1 0 1	1 7 12 44 25 41 9 17 4 6 22 33 13 2 3 11 23
2	1 0 0 0 1 0 1 0 1 0 1 1	1 7 31 2 11 3 9 36 17 4 22 6 18 72 5 10 19
2	1 0 0 0 0 1 1 1 0 1 0 1	1 21 11 50 39 13 6 4 14 16 25 26 3 2 7 8 27
6

かつて二週間近くかかった計算だが，このプログラムだと一瞬で答えが出てくる。

4. 実行結果

以下，ボールが129個以下において，解が存在しうると考えられるケース（ボールの個数が１＋素数または１＋素数のべき乗）において解の個数を調べた結果を以下に示す。素数 $p$ に対してボールの個数 $n = 1 + p^m$ とおくと，$m$ が大きくなるにつれて解の個数は減少する傾向があるように見える。

表３　ボールの個数が129個以下における解の個数

$n$	$n - 1$							解の個数	計算時間
	$p$	$p^2$	$p^3$	$p^4$	$p^5$	$p^6$	$p^7$
3	○							1	0.039
4	○							2	0.040
5		○						1	0.040
6	○							5	0.043
8	○							6	0.048
9			○					4	0.051
10		○						6	0.050
12	○							18	0.048
14	○							20	0.053
17				○				6	0.120
18	○							51	0.078
20	○							42	0.079
24	○							78	0.144
26		○						30	0.612
28			○					42	1.096
30	○							132	0.357
32	○							110	0.409
33					○			30	2.790
38	○							132	0.722
42	○							287	1.320
44	43							210	1.477
48	○							360	2.810
50		○						126	15.010
54	○							408	4.300
60	○							590	8.320
62	○							384	7.330
65						○		72	126.530
68	○							420	10.460
72	○							852	19.280
74	○							600	19.610
80	○							588	23.510
82				○				216	366.050
84	○							1,098	39.090
90	○							1,335	52.730
98	○							1,056	70.660
102	○							1,717	100.740
104	○							1,190	92.660
108	○							1,512	118.590
110	○							1,140	103.490
114	○							1,980	159.790
122		○						648	1343.320
126			○					828	2613.600
128	○							1,806	303.330
129							○	336	5628.350

ボールの個数と解の個数の関係を示す。ボールの個数 $n = 1 + p$ のとき，解の個数はボールの個数の二乗に比例するようである。一方，$n \ne 1 + p$ のときは明らかに解の個数が少ないが，それでもボールの個数の二乗に比例して解の個数が増加している傾向が見える。

◆$n = 1 + p$ のとき，●$n \ne 1 + p$ のとき

計算時間について，ボールの個数 $n = 1 + p$ のときはボールの個数 $n$ の五乗に比例 $O(n^5)$ の傾向により近づいている。一方，$n \ne 1 + p$ のときは計算時間が一桁多くなっている。

5. まとめ

前回，ボールの個数が１＋素数の場合におけるビリヤード問題の解を導出するプログラムを作成した。$n = 128$ のときでさえ計算時間はまだ５分程度であり，計算時間はボールの個数の五乗に比例するので，やる気と時間があればまだまだイケると思う。

今回，ボールの個数が１＋素数のべき乗の場合におけるビリヤード問題の解を導出するプログラムを作成した。ただし，べき乗数（指数）は７まで。ボールの個数が１＋素数の場合に比べて解の個数は一桁少ないが，逆に計算時間は一桁多くなっており，コチラのケースのほうがネックになっている。とくにべき乗数（指数）ごとにプログラムを作り分けているのはメンテナンスも大変で自分でも正直カッコ悪いと自覚しているが，腕に覚えのある人はプログラムの統合にチャレンジして欲しい。

最終回みたいな雰囲気だが，また次回に続く。

6. 参考文献

参考になったサイトを以下に示す。