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

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

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

原文 [Problem 94:Almost equilateral triangles](https://projecteuler.net/problem=94)

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

3辺の長さと面積が全て整数となる三角形を[ヘロンの三角形](https://ja.wikipedia.org/wiki/%E3%83%98%E3%83%AD%E3%83%B3%E3%81%AE%E4%B8%89%E8%A7%92%E5%BD%A2)と言うみたいですが、この問題では正三角形に近い$a=b=c\pm 1$のものを探す必要がありますが、これを以下のような手順で解いていきます。

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

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

まずヘロンの公式です。

```math
\begin{align}
S&=\sqrt{s(s-a)(s-b)(s-c)} \\
s&=\frac{a+b+c}{2} \\ \\
\end{align}
```
$s$も整数なのでcは偶数。したがって$c=2k$とおいて変形して行きます。

```math
\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で解く](https://qiita.com/masa0599/items/162c7e1989ae9a653a0c)」で紹介したようにsympyのdiop_DNを使う方法もありますが。今回は「[【Project Euler】Problem 66: ディオファントス方程式](https://qiita.com/masa0599/items/963b758c518e37716aab)」の中の「ペル方程式の解を平方根の連分数展開から求める」でやったように平方根の連分数から求める方法で作ったコードをもとにして作って行きます。

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

```python
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に収束する分数](https://qiita.com/masa0599/items/5daa44b06e25c514d0b6)」で作った連分数から収束する分数を求める関数**cf2frac**を変更しました。

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

```python
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$まで走らせた結果を表示しました。

```python
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
```