Project Euler 65,66

Problem 65

http://projecteuler.net/index.php?section=problems&id=65

連分数の計算法を適用する。最初に1を付加すると計算しやすい。

digits 0 = []
digits n = (mod n 10):(digits (div n 10))

a = 2:1:2:1:[ if n == 1 then 1 else n + 2 | n <- tail a ]
p = 1:(head a):[ a1 * p1 + p0 | (a1, p0, p1) <- zip3 (tail a) p (tail p) ]
main = print ((sum . digits) (p!!100))

Problem 66

http://projecteuler.net/index.php?section=problems&id=66

ペル方程式の初期解の解法を適用する。

int_square = floor . sqrt . fromIntegral
is_square n = n == m * m where m = int_square n

next n (m0,(a,b,c)) = (m, (div a2 d, div b2 d, div c2 d)) where
        b1 = b - m * c
        c2 = a * a * n - b1 * b1
        a2 = a * c
        b2 = -b1 * c
        d = gcd a2 (gcd b2 c2)
        m = floor (((sqrt (fromIntegral n)) * (fromIntegral a)
                            + (fromIntegral b)) / (fromIntegral c))

minimal_solution n =
        head [ x | (x, y) <- zip (tail p) (tail q), is_satisfied x y ] where
    a n = (0,(1,0,1)):[ next n x | x <- a n]
    r = map fst (tail (a n))
    p = 1:(head r):[ a * p1 + p0 | (a, p0, p1) <- zip3 (tail r) p (tail p) ]
    q = 0:1:[ a * q1 + q0 | (a, q0, q1) <- zip3 (tail r) q (tail q) ]
    is_satisfied x y = x^2 - n * (y^2) == 1

greater x y = if snd x >= snd y then x else y
a = [ (n, minimal_solution n) | n <- [2..1000], not (is_square n) ]
main = print (fst (foldr greater (0, 0) a))