ピタゴラス数の小ネタ

今までプログラマやってきて特に気にすることもなかったけど、project euler を始めてみるとピタゴラス数に関連する問題に出会うことがある。で、いろいろ調べてみると、ちょっと面白いネタというか Tips が見つかった。

■ 原始ピタゴラス数の無限リスト

このサイトで原始ピタゴラス数（primitive Pythagorean triple (PPT)）を次々に生成する行列が紹介されている。こんな行列。

（ -1　 -2　 2）
（ -2 　-1　 2）
（ -2　 -2　 3）

この行列に符号を3通りに変えたピタゴラス数を掛けると、子のピタゴラス数が3つ生成される。例えば(3,4,5)から始めると以下のようになる。

(-1 -2 2) (-3)   ( 5)
(-2 -1 2) ( 4) = (12) 
(-2 -2 3) ( 5)   (13)

(-1 -2 2) ( 3)   (15)
(-2 -1 2) (-4) = ( 8) 
(-2 -2 3) ( 5)   (17)

(-1 -2 2) (-3)   (21)
(-2 -1 2) (-4) = (20) 
(-2 -2 3) ( 5)   (29)

算出された子のピタゴラス数から、さらに３つずつの孫を作ることができ、これを再帰的に繰り返すとピタゴラス数の無限リストを作ることができる。Haskell で書くとこんな感じになる。

pyt1 (a,b,c) = ((-1)*a+(-2)*b+2*c,(-2)*a+(-1)*b+2*c,(-2)*a+(-2)*b+3*c)
pyt2 (a,b,c) = [pyt1(-a,b,c), pyt1(-a,-b,c), pyt1(a,-b,c)]
pyt3 ts = ts ++ pyt3 (concatMap pyt2 ts)

以下のような結果になる。

ghci< pyt3 [(3,4,5)]
[(3,4,5),(5,12,13),(21,20,29),(15,8,17),(7,24,25),(55,48,73),(45,28,53),(39,80,89),(119,120,169),(77,36,85),(33,56,65),(65,72,97),(35,12,37),(9,40,41),(105,88,137),(91,60,109),(105,208,233),(297,304,425),(187,84,205),(95,168,193),(207,224,305),(117,44,125),(57,176,185),(377,336,505),(299,180,349),(217,456,505),(697,696,985),(459,220,509),(175,288,337),(319,360,481),(165,52,173),(51,140,149),(275,252,373),(209,120,241),(115,252,277),(403,396,565),(273,136,305),(85,132,157),(133,156,205),(63,16,65),(11,60,61),(171,140,221),(153,104,185),(203,396,445),(555,572,797),・・・

■ 高さ＝底辺の二等辺三角形に漸近する漸化式

このサイトでは、ピタゴラス数を直角二等辺三角形に近づけていく漸化式が紹介されていた。直角を構成する辺 a と b の差を1に こんな式。

an+2 = 6*an+1 - an + 2  a1 = 3, a2 = 20
cn+2 = 6*bn+1 - bn - 2  b1 = 4, b2 = 21
cn+2 = 6*bn+1 - cn      c1 = 5, c2 = 29

で、ここから類推して高さ＝底辺の二等辺三角形（の片側の三角形）に漸近する漸化式もきっとあるだろうなと思って探してみた（きっかけは euler 138）。

まず直角を構成する2辺のうち長い辺と短い方を二倍したものとの差が1であるものを何個か探す。これは上述の無限リストに対して、条件に合うものをフィルタリングした。以下のようなものが見つかった（これ以上は数が大きくなりすぎて難しい）。

(15,8,17)
(273,136,305)
(4895,2448,5473)
(87841,43920,98209)
(1576239,788120,1762289)
(28284465,14142232,31622993)
(507544127,253772064,567451585)
(9107509825,4553754912,10182505537)

で見つかったピタゴラス数から連立方程式を立ててこれを解いた。手で計算するのは面倒なので、ネット上の連立方程式を解くプログラムを使った。例えば、斜辺ならこんな式を解けば良い。

305*x + 17y + z = 5473
5473*x + 305*y + z= 98209
98209*x + 5473*y + z =1762289

解）x=18, y=-1, z=0

で、c₁ = 17, c₂ = 305, c_n+2 = 18*c_n+1 - c_n という漸化式が得られる。Haskell で書くとこうなる。

l = 17:305:zipWith (\a b->18*b-a) l (tail l)

同様に、直角を構成する2辺についても漸化式が得られる。 a1 = 15, a2 = 273, a_n+3 = 17*a_n+2 + 17*a_n+1 - a_n (bも同様)

試しに20番目の要素を見てみると、a= 10151021471800938910964641, b = 5075510735900469455482320, c = 11349187026003431978487841)で、a² + b² = c², 1 = |2b-a| が確かに成立する。

==== ==== ==== ====

さっき problem 138 を解いて、euler project を150問解いた事になる。簡単な問題を選んでやってきたが、そろそろ難しくなってきた。かなり頑張ったけど、まだ半分にも至ってないのか…