How to solve this solvable 8th-degree algebraic equation by radicals?

Solve the following equation in radicals.

$$x^8-8x^7+8x^6+40x^5-14x^4-232x^3+488x^2-568x+1=0$$

I use Magma to verify that its Galois group is a solvable group.

R := RationalField(); 
R < x > := PolynomialRing(R); 
f := x^8-8*x^7+8*x^6+40*x^5-14*x^4-232*x^3+488*x^2-568*x+1; 
G := GaloisGroup(f); 
print G;
GroupName(G: TeX:=true);
IsSolvable(G);

The output of Magma(Online) is:

Permutation group G acting on a set of cardinality 8
Order = 16 = 2^4
    (2, 4)(6, 8)
    (1, 2, 3, 4)(5, 6, 7, 8)
    (1, 5)(2, 8)(3, 7)(4, 6)
C_2\times D_4
true

I also tried to calculate with PARI/GP(64-bit)v_2.13.3+GAP(64-bit)v_4.11.1, but failed.

gap> LoadPackage("radiroot");
true
gap> x := Indeterminate(Rationals,"x");;
gap> g := UnivariatePolynomial( Rationals, [1,-8,8,40,-14,-232,488,-568,1]);
x^8-8x^7+8x^6+40x^5-14x^4-232x^3+488x^2-568x+1
gap>  RootsOfPolynomialAsRadicals(g, "latex");
"/tmp/tmp.sfoZ6C/Nst.tex"
Error，AL_EXECUTABLE，the executable for PARI/GP，has to be set at /proc/
cygdrive/C/gap-4.11.1/pkg/aInuth-3.1.2/gap/kantin.gi : 205 called from

Well, in honor of an old cartoon I'll say a miracle occurs. But can we get behind the curtain to see how the special effects are made?

If you take the square root of, let us say, $2358$ by the standard "long division" method, you get $48$ with a remainder of $54$, which may be interpreted as the equation

$2358=48^2+54.$

We can adapt this method to determining the square root of a polynomial, and for the one given in this problem we end with this:

$x^8-8x^7+8x^6+40x^5-14x^4-232x^3+488x^2-568x+1=(x^4-4x^3-4x^2+4x+1)^2+(-192x^3+480x^2-576x)$

If the remainder were a constant times a square then we would be able to render our octic polynomial in the form $a^2-b^2=(a+b)(a-b)$ or perhaps $a^2+b^2=(a+bi)(a-bi)$, getting a pair of quartic factors which would then be solvable by radicals in the usual way. Sadly, we can't do that because the remainder is a cubic polynomial. Nonetheless, the fact that the coefficients of this remainder have a common factor makes one go "hmmm...". What if there were a way to modify the remainder so that it has an even degree and could be a square quantity (or next best, a constant times one)?

I started by noting that the square root of $2358$ as determined by the standard method comes out as $48$ with a remainder of $54$. But did I really have to render the "quotient" as $48$? If I allow a negative remainder in the final stage maybe I could render the root as $49$ instead, in which case the remainder is indeed negative and we get an expression equally valid as the first one I quoted:

$2358=48^2+54$ but also

$2358=49^2-43.$

We might even say that the second form is superior because, with the absolutely smaller remainder, it renders the rounded value of $\sqrt{2358}$ (correctly) as $49$ instead of $48$.

Now what can we do with our polynomial square root? Let us say that, just as we rendered the last digit of the root as $9$ instead of $8$ when we extracted the square root of $2358$, we leave the constant term in our quartic expression as something other than $1$. We get

$x^8-8x^7+8x^6+40x^5-14x^4-232x^3+488x^2-568x+1=(x^4-4x^3-4x^2+4x+h)^2+[(-2h+2)x^4+(8h-200)x^3+(8h+472)x^2+(-8h-568)x+(-h^2+1)]$

Can this remainder be a squared quantity, perhaps multiplied by a constant, for some value of $h$, presumably rational?

A necessary condition for this to occur in the quartic expression $ax^4+bx^3+cx^2+dx+e$ is $a/e=(b/d)^2$. Here we require

$\dfrac{-2h+2}{-h^2+1}=\left(\dfrac{8h-200}{-8h-568}\right)^2$

$\dfrac{2}{h+1}=\left(\dfrac{h-25}{h+71}\right)^2$

