import Data.List (find)

rref :: Fractional a => [[a]] -> [[a]]
rref m = f m 0 [0 .. rows - 1]
  where rows = length m
        cols = length $ head m

        f m _    []              = m
        f m lead (r: rs)
            | indices == Nothing = m
            | otherwise          = f m' (lead' + 1) rs
          where indices = find p l
                p (col, row) = m !! row !! col /= 0
                l = [(col, row) |
                    col <- [lead .. cols - 1],
                    row <- [r .. rows - 1]]

                Just (lead', i) = indices
                newRow = map (/ m!! i!! lead') $ m !! i

                m' = zipWith g [0..] $
                    replace r newRow $
                    replace i (m!! r) m
                g n row
                    | n == r    = row
                    | otherwise = zipWith h newRow row
                  where h = subtract . (* row!! lead')

replace :: Int -> a -> [a] -> [a]

replace n e l = a ++ e: b
  where (a, _: b) = splitAt n l