import Control.Monad (mplus,liftM)

type Matrix a = [[a]]

inverse :: (Monad m, Fractional a) => Matrix a -> m (Matrix a)
inverse m
    | null m       = fail "empty matrix" 
    | any (/=w) ws = fail "matrix not rectangular"
    | h /= w       = fail "matrix not quadractic"
    | otherwise    = liftM (map (drop h)) . clean . gauss $ e
    where h    = length m
	  w:ws = map length m
	  e    = zipWith (++) m $ munit h

munit n = 
    let k = n-1
	rep = flip replicate 0
    in [0..k] >>= \ i -> [rep i ++ [1] ++ rep (k-i)]

gauss m = 
    flip (maybe m) (pivot id m) $ \ (f,row:rows) ->
    map f . (row:) . map (0:) . gauss . map tail $ rows

pivot adjust m
    | all null m = fail "no pivot at all"
    | otherwise  = mplus (pivot1 m >>= return . (,) adjust)
		         (pivot ((0:) . adjust) (map tail m))

pivot1 (row@(0:_):rows) = pivot1 rows >>= \ (r:rs) -> return (r:row:rs)
pivot1 (row@(p:_):rows) = return . (first:) . map add $ rows
    where first         = map (/p) row
          add r@(x:_)   = sub r x first
pivot1 _                = fail "no pivot in first colum"

sub v c = zipWith (-) v . map (*c)

clean (row@(1:rt):rows) = 
    liftM (map (0:)) (clean . map tail $ rows) >>= \ crows ->
    return $ (foldl ( \ v (a,w) -> sub v a w ) row $ zip rt crows) : crows
clean (    (_: _):_   ) = fail "matrix singular"
clean m                 = return m