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

これはナイーブな方法では無理ですね。

理論的背景がよく分からない。ので、勉強中。

また関係ないことに興味が出る。