ScalaでProject Euler（75）

Problem 45

三角数は六角数を含むので、三角数は考えなくてもよいことになります。
単に五角数と六角数を昇順に並べてマージ法を使うのが簡単でしょう。

// Polygonal number
def P(p :Int) :Stream[Long] =
    Stream.from(1).map(n => n * ((p - 2L) * n - p + 4) / 2)

def coincident(s1 :Stream[Long], s2 :Stream[Long]) :Stream[Long] =
                                                    (s1, s2) match {
    case (h1 #:: t1, h2 #:: t2) if h1 < h2 => coincident(t1, s2)
    case (h1 #:: t1, h2 #:: t2) if h1 > h2 => coincident(s1, t2)
    case (h1 #:: t1, h2 #:: t2)            => h1 #:: coincident(t1, s2)
    case _ => Stream()
}

println (coincident(P(5), P(6)).take(3).toList)

しかし、数学を使うと速くなります。

P_m = H_n
m(3m - 1) / 2 = n(2n - 1)
36m² - 12m = 3(16n² - 8n)
(6m - 1)² - 1 = 3(4n - 1)² - 3
(6m - 1)² - 3(4n - 1)² = -2

この式で、

x = 6m - 1
y = 4n - 1

と置くと、

x² - 3y² = -2

となり、ほぼペル方程式で、

x_k + y_k√3 = (2 + √3)^k(1 + √3)

と計算できます。このうち、

x_k ≡ -1(mod 6)
y_k ≡ -1(mod 4)

を満たすものが解です。

def gen_Pell() = {
    // (2 + √3)^n(1 + √3)
    def next(xy :(Int,Int)) = {
        val (x, y) = xy
        (2 * x + 3 * y, x + 2 * y)
    }
    Iterator.iterate((1, 1))(next)
}

val it_sol = for((x, y) <- gen_Pell if x % 6 == 5 && y % 4 == 3;
                    m = (x + 1) / 6)
                yield m * (3L * m - 1) / 2

println (it_sol.take(3).toList)