Project Euler Problem 100

import Data.Ratio
import Data.List
import Data.Maybe

-- チャンドラグプタの方法
gupta1 :: Num t => t -> (t, t) -> (t, t)
gupta1 n (a,b) = (a*a+n*b*b,2*a*b) -- gupta identity self application
gupta2 :: Num t => t -> (t, t) -> (t, t) -> (t, t)
gupta2 n (a,b) (c,d) = (a*c+n*b*d,a*d+b*c)
lis :: [(Integer, Integer)]
lis = iterate (gupta2 2 (3,1)) (3,1)

-- chakravala methodの実装
chak :: Integral t => t -> [(t, t)]
chak n = iterate (gupta2 n (numerator a',numerator b')) (numerator a',numerator b')
  where
    (a',b',c') = imp seed
    lim = fromJust $ find (\x -> (x^2)>n) [1..]
    seed = findfirst n
    gupta (a,b,k) = (a*a+(n%1)*b*b,2*a*b,k*k) -- gupta identity self application
    imp x@(a,b,k) | k == 1%1 = (a,b,k)
                  | k == (-1)%1 = gupta x
                  | otherwise = imp $ findabkm n x
    -- 自明な初期値を求める。b=1固定でa,kを設定
    findfirst n = (a%1,1%1,k%1) 
      where 
        a = snd $ minimum [(abs(x^2 - n),x)|x<-[1..(lim+1)]]
        k = a^2 - n
    -- chakravala method繰り返し部分
    findabkm n (a,b,k) = (aa,bb,kk)
      where
        list = filter (\x -> (numerator (a+b*(x%1)))`rem`(numerator k)==0) [1..(lim+1)]
        m = snd $ minimum [(abs(x^2 - n),x%1)|x<-list]
        -- bhaskara
        aa = (a*m+(n%1)*b)/(abs k)
        bb = (a+b*m)/(abs k)
        kk = (m*m-(n%1))/k
--(3,1)が初期解のペル値を調べる
pells :: [((Integer, Integer), Integer)]
pells = takeWhile ((/=(3,1)).fst) $ map (\(l,r) -> (head (chak l),r)) $ filter (\(x,y) -> x /= (round $ sqrt $ fromIntegral x)^2) (zip [2..] [2..])
-- n=8となる
-- 後はrが求まったときにbを求める関係式を求めれば終わり
lis2 :: [(Integer, Integer)]
lis2 = iterate (gupta2 8 (3,1)) (3,1)
-- (2086968960817,313506783024) = (blue,red)
blues :: [Integer]
blues = tail $ map (\(t,r) -> (t - 2*r + 1)`quot`2) lis2 
anslist :: [(Integer, Integer)]
anslist = zip blues (map snd lis2)
ans1imp :: (Integer, Integer)
ans1imp = anslist !! ((fromJust $ findIndex (\(l,r) -> l+r>10^12) anslist))
ans1 :: Integer
ans1 = let (a,_) =ans1imp in a
-- 756872327473
main :: IO ()
main = print ans1