Zero solution does not result in zero

[8d15854a-f726-4f6b-88e7-82ec1970fbba]

In the following case, the solve function of sage 5.8 will yield a list of incorrect solutions, whereas the correct solution is missing. The computation involves complex numbers, symbolic expressions and trigonometric functions.

----------------------------------------------------------------------
| Sage Version 5.8, Release Date: 2013-03-15                         |
| Type "notebook()" for the browser-based notebook interface.        |
| Type "help()" for help.                                            |
----------------------------------------------------------------------
sage: phi = var('phi', domain='real')
sage: M = Matrix([
....: [-sqrt(-I)/(exp(I*phi)-1), sqrt(-I)/(exp(I*phi)-1)],
....: [-1/(exp(I*phi)-1)+exp(-I*phi), 1/(exp(I*phi)-1)]])
sage: z = M*vector((-I, 1))
sage: z = z[0]/z[1]
sage: zi = z.imag_part()
sage: s1 = [s.rhs() for s in solve(zi == 0, phi)]
sage: s1
[0, pi - arccos(1/2*sqrt(5) - 1/2), arccos(1/2*sqrt(5) + 1/2)]
sage: N(zi.subs(phi=s1[1]))
0.244669085661435
sage: zi.subs(phi=pi/2)
0

The correct solution was found from an easier form of zi, which is -1/2*sqrt(2)*cos(phi)/(sin(phi) + 1). I first reported this at asksage, and tmonteil confirmed that this is a bug, not my mistake.

Upstream: Reported upstream. No feedback yet.

CC: @kcrisman

Component: symbolics

Stopgaps: wrongAnswerMarker

Issue created by migration from https://trac.sagemath.org/ticket/14628

[kcrisman]

comment:2

One can confirm this graphically, too.

sage: plot(zi,(phi,0,pi),ticks=pi/4,tick_formatter=pi)
sage: plot(-1/2*sqrt(2)*cos(phi)/(sin(phi) + 1),(phi,0,pi),ticks=pi/4,tick_formatter=pi)

Note that

sage: solve(zi,phi,to_poly_solve='force')
[phi == 1/2*pi + 2*pi*z70]

works. So it is the native Maxima solve command that is causing the problem.

[kcrisman]

comment:3

Unfortunately, I can't debug this more now because of the annoying thing that Maxima doesn't allow copying and pasting of more than a certain number of characters, and I can't remember the workaround right now. What I can say is that it is something about our form of the imaginary part of z:

(%i1) display2d:false;
(%o1) false
(%i2) declare(phi,real);
(%o2) done
(%i3) f:(%i+1)*sqrt(-%i)/((%e^(%i*phi)-1)*((%i+1)/(%e^(%i*phi)-1)-%i*%e^(-%i*phi)));

(%o3) sqrt(-%i)*(%i+1)/(((%i+1)/(%e^(%i*phi)-1)-%i*%e^-(%i*phi))
                       *(%e^(%i*phi)-1))
(%i4) solve(imagpart(f),phi);
(%o4) [sin(phi) = -sqrt(-cos(phi)^2+cos(phi)+1),
       sin(phi) = sqrt(-cos(phi)^2+cos(phi)+1)]

on which to_poly_solve also works correctly.

(%i5) load(to_poly_solve);
Loading maxima-grobner $Revision: 1.6 $ $Date: 2009-06-02 07:49:49 $
(%o5) "/Users/.../sage-5.9.rc1/local/share/maxima/5.29.1/share/to_poly_solve/to_poly_solve.mac"
(%i6) %solve(imagpart(f),phi);
(%o6) %union([phi = 2*%pi*%z8+%pi/2])

It is strange that the error is so dramatically wrong.

[kcrisman]

comment:4

Ah, got it. The problem is related to #6862 in that phi wasn't actually real. Unfortunately, it doesn't seem to be quite that simple; somehow the imaginary part in Maxima is different from in Pynac.

