#3208 Can't solve this simple equation..

[nightwing]

solve([x^2+y^2=1, (x-2)^2+(y-3)^2=9],[x,y]);
then it returns []
why?

[David Billinghurst]

This is an algsys failure.

algsys([x^2+y^2=1,(x-2)^2+(y-3)^2=9],[x,y]);
[]

Duplicate of [bugs:#3150] and [bugs:#2736]
Probably related to [bugs:#1106] [bugs:#1266] [bugs:#1302] [bugs:#1369] [bugs:#1370] [bugs:#1469] [bugs:#1647] [bugs:#2038]

The ev(...,radcan) bug found below is [bugs:#3223]

and from the mailing list:

1. http://def.fe.up.pt/pipermail/maxima-discuss/2012/041218.html - with long discussion and interesting comments
2. http://def.fe.up.pt/pipermail/maxima-discuss/2012/041342.html

Also consider [bugs:#196] after we fix above.

[David Billinghurst]

I have started to explore this. It interesting to compare

algsys([x^2+y^2=1,(x-2)^2+(y-3)^2=9],[x,y]);
[]

and

algsys([x^2+y^2=1, (x-2)^2+(y-3)^2=8],[x,y]);
[[x = 12/13,y = 5/13],[x = 0,y = 1]]

Tracing the higher level functions in algsys

?trace(?algsys);
?trace(?algsys0);
?trace(?algsys1);
?trace(?condensesolnl);
?trace(?distrep);
?trace(?lofactors);
?trace(?findleastvar);
?trace(?bakalevel);
?trace(?presultant);

I can see that what is attempted is equivalent to (but not using sqrtdenest)

(%i4) load(sqdnst);
(%o4) "/usr/share/maxima/5.38.0/share/simplification/sqdnst.mac"
(%i5) eq1:y^2+x^2-1;
(%o5) y^2+x^2-1
(%i6) eq2:y^2-6*y+x^2-4*x+4;
(%o6) y^2-6*y+x^2-4*x+4
(%i7) r:resultant(eq1,eq2,y);
(%o7) 52*x^2-40*x-11
(%i8) X:solve(r,x);
(%o8) [x = -(3^(5/2)-10)/26,x = (3^(5/2)+10)/26]
(%i9) Y1:solve(subst(X[1],eq2),y);
(%o9) [y = -(sqrt(4077-28*3^(7/2))-78)/26,
       y = (sqrt(4077-28*3^(7/2))+78)/26]
(%i10) Y2:solve(subst(X[2],eq2),y);
(%o10) [y = -(sqrt(28*3^(7/2)+4077)-78)/26,
        y = (sqrt(28*3^(7/2)+4077)+78)/26]
(%i11) subst([X[1],Y1[1]],[eq1,eq2])$
(%i12) expand(sqrtdenest(%));
(%o12) [0,0]
(%i13) subst([X[1],Y1[2]],[eq1,eq2])$
(%i14) expand(sqrtdenest(%));
(%o14) [378/13-(4*3^(5/2))/13,0]
(%i15) subst([X[2],Y2[1]],[eq1,eq2])$
(%i16) expand(sqrtdenest(%));
(%o16) [0,0]
(%i17) subst([X[2],Y2[2]],[eq1,eq2])$
(%i18) expand(sqrtdenest(%));
(%o18) [(4*3^(5/2))/13+378/13,0]

This gives solutions [x = -(3^(5/2)-10)/26, y = -(sqrt(4077-28*3^(7/2))-78)/26] and [x = (3^(5/2)+10)/26, y = -(sqrt(28*3^(7/2)+4077)-78)/26]

Comparing the traces in the two cases above, I think that the failure occurs in the back-substitution of the solutions for x and y into eq1. I'm guessing that algsys doesn't simplify the expression to 0 and therefore rejects the correct solutions. More digging required.

[Aleksas]

--> eq1:x^2+y^2=1$ eq2:(x-2)^2+(y-3)^2=9$
solve([eq1,eq2-eq1],[x,y]);
(%o3) [[x=-(3^(5/2)-10)/26,y=(23^(5/2)+45)/78],[x=(3^(5/2)+10)/26,y=-(23^(5/2)-45)/78]]

[Hisham Sueyllam]

Also you can do the following:load(grobner)$poly_reduced_grobner([x^2+y^2-1, (x-2)^2+(y-3)^2-9], [x,y])$algsys(%,[x,y])$
to get your answer without having to figure out what to make out of the 2 eqs.
 hisham...

On Thursday, September 22, 2016 6:59 AM, Aleksas <aleksasd@users.sf.net> wrote:

--> eq1:x^2+y^2=1$ eq2:(x-2)^2+(y-3)^2=9$
solve([eq1,eq2-eq1],[x,y]);
(%o3) [[x=-(3^(5/2)-10)/26,y=(23^(5/2)+45)/78],[x=(3^(5/2)+10)/26,y=-(23^(5/2)-45)/78]] [bugs:#3208] Can't solve this simple equation..Status: open
Group: None
Labels: algsys solve
Created: Sun Sep 04, 2016 06:12 AM UTC by nightwing
Last Updated: Thu Sep 22, 2016 02:06 AM UTC
Owner: nobodysolve([x^2+y^2=1, (x-2)^2+(y-3)^2=9],[x,y]);
then it returns []
why?Sent from sourceforge.net because you indicated interest in https://sourceforge.net/p/maxima/bugs/3208/To unsubscribe from further messages, please visit https://sourceforge.net/auth/subscriptions/

[David Billinghurst]

Thanks. I am looking to fix the bugs in algsys, rather than a work around for this particular problem. As you can see above, there are many bug reports for similar simple problems.

[David Billinghurst]

More tracing. algsys is recursive. Annotated trace attached.

The test case above can be reduced to the following

display2d:false$
eq1:y^2+x^2-1$
eq2:y^2-6*y+x^2-4*x+4$
X2:x = (3^(5/2)+10)/26$

eq3:expand(676*subst(X2,eq2));
676*y^2-4056*y-28*3^(7/2)+2007
eq4:expand(676*subst(X2,eq1));
676*y^2+20*3^(5/2)-333

algsys([eq3,eq4],[y]);
[]

This has a solution

load(sqdnst)$
Y2:y = -(sqrt(28*3^(7/2)+4077)-78)/26$
expand(sqrtdenest(subst(Y2,eq3)));
 0
expand(sqrtdenest(subst(Y2,eq4)));
 0

[David Billinghurst]

Getting closer now.

The problem seems to be the call to lisp function presultant in algsys.lisp.

(defun presultant (p1 p2 var)
(cadr (ratf ($resultant (pdis p1) (pdis p2) (pdis (list var 1 1))))))

The resultant(eq3,eq4,y) above is zero. However, when we trace the lisp calls we get a non-zero result from presultant. Don't understand why.

(%i1) display2d:false$

(%i2) eq3:676*y^2-4056*y-28*3^(7/2)+2007$

(%i3) eq4:676*y^2+20*3^(5/2)-333$

(%i4) resultant(eq3,eq4,y);

(%o4) 0
(%i5) (?trace(?algsys),
?trace(?algsys0),
?trace(?algsys1),
?trace(?condensesolnl),
?trace(?distrep),
?trace(?lofactors),
?trace(?findleastvar),
?trace(?presultant),
?trace(?bakalevel),
?trace(?bakalevel1),
?trace(?mergesoln),
?trace(?callsolve),
?trace(?baksubst),
?trace(?ebaksubst))$

(%i6) algsys([eq3,eq4],[y]);
  0: (ALGSYS
      ((#:Y431 2 676 1 -4056 0 (#:|3^(1/2)432| 7 -28 0 2007))
       (#:Y431 2 676 0 (#:|3^(1/2)432| 5 20 0 -333))))
    1: (LOFACTORS (#:Y431 2 676 1 -4056 0 (#:|3^(1/2)432| 7 -28 0 2007)))
    1: LOFACTORS returned
         ((#:Y431 2 676 1 -4056 0 (#:|3^(1/2)432| 7 -28 0 2007)))
    1: (LOFACTORS (#:Y431 2 676 0 (#:|3^(1/2)432| 5 20 0 -333)))
    1: LOFACTORS returned ((#:Y431 2 676 0 (#:|3^(1/2)432| 5 20 0 -333)))
    1: (DISTREP
        (((#:Y431 2 676 1 -4056 0 (#:|3^(1/2)432| 7 -28 0 2007)))
         ((#:Y431 2 676 0 (#:|3^(1/2)432| 5 20 0 -333)))))
      2: (CONDENSESOLNL
          (((#:Y431 2 676 0 (#:|3^(1/2)432| 5 20 0 -333))
            (#:Y431 2 676 1 -4056 0 (#:|3^(1/2)432| 7 -28 0 2007)))))
      2: CONDENSESOLNL returned
           (((#:Y431 2 676 0 (#:|3^(1/2)432| 5 20 0 -333))
             (#:Y431 2 676 1 -4056 0 (#:|3^(1/2)432| 7 -28 0 2007))))
    1: DISTREP returned
         (((#:Y431 2 676 0 (#:|3^(1/2)432| 5 20 0 -333))
           (#:Y431 2 676 1 -4056 0 (#:|3^(1/2)432| 7 -28 0 2007))))
    1: (ALGSYS0
        ((#:Y431 2 676 0 (#:|3^(1/2)432| 5 20 0 -333))
         (#:Y431 2 676 1 -4056 0 (#:|3^(1/2)432| 7 -28 0 2007))))
      2: (ALGSYS1
          ((#:Y431 2 676 0 (#:|3^(1/2)432| 5 20 0 -333))
           (#:Y431 2 676 1 -4056 0 (#:|3^(1/2)432| 7 -28 0 2007))))
        3: (FINDLEASTVAR
            ((#:Y431 2 676 0 (#:|3^(1/2)432| 5 20 0 -333))
             (#:Y431 2 676 1 -4056 0 (#:|3^(1/2)432| 7 -28 0 2007))))
        3: FINDLEASTVAR returned
             ((#:Y431 2 676 1 -4056 0 (#:|3^(1/2)432| 7 -28 0 2007)) . #:Y431)
        3: (PRESULTANT (#:Y431 2 676 0 (#:|3^(1/2)432| 5 20 0 -333))
                       (#:Y431 2 676 1 -4056 0 (#:|3^(1/2)432| 7 -28 0 2007))
                       #:Y431)
        3: PRESULTANT returned
             (#:|3^(1/2)432| 14 358269184 12 511813120 10 182790400 7
              -59882135040 5 179646405120 0 -1201064537088)

[David Billinghurst]

now

(%i3) ?trace(?$resultant)$

(%i4) resultant(eq3,eq4,y);
  0: ($RESULTANT
      ((MPLUS SIMP) 2007 ((MTIMES SIMP) -28 ((MEXPT SIMP) 3 ((RAT SIMP) 7 2)))
       ((MTIMES SIMP) -4056 $Y) ((MTIMES SIMP) 676 ((MEXPT SIMP) $Y 2)))
      ((MPLUS SIMP) -333 ((MTIMES SIMP) 20 ((MEXPT SIMP) 3 ((RAT SIMP) 5 2)))
       ((MTIMES SIMP) 676 ((MEXPT SIMP) $Y 2)))
      $Y)
  0: $RESULTANT returned
       ((MTIMES RATSIMP) 7311616
        ((MPLUS RATSIMP)
         ((MTIMES RATSIMP) 49 ((MEXPT RATSIMP) ((MEXPT) 3 ((RAT) 1 2)) 14))
         ((MTIMES RATSIMP) 70 ((MEXPT RATSIMP) ((MEXPT) 3 ((RAT) 1 2)) 12))
         ((MTIMES RATSIMP) 25 ((MEXPT RATSIMP) ((MEXPT) 3 ((RAT) 1 2)) 10))
         ((MTIMES RATSIMP) -8190 ((MEXPT RATSIMP) ((MEXPT) 3 ((RAT) 1 2)) 7))
         ((MTIMES RATSIMP) 24570 ((MEXPT RATSIMP) ((MEXPT) 3 ((RAT) 1 2)) 5))
         -164268))
(%o4)                                  0

mmm. Add a resimplify to presultant. (not sure if this is a good idea. try anyway).

(defun presultant (p1 p2 var)
  (cadr (ratf (resimplify ($resultant (pdis p1) (pdis p2) (pdis (list var 1 1)))))))

Now the offending call to presultant returns 0. algsys still returns [].

[David Billinghurst]

Next problem is in ebaksubst when we backsubstitute solution Y2 into eq4 and the result isn't simplified to 0.

(%i1) load(sqdnst)$

(%i2) eq4:676*y^2+20*3^(5/2)-333$

(%i3) Y2:y = -(sqrt(28*3^(7/2)+4077)-78)/26$

(%i4) subst(Y2,eq4);

(%o4) (78-sqrt(28*3^(7/2)+4077))^2+20*3^(5/2)-333
(%i5) expand(sqrtdenest(%));

(%o5) 0

Calling sqrtdenest within ebaksubst works for this simple case, and algsys returns the correct solution for this case. It still fails on others.

(defun sqrtdenest (e) (meval `(($sqrtdenest) ,e)))

(defun ebaksubst (solnl lhsl)
  (mapcar #'(lambda (q) (cadr (ratf (sqrtdenest (what-the-$ev (pdis q)
                          (cons '(mlist) solnl)
                          '$radcan)))))
      lhsl))

We have identified two three problems:

1. the maxima function resultant returns an expression that can be simplified to 0, but this isn't done when it is called from lisp. An apparently non-zero resultant causes algsys to reject genuine solutions. This may also cause expression swell within algsys. Fixed with commit [5fd2f5]
2. ebaksubst returns expressions that can be simplified to zero by sqrtdenest, and need to be to detect solutions. These expressions occur for simple problems such as the intersection of two circles. Fix for simple cases in commit [2986fc]. More in [640ca7]
3. ev(...,radcan) introduces spurious complex numbers. Example below, also [bugs:#3223]. Fixed by avoiding inappropriate use of radcan [640ca7]

[Aleksas]

(%i2) eq1:x^2+y^2=1$ eq2:(x-2)^2+(y-3)^2=9$
(%i3) solve([eq1,eq1-eq2],[x,y]);
(%o3)
[[x=-(3^(5/2)-10)/26,y=(23^(5/2)+45)/78],[x=(3^(5/2)+10)/26,y=-(23^(5/2)-45)/78]]

2016-09-23 16:28 GMT+03:00 David Billinghurst billingd@users.sf.net:

Next problem is in ebaksubst when we backsubstitute solution Y2 into eq4
and the result isn't simplfied to 0.

(%i1) load(sqdnst)$

(%i2) eq4:676y^2+203^(5/2)-333$

(%i3) Y2:y = -(sqrt(28*3^(7/2)+4077)-78)/26$

(%i4) subst(Y2,eq4);

(%o4) (78-sqrt(283^(7/2)+4077))^2+203^(5/2)-333
(%i5) expand(sqrtdenest(%));

(%o5) 0

Calling sqdenest within ebaksubst works for this simple case, and algsys
returns the correct solution.

(defun sqdenest (e) (meval `(($sqrtdenest) ,e)))
(defun ebaksubst (solnl lhsl) (mapcar #'(lambda (q) (cadr (ratf (sqdenest (what-the-$ev (pdis q) (cons '(mlist) solnl) '$radcan))))) lhsl))

---

• [bugs:#3208] https://sourceforge.net/p/maxima/bugs/3208/ Can't solve
  this simple equation..*

Status: open
Group: None
Labels: algsys solve
Created: Sun Sep 04, 2016 06:12 AM UTC by nightwing
Last Updated: Fri Sep 23, 2016 12:23 PM UTC
Owner: nobody

solve([x^2+y^2=1, (x-2)^2+(y-3)^2=9],[x,y]);
then it returns []
why?

---

Sent from sourceforge.net because you indicated interest in
https://sourceforge.net/p/maxima/bugs/3208/

To unsubscribe from further messages, please visit
https://sourceforge.net/auth/subscriptions/

[David Billinghurst]

With some hacks in place to resimplify the resultant and apply sqrtdenest in ebaksubst, I find another problem.

The first reduced test case works

eq11:676*y^2-20*3^(5/2)-333$
eq21:676*y^2-4056*y+28*3^(7/2)+2007$
algsys([eq11,eq21],y);
[[y = sqrt(20*3^(5/2)+333)/26]]

but the second reduced test case still returns []

eq12:676*y^2+20*3^(5/2)-333$
eq22:676*y^2-4056*y-28*3^(7/2)+2007$
algsys([eq12,eq22],[y]);
[]

Some tracing shows that algsys is attempting

(%i64) Y12:solve(eq12,y);
(%o64) [y = -sqrt(333-20*3^(5/2))/26,y = sqrt(333-20*3^(5/2))/26]

(%i65) ev(eq22,Y12[1],radcan);
(%o65) 468*sqrt(20*sqrt(3)-37)*%i-104*3^(5/2)+2340

(%i66) float(%);
(%o66) -2.273736754432321e-13

This looks wrong as

(%i82) float(Y12[1]);
(%o82) y = -0.1772190444071826

(%i83) float(eq22),%;
(%o83) 1437.600888231064

(%i84) float(Y12[2]);
(%o84) y = 0.1772190444071826

(%i85) float(eq22),%;
(%o85) -9.094947017729282e-13

Something wrong elsewhere in maxima? Reported as [bugs:#3223]

[David Billinghurst]

Fixed by commit [2986fc]

(%i1) solve([x^2+y^2=1, (x-2)^2+(y-3)^2=9],[x,y]);

               5/2                   3/2
              3    - 10        (- 2 3   ) - 15
(%o1) [[x = - ---------, y = - ---------------],
                 26                  26
                                                  5/2                3/2
                                                 3    + 10        2 3    - 15
                                            [x = ---------, y = - -----------]]
                                                    26                26

[nightwing]

Wow! Where is the bug?

发自我的小米手机

在 David Billinghurst <billingd@users.sf.net>，2016年10月2日 下午2:45写道：

• status: open --> closed
• Comment:

Fixed by commit [2986fc]

(%i1) solve([x^2+y^2=1, (x-2)^2+(y-3)^2=9],[x,y]);

               5/2                   3/2
              3    - 10        (- 2 3   ) - 15
(%o1) [[x = - ---------, y = - ---------------],
                 26                  26
                                                  5/2                3/2
                                                 3    + 10        2 3    - 15
                                            [x = ---------, y = - -----------]]
                                                    26                26

---

[bugs:#3208] Can't solve this simple equation..

Status: closed
Group: None
Labels: algsys solve
Created: Sun Sep 04, 2016 06:12 AM UTC by nightwing
Last Updated: Sun Sep 25, 2016 02:28 PM UTC
Owner: nobody

solve([x^2+y^2=1, (x-2)^2+(y-3)^2=9],[x,y]);
then it returns []
why?

---

Sent from sourceforge.net because you indicated interest in https://sourceforge.net/p/maxima/bugs/3208/

To unsubscribe from further messages, please visit https://sourceforge.net/auth/subscriptions/

<link href="https://sourceforge.net/p/maxima/bugs/3208/"> <meta content="View">

<meta content="View">