We turn this to a cubic polynomial equation for $h$, seek rational roots and discover $h=49$. We again go "hmmm...", for not only did we hit on a rational root but we incremented $h$ from its earlier value ($1$) by half the common factor of $96$ we saw in the earlier remainder.

We insert $h=49$ and obtain

$x^8-8x^7+8x^6+40x^5-14x^4-232x^3+488x^2-568x+1=(x^4-4x^3-4x^2+4x+49)^2-96[x^4-2x^3-9x^2+10x+25]$

If the bracketed quantity were to be a square, it would be $(x^2-x-5)^2$ to match the degree 4, degree 3, degree 1 and degree 0 terms (which our equation for $h$ was designed to do). But do we get the proper degree 2 term? In fact:

$(x^2-x-5)^2=x^4-2x^3-9x^2+10x+25$

and we have hit on our squared remainder!

So now we just factor the octic polynomial as a difference of squares whose roots contain $\sqrt{96}$ or equivalently $\sqrt{6}$:

$x^8-8x^7+8x^6+40x^5-14x^4-232x^3+488x^2-568x+1=[(x^4-4x^3-4x^2+4x+49)+4\sqrt6(x^2-x-5)][(x^4-4x^3-4x^2+4x+49)-4\sqrt6(x^2-x-5)]$

and we then solve each quartic factor by the usual method.

The roots, with all radicals defined as nonnegative real numbers, are

$1+\sqrt2+\sqrt3+\sqrt[4]{3}$

$1-\sqrt2+\sqrt3+\sqrt[4]{3}$

$1+\sqrt2+\sqrt3-\sqrt[4]{3}$

$1-\sqrt2+\sqrt3-\sqrt[4]{3}$

$1+\sqrt2-\sqrt3\pm i\sqrt[4]{3}$

$1-\sqrt2-\sqrt3\pm i\sqrt[4]{3}$

This set of roots conforms with the $C_2×D_4$ symmetry from the Galois group calculation.

Just to explain why the attempt to calculate in GAP failed. The error message was produced by the Alnuth GAP package on which RadiRoot depends, and Alnuth in its turn requires PARI/GP. If everything is installed, the code in question worked (also note the correct order of the coefficients):

gap> LoadPackage("radiroot");
true
gap> x := Indeterminate(Rationals,"x");;
gap> g := UnivariatePolynomial( Rationals, [1,-568,488,-232,-14,40,8,-8,1]);
x^8-8*x^7+8*x^6+40*x^5-14*x^4-232*x^3+488*x^2-568*x+1
gap> RootsOfPolynomialAsRadicals(g, "latex");
"/var/folders/dt/some_random_path/tempfilename.tex"

It wrote a temporary file and displayed its name. The file contains all commands needed to compile it with LaTeX, and I just show here the crucial part of the output (without modifying the content, so it has some redundant details) which says the following:

An expression by radicals for the roots of the polynomial $$x^{8} - 8x^{7} + 8x^{6} + 40x^{5} - 14x^{4} - 232x^{3} + 488x^{2} - 568x + 1$$ with the $n$-th root of unity $\zeta_n$ and $$\omega_1 = \sqrt[2]{2},$$ $$\omega_2 = \sqrt[2]{3},$$ $$\omega_3 = \sqrt[2]{-\omega_2},$$ $$\omega_4 = \sqrt[2]{\omega_2},$$ is: $$1-\omega_1-\omega_2-\omega_3$$

For other sources of help with GAP, especially with specific GAP packages, please see the GAP tag description.

A Pari computation shows that the discriminant of the degree $8$ number field $F$ generated by a root of the polynomial $f=𝑥^8−8𝑥^7+8𝑥^6+40𝑥^5−14𝑥^4−232𝑥^3+488𝑥^2−568𝑥+1$ is the product of a power of $2$ and a power of $3$. Therefore any quadratic subfield of $F$ is of the form ${\bf Q}(\sqrt{\pm d})$ with $d$ a divisor of $6$. A Pari computation shows that $f$ factors into a product of two degree $4$ irreducible factors over ${\bf Q}(\sqrt{2})$ and over ${\bf Q}(\sqrt{3})$. This implies that we have ${\bf Q}(\sqrt{2},\sqrt{3})\subset F$. A final Pari computation shows that $f$ is a product of four quadratic polynomials in ${\bf Q}(\sqrt{2},\sqrt{3})[X]$. Their roots are equal to $1+\sqrt{2}+{\root 4 \of 3}(1+{\root 4 \of 3})$ and its conjugates.