sage: phi = var('phi')
sage: z = (I + 1)*sqrt(-I)/((e^(I*phi) - 1)*((I + 1)/(e^(I*phi) - 1) - I*e^(-I*phi)))
sage: z.maxima_methods().imagpart()
-1/2*(((cos(phi) - 1)*((sin(phi) + cos(phi) - 1)/((cos(phi) - 1)^2 + sin(phi)^2) - sin(phi)) - 
<more stuff>
sage: z.imag()
-sqrt(2)*e^imag_part(phi)*cos(-real_part(phi))/((e^(-2*imag_part(phi))*sin(real_part(phi))^2 + 
<a whole lot of stuff with imaginary and real part>
sage: assume(phi,'real')
sage: z.imag().simplify()
sqrt(2)*sin(phi)^2/((sin(phi)^2 + 
<a whole bunch of stuff without imag part and real part, correctly>
+ 2*cos(phi)^2/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1)^2 - 4*cos(phi)/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1)^2 + 2/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1)^2))

But it isn't evident that the last thing and the maxima_methods() thing are the same (unless fully simplified, but that's not what we're talking about here, as see the post mentioned it the description).

[8d15854a-f726-4f6b-88e7-82ec1970fbba]

comment:5

One way to deal with the line length limitation would be splitting the expression into several parts:

x1:    - sqrt(2)*sin(-phi)*sin(phi)/((sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1)*(cos(-phi)^2 + 2*sin(-phi)*sin(phi)/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1) + 2*sin(-phi)*cos(phi)/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1) + 2*sin(phi)*cos(-phi)/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1) + 2/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1)^2 + sin(-phi)^2 - 2*cos(-phi)*cos(phi)/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1) - 2*sin(-phi)/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1) + 2*cos(-phi)/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1) + 2*sin(phi)^2/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1)^2 + 2*cos(phi)^2/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1)^2 - 4*cos(phi)/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1)^2));
x2: x1 + sqrt(2)*cos(-phi)*cos(phi)/((sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1)*(cos(-phi)^2 + 2*sin(-phi)*sin(phi)/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1) + 2*sin(-phi)*cos(phi)/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1) + 2*sin(phi)*cos(-phi)/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1) + 2/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1)^2 + sin(-phi)^2 - 2*cos(-phi)*cos(phi)/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1) - 2*sin(-phi)/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1) + 2*cos(-phi)/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1) + 2*sin(phi)^2/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1)^2 + 2*cos(phi)^2/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1)^2 - 4*cos(phi)/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1)^2));
x3: x2 - sqrt(2)*cos(-phi)/((sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1)*(cos(-phi)^2 + 2*sin(-phi)*sin(phi)/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1) + 2*sin(-phi)*cos(phi)/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1) + 2*sin(phi)*cos(-phi)/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1) + 2/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1)^2 + sin(-phi)^2 - 2*cos(-phi)*cos(phi)/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1) - 2*sin(-phi)/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1) + 2*cos(-phi)/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1) + 2*sin(phi)^2/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1)^2 + 2*cos(phi)^2/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1)^2 - 4*cos(phi)/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1)^2));
x4: x3 - sqrt(2)*sin(phi)^2/((sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1)^2*(cos(-phi)^2 + 2*sin(-phi)*sin(phi)/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1) + 2*sin(-phi)*cos(phi)/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1) + 2*sin(phi)*cos(-phi)/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1) + 2/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1)^2 + sin(-phi)^2 - 2*cos(-phi)*cos(phi)/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1) - 2*sin(-phi)/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1) + 2*cos(-phi)/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1) + 2*sin(phi)^2/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1)^2 + 2*cos(phi)^2/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1)^2 - 4*cos(phi)/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1)^2));
x5: x4 - sqrt(2)*cos(phi)^2/((sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1)^2*(cos(-phi)^2 + 2*sin(-phi)*sin(phi)/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1) + 2*sin(-phi)*cos(phi)/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1) + 2*sin(phi)*cos(-phi)/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1) + 2/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1)^2 + sin(-phi)^2 - 2*cos(-phi)*cos(phi)/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1) - 2*sin(-phi)/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1) + 2*cos(-phi)/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1) + 2*sin(phi)^2/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1)^2 + 2*cos(phi)^2/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1)^2 - 4*cos(phi)/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1)^2));
x6: x5 + 2*sqrt(2)*cos(phi)/((sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1)^2*(cos(-phi)^2 + 2*sin(-phi)*sin(phi)/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1) + 2*sin(-phi)*cos(phi)/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1) + 2*sin(phi)*cos(-phi)/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1) + 2/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1)^2 + sin(-phi)^2 - 2*cos(-phi)*cos(phi)/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1) - 2*sin(-phi)/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1) + 2*cos(-phi)/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1) + 2*sin(phi)^2/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1)^2 + 2*cos(phi)^2/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1)^2 - 4*cos(phi)/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1)^2));
x7: x6 - sqrt(2)/((sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1)^2*(cos(-phi)^2 + 2*sin(-phi)*sin(phi)/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1) + 2*sin(-phi)*cos(phi)/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1) + 2*sin(phi)*cos(-phi)/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1) + 2/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1)^2 + sin(-phi)^2 - 2*cos(-phi)*cos(phi)/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1) - 2*sin(-phi)/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1) + 2*cos(-phi)/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1) + 2*sin(phi)^2/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1)^2 + 2*cos(phi)^2/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1)^2 - 4*cos(phi)/(sin(phi)^2 + cos(phi)^2 - 2*cos(phi) + 1)^2));

This together with the %i1 and %i2 from your example will reproduce the issue:

[phi = 0,phi = %pi-acos(sqrt(5)/2-1/2),phi = acos(sqrt(5)/2+1/2)]

But now I see that you posted a new comment while I was logging in… Reading that, it doesn't seem like a final solution, so I'll post mine anyway, in case it might be useful at some point.

[kcrisman]

comment:6

Replying to @gagern:

One way to deal with the line length limitation would be splitting the expression into several parts:

True, though that wasn't the workaround I was trying to remember.

But now I see that you posted a new comment while I was logging in… Reading that, it doesn't seem like a final solution, so I'll post mine anyway, in case it might be useful at some point.

Absolutely. Hmm, that is actually very useful to know. I wonder why there is that huge difference? Note the Maxima warning...

solve: using arc-trig functions to get a solution.
Some solutions will be lost.

Yeah, I'll say! This is still right, as we expect.

(%i16) %solve(x7,phi);

(%o16) %union([phi = 2*%pi*%z21+%pi/2])

This is something I can report to the Maxima bug tracker, and I'll do so now.

[kcrisman]

Upstream: Reported upstream. No feedback yet.

[kcrisman]

comment:7

This is now Maxima bug 2592.

[ea1d0bf8-c27a-4548-8cb7-de0b1d02441a]

Stopgaps: wrongAnswerMarker