mdeterminant   xtlang

(bind-func mdeterminant
  "matrix determinant - matrix must be square (= nrows ncols)"
  (lambda (m1:double* nrows:i64)
    (if (= nrows 1) (pref m1 0)
        (let ((det:double 0.0) (s:double 1.0)
              (c 0) (i 0) (j 0) (m 0) (n 0)
              (mb:double* (salloc (* nrows nrows))))
          (dotimes (c nrows)
            (set! m 0) (set! n 0)
            (dotimes (i nrows)
              (dotimes (j nrows)
                (pset! mb (+ (* i nrows) j) 0.0)
                (if (and (<> i 0) (<> j c))
                    (begin
                      (pset! mb (+ (* m (- nrows 1)) n)
                             (pref m1 (+ (* i nrows) j)))
                      (if (< n (- nrows 2))
                          (set! n (+ n 1))
                          (begin
                            (set! n 0)
                            (set! m (+ m 1)))))
                    (begin 1))))
            (set! det (+ det (* s (* (pref m1 c)
                                     (mdeterminant mb (- nrows 1))))))
            (set! s (* -1.0 s)))
          det))))