Solveset returns empty set for cos(2x) = 1/2

[tarjeiba]

Running in https://live.sympy.org/, where the version seems to be 1.5.1 you get the following:

>>> solveset(cos(2*x) - 1/2, x, Interval(0, 2*pi))
{π6,5π6,7π6,11π6}

However, running 1.7.1, you get:

>>> solveset(cos(2*x) - 1/2, x, Interval(0, 2*pi))
∅

Any idea where to begin looking for a solution?

[oscarbenjamin]

Master gives the same as 1.7.1:

In [1]: solveset(cos(2*x) - 1/2, x, Interval(0, 2*pi))
Out[1]: ∅

Any idea where to begin looking for a solution?

I would use git bisect to identify the commit that broke this.

[tarjeiba]

You wouldn't happen to know what commit live.sympy.com is running? Checking out the commit tagged sympy-1.5 Head: 7f92c96417 Merge pull request #18015 from oscarbenjamin/bump_version_15 still returns the empty set.

Edit: I am testing via bin/isympy in the local repo folder, sympy.__version__ is 1.5. However, should I uninstall my global sympy installation to ensure that I'm running the code I want?

[oscarbenjamin]

I see now there are a couple of compounding issues. SymPy live is running the commit with tag sympy-1.5.1 I think. However it runs that under Python 2.7. In Python 2.7 the expression 1/2 gives the int 0 whereas in Python 3.x 1/2 gives float 0.5. However in sympy live 1/2 is translated to Rational(1, 2) which won't happen if you run locally under any versions of Python or SymPy.

Actually the equation is correctly solved by solveset on current master if you use Rational(1, 2) or S(1)/2:

In [1]: solveset(cos(2*x) - S(1)/2, x, Interval(0, 2*pi))
Out[1]: 
⎧π  5⋅π  7⋅π  11⋅π⎫
⎨─, ───, ───, ────⎬
⎩6   6    6    6  ⎭

[tarjeiba]

Ah. I see. That explains it. Thank you for looking into it.

I assume the casting of 1/2 into Rational(1, 2) is intentional? Do you know why it's set only for sympy live?

[oscarbenjamin]

It isn't possible for sympy to control how 1/2 is evaluated in Python. I guess that because sympy live is made by sympy someone has added something to parse the input expressions so that they become rational.

The issue here then is why solveset does not work for a float:

In [2]: solveset(cos(2*x) - 0.5, x, Interval(0, 2*pi))
Out[2]: ∅

[hyadav2k]

sympy/sympy/solvers/solveset.py

Lines 626 to 631 in 9e8f62e

	if isinstance(solns, FiniteSet):
	if any(isinstance(s, RootOf) for s in solns):
	raise _SolveTrig1Error("polynomial results in RootOf object")
	# revert the change of variable
	cov = cov.subs(x, symbol/mu)
	result = Union(*[inverter(cov, s, symbol)[1] for s in solns])

The issue lies here, for solveset(cos(2*x) - 0.5, x, Interval(0, 2*pi)) and solveset(2*cos(2*x) - 1, x, Interval(0, 2*pi)), the result value varies, specifically for solveset(cos(2*x) - 0.5, x, Interval(0, 2*pi)), it comes in complex domain, while for solveset(2*cos(2*x) - 1, x, Interval(0, 2*pi)), it comes correct i.e. the correct solution.

sympy/sympy/solvers/solveset.py

Line 650 in 9e8f62e

return ConditionSet(symbol, cond, Intersection(result, domain))

Ultimately, Intersection(result, domain) gives us EmptySet.

[oscarbenjamin]

The problem also happens when solving over the reals:

In [27]: solveset(cos(2*x) - 0.5, x, Reals)
Out[27]: ∅

This is because solveset generates a solution over the Complexes and then tries to extract the real solutions:

In [28]: solveset(cos(2*x) - 0.5, x, Complexes)
Out[28]: 
⎧-ⅈ⋅(ⅈ⋅(2⋅n⋅π - 1.0471975511966) - 1.11022302462516e-16)         ⎫   ⎧-ⅈ⋅(ⅈ⋅(2⋅n⋅π + 1.0471975511966) - 1.11022302462516e-16)         ⎫
⎨──────────────────────────────────────────────────────── | n ∊ ℤ⎬ ∪ ⎨──────────────────────────────────────────────────────── | n ∊ ℤ⎬
⎩                           2                                    ⎭   ⎩                           2                                    ⎭

However because of the use of Float the expression for the solution has a small imaginary part coming from a rounding error:

In [29]: solveset(cos(2*x) - 0.5, x, Complexes) & Reals
Out[29]: ∅

In [31]: solveset(cos(2*x) - 0.5, x, Complexes).args[0].args[0].args[1]
Out[31]: 
-ⅈ⋅(ⅈ⋅(2⋅n⋅π - 1.0471975511966) - 1.11022302462516e-16) 
────────────────────────────────────────────────────────
                           2                            

