From: Mattias Engdegård <mattiase@acm.org>
Date: Wed, 13 May 2020 13:17:10 +0000 (+0200)
Subject: Calc: fix LU decomposition for non-numeric matrices (bug#41223)
X-Git-Tag: emacs-28.0.90~7364
X-Git-Url: http://git.eshelyaron.com/gitweb/?a=commitdiff_plain;h=1d559581b3bcd91644e44b1e3a3788614d99924f;p=emacs.git

Calc: fix LU decomposition for non-numeric matrices (bug#41223)

Computing determinant and inverse for on some matrices containing
non-numeric elements failed or gave the wrong result.
Reported by Mauro Aranda.

* lisp/calc/calc-mtx.el (math-do-matrix-lud): Don't use zero as pivot.
* test/lisp/calc/calc-tests.el (calc-matrix-determinant): New test.
---

diff --git a/lisp/calc/calc-mtx.el b/lisp/calc/calc-mtx.el
index fe241b57c60..2850b33721b 100644
--- a/lisp/calc/calc-mtx.el
+++ b/lisp/calc/calc-mtx.el
@@ -275,7 +275,7 @@ in LUD decomposition."
 		k (1+ k)))
 	(setcar (nthcdr j (nth i lu)) sum)
 	(let ((dum (math-lud-pivot-check sum)))
-	  (if (Math-lessp big dum)
+	  (if (or (math-zerop big) (Math-lessp big dum))
 	      (setq big dum
 		    imax i)))
 	(setq i (1+ i)))
diff --git a/test/lisp/calc/calc-tests.el b/test/lisp/calc/calc-tests.el
index 6db5426ff6d..9e36d91ac3e 100644
--- a/test/lisp/calc/calc-tests.el
+++ b/test/lisp/calc/calc-tests.el
@@ -345,6 +345,30 @@ An existing calc stack is reused, otherwise a new one is created."
   (should (Math-num-integerp '(float 1 0)))
   (should-not (Math-num-integerp nil)))
 
+(ert-deftest calc-matrix-determinant ()
+  (should (equal (calcFunc-det '(vec (vec 3)))
+                 3))
+  (should (equal (calcFunc-det '(vec (vec 2 3) (vec 6 7)))
+                 -4))
+  (should (equal (calcFunc-det '(vec (vec 1 2 3) (vec 4 5 7) (vec 9 6 2)))
+                 15))
+  (should (equal (calcFunc-det '(vec (vec 0 5 7 3)
+                                     (vec 0 0 2 0)
+                                     (vec 1 2 3 4)
+                                     (vec 0 0 0 3)))
+                 30))
+  (should (equal (calcFunc-det '(vec (vec (var a var-a))))
+                 '(var a var-a)))
+  (should (equal (calcFunc-det '(vec (vec 2 (var a var-a))
+                                     (vec 7 (var a var-a))))
+                 '(* -5 (var a var-a))))
+  (should (equal (calcFunc-det '(vec (vec 1 0 0 0)
+                                     (vec 0 1 0 0)
+                                     (vec 0 0 0 1)
+                                     (vec 0 0 (var a var-a) 0)))
+                 '(neg (var a var-a)))))
+
+
 (provide 'calc-tests)
 ;;; calc-tests.el ends here