CubicResolvent := proc ( f::polynom )
	local a, b, c, d, p, q, r, s, t, g, discr, sqdr;
	if irreduc(f)
 		then 
		a := coeff(f,x,3):
		b := coeff(f,x,2):
		c := coeff(f,x,1):
		d := coeff(f,x,0):
		p := (-3*a^2)/8 + b:
		q := a^3/8 - (a*b)/2 + c:
		r := (-3*a^4)/256 + (a^2*b)/16 - (a*c)/4 + d:
		s := p^2-4*r:
		t := -q^2:
		g := sort(x^3 + 2*p*x^2 + s*x + t):
		print(` The cubic resolvent is`, g);
		discr := 4*p^2*s^2 + 36*p*s*t - 32*p^3*t - 4*s^3 - 27*t^2;
		print(`The discriminant is` , discr);
		sqdr := sqrt(discr):
		print(`The square root of the discriminant is `, sqdr);
		if irreduc(g)
			then print(`The cubic resolvent does not factor.`);
			else print(`The cubic resolvent factors.`)
			fi;
		galoisgroup( typematch(sqdr,integer), irreduc(g) );
	else print(`The quartic is reducible.`);
	fi;
end:

galoisgroup := proc ( v::boolean, w::boolean )
	if not v and w
		then print(`So the Galois Group is S4.`);
	elif v and w
		then print(`So the Galois Group is A4.`);
	elif v and not w
		then print(`So the Galois Group is V.`);
	else print(`So the Galois Group is D4 or C4.`)
	fi;
end:

coeffsredquartic := proc( p::polynom )
	local a,q;
	a := coeff( p, x^3);
	q := expand( subs( x=x-a/4, p ) );
	[ coeff(q,x,0), coeff(q,x,1), coeff(q,x,2) ];
end:

ellipse := proc (a, b, c, t)
	local e,f,g,h,s;
	e := evalf( sqrt( b ^ 2 / (4*t) + (c - t) ^ 2 / 4 - a ) );
	f := evalf( sqrt(t) );
	g := b/(2*f);
	h := (c-t)/2;
	plot( [ (e*sin(s) - g)/f, e*cos(s) - h, s=0..2*Pi ], colour=black, axes=none );
end:

hyperbola := proc(a, b, c, t)
	local e,f,g,h,k,s;
	e := evalf( sqrt(-t) );
	f := -b /(2*sqrt(-t));
	g := (c - t)/2;
	h := b^2/(4*t) + (c - t)^2/4 - a;
	if h < 0 
		then k := evalf( sqrt(-h) );
		plot( {[ (k*cosh(s) - f)/e, k*sinh(s) - g, s=-5..5 ], [ -(k*cosh(s) + f)/e, -k*sinh(s) - g, s=-5..5 ]},
		colour=[black,black], axes=none );
		else k := evalf( sqrt(h) );
		plot( {[ (k*sinh(s) - f)/e, k*cosh(s) - g , s=-5..5 ],  [ -(k*sinh(s) + f)/e, -k*cosh(s) - g, s=-5..5 ]},
		colour=[black,black], axes=none );
	fi;
end:

parabola := proc( a, b, c, xmin, xmax)
	local s;
	if b <> 0 then
		plot( -(s^2 + c*s + a)/b, s=xmin..xmax, axes=none);
	else
	fi;
end:

line := proc( q::list,r::list )
	plot( (x-q[1])*(r[2]-q[2])/(r[1]-q[1]) + q[2], x=xmin..xmax );
end:

QuarticPlot := proc( p::polynom )
	local s,dx,c;
	global xmin,xmax;
	s := [fsolve(p)];
	if nops(s) < 4 then
		print(`The quartic does not have four real roots!`);
	else
		dx := abs( max(op(s)) - min(op(s)) );
		xmin := min(op(s)) - dx;
		xmax := max(op(s)) + dx;
		c := coeffsredquartic(p);
		display( { seq( hyperbola(op(1,c),op(2,c),op(3,c),t), t=-10..-1 ), seq( ellipse(op(1,c),op(2,c),op(3,c),t), t=1..10 ),
			parabola(op(1,c),op(2,c),op(3,c),xmin,xmax), line( [s[1],s[1]^2], [s[2],s[2]^2] ), line( [s[3],s[3]^2], [s[4],s[4]^2] ),
			line( [s[1],s[1]^2], [s[3],s[3]^2] ), line( [s[2],s[2]^2], [s[4],s[4]^2] ),
			line( [s[1],s[1]^2], [s[4],s[4]^2] ), line( [s[2],s[2]^2], [s[3],s[3]^2] ) } , 
			view=[-6..6,-5..15], scaling=CONSTRAINED );
	fi;
end: