原始ピタゴラス数の生成
「ピタゴラスの定理(Wikipedia)」等にあるように、原始ピタゴラス数の生成は以下の式が有名ですが、今回はMatrixを用いて原始ピタゴラス数の木を作る方法を紹介します。
(a, b, c) = (m^2 − n^2, 2mn, m^2 + n^2)
原始ピタゴラス数の木
原始ピタゴラス数の木（高校数学の美しい物語） に定理と証明が記載されています。定理は以下のようなものです。
\begin{align}
&すべての原始ピタゴラス数は \\
&(3,4,5)に以下の行列A,B,Cの任意の組合せの積をM^nとすると \\
&(a,b,c)=(3,4,5)M^nで表せる \\ \\
&A=\left(
\begin{array}{rr}
1 & 2 & 2 \\
-2 & -1 & -2 \\
2 & 2 & 3 \\
\end{array}
\right) \\
&B=\left(
\begin{array}{rr}
1 & 2 & 2 \\
2 & 1 & 2 \\
2 & 2 & 3 \\
\end{array}
\right) \\
&C=\left(
\begin{array}{rr}
-1 & -2 & -2 \\
2 & 1 & 2 \\
2 & 2 & 3 \\
\end{array}
\right) \\
\end{align}
これを図に表すと以下のようになります。レベル4までが求められています。
幅優先探索のプログラムにすると、以下のようになります。
import numpy as np
A = np.matrix([[ 1, 2, 2],[-2,-1,-2],[ 2, 2, 3]])
B = np.matrix([[ 1, 2, 2],[ 2, 1, 2],[ 2, 2, 3]])
C = np.matrix([[-1,-2,-2],[ 2, 1, 2],[ 2, 2, 3]])
ABC = [A,B,C]
pt = np.matrix([3,4,5])
L = [[pt], []]
lp1,lp2,count = 0,1,1
for d in range(1,3+1):
print(d, L[lp1])
L[lp2] = []
for pt in L[lp1]:
for M in ABC:
L[lp2].append(pt @ M)
lp1,lp2 = lp2,lp1
# 1 [matrix([[3, 4, 5]])]
# 2 [matrix([[ 5, 12, 13]]), matrix([[21, 20, 29]]), matrix([[15, 8, 17]])]
# 3 [matrix([[ 7, 24, 25]]), matrix([[55, 48, 73]]), matrix([[45, 28, 53]]), matrix([[39, 80, 89]]), matrix([[119, 120, 169]]), matrix([[77, 36, 85]]), matrix([[33, 56, 65]]), matrix([[65, 72, 97]]), matrix([[35, 12, 37]])]
原始ピタゴラス数の木アルゴリズムの応用
このアルゴリズムは面白いのですが、既にピタゴラス数を効率的に生成するアルゴリズムはあるので、その実用性はあまりないのかなと思われますが。実はあるのです。
この定理は、$a^2+b^2=c^2$だけではなくて以下の式でも成り立つからです。
\large a^2+b^2=c^2+k
これを確かめるために以下の問題を考えてみます。
問題：$a^2+b^2=c^2-1$となる (a,b,c)の組の数をを$(a+b+c) \le 100$で求めよ
上のプログラムのMatrixの初期値ptsを(2,2,3)に変えることで簡単に求まります。ただし$a$と$b$を交換した答えが両方生成されるので、$a \le b$の物だけカウントしています。
pts = [np.matrix([2,2,3])]
print(0, pts[0], np.sum(pts[0]))
N = 100
L = [pts, []]
lp1,lp2,count = 0,1,1
for d in itertools.count(1):
L[lp2] = []
for pt in L[lp1]:
for M in ABC:
pt1 = pt @ M
if np.sum(pt1) <= N:
L[lp2].append(pt1)
if pt1[0,0] <= pt1[0,1]: # a <= b
print(d, pt1, np.sum(pt1))
count += 1
if len(L[lp2]) == 0: break
lp1,lp2 = lp2,lp1
print(f"N = {N}, count={count}")
# 0 [[2 2 3]] 7
# 1 [[4 8 9]] 21
# 1 [[12 12 17]] 41
# 2 [[ 6 18 19]] 43
# 2 [[18 30 35]] 83
# 3 [[ 8 32 33]] 73
N = 100, count=6
このアルゴリズムはEuler Project: Problem 223, 224を解くのに役に立ちます。
（開発環境：Google Colab）