This is a further answer to the question, although fidbc and Emmanuel Charpentier did already said all.

But the given polynomial is a special one, it is a **reciprocal polynomial in $q^2$**, so we may go some steps further. (But unfortunately i could not find an explicit form either.)

So it is natural to make the substitution
$$ t = q^2+\frac 1 {q^2}=q^2+r^2\ , r =\frac 1q\ .$$
I am doing this substitutions my way...

```
R.<q,r, t> = PolynomialRing(QQ)
W = q^32 - q^30 + 3*q^28 - 3*q^26 + 6*q^24 - 6*q^22 + 9*q^20 - 9*q^18 + 12*q^16 - 9*q^14 + 9*q^12 - 6*q^10 + 6*q^8 - 3*q^6 + 3*q^4 - q^2 + 1
J = R.ideal( [ W, q^2 + r^2 - t, q*r-1 ] )
K = J.elimination_ideal( [ q,r ] )
P = K.gens()[0]
print P
```

and the polynomial in $t$ is after substitution (or rather elimination):

```
sage: print P
t^8 - t^7 - 5*t^6 + 4*t^5 + 8*t^4 - 5*t^3 - 4*t^2 + t + 2
```

(*The further lines are misleading unless we are searching for exact solutions...*)

The roots of this polynomial are:

```
sage: var( 't' );
sage: P = t^8 - t^7 - 5*t^6 + 4*t^5 + 8*t^4 - 5*t^3 - 4*t^2 + t + 2
sage: for root in P.roots( ring=QQbar, multiplicities=False ):
....: print root
....:
1.514032200639550?
1.791264905617066?
-1.480484671919947? - 0.2566820580880113?*I
-1.480484671919947? + 0.2566820580880113?*I
-0.4934011167978715? - 0.4011241207989847?*I
-0.4934011167978715? + 0.4011241207989847?*I
0.8212372355895106? - 0.3652199571071136?*I
0.8212372355895106? + 0.3652199571071136?*I
```

and the two real $t$-roots correspond to the eight roots of absolute value one of the initial polynomial:

```
sage: R.<q> = PolynomialRing(QQ)
....: W = q^32 - q^30 + 3*q^28 - 3*q^26 + 6*q^24 - 6*q^22 + 9*q^20 - 9*q^18 + 12*q^16 - 9*q^14 + 9*q^12 - 6*q^10 + 6*q^8 - 3*q^6 + 3*q^4 - q^2 + 1
sage: for w in W.roots( ring=QQbar, multiplicities=False ):
....: if abs(w) != 1: continue
....: print "w = %s |w| = %.2f w^2 + w^-2 = %.8f" % ( w, abs(w), w^2+w^-2 )
....:
w = -0.9735585377388802? - 0.2284376798948314?*I |w| = 1.00 w^2 + w^-2 = 1.79126491
w = -0.9735585377388802? + 0.2284376798948314?*I |w| = 1.00 w^2 + w^-2 = 1.79126491
w = -0.9372876026918778? - 0.3485569535099145?*I |w| = 1.00 w^2 + w^-2 = 1.51403220
w = -0.9372876026918778? + 0.3485569535099145?*I |w| = 1.00 w^2 + w^-2 = 1.51403220
w = 0.9372876026918778? - 0.3485569535099145?*I |w| = 1.00 w^2 + w^-2 = 1.51403220
w = 0.9372876026918778? + 0.3485569535099145?*I |w| = 1.00 w^2 + w^-2 = 1.51403220
w = 0.9735585377388802? - 0.2284376798948314?*I |w| = 1.00 w^2 + w^-2 = 1.79126491
w = 0.9735585377388802? + 0.2284376798948314?*I |w| = 1.00 w^2 + w^-2 = 1.79126491
```

So it is enough to understand the roots of
$$ P = t^{8} - t^{7} - 5 t^{6} + 4 t^{5} + 8 t^{4} - 5 t^{3} - 4 t^{2} + t + 2\ . $$
Now - sometimes, but only sometimes - there is an implemented luck that allows to factorize a polynomial over some smaller field. In our case there is no obvious choice, since we cannot find easy subfields of the splitting field of $P$.

```
sage: K.<a> = NumberField( P )
sage: for L, _, _ in K.subfields():
....: print "Subfield of degree %s and discriminant %s" % ( L.degree(), L.disc().factor() )
Subfield of degree 1 and discriminant 1
Subfield of degree 8 and discriminant -1 * 3^9 * 21157
```

Now i would need some more mathematical insight about the initial `W`

, this in the hope that we can still solve by radicals. Sometimes, after adjoining a special element there appear subfields, and one can split over them. For instance, starting with the simple polynomial $Q=x^3-x-1$, there are no non-trivial subfields of `NumberField(Q)`

. But if we adjoin $\sqrt{-23}$...

```
sage: _.<x> = QQ[]
sage: Q = x^3 -x -1
sage: K.<b> = NumberField( Q )
sage: for L, _, _ in K.subfields(): print L.degree(), L.disc()
1 1
3 -23
K.<b,c> = NumberField( Q, x^2+23 )
sage: for L, _, _ in K.subfields(): print L.degree(), L.disc()
1 1
2 -23
3 -23
3 -23
3 -23
6 -12167
```

there is a field of degree $2$, that is helping us to "solve the cubic". Of course, things would be clear, if sage could deliver in real time the Galois closure, for instance...

```
R.<t> = PolynomialRing( QQ )
P = t^8 - t^7 - 5*t^6 + 4*t^5 + 8*t^4 - 5*t^3 - 4*t^2 + t + 2
K.<a> = NumberField( P )
L.<b> = K.galois_closure()
```

(Had to abort the computation after a while.) So is there geometric insight for the initial polynomial `W`

. (E.g. do we count points of some variety...?!)

also for approx solutions, try with:

`solve(f, x, to_poly_solve=True)`