wrong linsolve solution to system with parameters

[rtoy]

Imported from SourceForge on 2024-07-08 01:54:28 Created by cypgay on 2020-08-23 21:04:05 Original: https://sourceforge.net/p/maxima/bugs/3651

---

I tumbled on something bizarre and boiled it down to the following 2 equation system with 2 unknowns and two parameters. The result provided by linsolve is generally incorrect. The code below provides the correct solution worked out by matrix inversion for comparison. Thanks for your help.

/* unexpected behaviour of linsolve
for a system of two equations with parameters */

/* no specific value given to K : the error remains */
/* for two obvious values of K, the system is singular : */
/* K: -sqrt(3); */
/* K: -1; */ 
/* for all other values of K the error disappears : */
/* K: sqrt(3); */
/* K: 1; */

/* system of two equations eq1,eq2 and two unknowns x,y
with two parameters K,e */
eq1: (K+sqrt(3))*(x-y);
eq2: K*e+(K+1)*(y+x);

/* extracting the matrix of the system */
A: diff(eq1,x);
B: diff(eq1,y);
E: -diff(eq1,e);
C: diff(eq2,x);
D: diff(eq2,y);
F: -diff(eq2,e);

/* inverting the matrix */
determ: factor(A*D-B*C);
Ainv: D/determ;
Dinv: A/determ;
Binv: -B/determ;
Cinv: -C/determ;

/* solution from inverse matrix */
xok: factor((Ainv*E+Binv*F)*e);
yok: factor((Cinv*E+Dinv*F)*e);

/* checking the matrix solution */
check_eq1: factor(subst([x=xok,y=yok],eq1));
check_eq2: factor(subst([x=xok,y=yok],eq2));

/* no specific value given to e : the error remains */
/* error only for some values of e */
/* e: sqrt(3) */ /* the error remains */
/* e: 1 */ /* the error disappears */

/* applying the specific value before solving */
eq1: ev(eq1);
eq2: ev(eq2);

/* solving using linsolve */
sol:linsolve([factor(eq1),factor(eq2)],[x,y]);

/* checking linsolve solution : these two lines should yield zero */
factor(subst(sol,eq1));
factor(subst(sol,eq2));

[rtoy]

Imported from SourceForge on 2024-07-08 01:54:29 Created by cypgay on 2020-08-23 21:09:50 Original: https://sourceforge.net/p/maxima/bugs/3651/#670f

---

Forgot to say : the sqrt(3) is important : it seems to work if that becomes an integer or a fraction. Version used from the command line : Maxima 5.41.0 (Ubntu 18.04.5 package).

[rtoy]

Imported from SourceForge on 2024-07-08 01:54:32 Created by robert_dodier on 2020-09-04 06:14:40 Original: https://sourceforge.net/p/maxima/bugs/3651/#df98

---

Indeed, sqrt(3) triggers the bug, likewise sqrt(5), sqrt(7), etc. However, sqrt(z) (with z being a symbol), or any simple symbol or number, or sqrl(3) or sqrt3 (i.e. similar to but not sqrt(3)) all work OK.

I find that the following patch is enough to get the correct result.

diff --git a/src/mat.lisp b/src/mat.lisp
index 772509c0a..a87d7724d 100644
--- a/src/mat.lisp
+++ b/src/mat.lisp
@@ -241,9 +241,9 @@
           (mess2 (aref ax (aref *row* i) (aref *col* i))  ))
       (cond ((equal 0 mess1) 0)
         ((equal 0 mess2) 0)
-        (t   ;;   (pquotient mess1 mess2) ; fixed by line below. RJF 1/12/2017
+        (t   (pquotient mess1 mess2) ; fixed by line below. RJF 1/12/2017

-         (car (ratreduce mess1 mess2))
+         ;;   (car (ratreduce mess1 mess2))
          )
         ))))
     (do ((l (1+ i) (1+ l)))

The call to RATREDUCE was introduced in commit b6dd08c which was to fix a bug which was reported to the mailing list, see: https://sourceforge.net/p/maxima/mailman/message/35594033/

I find that changing it back to PQUOTIENT causes it to fail on the example mentioned in the bug report, so I guess I can't just change it back. I'll investigate some more. I'm not sure what RATREDUCE is doing, I guess that's the first thing to look at.