In [32]: solveset(cos(2*x) - 0.5, x, Complexes).args[0].args[0].args[1].expand()
Out[32]: n⋅π - 0.523598775598299 + 5.55111512312578e-17⋅ⅈ

There are several ways that we could address this:

1. Make intersection cleverer about Float. There are other issues arising from this (e.g. RootSet and .evalf for sets #20356).
2. Add specific trig handlers for the real case in solveset.
3. Make solveset return a better form for the solution. It isn't clear to me why solveset has generated an expression with -I*I in it. (Maybe it translates cos to exp internally).

I don't know which is the best approach to take in solveset here.

@gschintgen you have more experience of the solveset code than me. Do you have an opinion on what would be the best approach?

[gschintgen]

A few thoughts:

1. The proper solution will be to fix/implement solveset: current state of _solve_trig and _invert #17334. (I have rough proof of concept code but not enough time to follow through. Stumbling blocks I encountered: bugs in sets code (most of them should be fixed by now), no proper way to handle conditions as e.g. PiecewiseSet, see Improve _solve_trig1 (handle rational & symbolic coefficients) #19507)
2. In the example given here, the main trig solver _solve_trig1 is used. It's based on .rewrite(exp) as conjectured by @oscarbenjamin. It will probably fail quite systematically for Floats: going back to the initial variable requires solving exp(I*x) = some-complex-float-lying-on-unit-circle and the chances of its solutions landing exactly on the real line are slim. And those erroneously non-real solutions will then be discarded. The fix would require to detect these small residual imaginary parts and chop them. But what threshold to use? This type of issue must have arisen in other parts of the codebase? How is it handled there?
3. The -I*I should be rather easy to fix. Some kind of shallow expand should do the trick. Either here

   sympy/sympy/solvers/solveset.py

   Line 654 in 95fbdbc

   result = Union(*[inverter(cov, s, symbol)[1] for s in solns])

   to fix up the return value or, preferably, here:

   sympy/sympy/solvers/solveset.py

   Line 343 in 95fbdbc

   exp_invs = Union(*[imageset(Lambda(n, I*(2*n*pi + arg(g_y)) +

[oscarbenjamin]

Thanks @gschintgen that's very useful.

1. The proper solution will be to fix/implement solveset: current state of _solve_trig and _invert #17334.

Yes, I guess that's sort of what I meant by point 2 "Add specific trig handlers" rather than solving with exp over the complex numbers.

2. The fix would require to detect these small residual imaginary parts and chop them. But what threshold to use? This type of issue must have arisen in other parts of the codebase? How is it handled there?

The only other places that handle trig equations are nonlinsolve and solve. It seems that nonlinsolve returns the same as solveset (I assume it uses solveset internally) whereas solve gives:

In [4]: solve(cos(2*x) - 0.5)
Out[4]: [0.523598775598299, 2.61799387799149]

I think that solve here is just using a more direct solver for trig functions. In solve hyperbolic trig functions are rewritten as exp but normal trig functions are handled directly.

The solve function has a keyword argument "rational" which defaults to True and results in all Floats being converted to Rational on input. This is a general workaround for the difficulties of handling floats.

Other similar cases are Poly which generates both real and non-real roots but has a whole load of code for carefully isolating and extracting real roots as separate from non-real roots. The full formulas from roots when numerically evaluated will give small imaginary parts though:

In [38]: p = 4*x**3 - 3*x - S.Half

In [39]: Poly(p).nroots()
Out[39]: [-0.766044443118978, -0.17364817766693, 0.939692620785908]

In [40]: real_roots(p)
Out[40]: 
⎡       ⎛   3             ⎞         ⎛   3             ⎞         ⎛   3             ⎞⎤
⎣CRootOf⎝8⋅x  - 6⋅x - 1, 0⎠, CRootOf⎝8⋅x  - 6⋅x - 1, 1⎠, CRootOf⎝8⋅x  - 6⋅x - 1, 2⎠⎦

In [42]: [r.n() for r in real_roots(p)]
Out[42]: [-0.766044443118978, -0.17364817766693, 0.939692620785908]

In [41]: [r.n() for r in roots(p)]
Out[41]: [0.939692620785908 - 0.e-22⋅ⅈ, -0.766044443118978 - 0.e-23⋅ⅈ, -0.17364817766693 + 0.e-21⋅ⅈ]

The solution then is to use real_roots if you want real roots symbolically and nroots if you want numerical approximations of the roots. Going via the radical form makes everything more difficult which I guess is analogous to the way that solveset rewrites trig functions as exp, works over the complexes and then tries to figure out what is real at the end.

I have rough proof of concept code but not enough time to follow through

If you could push that up as it is with a short explanation then I think that would be very helpful if someone else was looking to take this on. It could be a PR or it could just be an issue with the diff posted in or something.

A proof of concept is still extremely useful as a starting point if it represents a good approach as opined by someone who has studied the existing code and its limitations carefully. If you don't have time to finish it yourself then the most useful contribution (if you had any time) would be being able to review a PR from someone else if they did manage to finish it off.