Problem 66

http://projecteuler.net/index.php?section=problems&id=66
とりあえず、D=61でつまった。
考える。考える。考える。
それでもD=61の解がでない。
そこで、禁断のgoogle先生。
どうやらこの問題は Pell Equationというらしい。
でmathworld。
http://mathworld.wolfram.com/PellEquation.html

とりあえず、アルゴリズムは分かった。前の問題の関数が利用できる。
ということで

import Data.List 
import Data.Ratio
import Data.Ord

next n (_,(b,c)) = 
    let c' = (n-b*b) `div` c
        n' = floor.sqrt.fromIntegral $ n
        (a,r) = (n'+b)`divMod`c'
    in (a,(n'-r,c'))

expand n = expand' . iterate (next n) $ (n',(n',1))
    where n' = floor.sqrt.fromIntegral $ n
          expand' (x:y:zs) = (fst x, map fst.(y:).takeWhile (y/=) $ zs) 

dio d | odd.length $ b = (f $ a:b++b,d)
      | otherwise = (f $ a:b,d)
      where (a,b) = expand d
            f = numerator.foldr1(\x y->x+1/y).init.map fromIntegral

p066 = snd.maximumBy (comparing fst).map dio $ [1..1000] \\ map (^2) [1..31]

最大のxはなんと
16421658242965910275055840472270471049
これはナイーブな方法では無理ですね。
理論的背景がよく分からない。ので、勉強中。
また関係ないことに興味が出る。