漸く森博嗣「笑わない数学者」のビリヤードの問題／魔円陣を解いたものの要領を良くする

1. そもそもどういう計算だったか

　 $p$ を素数としての有限体 $GF(p^n)$ を3次の既約多項式から3次拡大して ${GF(p^n)^3}$ を拵えてそこから有限射影平面を構成して, その巡回平面を求めることで魔円陣／(巡回)完全差集合を求める事ができるという話なのだけど, 「既約多項式とそこからの生成元から射影平面を導いて巡回平面を求める」過程が少しサイズが大きくなると大変なことになる。
　これを何とか要領良くできないかという話。
　実際に Julia で Nemo という Pkg を用いて, 有限体を処理して実行するプログラムを書いてみて出来上がったのだけど, $p=2,n=10$ で試してみたところ, いつまでも終わらないので「この程度で終わらないと駄目やん」となったというわけ。

2. 具体的な例での確認

例えば, $p=3$ の場合で確認してみよう。
　 $n=1$ であれば, 未だ少しは簡単だけど, それでも $GF(3)$ を3次拡大する際には3次の既約多項式を用いたりする。例えば, $x^3 + 2x + 1$ で良いらしい。
　先ずこの $x^3 + 2x + 1$ が既約多項式であることはどうすれば確認できるかというと, $GF(3)=\big\{0,1,2\big\}$ だから, $x=0$ とすれば, $x^3 + 2x + 1=1\ne 0$ , $x=1$ とすれば, $x^3 + 2x + 1=1+2+1\equiv 1\ne 0$ , $x=2$ とすれば, $x^3 + 2x + 1=8+4+1\equiv 1\ne 0$ となり, $GF(3)$ の元はいずれも解でないことと, 3次の多項式が因数分解できるとすれば, 必ず1次式を含むことから分かる。
　そうなると当然 $x^3 + 2x + 1=0$ の解は $GF(3)$ には無いので, 新たな数として $\alpha$ を用意することとなる。 $\alpha^3+2\alpha+1=0$ だから $\alpha^3=-2\alpha-1\equiv \alpha+2$ となることに注意すれば,
　$\alpha^1=\alpha$
　$\alpha^2=\alpha^2$
　$\alpha^3=\alpha+2$
　$\alpha^4=\alpha^2+2\alpha$
　$\alpha^5=\alpha^3+2\alpha^2=2\alpha^2+\alpha+2$
　$\alpha^6=2\alpha^3+\alpha^2+2\alpha=2\alpha+4+\alpha^2+2\alpha\equiv\alpha^2+\alpha+1$
　$\alpha^7=\alpha^3+\alpha^2+\alpha=\alpha+2+\alpha^2+\alpha=\alpha^2+2\alpha+2$
　$\alpha^8=\alpha^3+2\alpha^2+2\alpha=\alpha+2+2\alpha^2+2\alpha\equiv2\alpha^2+2$
　$\alpha^9=2\alpha^3+2\alpha=2\alpha+4+2\alpha\equiv \alpha+1$
　$\alpha^{10}=\alpha^2+\alpha$
　$\alpha^{11}=\alpha^3+\alpha^2=\alpha^2+\alpha+2$
　$\alpha^{12}=\alpha^3+\alpha^2+2\alpha=\alpha+2+\alpha^2+2\alpha\equiv\alpha^2+2$
　$\alpha^{13}=\alpha^3+2\alpha=\alpha+2+2\alpha\equiv 2$
　$\alpha^{14}=2\alpha$
　$\phantom{\alpha{15}}\vdots$
　と $\alpha$ の累乗で次々と異なる元が作成されることが分かる。 $\alpha^{14}=2$ となるので, このあとは $\alpha^1$ から $\alpha^{13}$ の2倍の元が現れ, $\alpha^{26}=2\alpha^{13}=4\equiv 1$ となり, $\alpha^{27}=\alpha$ となり巡回する。
　よく知られる通り, 有限体は積に関して巡回群になっているというわけだ。
　書かれたものを見ている限り, 何だ簡単な計算だなぁと思うかも知れないが, 計算間違いは容易に起こるし, 面倒な計算では無いだろうか。それでもコンピュータに任せれば「いとも簡単に」計算結果を見せてはくれるが, 実際には結構な計算リソースを使ってそうだ。
　よくよく見て考えれば, $\alpha^2$ と $\alpha$ と定数項の3つしかないのだし, 「何とか楽できないものか」と考えるのは至極当然だろう。
　 $\alpha^i=a_i\alpha^2+b_i\alpha+c_i$ とし, 面倒なので $(a_i,b_i,c_i)$ と表してみれば,
　$\alpha^0\ \leftrightarrow\ (0,0,1)$
　$\alpha^1\ \leftrightarrow\ (0,1,0)$
　$\alpha^2\ \leftrightarrow\ (1,0,0)$
　$\alpha^3\ \leftrightarrow\ (0,1,2)$
　$\alpha^4\ \leftrightarrow\ (1,2,0)$
　$\alpha^5\ \leftrightarrow\ (2,1,2)$
　$\phantom{\alpha{15}}\vdots$
のようになることが分かる。 $\alpha$ を掛けると, ベクトル表記は成分が1つ左にずれて, 溢れた分は, $\alpha^3$ だから $\alpha+2$ となり, 成分に加わる。
　つまり, $(a_{i+1},b_{i+1},c_{i+1})=(b_i,c_i+a_i,2a_i)$ となるわけだ。これだと, 「機械的に計算する」といっても随分と要領良く楽に計算できることになる。

