ScalaでProjectEuler39を解いてみた. 後で自分のコードを見て理解できる気がしないので メモ.

問題をPukiWiki からコピー.

辺の長さが \({a,b,c}\) と整数の3つ組である直角三角形を考え, その周囲の長さを \(p\) とする. \(p = 120\) のときには3つの解が存在する:

{20,48,52}, {24,45,51}, {30,40,50}

\(p \le 1000\) のとき解の数が最大になる \(p\) はいくつか?

解法

$p_{\mathrm{max}} = 1000$と置く.

\({a, b, c}\)はピタゴラスの三角形の三組なので 一意的な自然数 \(u, v, g\)があって

  • $u > v$が成り立つ.
  • \(u\)と\(v\)の偶奇が異なる.
  • \(u\)と\(v\)は互いに素である.
  • $ a = g(u^2 - v^2), b = 2guv, c = g(u^2 + v^2) $が成り立つ.

この$u, v, g$を使うと

\[\text{(周囲の長さ)} = a + b + c = 2gu(u + v)\]

が成り立つ.

$u, v$及び$g$の値を評価したい.

\[4gv^2 = 2gv(v + v) < 2gu(u + v) \le p_{\mathrm{max}}\]

すなわち

\[v < \frac{1}{2} \sqrt{\frac{p_{\mathrm{max}}}{g}} \le \frac{\sqrt{p_{\mathrm{max}}}}{2}.\]

ここで $p_{\mathrm{max}} = 1000$を代入すると

\[v < 5 \sqrt{10} \fallingdotseq 15.8114\]

となるので $v \le 15$ となる. 15は十分小さいので 十分評価されていると考えられる.

次に$u$を評価する.

\[2u^2 + 2vu \le 2gu(u + v) \le p_{\mathrm{max}}\]

より

\[2u^2 + 2vu - p_{\mathrm{max}} \le 0.\]

これより

\[\begin{align} u & \le \frac{-v + \sqrt{v^2 + 2p_{\mathrm{max}}}}{2} \\ & < \frac{-v + \sqrt{(v + p_{\mathrm{max}})^2}}{2} \\ & = \frac{-v + (v + p_{\mathrm{max}})}{2} = \frac{p_{\mathrm{max}}}{2} \end{align}\]

が成り立つ.

最後に

\[2gu(u + v) \le p_{\mathrm{max}}\]

より

\[g \le \frac{p_{\mathrm{max}}}{2u(u + v)}.\]

方針

  • \(v \le \frac{\sqrt{p_{\mathrm{max}}}}{2}\) のそれぞれに対して
    • $ v + 1 \le u \le \frac{p_{\mathrm{max}}}{2} $を
      • $u$ と $v$ が互いに素
      • $u$と$v$の偶奇が異なる
    • ようにとる(それらを$u_s$とおく)
    • $u_s$それぞれに対して $g_s = (1 \text{ to } \frac{p_{\mathrm{max}}}{2u(u + v)})$とおく
  • $ a + b + c \le p $ となっているか確認する

ソースコード

package projecteuler.exercises

import projecteuler.Answer

object Ex039 extends Answer {
  def gcd(a: Int, b: Int): Int = if(b ==0) a else gcd(b, a%b)
  def relativelyPrime(a: Int, b: Int): Boolean = gcd(a, b) == 1

  def result(pMax: Int) =
    (
    for (
      v <- (1 to ((Math.sqrt(pMax) / 2)).toInt);
      u <- (v + 1 to pMax / 2).by(2) if relativelyPrime(u, v);
      g <- (1 to pMax / (2 * u * (u + v))))
        yield (g * (u * u - v * v), 2 * g * u * v, g * (u * u + v * v))
    )
    .filter(t => t._1 + t._2 + t._3 <= pMax)
    .groupBy(t => t._1 + t._2 + t._3)
    .maxBy(e => e._2.size)

  def run(): Any = {
    result(1000)._1 // => 840
  }
}