【Project Euler】Problem 94: 正三角形に近いヘロンの三角形

• 本記事はProjectEulerの「100番以下の問題の説明は記載可能」という規定に基づいて回答のヒントが書かれていますので、自分である程度考えてみてから読まれることをお勧めします。

問題： 正三角形に近いヘロンの三角形

原文 Problem 94:Almost equilateral triangles

問題の要約：　三辺の長さが整数$a,b,c$で$a=b=c\pm 1$で、面積も整数になる三角形で、三辺の和が$10^9$を超えないもの三辺の和の合計を求めよ

3辺の長さと面積が全て整数となる三角形をヘロンの三角形と言うみたいですが、この問題では正三角形に近い$a=b=c\pm 1$のものを探す必要がありますが、これを以下のような手順で解いていきます。

1. ヘロンの公式（三辺から面積を求める）使って方程式（ペル方程式）を作る
2. 【ペル方程式を解く（その１）】平方根の連分数展開
3. 【ペル方程式を解く（その２）】連分数から解を発生させる
4. ペル方程式の解から元のa,b,cを求める

ヘロンの公式（三辺から面積を求める）使って方程式を作る

まずヘロンの公式です。

\begin{align}
S&=\sqrt{s(s-a)(s-b)(s-c)} \\
s&=\frac{a+b+c}{2} \\ \\
\end{align}

$s$も整数なのでcは偶数。したがって$c=2k$とおいて変形して行きます。

\begin{align}
&c = 2k、a=b=2k \pm 1\\
&したがって\\
&s=3k\pm 1 \\
S&=\sqrt{(3k\pm 1)\times k \times k \times (k \pm 1)} \\
&= k \times \sqrt{3k^2 \pm 4k + 1} \\
&Sが整数なので、ルートの中が平方数\\
&3k^2 \pm 4k + 1=y^2 \ \ とおく \\
&これを変形すると　\\
&(3k \pm 2)^2-1=3y^2 \\
&すると以下のようにペル方程式となります。\\
&x^2-3y^2=1　(x = 3k \pm 2)
\end{align}

このペル方程式を満たす整数$(x,y)$を求め、$x$から$k$そして周囲の長さを求めます。

【ペル方程式を解く（その１）】平方根の連分数展開

ペル方程式の解を生成していく方法には「ペル型方程式をsympyのdiop_DNで解く」で紹介したようにsympyのdiop_DNを使う方法もありますが。今回は「【Project Euler】Problem 66: ディオファントス方程式」の中の「ペル方程式の解を平方根の連分数展開から求める」でやったように平方根の連分数から求める方法で作ったコードをもとにして作って行きます。

まず平方根を連分数展開する関数cfsqrtで3の平方根の連分数$[1, [1, 2]]$が求まります。

def cfsqrt(m):
  n0 = int(m**(1/2))
  if n0*n0 == m: return [n0]
  n,a,b,cf = n0,1,0,[]
  while True:
    b = n*a -b
    a = (m-b*b)//a
    n = (n0+b)//a
    cf.append(n)
    if a == 1: break
  return [n0,cf]

cfp = cfsqrt(3)
print(cfp)
# [1, [1, 2]]

【ペル方程式を解く（その２）】連分数からペル方程式の解を発生させる

次にこの連分数からペル方程式の解を発生させる関数cfp2pellを作ります。このコードは「【Project Euler】Problem 65: ネイピア数eに収束する分数」で作った連分数から収束する分数を求める関数cf2fracを変更しました。

循環部分（今回は[1,2]の部分）の長さが偶数か奇数かで処理が変わってきます。奇数の場合には変数cntを使って2回に1回出力するようにしました。yieldを使っているのでnext関数を呼ぶごとに次のペル方程式の解が出てくる仕様です。

def cfp2pell(cfp):
  if len(cfp) == 1: yield (cfp[0],1)
  p, prd = 0, cfp[1] 
  n0, l = cfp[0], len(prd)
  cnt = 0 if l>1 else 1   # alternate counter for l is odd case
  pn_2, pn_1, qn_2, qn_1, an_1 = 1, n0, 0, 1, prd[p]

  while True:
    pn, qn = an_1 * pn_1 + pn_2, an_1 * qn_1 + qn_2
    #print(p, cnt, pn, qn)
    if p==max(l-2,0) : # if p = kl-1
      if l%2==1: cnt = 1-cnt  # anternate count
      if cnt == 0: yield (pn, qn)
      
    p = (p+1)%l 
    an_1, pn_1, pn_2, qn_1, qn_2 =  prd[p], pn, pn_1, qn, qn_1

ペル方程式の解から元のa,b,cを求める

2つの関数cfsqrtとcfp2pellを使ってメインを作ります。cfp2pellから生成されたペル方程式の解$(x,y)$を$k$に変換します。$\pm$があるので二つの値の各々の整数チェックをして周囲の長さ（$6k \pm 2$）が上限になるまで足して行きます。
以下のプログラムでは$10^3$まで走らせた結果を表示しました。

cfp = cfsqrt(3)
xyPell = cfp2pell(cfp)
peri_sum, cont = 0, True
while cont:
  (x,y) = next(xyPell)
  for i in [-1,1]:
    k = x+2*i
    if k%3==0:  # check k is multiple of 3
      k = k // 3
      if k == 0: continue
      S, peri = k*y, 6*k-2*i 
      if peri > 10**3: 
        cont = False
        break
      print(f"a,b={2*k-i}, c={2*k}, peri={peri}, S={S}")
      peri_sum += peri
print(F"answer = {peri_sum}")
# a,b=5, c=6, peri=16, S=12
# a,b=17, c=16, peri=50, S=120
# a,b=65, c=66, peri=196, S=1848
# a,b=241, c=240, peri=722, S=25080
# answer = 984