3. 一般化

　$GF(p)$ の3次の既約多項式を $x^3+Ax^2+Bx+C$ とし, $x^3+Ax^2+Bx+C=0$ の解を $\alpha$ とすれば, $\alpha^3=-A\alpha^2-B\alpha-C$ となる。
　同じように $\alpha^i$ を $a_i\alpha^2+b_i\alpha+c_i$ とし, $(a_i,b_i,c_i)$で扱うとすると,

\begin{align}
　\alpha^i=a_i\alpha^2+b_i\alpha+c_i
\end{align}

だから,

\begin{align}
\alpha^{i+1}&=a_i\alpha^3+b_i\alpha^2+c_i\alpha\\
&=a_i(-A\alpha^2-B\alpha-C)+b_i\alpha^2+c_i\alpha\\
&=(-Aa_i+b_i)\alpha^2+(-Ba_i+c_i)\alpha-a_iC\\
&=a_{i+1}\alpha^2+b_{i+1}\alpha+c_{i+1}
\end{align}

つまり,

\begin{align}
a_{i+1}&=-Aa_i+b_i\\
b_{i+1}&=-Ba_i+c_i\\
c_{i+1}&=-a_iC
\end{align}

という漸化式で計算できるということになります。※もちろん計算はすべてGF(p)上で行われるので, 所謂$p$の剰余系な計算となります。

4. 問題はないのか

　当然問題があって, 「$GF(p)$に対する$A,B,C$はどうやって設定するか」ということです。でもまぁ計算の要領は良くなった筈だし, どのみち「既約多項式か？」と「原始元／生成元か？」は地道に確認するしかなさそうなので, 「地道に確認する」方向で話を進めましょう。
　$A,B,C$の候補は, $GF(p)$を元にしている場合は, $0,1,2,...,p-1$です。ただし, $C\ne 0$です。これは$C=0$だと明らかに

\begin{align}
x^3+Ax^2+Bx+C=x^3+Ax^2+Bx=x(x^2+Ax+B)
\end{align}

と因数分解できて, 既約多項式にならないからです( $x=0$ のとき, $x^3+Ax^2+Bx+C=C=0$ でも分かります)。
　あと $A=B=0$ の場合も除外して良いでしょう。この場合は $x^3=-C$ となり,

\begin{align}
\alpha^1&=\alpha\\
\alpha^2&=\alpha^2\\
\alpha^3&=-C\\
\alpha^4&=-C\alpha\\
&\vdots
\end{align}

となり, $\alpha$ が原始元となりそうにないからです。
　それ以外の$A,B,C$について, 「 $x^3+Ax^2+Bx+C$ が $1,2,\cdots,p-1$ に対して $0$ にならないこと」(既約多項式かどうかのチェック)と $(a_i,b_i,c_i)$ の列を確認して, 「 $i=1,2,\cdots,p^2+p$ で$a_i=b_i=0$ となる事がないこと」が確認できれば, 採用して計算を進める事ができそうです。

5. 具体的に試してみる

　$p=3$の場合で具体的に試してみましょう。
　$A=0,\ B=1,\ C=2$としてみます。 $x^3+x+2$ という事になります。この場合,

\begin{align}
x=0\ とすれば,\ &x^3+x+2=2\ne 0\\
x=1\ とすれば,\ &x^3+x+2=1+1+2=4\equiv 1\ne 0\\
x=2\ とすれば,\ &x^3+x+2=8+2+2=12\equiv 0=0\\
\end{align}

