Sympy offers us a solution, but sage is not (yet) ready to parse it totally. But with a little Sympy's doc reading and a little manual unscrewing, one can get it:

```
sage: [u.args[0].args[1]._sage_() for u in c_p.solve(theta, algorithm="sympy")]
[2*pi*n + arctan((2*pi - 5)/(pi*sqrt(20/pi - 25/pi^2 + 12))),
pi + 2*pi*n - arctan((2*pi - 5)/(pi*sqrt(20/pi - 25/pi^2 + 12))),
-pi + 2*pi*n + arctan(1/2*(2*pi + 5)/(pi*sqrt(-5/pi - 25/4/pi^2 + 3))),
2*pi + 2*pi*n - arctan(1/2*(2*pi + 5)/(pi*sqrt(-5/pi - 25/4/pi^2 + 3)))]
```

A bit of fooling around with `latex`

and string splicing allows us to get the $\LaTeX$ form:

$$2 \pi n + \arctan\left(\frac{2 \pi - 5}{\pi \sqrt{\frac{20}{\pi} - \frac{25}{\pi^{2}} + 12}}\right)$$
$$\pi + 2 \pi n - \arctan\left(\frac{2 \pi - 5}{\pi \sqrt{\frac{20}{\pi} - \frac{25}{\pi^{2}} + 12}}\right)$$
$$-\pi + 2 \pi n + \arctan\left(\frac{2 \pi + 5}{2 \pi \sqrt{-\frac{5}{\pi} - \frac{25}{4 \pi^{2}} + 3}}\right)$$
$$2 \pi + 2 \pi n - \arctan\left(\frac{2 \pi + 5}{2 \pi \sqrt{-\frac{5}{\pi} - \frac{25}{4 \pi^{2}} + 3}}\right)$$

Note: `Maxima`

could also get it (via `to_poly_solve`

), but gives pompous pedantic nonsense involving logs of complex expressions, probably due to the `domain:complex`

default...

Note 2: Substituting these values in `c_p`

gets rid of the trig functions easily (via `trig_reduce()`

after declaring `n`

integer...). However, the resulting expressions involve sums ad products of scalar and square roots values ; proving them to be 0 is another possibly nontrivial problem... It can be "checked" numerically, and formally proved for the last two:

```
sage: L2=[u.args[0].args[1]._sage_() for u in c_p.solve(theta, algorithm="sympy")]
sage: [c_p.subs(theta==u).trig_reduce().n() for u in L2]
[2.22044604925031e-16,
2.22044604925031e-16,
2.22044604925031e-16,
2.22044604925031e-16]
sage: [bool(c_p.subs(theta==u).trig_reduce()==0) for u in L2]
[False, False, True, True]
```

**EDIT :** A few notes for the terminally curious:

`fricas`

and `giac`

give the same results as `maxima`

(i. e. Sage).

I have been unable to obtain no-nonsensical results from `to_poly_solve`

.

One can also screw around `Mathematica`

's answers to get another form of the results (closer to what Sage/Maxima gives) :

```
sage: var("k", domain="integer")
k
sage: LM=[u[1][2][1].sage(locals={"C":lambda x:k}) for u in mathematica.Solve(c_p
....: ==0, theta)]
sage: LM
[pi + 2*pi*k + arcsin(1/4*(2*pi + 5)/pi),
2*pi*k - arcsin(1/4*(2*pi + 5)/pi),
pi + 2*pi*k - arcsin(1/4*(2*pi - 5)/pi),
2*pi*k + arcsin(1/4*(2*pi - 5)/pi)]
```

which is:

$$\pi + 2 \pi k + \arcsin\left(\frac{2 \pi + 5}{4 \pi}\right)$$
$$2 \pi k - \arcsin\left(\frac{2 \pi + 5}{4 \pi}\right)$$
$$\pi + 2 \pi k - \arcsin\left(\frac{2 \pi - 5}{4 \pi}\right)$$
$$2 \pi k + \arcsin\left(\frac{2 \pi - 5}{4 \pi}\right)$$

which turns out to be somewhat easier to check/prove:

```
sage: [c_p.subs(theta==u).trig_reduce().expand() for u in LM]
[0, 0, 0, 0]
```