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]
