Tschirnhaus' transformation

Calculating power sums with Newton's relations

We first need to implement a function for the sum s(n) of the n-th powers of all of the roots of a quintic equation a0 x^5 + a1 x^4 + a2 x^3 + a3 x^2 + a4 x + a5 = 0. In general Newton's relations are given recursively by s(n) = - (a1 s(n-1) + a2 s(n-2) + ... + an-1 s(1) + an n)/a0.
We implement this in the function:

Psi[q_, x_, n_Integer] := Psi[q, x, n] =
-(n Coefficient[q, x, 5-n] + 
     Sum[Psi[q, x, n-j] Coefficient[q, x, 5-j],
         {j, n-1}]) / Coefficient[q, x, 5]

It is interesting that it is not necessary to solve the quintic to find the roots explicitly; the coefficients of the quintic have all the information we need.
Let us give an example. Here is a quintic equation:

r = (2. - 2.x + 4.x^2 + x^3 + 5.x^4 - 7x^5 == 0)

Here we calculate the first ten power sums with Psi:

Table[Psi[r[[1]], x, i], {i, 10}]

For comparison, we find the roots numerically and calculate the same power sums by brute force:

roots = x /. NSolve[r, x]

Table[Plus @@ (roots^i), {i, 10}] // Chop

Transformation to a principal quintic

We will show how to transform a quintic a0 x^5 + a1 x^4 + a2 x^3 + a3 x^2+ a4 x + a5 = 0 to one of the form y^5 + b3 y^2 + b4 y + b5 = 0. Newton's relations for the yj in terms of the bj is a linear system in the bj; solving for the bj expresses them in terms of the power sums 
s(n, yj). We will show that these power sums can be expressed in terms of the aj, so that the bj will be expressed in terms of the aj.
For a quintic to have no quartic or cubic term it is necessary that the sums of the roots and the sums of the squares of the roots vanish: s(1, yj) = 0, s(2, yj) = 0 (*). 
Assume that the roots yj of the new quintic are related to the roots xj of the original quintic by yj = xj^2 + alpha xj + beta. Substituting this into (*) yields two equations for alpha and beta which can be multiplied out, simplified by using Newton's relations for the power sums in the xj, and finally solved, so that alpha and beta can be expressed using radicals in terms of the coefficients aj. 
Again by substitution into (*), we can calculate  s(3, yj), s(4, yj), s(5, yj) in terms of alpha and beta and the xj; by the previous solution for alpha and beta and again by using Newton's relations for the power sums in the xj we can ultimately express these power sums in terms of the aj. 
We implement this in the following definition. The result is a list with two elements; the first element is the transformation of the roots as a pure function and the second element is the new quintic.

