type Vector a = [a]
type Matrix a = [[a]]
 
gauss :: Fractional a => Matrix a -> Vector a -> Vector a
gauss a b =
    let m = zipWith join a b
          where join xs y = xs ++ [y]
        n = length b
        m2 = foldl deal_with_jth_row m [1..n-1]
        x = accum_scanr1 solve_ith_eqn m2
    in 
        x
 
deal_with_jth_row :: Fractional a => Matrix a -> Int -> Matrix a
deal_with_jth_row m j = first_j_rows ++ (map put_0_in_jth_pos other_rows)
         where
           first_j_rows = take j m
           other_rows   = drop j m
           put_0_in_jth_pos row =  zipWith join row (m !! (j-1))
               where
                 join x y = mult * y + x
                 mult = -(row !! (j-1) / m !! (j-1) !! (j-1))
 
 
solve_ith_eqn :: Fractional a => Vector a -> Vector a -> a
solve_ith_eqn row zs = 
    (foldl (+) (row!!n) (zipWith f (drop i row) zs)) / (row!!(i-1))
   where n = (length row) - 1
         i = n - length zs
         f x y = -(x * y)
 
 
accum_scanr1 :: (a -> [b] -> b) -> [a] -> [b]
accum_scanr1 f zs = map head (scanr' g [] zs)
   where 
        g x ys = (f x ys):ys
        scanr' f a xs = init(scanr f a xs)