(inv-mat-fn size)
Produces optimized code for a function that inverts matrices
of the given size. The resulting code has no loops, branches, or
recursive calls. The inverse is computed using Cramer's rule and
the required determinants are computed using recursive cofactor
expansion along the first column. All intermediate determinants
are bound to local variables since many of them are used multiple
times.
Source
(defmacro inv-mat-fn
"Produces optimized code for a function that inverts matrices
of the given size. The resulting code has no loops, branches, or
recursive calls. The inverse is computed using Cramer's rule and
the required determinants are computed using recursive cofactor
expansion along the first column. All intermediate determinants
are bound to local variables since many of them are used multiple
times."
[size]
(let [sub-det-syms (memoize (fn [rows cols]
(with-meta (gensym)
;; we only tag the matrix entries
(when-not (next rows) {:tag 'double}))))
m (vec (for [i (range size) j (range size)]
(sub-det-syms `(~i) `(~j))))]
`(fn ~[m]
(let ~(mat-inv-fn-binds size sub-det-syms)
~(vec (for [i (range size) j (range size)]
(let [detA (sub-det-syms (range size) (range size))
detA_ij (sub-det-syms (remove #(= j %) (range size))
(remove #(= i %) (range size)))]
(if (even? (+ i j))
`(/ ~detA_ij ~detA)
`(- (/ ~detA_ij ~detA))))))))))