(defun mmul (A B)
  (let* ((m (car (array-dimensions A)))
         (n (cadr (array-dimensions A)))
         (l (cadr (array-dimensions B)))
         (C (make-array `(,m ,l) :initial-element 0)))
    (loop for i from 0 to (- m 1) do
              (loop for k from 0 to (- l 1) do
                    (setf (aref C i k)
                          (loop for j from 0 to (- n 1)
                                sum (* (aref A i j)
                                       (aref B j k))))))
    C))

;; Cofactor i j
(defun matrix-cofactor (A i j)
    (let* (
        (tmp 0)
        (n (car (array-dimensions A)))
        (ans (make-array `(,(- n 1) ,(- n 1)) :initial-element 0)))
    (loop for ii from 0 to (- i 1) do 
        (loop for jj from 0 to (- j 1) do (setf (aref ans ii jj) (aref A ii jj)))
        (loop for jj from (+ j 1) to (- n 1) do (setf (aref ans ii (- jj 1)) (aref A ii jj))))
    (loop for ii from (+ i 1) to (- n 1) do 
        (loop for jj from 0 to (- j 1) do (setf (aref ans (- ii 1) jj) (aref A ii jj)))
        (loop for jj from (+ j 1) to (- n 1) do (setf (aref ans (- ii 1) (- jj 1)) (aref A ii jj))))
    ans))

;; Calculate det value of matrix A
(defun matrix-det (A)
    (let* ((n (car (array-dimensions A))) (tmp 0))
    (if (= n 1)
        (aref A 0 0)
        (let ((tmp 0) (MA 0))
            (loop for j from 0 to (- n 1) do
                (setf MA (matrix-cofactor A 0 j))
                (if (= (mod j 2) 0)
                    (setf tmp (+ tmp (* (aref A 0 j) (matrix-det MA))))
                    (setf tmp (- tmp (* (aref A 0 j) (matrix-det MA))))))
            tmp))))

;; Calculate inverse matrix of matrix A
(defun matrix-inverse (A)
    (let* ((n (car (array-dimensions A))) (tmp 0) (inv-A (make-array `(,n ,n) :initial-element 0)) (det-val 0))
        (setf det-val (matrix-det A))
        (loop for i from 0 to (- n 1) do
            (loop for j from 0 to (- n 1) do
                (setf (aref inv-A j i) (/ (matrix-det (matrix-cofactor A i j)) det-val))))
        inv-A))

(let ((g 0) (g0 0))
  (setf g (make-array '(3 3) :initial-contents '((1 0 0) (2 2 0) (3 4 5))))
  (setf g0 (make-array '(1 1) :initial-contents '((2))))
  ; (matrix-cofactor g 1 1)
  (format t "~A~%" (matrix-inverse g)))