ですから, 既約多項式ではないので駄目です。

\begin{align}
x^3+x+2=x^3+1+x+1=(x+1)(x^2-x+1+1)\equiv(x+1)(x^2+2x+2)
\end{align}

と因数分解できたのですね。（式を見ていても気が付かないのが普通かも...）
　$A=0,\ B=1,\ C=1$だとどうでしょう。$x^3+x+1$ という事になります。この場合は...

\begin{align}
x=0\ とすれば,\ &x^3+x+1=1\ne 0\\
x=1\ とすれば,\ &x^3+x+1=1+1+1=3\equiv 0\\
\end{align}

駄目ですね。
　$A=0,\ B=2,\ C=1$だとどうでしょう。$x^3+2x+1$ という事になります。この場合は...

\begin{align}
x=0\ とすれば,\ &x^3+2x+1=1\ne 0\\
x=1\ とすれば,\ &x^3+2x+1=1+2+1=4\equiv 1\ne 0\\
x=2\ とすれば,\ &x^3+2x+1=8+4+1=13\equiv 1\ne 0\\
\end{align}

既約多項式チェックはクリアです。
　では$(a_i,b_i,c_i)$の列をチェックしてみましょう。

\begin{align}
a_{i+1}&=-Aa_i+b_i=b_i\\
b_{i+1}&=-Ba_i+c_i=-2a_i+c_i\equiv a_i+c_i\\
c_{i+1}&=-a_iC=-a_i\equiv 2a_i
\end{align}

となりますね。そうすると, $1=(a_0,b_0,c_0)=(0,0,1)$から始めてみれば,

\begin{align}
(a_0,b_0,c_0)&=(0,0,1)\\
(a_1,b_1,c_1)&=(0,1,0)\\
(a_2,b_2,c_2)&=(1,0,0)\\
(a_3,b_3,c_3)&=(0,1,2)\\
(a_4,b_4,c_4)&=(1,2,0)\\
(a_5,b_5,c_5)&=(2,1,2)\\
(a_6,b_6,c_6)&=(1,4,4)\equiv(1,1,1)\\
(a_7,b_7,c_7)&=(1,2,2)\\
(a_8,b_8,c_8)&=(2,3,2)\equiv(2,0,2)\\
(a_{9},b_{9},c_{9})&=(0,4,4)\equiv(0,1,1)\\
(a_{10},b_{10},c_{10})&=(1,1,0)\\
(a_{11},b_{11},c_{11})&=(1,1,2)\\
(a_{12},b_{12},c_{12})&=(1,3,2)\equiv(1,0,2)\\
(a_{13},b_{13},c_{13})&=(0,3,2)\equiv(0,0,2)
\end{align}

となります。$p^2+p+1=3^2+3+1=13$で丁度$a_{13}=b_{13}=0$となり, それに至る間では$a_i=b_i=0$となることはないので, 3次拡大が無事にできた上に有限射影平面も構築できている事が分かります。
　これらの $(a_i,b_i,c_i)$ のうち, $a_i=0$ となるのは, $i=0,1,3,9,13$ のときで, これが(巡回)完全差集合となっていて, その隣接2項間の差を取った $1,2,6,4$ が魔円陣となっています。

6. コードにする

　これを Julia のコードにすれば, まだ $GF(p)$ の場合だけだけど, 少しは要領のよい計算プログラムができそうだ。
　具体的に試してみた通り, 素数$p$に対して, $A,B,C$を逐次試して,

$C\ne 0$や$A=B=0$でない
$x=1,...,p-1$を代入し $x^3+Ax^2+Bx+C$が$0$にならない事(既約多項式であること)を確認
$(a_0,b_0,c_0)=(0,0,1)$ から始めて, $(a_i,b_i,c_i)$ を漸化式から順に計算し, $i=p^2+p+1$のとき始めて$a_i=b_i=0$になることを確認する

という段取りを組めれば良い。
　できあがった点群から, 例えば $a_i=0$ となるものを選べば, 巡回平面の一つとなっており, そこから(巡回)完全差集合が得られて, 魔円陣を求める事ができる。
　やってみた。

　$A,B,C$の組を真面目に逐次試すと矢鱈と時間が掛かることが分かったので, 途中で, ランダムに$A,B,C$の組を設定して条件に合うかなどチェックしてというプログラムに変えた。
　こ、これは...強烈に早い...$GF(1009)$の3次拡大からの魔円陣が, 0.628秒だった。(乱数の具合にも拠るのだろうけどだ。)