PrincipalTransform[p_ == 0, x_, y_]:=
Module[{alpha, beta, xi},
{alpha, beta} = {alpha, beta} /. Last[Solve[
  {5 (xi^2+ alpha xi + beta) == 0, 
   Expand[5 (xi^2 + alpha xi + beta)^2] == 0} /.
                xi^n_.->1/5 Psi[p, x, n],{alpha, beta}]];
  {
   Evaluate[#^2 + alpha # + beta]&, 
   (y^5 - Sum[y^(5-j)/j Collect[(xi^2 + alpha xi + beta)^j + 
   		4 beta^j, xi] /.
    		xi^n_. -> Psi[p, x, n], {j, 3, 5}]) == 0
   }
   ] /; MatchQ[CoefficientList[p, x],
   		{_, _, _, _, _?(# =!= 0 &), _}]

As an example we use the same equation r as above. The principal quintic of r is:

p = PrincipalTransform[r, x, y]

Applying the pure function (the first element of this result) to the roots xi of the old equation gives the roots of the transformed equation yi:

(First[p] /@ roots) // Sort

Solving numerically gives us a check:

NSolve[Last[p], y]

Transformation to a Bring-Jerrard quintic

We can also use a Tschirnhaus transformation to take a quintic y^5 + b3 y^2 + b4 y + b5 = 0 to one of the form z^5 + c4 z + c5 = 0. Now the sum of the third powers of the roots has to vanish as well: s(3, yj) = 0.
We assume that the roots zi of the new Bring-Jerrard quintic are related to the roots yi of the principal quintic by zi = alpha yi^4 + beta yi^3 + gamma yi^2 + delta yi  + epsilon. In a similar way to the previous reduction we can express the coefficients cj in terms of the bj. 
This is implemented in the following definition. The result is a list with two elements; the first element is the transformation of the roots as a pure function and the second element is the new quintic.

BringJerrardTransform[p_ == 0, y_, z_] :=
Module[{alpha, beta, gamma, delta, epsilon, kappa, lambda,
        mu, nu, psi, xi, zeta, a, b, c, g, h},
psi[t_] := Expand[5t]/.xi^n_. -> Psi[p, y, n]/5;
{a, b, c} = Map[Psi[p, y, #]&, {3, 4, 5}];
g=5a xi^3 - 5 b xi^2 - a^2; h = 5 a xi^4 - 5 c xi^2 - a b;
{lambda, mu, nu} = Map[psi, {g^2, 2 g h, h^2}];
kappa = -mu/(2 lambda) + Sqrt[mu^2/(4 lambda^2) - nu/lambda];
delta = Solve[psi[(zeta xi+ kappa g + h)^3] == 0, 
     zeta][[1, 1, 2]];
alpha = 5a; beta = 5a kappa; gamma = -5b kappa-5 c; 
epsilon = -a^2 kappa- a b;
 {
 Evaluate[alpha #^4 + beta #^3 + 
     gamma #^2 + delta # + epsilon]&,
 z^5 - Sum[z^(5-j)/j Collect[psi[(delta xi + 
     kappa g + h)^j],xi],{j,4,5}] == 0
  }
 ]/; MatchQ[CoefficientList[p, y],{_, _, _, 0, 0, _}]

We continue our example:

b = BringJerrardTransform[p[[2]], y, z]

Canonical quintic

A Bring-Jerrard quintic has the form z^5 + c4 z + c5 = 0. It is possible to write such a quintic by a simple rescaling in the form t^5 - t - r = 0, a canonical quintic. (It is not possible to transform a sextic to another that depends on only one parameter.)

CanonicalTransform[z_^5 + e_. z_ + f_ == 0, z_, t_] :=
          {#/(-e)^(1/4)&, t^5 - t + f/(-e)^(5/4) == 0}

Here is our example reduced to a canonical quintic:

c = CanonicalTransform[b[[2]], z, t]

Transforming back

After having solved a canonical quintic via Hermite's method or with hypergeometric functions, or a principal quintic via Klein's approach, it is necessary to reverse the Tschirnhaus transformations to get back to the roots of the original quintic. Since these polynomial transformations are nonlinear inverting them will give mutiple roots; the spurious ones can be  eliminated by seeing which one satisfies the original equation. The following function implements this. The original polynomial is p and s is the list of candidate solutions.

ChooseRoot[p_, s_] := 
         s[[Position[#, Min[#]][[1, 1]]&[Abs[p /. s]]]]

These are the solutions of t^5 -t + 3.24201 = 0:

tau = {-1.356724, 
     -0.317423 - 1.26515 I, 
     -0.317425 + 1.26515 I, 
      0.9957874 - 0.64192 I, 
      0.9957874 + 0.64192 I};

These are the solutions of the Bring-Jerrard quintic:

zeta = 
     Table[Solve[First[c][z] == tau[[i]], z][[1, 1, 2]],
          {i, 5}]

These are the solutions of the principal quintic:

eta = Table[ChooseRoot[p[[2, 1]],
          Solve[First[b][y] == zeta[[i]], y] ][[1, 2]],
              {i, 5}]

Finally, these are the solutions of the original quintic:

Sort[Table[ChooseRoot[r[[1]],
          Solve[First[p][x] == eta[[i]], x] ][[1, 2]],
              {i, 5}]]

This checks with the direct numerical solution of the original quintic:

x /. NSolve[r, x]