６６問目

http://projecteuler.net/index.php?section=problems&id=66
不定方程式 x² - Dy² = 1 について、D=13 のとき x の最小値は 649² - 13×180² = 1 となる。Dが平方数のとき、自然数の範囲では解がないと想定されている。D ≦ 7 で平方数でないDに対してxの最小値を考えると、D=5 のときに最大値9が得られる。D ≦ 1000 としたとき、xが最大となるようなDの値を求める。
1000以下ならbrute forceで行けるだろうとか思ったら、甘かった。題意の方程式はペル方程式というらしく、〜プログラム@Bal4u〜: ペル方程式の最小解　pellEquとかペル方程式の解の列挙方法 - まめめもとかペル方程式の解の構成とか見て、漸化式を立てて解きました。Webを探せば、答えそのものも出てきます。
基本的には上述のリンク先に従って組めばいいのですが、このコード特有の問題として、13, 14行目があります。このように as', xs' を明示的に定義しないと、リストだと判断されずに毎回関数の値を求めてしまい、メモ化が機能しないのです。

import Data.List

main = print $ euler066 1000

euler066 :: Integer -> Integer
euler066 m = snd $ maximum
    [ (last $ xs i, i) | i <- [1..m] \\ takeWhile (<= m) [ n^2 | n <- [1..] ] ]

xs :: Integer -> [Integer]
xs n = 1 : (head as') : zipWith3 (\ a x1 x0 -> a * x1 + x0) (tail as') (tail xs') xs'
    where as' = as n
          xs' = xs n

as :: Integer -> [Integer]
as n = (\ (s, t, a) -> a) $ unzip3 $ take (index * ((mod index 2) + 1)) rems'
    where index = (findIndices (\ (s, t, a) -> t == 1) rems') !! 1
          rems' = rems n

rems :: Integer -> [(Integer, Integer, Integer)]
rems n = iterate next (0, 1, r)
    where r = floor $ sqrt $ fromIntegral n
          next (s, t, a) = let s' = a * t - s
                               t' = div (n - s' ^ 2) t
                               a' = div (r + s') t'
                           in (s', t', a')

コメントで添削していただいたので、修正版を。主なポイントは↓です。

• (s, t, a) を (s, (t, a)) とすることで、不要な s を無視。
• index が偶数なら index まで、奇数なら index*2 までというルールを lcm 2 index で表現。
• x1 と x0 を (x1, x0) とセットで扱うことで、zipWith3 を scanl で簡潔に表現できる。さらに最後の要素しか見ていないので実は foldl でよい。

変なメモ化は要らないし、速度は当社比４倍です。

import Data.List

main = print $ euler066 1000

euler066 :: Integer -> Integer
euler066 m = snd $ maximum [ (f i, i) 
    | i <- takeWhile (<= m) $ concat [ [n^2+1..(n+1)^2-1] | n <- [1..] ] ]

f :: Integer -> Integer
f n = fst $ foldl (\ (x1, x0) a  -> (a * x1 + x0, x1)) (1, 0) $ as n

as :: Integer -> [Integer]
as n = snd $ unzip $ take (lcm 2 index) $ rems n
    where index = (!! 1) $ findIndices ((== 1) . fst) $ rems n

rems :: Integer -> [(Integer, Integer)]
rems n = snd $ unzip $ iterate next (0, (1, r))
    where r = floor $ sqrt $ fromIntegral n
          next (s, (t, a)) = let s' = a * t - s
                                 t' = div (n - s' ^ 2) t
                                 a' = div (r + s') t'
                             in (s', (t', a'))