Fix linear equation system solving in Calc (bug#35374)

* lisp/calc/calcalg2.el (math-try-solve-for):
To solve Ax^n=0 where A is a nonzero constant and x the variable to
solve for, solve x^n=0 instead of solving A=0 (which obviously fails)
or something equally stupid.
* test/lisp/calc/calc-tests.el (calc-test-solve-linear-system): New.

2 files changed, 109 insertions, 0 deletions

diff --git a/lisp/calc/calcalg2.el b/lisp/calc/calcalg2.el
index 18243bfc749..2a716633ae6 100644
--- a/lisp/calc/calcalg2.el
+++ b/lisp/calc/calcalg2.el
@@ -2417,6 +2417,12 @@
                             ((= (length math-t1) 2)
                              (apply 'math-solve-linear
                                     (car math-t2) math-try-solve-sign math-t1))
+                            ((= (length math-t1) 1)
+                             ;; Constant polynomial.
+                             (if (eql (nth 2 math-t2) 1)
+                                 nil    ; No possible solution.
+                               ;; Root of the factor, if any.
+                               (math-try-solve-for (nth 2 math-t2) 0 nil t)))
                             (math-solve-full
                              (math-poly-all-roots (car math-t2) math-t1))
                             (calc-symbolic-mode nil)
diff --git a/test/lisp/calc/calc-tests.el b/test/lisp/calc/calc-tests.el
index e1ee20b5d2d..36a81dc2b71 100644
--- a/test/lisp/calc/calc-tests.el
+++ b/test/lisp/calc/calc-tests.el
@@ -215,6 +215,109 @@ An existing calc stack is reused, otherwise a new one is created."
   (should (equal (math-absolute-from-julian-dt -101 3 1) -36832))
   (should (equal (math-absolute-from-julian-dt -4713 1 1) -1721425)))
 
+(ert-deftest calc-test-solve-linear-system ()
+  "Test linear system solving (bug#35374)."
+  ;;   x + y =   3
+  ;;  2x - 3y = -4
+  ;; with the unique solution x=1, y=2
+  (should (equal
+           (calcFunc-solve
+            '(vec
+              (calcFunc-eq (+ (var x var-x) (var y var-y)) 3)
+              (calcFunc-eq (- (* 2 (var x var-x)) (* 3 (var y var-y))) -4))
+            '(vec (var x var-x) (var y var-y)))
+           '(vec (calcFunc-eq (var x var-x) 1)
+                 (calcFunc-eq (var y var-y) 2))))
+
+  ;;  x + y = 1
+  ;;  x + y = 2
+  ;; has no solution
+  (should (equal
+           (calcFunc-solve
+            '(vec
+              (calcFunc-eq (+ (var x var-x) (var y var-y)) 1)
+              (calcFunc-eq (+ (var x var-x) (var y var-y)) 2))
+            '(vec (var x var-x) (var y var-y)))
+           '(calcFunc-solve
+             (vec
+              (calcFunc-eq (+ (var x var-x) (var y var-y)) 1)
+              (calcFunc-eq (+ (var x var-x) (var y var-y)) 2))
+             (vec (var x var-x) (var y var-y)))))
+  ;;   x - y = 1
+  ;;   x + y = 1
+  ;; with the unique solution x=1, y=0
+  (should (equal
+           (calcFunc-solve
+            '(vec
+              (calcFunc-eq (- (var x var-x) (var y var-y)) 1)
+              (calcFunc-eq (+ (var x var-x) (var y var-y)) 1))
+            '(vec (var x var-x) (var y var-y)))
+           '(vec (calcFunc-eq (var x var-x) 1)
+                 (calcFunc-eq (var y var-y) 0))))
+  ;;  2x - 3y +  z =  5
+  ;;   x +  y - 2z =  0
+  ;;  -x + 2y + 3z = -3
+  ;; with the unique solution x=1, y=-1, z=0
+  (should (equal
+           (calcFunc-solve
+            '(vec
+              (calcFunc-eq
+               (+ (- (* 2 (var x var-x)) (* 3 (var y var-y))) (var z var-z))
+               5)
+              (calcFunc-eq
+               (- (+ (var x var-x) (var y var-y)) (* 2 (var z var-z)))
+               0)
+              (calcFunc-eq
+               (+ (- (* 2 (var y var-y)) (var x var-x)) (* 3 (var z var-z)))
+               -3))
+            '(vec (var x var-x) (var y var-y) (var z var-z)))
+           ;; The `float' forms in the result are just artefacts of Calc's
+           ;; current solver; it should be fixed to produce exact (integral)
+           ;; results in this case.
+           '(vec (calcFunc-eq (var x var-x) (float 1 0))
+                 (calcFunc-eq (var y var-y) (float -1 0))
+                 (calcFunc-eq (var z var-z) 0))))
+  ;;   x = y + 1
+  ;;   x = y
+  ;; has no solution
+  (should (equal
+           (calcFunc-solve
+            '(vec
+              (calcFunc-eq (var x var-x) (+ (var y var-y) 1))
+              (calcFunc-eq (var x var-x) (var y var-y)))
+            '(vec (var x var-x) (var y var-y)))
+           '(calcFunc-solve
+             (vec
+              (calcFunc-eq (var x var-x) (+ (var y var-y) 1))
+              (calcFunc-eq (var x var-x) (var y var-y)))
+             (vec (var x var-x) (var y var-y)))))
+  ;;  x + y + z = 6
+  ;;  x + y     = 3
+  ;;  x - y     = 1
+  ;; with the unique solution x=2, y=1, z=3
+  (should (equal
+           (calcFunc-solve
+            '(vec
+              (calcFunc-eq (+ (+ (var x var-x) (var y var-y)) (var z var-z)) 6)
+              (calcFunc-eq (+ (var x var-x) (var y var-y)) 3)
+              (calcFunc-eq (- (var x var-x) (var y var-y)) 1))
+            '(vec (var x var-x) (var y var-y) (var z var-z)))
+           '(vec
+             (calcFunc-eq (var x var-x) 2)
+             (calcFunc-eq (var y var-y) 1)
+             (calcFunc-eq (var z var-z) 3))))
+  ;; x = 3
+  ;; x + 4y^2 = 3                   (ok, so this one isn't linear)
+  ;; with the unique (double) solution x=3, y=0
+  (should (equal
+           (calcFunc-solve
+            '(vec
+              (calcFunc-eq (var x var-x) 3)
+              (calcFunc-eq (+ (var x var-x) (* 4 (^ (var y var-y) 2))) 3))
+            '(vec (var x var-x) (var y var-y)))
+           '(vec (calcFunc-eq (var x var-x) 3)
+                 (calcFunc-eq (var y var-y) 0)))))
+
 (provide 'calc-tests)
 ;;; calc-tests.el ends here