bug#35374: Calc; Div by 0 in solving system

bug#35374: Calc; Div by 0 in solving system

[Christopher Howard]

[-- Attachment #1: Type: text/plain, Size: 662 bytes --]

Place this on stack with algebriac mode:

[2 * x - 3 * y + z = 5, x + y - 2 * z = 0, -x + 2 * y + 3 * z = -3]

run command 'a S x,y,z'

Receive error 'Division by zero: 1'

I enabled Debug on Error but could not find the backtrace

Emacs  : GNU Emacs 25.1.1 (x86_64-pc-linux-gnu, GTK+ Version 3.22.11)
 of 2017-09-14, modified by Debian
Package: Calc

-- 
Christopher Howard
p: +1 (907) 374-0257
w: https://librehacker.com
xmpp: creationist@member.fsf.org
otr: E9685B53 01F038DD D29281C9 30FDA71E BD0095D4
gnupg: 23FD5CC5 (keys.gnupg.net)
radio: KL1TL
featured: https://emailselfdefense.fsf.org/en/ id="-x-evo-selection-
start-marker">

[-- Attachment #2: Type: text/html, Size: 865 bytes --]

bug#35374: Calc; Div by 0 in solving system

[Mattias Engdegård]

[-- Attachment #1: Type: text/plain, Size: 343 bytes --]

tags 35374 patch
stop

> Maybe someone more knowledgeable about Calc can help further.

Much to our surprise, very few people volunteered their service. Now, despite being grossly under-qualified for the task, I did give it a go. After all, we can't have it fail on simple linear equation systems now can we.

Please try this patch.


[-- Attachment #2: 0001-Fix-linear-equation-system-solving-in-Calc-bug-35374.patch --]
[-- Type: application/octet-stream, Size: 6111 bytes --]

From 980a1f0b73d2aa58335a47ad4a68005a28307a4b Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Mattias=20Engdeg=C3=A5rd?= <mattiase@acm.org>
Date: Sun, 22 Sep 2019 15:03:02 +0200
Subject: [PATCH] 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.
---
 lisp/calc/calcalg2.el        |   6 ++
 test/lisp/calc/calc-tests.el | 103 +++++++++++++++++++++++++++++++++++
 2 files changed, 109 insertions(+)

diff --git a/lisp/calc/calcalg2.el b/lisp/calc/calcalg2.el
index 18243bfc74..2a716633ae 100644
--- a/lisp/calc/calcalg2.el
+++ b/lisp/calc/calcalg2.el
@@ -2417,6 +2417,12 @@ math-try-solve-for
 			    ((= (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 e1ee20b5d2..36a81dc2b7 100644
--- a/test/lisp/calc/calc-tests.el
+++ b/test/lisp/calc/calc-tests.el
@@ -215,6 +215,109 @@ calc-test-calendar
   (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
 
-- 
2.21.0 (Apple Git-122)

bug#35374: Calc; Div by 0 in solving system

[Christopher Howard]

I am now using Emacs 26.2 (through Guix) but am able to reproduce the
original error in 26.2. I patched just this function (math-try-solve-
for) and reloaded it. The error went away, and the correct solution was generated.

I did not try running the tests.

-- 
https://librehacker.com
Christopher Howard
p: +1 (907) 374-0257
social feed: https://gnusocial.club/librehacker
xmpp: creationist@member.fsf.org
otr: E9685B53 01F038DD D29281C9 30FDA71E BD0095D4
gnupg: 23FD5CC5 (keys.gnupg.net)
radio: KL1TL
featured: https://emailselfdefense.fsf.org/en/

On Sun, 2019-09-22 at 19:34 +0200, Mattias Engdegård wrote:
> tags 35374 patch
> stop
> 
> > Maybe someone more knowledgeable about Calc can help further.
> 
> Much to our surprise, very few people volunteered their service. Now,
> despite being grossly under-qualified for the task, I did give it a
> go. After all, we can't have it fail on simple linear equation
> systems now can we.
> 
> Please try this patch.
> 

bug#35374: Calc; Div by 0 in solving system

[Mattias Engdegård]

tags 35374 fixed
close 35374 26.4
stop

23 sep. 2019 kl. 07.28 skrev Christopher Howard <christopher.howard@qlfiles.net>:
> 
> I am now using Emacs 26.2 (through Guix) but am able to reproduce the
> original error in 26.2. I patched just this function (math-try-solve-
> for) and reloaded it. The error went away, and the correct solution was generated.

Thanks for verifying! Pushed to emacs-26.

bug#35374: Calc; Div by 0 in solving system

[Mattias Engdegård]

21 apr. 2019 kl. 17.37 skrev Christopher Howard <christopher.howard@qlfiles.net>:

> I enabled Debug on Error but could not find the backtrace

Reproduced in 26.1. Backtrace below.

I quickly got lost in calcalg2.el, but it looks like the equation

  2(z - 1.0) - (5 + 3(z - 1.0) - z) / 2 + 3z + 3 = 0,

isn't solved properly by math-try-solve-for: math-decompose-poly returns the correct decomposition of the LHS, (4) * z, but then the code decides to find a root of the polynomial (4) which of course doesn't exist and nothing good comes out of it. I'm not sure why the floating-point numbers (float 1 0) are there in the first place.

Maybe someone more knowledgeable about Calc can help further.

* math-reject-arg((float 1 0) "*Division by zero")
 math-div-by-zero((float 1 0) 0)
 math-div((float 1 0) 0)
 math-float((frac 1 0))
 math-mul((frac 1 0) (float 138629436111989 -14))
 math-mul-objects-fancy((frac 1 0) (cplx (float 138629436111989 -14) (float 314159265358979 -14)))
 math-mul((frac 1 0) (cplx (float 138629436111989 -14) (float 314159265358979 -14)))
 math-pow-fancy(-4 (frac 1 0))
 math-pow(-4 (frac 1 0))
 math-poly-laguerre-root((4) 0 t)
 math-poly-any-root((4) 0 t)
 math-try-solve-for((+ (+ (- (* 2 (- (var z var-z) (float 1 0))) (/ (- (+ 5 (* 3 (- (var z var-z) (float 1 0)))) (var z var-z)) 2)) (* 3 (var z var-z))) 3) 0 nil nil)
 math-try-solve-for((calcFunc-eq (+ (- (* 2 (- (var z var-z) (float 1 0))) (/ (- (+ 5 (* 3 (- (var z var-z) (float 1 0)))) (var z var-z)) 2)) (* 3 (var z var-z))) -3) 0 nil)
 math-solve-for((calcFunc-eq (+ (- (* 2 (- (var z var-z) (float 1 0))) (/ (- (+ 5 (* 3 (- (var z var-z) (float 1 0)))) (var z var-z)) 2)) (* 3 (var z var-z))) -3) 0 (var z var-z) nil)
 math-solve-system-rec((((calcFunc-eq (+ (- (* 2 (- (var z var-z) (float 1 0))) (/ (- (+ 5 (* 3 (- (var z var-z) (float 1 0)))) (var z var-z)) 2)) (* 3 (var z var-z))) -3))) ((var z var-z)) (((var y var-y) (- (var z var-z) (float 1 0))) ((var x var-x) (/ (- (+ 5 (* 3 (- (var z var-z) (float 1 0)))) (var z var-z)) 2))))
 math-solve-system-rec((((calcFunc-eq (- (+ (/ (- (+ 5 (* 3 (var y var-y))) (var z var-z)) 2) (var y var-y)) (* 2 (var z var-z))) 0)) ((calcFunc-eq (+ (- (* 2 (var y var-y)) (/ (- (+ 5 (* 3 (var y var-y))) (var z var-z)) 2)) (* 3 (var z var-z))) -3))) ((var y var-y) (var z var-z)) (((var x var-x) (/ (- (+ 5 (* 3 (var y var-y))) (var z var-z)) 2))))
 math-solve-system-rec((((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))) ((var x var-x) (var y var-y) (var z var-z)) nil)
 math-solve-system((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)) nil)
 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)))
 apply(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))))
 math-normalize((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))))
 math-simplify((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))))
 calc-normalize-fancy((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))))
 calc-normalize((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))))
 math-evaluate-expr((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))))
 mapcar(math-evaluate-expr ((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)))))
 #f(compiled-function () #<bytecode 0x4a71afd9>)()
 calc-do(#f(compiled-function () #<bytecode 0x4a71afd9>) 126)
 calc-evaluate(1)
 funcall-interactively(calc-evaluate 1)
 call-interactively(calc-evaluate nil nil)
 command-execute(calc-evaluate)