mmmul   xtlang

(bind-func mmmul
  (lambda (m1:double* r1:i64 c1:i64 m2:double* r2:i64 c2:i64 result:double*)
    (if (= c1 r2)
        (let ((r 0) (c 0) (k 0))
          (dotimes (r r1)
            (dotimes (c c2)
              (pset! result (+ c (* r c2)) 0.0)
              (dotimes (k c1)
                (pset! result (+ c (* r c2))
                       (+ (pref result (+ c (* r c2)))
                          (* (pref m1 (+ k (* r c1)))
                             (pref m2 (+ (* k c2) c))))))))
          result)
        (begin (printf "bad matrix multiplication [%lld,%lld]*[%lld,%lld]\n"
                       r1 c1 r2 c2)
               (cast null double*)))))