#3223 ev(...,radcan) issue

[David Billinghurst]

This is a subcase reduced from [bugs:#3208], which is a problem with algsys

We are looking for solutions of eq12 and eq22. We can solve eq12 for y, giving two solutions in list Y12. Numerically, we see that y = (+/-) 0.17721, and that both roots solve eq12 but only the positive root satisfies eq22. (With some messing around I can do this symbolically, too).

(%i1) eq12:676*y^2+20*3^(5/2)-333$
(%i2) eq22:676*y^2-4056*y-28*3^(7/2)+2007$

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

(%i4) yn1:float(Y12[1]);
(%o4) y = -0.1772190444071826

(%i5) float(eq12),yn1;
(%o5) 0.0

(%i6) float(eq22),yn1;
(%o6) 1437.600888231064

(%i7) yn2:float(Y12[2]);
(%o7) y = 0.1772190444071826

(%i8) float(eq12),yn2;
(%o8) 0.0

(%i9) float(eq22),yn2;
(%o9) -9.094947017729282e-13

The problem is that algsys substitutes the roots Y12 into eq22 using ev(eq22,y=..,radcan) and the result here is wrong.

(%i17) ev(eq22,Y12[1],radcan);
(%o17) 468*sqrt(20*sqrt(3)-37)*%i-104*3^(5/2)+2340
(%i18) float(%);
(%o18) -2.273736754432321e-13

(%i19)  ev(eq22,Y12[2],radcan);
(%o19) (-468*sqrt(20*sqrt(3)-37)*%i)-104*3^(5/2)+2340
(%i20) float(%);
(%o20) 1437.600888231063

[David Billinghurst]

I think this is a bug in radcan, Continuing from above:

(%i24) eqn:subst(Y12[1],eq22);
(%o24) (-104*3^(5/2))+156*sqrt(333-20*3^(5/2))+2340

(%i25) float(eqn);
(%o25) 1437.600888231064

(%i26) radcan(eqn);
(%o26) 468*sqrt(20*sqrt(3)-37)*%i-104*3^(5/2)+2340

(%i27) float(%);
(%o27) -2.273736754432321e-13

and looking at the subexpressions

(%i33) eqn2:sqrt(333-20*3^(5/2));
(%o33) sqrt(333-20*3^(5/2))

(%i34) float(%);
(%o34) 4.607695154586747

(%i35) radcan(eqn2);
(%o35) 3*sqrt(20*sqrt(3)-37)*%i

(%i36) float(%);
(%o36) -4.60769515458674

and even simpler

(%i45) eq3:sqrt(6-3^(3/2));
(%o45) sqrt(6-3^(3/2))

(%i46) float(%);
(%o46) 0.8965754721680539

(%i47) radcan(eq3);
(%o47) sqrt(sqrt(3)-2)*sqrt(3)*%i

(%i48) float(%);
(%o48) -0.8965754721680536

[Stavros Macrakis]

This isn't related to nested roots in particular. It has to do with radcan treating sqrt(EX) and -sqrt(EX) as, in some sense, interchangeable. This is working as designed, so it looks like radcan should not be called in this case.

[David Billinghurst]

This is fixed by commits [2986fc] and [640ca7], which avoids the calls to radcan.

[Stavros Macrakis]

By the way, ev(...,eq,radcan) is tempting fate because of the idiosyncracies of ev and the possible side-effects of evaluating (...), and should certainly never be used in code.

radcan(subst(eq,...)) is much better.