pell'sEquation :: Integer -> [(Integer, Integer)]
pell'sEquation d = iterate (f (x, y)) (x, y)
  where
    f (a1, b1) (a2, b2) = (a1 * a2 + d * b1 * b2, a1 * b2 + a2 * b1)
    (x, y) = loop 0 1 k0 1 0
      where
        k0 = (truncate . sqrt . fromIntegral) d
        loop g1 h1 k1 y0 y1
          | g1 == g2  = g $ f (x1, y1) (x1, y1)
          | h1 == h2  = let ns = g $ f (x1, y1) (x2, y2) in f ns ns
          | otherwise = loop g2 h2 k2 y1 y2
          where
            g2 = k1 * h1 - g1
            h2 = div (d - g2^2) h1
            k2 = div (k0 + g2) h2
            y2 = y0 + k1 * y1
            x1 = g1 * y1 + h1 * y0
            x2 = g2 * y2 + h2 * y1
            g (x, y) = (div x h1, div y h1)
pell'sEquation d = iterate (f (x, y)) (x, y)
  where
    f (a1, b1) (a2, b2) = (a1 * a2 + d * b1 * b2, a1 * b2 + a2 * b1)
    (x, y) = loop 0 1 k0 1 0
      where
        k0 = (truncate . sqrt . fromIntegral) d
        loop g1 h1 k1 y0 y1
          | g1 == g2  = g $ f (x1, y1) (x1, y1)
          | h1 == h2  = let ns = g $ f (x1, y1) (x2, y2) in f ns ns
          | otherwise = loop g2 h2 k2 y1 y2
          where
            g2 = k1 * h1 - g1
            h2 = div (d - g2^2) h1
            k2 = div (k0 + g2) h2
            y2 = y0 + k1 * y1
            x1 = g1 * y1 + h1 * y0
            x2 = g2 * y2 + h2 * y1
            g (x, y) = (div x h1, div y h1)