# Rational roots of x^5 + a x + b

I'm validating a formula from these article:

If irreducible rational equation $$x^5 + a x + b$$ has solutions, then it can be written as:

\begin{aligned} v_ 1=+\sqrt{d}+\sqrt{d-e\sqrt{d}}\\ v_ 2=-\sqrt{d}-\sqrt{d+e\sqrt{d}}\\ v_ 3=-\sqrt{d}+\sqrt{d+e\sqrt{d}}\\ v_ 4=+\sqrt{d}-\sqrt{d-e\sqrt{d}}\\ \end{aligned}

\begin{aligned} u_ 1=\sqrt[5]{\frac{p^5 v_ 1^2 v_ 3}{d^2}}\\ u_ 2=\sqrt[5]{\frac{p^5 v_ 3^2 v_ 4}{d^2}}\\ u_ 3=\sqrt[5]{\frac{p^5 v_ 2^2 v_ 1}{d^2}}\\ u_ 4=\sqrt[5]{\frac{p^5 v_ 4^2 v_ 2}{d^2}}\\ \end{aligned}

$$x_i = w^i u_1 + w^{2i} u_2+ w^{3i} u_3+ w^{4i} u_4, i=\{0,1,2,3,4\},$$

where $$w=e^{2\pi i / 5}$$.

I translated it into the following codes:

getPseudoRoot[a_, b_, paras_] := Module[
{u, v},
Subscript[v, 1] = + Sqrt[d] + Sqrt[d - e Sqrt[d]];
Subscript[v, 2] = - Sqrt[d] - Sqrt[d + e Sqrt[d]];
Subscript[v, 3] = - Sqrt[d] + Sqrt[d + e Sqrt[d]];
Subscript[v, 4] = + Sqrt[d] - Sqrt[d - e Sqrt[d]];
Subscript[u, 1] = p^5 Subscript[v, 1]^2 * Subscript[v, 3] / d^2;
Subscript[u, 2] = p^5 Subscript[v, 3]^2 * Subscript[v, 4] / d^2;
Subscript[u, 3] = p^5 Subscript[v, 2]^2 * Subscript[v, 1] / d^2;
Subscript[u, 4] = p^5 Subscript[v, 4]^2 * Subscript[v, 2] / d^2;
Table[Subscript[u, i], {i, 1, 4}] /. paras
];
getRational[a_, b_] := Solve[
{
a == 5p^4(3 - 4 e q) / (q^2 + 1),
b == -4p^5(11e + 2  q) / (q^2 + 1),
e^2 == 1,
d == q^2 + 1,
q >= 0,
p != 0
},
{p, e, q, d}, Rationals
];
getRoot[a_, b_] := Module[
{paras, pseudo, z},
paras = getRational[a, b];
If[Length@paras == 0, Return[Failure["NoSolution", <||>]]];
pseudo = getPseudoRoot[a, b, First@paras];
Table[Sum[pseudo[[i]]^(1 / 5)Exp[2Pi I / 5]^(i j), {i, 1, 4}], {j, 0, 4}]
];


However, the solution calculated by this formula is completely different from the numerical solution.

SortBy[getRoot[11, 44] // N, Im]
SortBy[x /. NSolve[x^5 + 11 x + 44 == 0, x], Im]


What's wrong with my code that I can't reproduce the formula?

### Update 20220517

On further research, I found a new case where \zeta can be an imaginary number:

eq = x^5 + x^4 - 4 x^3 - 3 x^2 + 3 x + 1;
{a, b} = {
-(1 / 5),
{
Root[25937424601 + 492715403125 # + 4562314453125 #^2 + 29876708984375 #^3 + 95367431640625 #^4& , 1, 0],
Root[25937424601 + 492715403125 # + 4562314453125 #^2 + 29876708984375 #^3 + 95367431640625 #^4& , 2, 0],
Root[25937424601 + 492715403125 # + 4562314453125 #^2 + 29876708984375 #^3 + 95367431640625 #^4& , 3, 0],
Root[25937424601 + 492715403125 # + 4562314453125 #^2 + 29876708984375 #^3 + 95367431640625 #^4& , 4, 0]
}
};
Solve[Factor[eq, Extension -> ToRadicals@First[b]^(1/5)] == 0]


Which case none of the following methods give the correct root

x /. NSolve[eq == 0, WorkingPrecision -> 50]
N[Tr@Surd[b, 5] + a, 50]
N[Tr[b^(1 / 5)] + a, 50]
N[Tr[Sign[b] * Abs[b]^(1 / 5)] + a, 50]


How should I get the correct root in this case?

• Aren't you missing the square roots, e.g. Subscript[u, 1] = Sqrt[p^5 + ...? Commented May 7, 2022 at 8:07
• @ChrisDegnen, it was in pseudo[[i]]^(1 / 5), I put it outside because it affects FullSimplify Commented May 7, 2022 at 8:34
• Programming note: The localization of u and v is subverted by the use of Subscript. Best not to use Subscript as a variable to store values. Commented May 7, 2022 at 17:30

You get the array pseudo

(* {-0.510481,-0.143328,-6.43215,1.80596} *)


These are real numbers, but if any of these numbers is negative, you get a complex value, when exponentiated to one-fifth.

Change last line

Table[Sum[Sign[pseudo[[i]]]*Abs[pseudo[[i]]]^(1/5) * Exp[2Pi I / 5]^(i j), {i, 1, 4}], {j, 0, 4}]


Now I get the same numerical result as when using NSolve

(* {-0.861241 - 1.9105 I, 1.80012 - 1.44746 I, -1.87775, 1.80012 + 1.44746 I, -0.861241 + 1.9105 I} *)

• Seems Surd[x, 5] also works, I'm looking for more examples to verify. Commented May 7, 2022 at 15:40
• Yes, also possible. Commented May 7, 2022 at 15:52
• It's interesting that Solve[x^5 + 11 x + 44 == 0, x] // ToRadicals is unable to express roots in radicals. Commented May 7, 2022 at 16:01
• @VaclavKotesovec Solve[Factor[x^5 + 11 x + 44, Extension -> First@ getPseudoRoot[11, 44, {p -> -1, e -> 1, q -> 2/11, d -> 125/121}]^(1/5)] == 0] gives the roots in terms of radicals. Commented May 7, 2022 at 18:26