Pell's equation solver

pell.hs

import Data.List

-- ペル方程式 x^2 - d * y^2 = ±1について、
-- 指定されたdに対して最小の解(x,y)を返す
-- 参考：http://www004.upp.so-net.ne.jp/s_honma/pell/pell.htm
minSolutionForPellEquation d
  | even n    = ((x_m0^2+d*y_m0^2) `div` h_m,
                 (2*x_m0*y_m0)     `div` h_m)
  | otherwise = ((x_m0*x_m1+d*y_m0*y_m1) `div` h_m,
                 (x_m0*y_m1+x_m1*y_m0)   `div` h_m)
  where
    gs = 0 : zipWith3 gEquation gs ks hs
    gEquation g k h = - g + k * h
    hs = 1 : zipWith hEquation (tail gs) hs
    hEquation g h = (d - g^2) `div` h
    k0 = floor . sqrt $ fromIntegral d
    ks = k0 : zipWith kEquation (tail gs) (tail hs)
    kEquation g h = (k0+g) `div` h
    ys = 0 : 1 : zipWith3 yEquation ys (tail ks) (tail ys)
    yEquation y1 k y2 = y1 + k*y2
    xs = 1 : zipWith4 xEquation (tail gs) (tail ys) (tail hs) ys
    xEquation g y2 h y1 = g*y2 + h*y1
    n = moirIndex gs hs
    m = n `div` 2
    h_m = hs !! m
    x_m0 = xs !! m
    x_m1 = xs !! (m+1)
    y_m0 = ys !! m
    y_m1 = ys !! (m+1)

-- gとhの数列を受け取って、h_n = 1となる最小のnを返す。
-- Moirの定理より。
-- 参考：http://www004.upp.so-net.ne.jp/s_honma/pell/pell.htm
moirIndex gs hs = moirIndex' gs hs 0
  where
    moirIndex' gs@(g1:g2:_) hs@(h1:h2:_) m
      | g1 == g2 = 2*m
      | h1 == h2 = 2*m+1
      | otherwise = moirIndex' (tail gs) (tail hs) (m+1)