Solutions in radicals

This is an implementation of an algorithm for solving quintics in radicals, published by David Dummit : "Solving Solvable Quintics",  Math. of Comp. 57, 387 (1991).

SolvableQ::usage =
"SolvableQ[poly, x] yields True if poly, a quintic in 
the variable x, is solvable in radicals, and yields False 
otherwise.";

SolveQuintic::usage = 
"RadicalQuinticSolve[lhs == rhs, x] attempts to solve the 
reduced quintic lhs == rhs in the variable x in terms of 
radicals.";

SolveQuintic::rad =
"This quintic is unsolvable in radicals.";

SolvableQ[e_, x_] := QuinticByRadicalsQ[e, x]//First

QuinticByRadicalsQ[expr_, x_] :=
	Module[ {coefs, factor},
		If[ Head[factor=Factor[expr]] === Times &&
		    Length[Complement[Exponent[#, x]&/@(List@@factor),
		    		{0,5}]] > 0
		    ,
		    Return[{True, $Failed}]
		];
		coefs = CoefficientList[expr, x];
		If[ And@@((# === Integer || 
		    # === Rational)&/@(Head/@coefs)) 
		    ,
		    If[ Last[coefs] =!= 1
		        ,
			coefs = coefs/Last[coefs]
		    ];
		    If[ coefs[[-2]] =!= 0
		        ,
		        Return[ QuinticByRadicalsQ[expr/.
		        x->x-coefs[[-2]]/(5 coefs[[-1]]), x] ]
		    ];
		    factor = Factor[resolvent @@ 
		    	Append[Drop[coefs, -2]/Last[coefs], x]];
		    If[ Head[factor] === Times
		        ,
		        If[ MemberQ[Exponent[#, x]& /@ 
		            (List@@factor), 1]
		            ,
			        {True, (Cases[factor, p_. x + r_.][[1]])/.
			    	      {p_. x + r_. :> -r/p} }
			        ,
			        {False, $Failed}
		        ]
		        ,
		        {False, $Failed}
		    ]
		    ,
		    {False, $Failed}
		]
	]/; PolynomialQ[expr, x] && Exponent[expr, x] == 5

resolvent[s_,r_,q_,p_,x_] := 
x^6+8r x^5+(2p q^2-6p^2 r+40r^2-50q s) x^4 +
(-2q^4+21p q^2 r-40p^2 r^2 +160r^3-15p^2 q s-400q r s+
125p s^2) x^3+(p^2 q^4-6p^3 q^2 r-8 q^4 r+9p^4 r^2+
76p q^2 r^2-136p^2 r^3+400 r^4-50p q^3 s+90p^2 q r s-
1400q r^2 s+625q^2 s^2+500p r s^2) x^2+(-2p q^6+19p^2 q^4 r-
51p^3 q^2 r^2+3q^4 r^2+32p^4 r^3+76p q^2 r^3-256p^2 r^4+
512r^5-31p^3 q^3 s-58q^5 s+117p^4 q r s+105p q^3 r s+
260p^2 q r^2 s-2400q r^3 s-108p^5 s^2-325p^2 q^2 s^2+
525p^3 r s^2+2750q^2 r s^2-500p r^2 s^2+625 p q s^3-
3125s^4) x +(q^8-13p q^6 r+p^5 q^2 r^2+65 p^2 q^4 r^2-
4p^6 r^3-128 p^3 q^2 r^3 + 17 q^4 r^3+48 p^4 r^4-16p q^2 r^4-
192 p^2 r^5+256r^6-4p^5 q^3 s- 12p^2 q^5 s+18 p^6 q r s+
12p^3 q^3 r s-124 q^5 r s+196 p^4 q r^2 s +590p q^3 r^2 s-
160 p^2 q r^3 s-1600q r^4 s-27p^7 s^2-150p^4 q^2 s^2 - 
125p q^4 s^2 - 99 p^5 r s^2 - 725 p^2 q^2 r s^2 +
1200 p^3 r^2 s^2 +3250q^2 r^2 s^2-2000p r^3 s^2-
1250p q r s^3+3125p^2 s^4-9375 r s^4);

Discriminant[poly_, var_]:=
	Resultant[poly,D[poly,var],var]/Coefficient[poly,var,5]

SolveQuintic[a_ == b_, x_] :=
	Module[ {temp},
		temp = QuinticByRadicalsQ[a-b, x];
		If[ temp[[1]]
		    ,
		    If[ temp[[2]] === $Failed
		        ,
			    Return[ Solve[ a == b, x] ]
			    ,
			    If[ (# === Integer || # === Rational)& 
			        [Head[temp[[2]]]]
			        ,
			        tosolve[a-b, x, temp[[2]]]
			        ,
			        Message[Solve::rad];
			        Return[$Failed]
			    ]
		    ]
		    ,
		    Message[Solve::rad];
		    Return[$Failed]
		]
	]

tosolve[pol_, x_, r_] :=
	Module[ {coefs},
		coefs = CoefficientList[pol, x];
		If[ Last[coefs] =!= 1
		    ,
		    coefs = coefs/Last[coefs]
		];
		If[ coefs[[-2]] =!= 0
		    ,
		    coefs = CoefficientList[
		    	pol/. x->x-coefs[[-2]]/(5 coefs[[-1]]), x]
		];
		If[ coefs[[3]] == coefs[[4]] == coefs[[5]] == 0
		    ,
		    canonicalQuintic[coefs[[2]], coefs[[1]], x, r,
		    		Sqrt[Discriminant[pol, x]]]
		    ,
		    $Failed
		]
	]

canonicalQuintic[a_, b_, x_, r_, dis_] :=
	Module[ {temp, l0, l1, l2, l3, l4},
		temp = Roots[x^2 + (T1[a,b,r] + T2[a,b,r] dis) x + 
			   (T3[a,b,r] + T4[a,b,r] dis) == 0, x];
		l1 = temp[[1,2]]/.Sqrt[w_] :> I Sqrt[-w]/; Negative[w];
		l4 = temp[[2,2]]/.Sqrt[w_] :> I Sqrt[-w]/; Negative[w];
		temp = Roots[x^2+(T1[a,b,r] - T2[a,b,r] dis) x + 
			   (T3[a,b,r] - T4[a,b,r] dis) == 0, x];
		l2 = temp[[1,2]]/.Sqrt[w_] :> I Sqrt[-w]/; Negative[w];
		l3 = temp[[2,2]]/.Sqrt[w_] :> I Sqrt[-w]/; Negative[w];
		l0 = Expand[-(l1+l2+l3+l4)];
		If[ (Expand[(l1-l4)(l2-l3) - V[a,b,r] dis]/.
			Sqrt[w_] :> Sqrt[w//Expand]) =!= 0,
		    temp = l1;
		    l1 = l4;
		    l4 = temp
		];
		Inner[List[Rule[#1,#2]]&,Table[x,{5}],
			LagrangeResolvent[{l0,l1,l2,l3,l4}],List]
	]

LagrangeResolvent[{l0_,l1_,l2_,l3_,l4_}] :=
	Module[{prim1,prim2,prim3,prim4,r1,r2,r3,r4,p,q},
		prim1 = (-q^2/Sqrt[5] + 2^(1/2)*I*p)/4;
		prim2 = (p^2/Sqrt[5] - 2^(1/2)*I*q)/4;
		prim3 = (p^2/Sqrt[5] + 2^(1/2)*I*q)/4;
		prim4 = (-q^2/Sqrt[5] - 2^(1/2)*I*p)/4;
		r1 = simp[{l0,l1,l2,l3,l4} . 
		          {1,prim1,prim2,prim3,prim4},p,q];	
		r2 = simp[{l0,l3,l1,l4,l2} . 
	              {1,prim1,prim2,prim3,prim4},p,q];
		r3 = simp[{l0,l2,l4,l1,l3} . 
		          {1,prim1,prim2,prim3,prim4},p,q];
		r4 = simp[{l0,l4,l3,l2,l1} . 
		          {1,prim1,prim2,prim3,prim4},p,q];
		r1 = If[Positive[r1], r1^(1/5), -(-r1)^(1/5)];
		r2 = If[Positive[r2], r2^(1/5), -(-r2)^(1/5)];
		r3 = If[Positive[r3], r3^(1/5), -(-r3)^(1/5)];
		r4 = If[Positive[r4], r4^(1/5), -(-r4)^(1/5)];
		prim1 = (-1-5^(1/2)+2^(1/2)*I*(5-5^(1/2))^(1/2))/4;
		prim2 = (-1+5^(1/2)-2^(1/2)*I*(5+5^(1/2))^(1/2))/4;
		prim3 = (-1+5^(1/2)+2^(1/2)*I*(5+5^(1/2))^(1/2))/4;
		prim4 = (-1-5^(1/2)-2^(1/2)*I*(5-5^(1/2))^(1/2))/4;
		{{r1,r2,r3,r4} . {1,1,1,1},	
		 {r1,r2,r3,r4} . {prim4,prim3,prim2,prim1},
		 {r1,r2,r3,r4} . {prim3,prim1,prim4,prim2},
		 {r1,r2,r3,r4} . {prim2,prim4,prim1,prim3},
		 {r1,r2,r3,r4} . {prim1,prim2,prim3,prim4} }/5
	]

simp[expr_, a_,b_] :=
	Collect[expr, {a, b}]/.{
		w_ a^2 :> Expand[w (5 - 5^(1/2))],
		w_ b^2 :> Expand[w (5 + 5^(1/2))]
		} /.{
		a -> (5 - 5^(1/2))^(1/2),
		b -> (5 + 5^(1/2))^(1/2)
		} /.{Sqrt[w_] :> Sqrt[w//Expand]}

T1[a_, b_, r_] := 
	(512 a^5 - 15625 b^4 + 768 a^4 r + 416a^3 r^2 + 
	112 a^2 r^3 + 24 a r^4 + 4 r^5)/(50 b^3);
T2[a_, b_, r_] := 
	(3840 a^5 - 78125b^4+4480a^4 r +2480a^3 r^2 +
	760a^2 r^3+ 140a r^4+30r^5)/(512a^5 b + 6250 b^5);
T3[a_, b_, r_] := 
	(-18880 a^5 + 781250b^4-34240a^4 r-21260a^3 r^2-
	5980a^2 r^3-1255 a r^4-240r^5)/(2 b^2);
T4[a_, b_, r_] := 
	(68800 a^5+25000a^4 r + 11500a^3 r^2+3250a^2 r^3+
	375 a r^4+100 r^5)/(512 a^5+6250 b^4);
V[a_, b_, r_] := 
	(-1036800 a^5+48828125 b^4-2280000a^4 r-
	1291500a^3 r^2-399500a^2 r^3-76625a r^4-16100r^5)/
	(256a^5+3125b^4);

Example

SolvableQ[x^5 + 20 x + 32, x]

??SolvableQ

SolveQuintic[x^5 + 20 x + 32 == 0, x]

N[%]

NRoots[x^5 + 20 x + 32 == 0, x]

SolveQuintic[x^5 + 20 x + 31 == 0, x]