Galois groups of quintics

GaloisGroup::usage = 
"GaloisGroup[quintic] computes the Galois group of the given
quintic.";

GaloisGroup::red = "The quintic is reducible.";

GaloisGroup[expr_] :=
	Module[ {var},
		var = Union[Select[Level[expr, {-1}], 
			  Head[#]===Symbol&]];
		If[ Length[var] == 1 && Exponent[expr,var[[1]]] == 5
		    , 
		    galois[expr, var[[1]]]
		    ,
		    Return[$Failed]
		]
	]
	
galois[expr_, var_] :=
	If[Exponent[expr, var] == 5 &&
	   And@@((# === Integer || # === Rational)&/@(
	   Head/@CoefficientList[expr,var]))
	   ,    
	   Module[ {factor},
		    If[ Head[factor=Factor[expr]] === Times &&
		    Length[Complement[Exponent[#, var]&/@(List@@factor),
		    		{0,5}]] > 0 
		    ,
		    Message[GaloisGroup::red];
		    Select[
		    	Irreducible[#,var,Exponent[#,var]]& /@ 
				(List@@factor), # =!= NULL &]
		    ,
		    If[ Head[factor=Factor[expr]] === Power
		        ,
		        Message[GaloisGroup::red];
			    Table[SymmetricGroup[1],{i,Exponent[expr,var]}]
			    ,
		        {Irreducible[expr, var, Exponent[expr,var]]}
		    ]
		]
	  ] 
	  ,
	  $Failed
	]

Irreducible[expr_, var_, n_] :=
	Module[ {isfactorable, introot, roots},
		{isfactorable, introot, roots} = 
			Resolvent[expr/Coefficient[expr,var,n], var];
		If[ isfactorable
	    	,
		    SolvableGroup[Sqrt[Discriminant[expr, var,5]], 
		    	expr, var, introot, roots]
		    ,
		    If[ IntegerQ[Sqrt[Discriminant[expr, var,5]]]
			    ,
			    AlternatingGroup[n] 
			    ,
			    SymmetricGroup[n] 
		    ]
		    ,
		    $Failed
		]
	] /; n == 5

Irreducible[__] := NULL
	
Discriminant[poly_, var_, n_]:=
	Resultant[poly,D[poly,var],var]/Coefficient[poly,var,n] *
	(-1)^(n (n-1)/2)
      

Resolvent[expr_, var_] :=
	Module[ {roots, integralroot, temp},

		roots = (List @@ NRoots[expr==0, var,20])/.
					Equal[_,p_] :> p;
		temp = Stabilizer/@ Representatives[roots];
		
		integralroot = 
			Position[ Chop/@((Round[#] - #)&/@temp), _Integer];

		If[ Length[integralroot] >= 1
		    ,
		    {True, temp[[integralroot[[1,1]]]]//Round, roots}
		    ,
		    {False, $Failed, $Failed}
		]
	]


Representatives[{x1_,x2_,x3_,x4_,x5_}] :=
	{{x1,x2,x3,x4,x5},{x2,x1,x3,x4,x5},{x3,x2,x1,x4,x5},
	 {x4,x2,x3,x1,x5},{x5,x2,x3,x4,x1},{x1,x5,x3,x4,x2}}

Stabilizer[{x1_,x2_,x3_,x4_,x5_}] :=
	x1^2 (x2 x5 + x3 x4) + x2^2(x1 x3 + x4 x5) + 
	x3^2 (x1 x5 + x2 x4) + x4^2 (x1 x2 + x3 x5) + 
	x5^2 (x1 x4 + x2 x3)

SolvableGroup[disc_, r__] :=
	If[ !IntegerQ[disc]
	    ,
	    MetacyclicGroup[20] 
	    ,
	    DihedralOrCyclic[r]
	]
	
DihedralOrCyclic[expr_, var_, introot_, roots_List] := 
	Module[ {disc, p},
		
	Fnew[x1_,x2_,x3_,x4_,x5_] := 
		x1 x2^2 + x2 x3^2 + x3 x4^2 + x4 x5^2 + x5 x1^2;
	
	NewResolvent[{x1_,x2_,x3_,x4_,x5_}] :=
		(var - Fnew[x1,x2,x3,x4,x5]) *
		(var - Fnew[x2,x1,x5,x4,x3]);
		
		p = Position[Round /@ (Stabilizer /@ 
				Representatives[roots]), introot][[1,1]];
		disc = Discriminant[NewResolvent[
			Representatives[roots][[p]] ], var, 2];
		
		disc = disc /. {e_Real :> Round[e], 
		                e_Complex :> Round[e]};
		
		If[ disc == 0, 
		    Return[SingularCase[expr, var, roots, 5]] 
		];
		
		If[ IntegerQ[Sqrt[disc]]
		    ,
		    CyclicGroup[5]  
		    ,
		    DihedralGroup[10] 
		]
	]

SingularCase[expr_, var_, roots_, n_] :=
	Module[ {newp, newroots,stabs,realr,p,disc},
		newp = Tschirnhaus[expr,var,n];
		newroots = newp /. ({var->#}&/@roots);
		stabs = Stabilizer /@ N[Representatives[newroots],32];
		realr = Select[stabs//Round, IntegerQ];
		If[ Length[realr] == 0
		    ,
		    Return[$Failed]
		];
		
		p = Position[stabs//Round, First[realr]] [[1,1]];
		
		disc = Discriminant[NewResolvent[
			Representatives[newroots][[p]] ], var, 2];
			
		If[ IntegerQ[Sqrt[disc]]
		    ,
		    CyclicGroup[5]  
		    ,
		    DihedralGroup[10] 
		]
	]
	
Tschirnhaus[expr_, var_, n_] :=
	Module[{z, randp, newp},
		randp = Sum[ Random[Integer,{-1,0}] z^(n-k-2), 
					 {k,0,n-2}];
		newp = Resultant[expr/.var->z, var - randp, z];

		If[ FreeQ[PolynomialGCD[newp, D[newp, var]], var]
		    ,
		    newp
		    ,
		    Tschirnhaus[expr,var,n]
		]
	]

Examples.

GaloisGroup[x^5 + 20 x + 32]

GaloisGroup[x^5 + x^4 - 4 x^3 - 3 x^2 + 3 x + 1]

GaloisGroup[x^5 + x^4 - 12 x^3 - 21 x^2 + x + 5]

GaloisGroup[x^5 + 2]

GaloisGroup[x^5 + 20 x + 16]

GaloisGroup[x^5 - 5 x + 12]

GaloisGroup[x^5 + 15 x + 12]

GaloisGroup[x^5 + 5 x^3 + 10 x^2 + 5 x + 4]

GaloisGroup[y^5 - 33826005 y - 4140303012 ]