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 fe241b57c60fb10d57b8619352c7464b420c7c64..2850b33721b6c5c0d6e4ebfd1bb0f6baa728707b 100644 (file)
--- 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 6db5426ff6dd4fef23f29cf15fde3732905b5ed1..9e36d91ac3eeb7ff9183acaf26aeffbbdfd9dff2 100644 (file)
--- 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