3×3平方数魔方陣を探索してみた話
何の気なしにnoteの数学タグを探索していると、次の記事が目に入りました。
どうやら平方数のみで構成された3×3の魔方陣を作る試みのよう。さらに調べてみるとどうやらこのような魔方陣が存在するかは未解決問題のようです。
これに興味を持った私が、自分なりにこの問題について探索してみたというのがこの記事の内容です。
ただ、先に言っておくと、実は私は存在しないことの証明には興味がありません。というのも、長年未解決問題のままということは一朝一夕でひらめくアイデア程度では解けないと思ったからです(ラマヌジャンレベルなら別ですが)。なので、あくまで数値計算で魔方陣を探索してみよう、ということです。
魔方陣のパラメータの構成
まず、平方数魔方陣とはどういうものかをしっかり確認しましょう。表を作るのがめんどくさいので書きやすい行列で魔方陣を代用することにします。魔方陣の各成分を
$$
\begin{pmatrix}
a & b & c \\ d & e & f \\ g & h & i
\end{pmatrix}
$$
とし、それぞれの各成分に対して、
$$
a=A^2, b=B^2, \dots
$$
となる自然数$${A,B,C,D,E,F,G,H,I}$$が存在するとしましょう。(これらは互いに異なる数だとします。)この時、魔方陣が成立するということは、上の9個の変数に加えてある自然数$${S}$$を含めた10個の自然数が
$$
A^2+B^2+C^2=S \\
D^2+E^2+F^2=S \\
G^2+H^2+I^2=S \\
A^2+D^2+G^2=S \\
B^2+E^2+H^2=S \\
C^2+F^2+I^2=S \\
A^2+E^2+F^2=S \\
C^2+E^2+G^2=S
$$
の8個の条件式を満たすということです(実際には6つ目の式は上5つの式
から作れるので、独立な条件式は7つです)。
よくある普通のルールですね。
余談ですが、最近ルールが改変されたゲキムズ理詰めマインスイーパーやっているので、普通のルールという言葉に安心を覚えてしまいます。
パラメータの条件
条件の整理
さて、これらの条件式を整理すると、
$$
S=3E^2 \\
B^2+H^2=2 E^2\\
D^2+F^2=2 E^2\\
2A^2=F^2+H^2\\
2C^2=H^2+D^2\\
2I^2=D^2+B^2\\
2G^2=B^2+F^2
$$
となります。つまり、魔方陣は
1. $${E}$$を決定し、
2. 第2式および第3式を満たすように$${(B,H),(D,F)}$$の組たちを決めれば
あとは自動的に決まります。
Eへの制限
この時点で、もし$${2E^2}$$が$${2E^2=x^2+y^2}$$のように、二つの平方数に分解できなければ、あるいはできてもその方法が1通りしかなければその$${E}$$を使って平方数魔方陣を作ることはできないことになります(今回同じ数が含まれる魔方陣は考えていないので、$${2E^2=E^2+E^2}$$という自明な分解は無視しています)。
さらに、ナナメも考えると$${2E^2}$$は自明な分解を除いて少なくとも4通りの異なる分解を持たなければなりません。これが一つ目の重要なヒントですね。
分解のパターンの数というのはヤコビの二平方定理で簡単に(?)計算できます。詳しくは次のページを見てみてください。
ところで、分解のパターンの数が4つ以上でないといけないという話をしましたが、実は1つ以上か0かというのだけは(素因数分解できれば)簡単に求めることができます。実際に、Wikipediaで泳ぐだけでも次の記事が見つかります。
これらをまとめると、
素数$${p}$$が2つの平方数に分解できる必要十分条件は、$${p=2}$$、あるいは$${p\equiv 1\;(\text{mod} \,4)}$$である。また、その分解の方法は1通りである。
合成数が高々2つの平方数に分解できる必要十分条件は、その素因数分解は$${ q\equiv 3\;(\text{mod} \,4)}$$を満たすすべての素数$${q}$$に対して$${q^{2k-1}}$$の形を含まないことである。
となります(高々、と言っているのが少々面倒で、$${9=0^2+3^2}$$という分解も「高々2つの平方数に分解」に含まれるので、2つ目の定理を今回使う分には少し注意が必要です)。
ちなみに4で割った余りが1の素数はピタゴラス素数と名付けられており、平方数への分解で重要になる素数です。ヤコビの二平方定理内でも出没してますね。
しかし、もしかしたら「$${2E^2}$$は常に合成数でしかも常に$${q^{2k-1}}$$の形を含まないのだから、上の条件は今回の話には役に立たないのではないか」と思われるかもしれません。実際それはごもっともでこのままでは役に立ちません。「何のためにこの話をしているのか」と思われるかもしれませんね。
しかし、上の定理を発展させて次が言えるとしたらどうでしょうか。
定理?
自然数$${E}$$に対して$${2E^2=x^2+y^2}$$を満たす自然数$${x,y}$$が存在するとき、
$${E\equiv 0 \;(\text{mod} \;q) \implies x,y \equiv 0 \;(\text{mod} \;q)}$$
が成り立つ。ここで、$${q}$$は$${q\equiv 3\;(\text{mod}\,4)}$$を満たす素数である。
例えば$${E=35}$$とすると、これは$${35 \equiv 0 \;(\text{mod} \,7)}$$であるので$${2\times 35^2=2450}$$の分解は7の倍数のみで構成されているというのがこの定理の主張です。実際、$${2450}$$の分解は$${2450=7^2+49^2=35^2+35^2}$$のみなので成り立っていますね。
注意点として、これは証明できていないし、過去に先人が証明していたかは知らないのでもしかしたら反例が存在するかもしれないということに言及しておきます。お手数ですが見つけたら教えてください!
さて、上の定理が成り立っているとして、果たして今回の問題にどのようにかかわってくるのでしょうか。
試しに、Eが$${q}$$の倍数だとしましょう。この時定理から$${B,H,D,F}$$も$${q}$$の倍数となります。そして、$${2 A^2=B^2+H^2}$$などから$${A,C,G,I}$$も$${q}$$の倍数となります。すると…、なんと魔方陣の数字$${A^2,\dots,I^2}$$すべてが$${q^2}$$の倍数であることがわかります。
つまり、Eが$${q}$$の倍数なら全体も$${q^2}$$の倍数になるので魔方陣全体を$${q^2}$$で割ってより小さい数字の魔方陣を作ることができるのです。したがって、魔方陣を探索するうえで$${E}$$が$${q}$$の倍数の場合は考えなくてよいことになります。
実は上の定理は$${q}$$を2に置き換えても成り立っています(これは平方剰余を考えれば簡単に証明できます)。したがって2の倍数も考える必要もない考えるべき$${E}$$はさらに減らせます。
以上から、魔方陣を探索するうえでは(上の定理が成り立っているなら)$${E}$$の素因数がピタゴラス素数のみであるものだけを探索すれば良いことがわかります。このような数は大体全体の1割ぐらいなので考えるのが楽になる…かもしれません(実際はメモリの都合で現実的ではないです)。
実際にやってみた
プログラムをmathematicaで組んで実際に探索してみました。mathematicaにはPowersRepresentationsという平方数への分解を作ってくれる関数が標準で搭載されているのでコードを組む上ではそこまで難しいことはありません。あとは実際に探索するだけ…ではありますが単純に時間がかかるのでもっと制限を見つけるか効率的な探し方をしないと発見は厳しそうです。
今回私は$${E=2^{25}=33554432}$$まで探索して見つかりませんでした。最後にコードを載せておくのでもしmathematicaを購入されている方は試してみて下さい。
コード(mathematica)
listPrime =
ParallelMap[
If[Mod[#, 4] =!= 1, #, Nothing] &, (Prime /@ Range[2^20])];
divCheck[e_] := NoneTrue[listPrime, Divisible[e, #] &];
ParallelDo[
If[divCheck[e] && SquaresR[2, 2 e^2] >= 4 + 8*4,
With[{y = PowersRepresentations[2 e^2, 2, 2][[;; -2]]},
Table[If[
AllTrue[Tuples[{y[[i]], y[[j]]}],
IntegerQ[Sqrt[(#[[1]]^2 + #[[2]]^2)/2]] &],
PrintTemporary["Find! e = " <> ToString[e]];
Print[{{Sqrt[(y[[i, 2]]^2 + y[[j, 2]]^2)/2], y[[i, 1]],
Sqrt[(y[[i, 2]]^2 + y[[j, 1]]^2)/2]}, {y[[j, 1]], e,
y[[j, 2]]}, {Sqrt[(y[[i, 1]]^2 + y[[j, 2]]^2)/2], y[[i, 2]],
Sqrt[(y[[i, 1]]^2 + y[[j, 1]]^2)/2]}}]], {i, Length[y]}, {j,
i + 1, Length[y]}]]], {e, 1, 2^10}]
探索する範囲を変えたい場合は最後の2^10を別の数字に変えて下さい。
ちなみに、$${2^{10}}$$までだと10秒くらいで終わります。
コードの中身を説明すると、
まずピタゴラス素数ではない小さめの素数を生成し、それらで割れないときに真を返す関数divCheckを定義しています。
あとは$${E}$$を順番に当てはめてdivCheckを通過し、平方数への分割が自明なのを除いて4パターン以上あるときのみ具体的に構成してみて成功したらPrintで出力するというものになります。