findQuarticRoots<-function(a) {
if (is.nul(a[5]))
	return(findCubicRoots(a[1:4]))
a3<-a[4]/a[5]; a2<-a[3]/a[5]; a1<-a[2]/a[5]; a0<-a[1]/a[5]
p<-a2-(3*(a3^2)/8); q<-((a3^3)/8)-(a2*a3/2)+a1
r<-a0+(a2*(a3^2)/16)-((a3*a1)/4)-(3*(a3^4)/256); sb<-a3/4
if (is.nul(r)) {
	cr<-findCubicRoots(c(q,p,0,1))
	return(lapply(c(list(0),cr),function(x) x-sb))
	}
else {
	cr<-findCubicRoots(c(((r*p)/2)-((q^2)/8),-r,-p/2,1))
	z<-cr[[1]]; u<-(z^2)-r; v=(2*z)-p
	if (is.nul(u)) u<-0 else if (is.pos(u)) u<-sqrt(u) else return(list())
	if (is.nul(v)) v<-0 else if (is.pos(v)) v<-sqrt(v) else return(list())
	y1<-findQuadraticRoots(c(z-u,v,1))
	y2<-findQuadraticRoots(c(z+u,-v,1))
	return(lapply(c(y1,y2),function(x) x-sb))
	}
}

findCubicRoots<-function(a) {
if (is.nul(a[4])) 
	return(findQuadraticRoots(a[1:3]))
p<-a[3]/a[4]; q<-a[2]/a[4]; r<-a[1]/a[4]
p3<-p/3; r3<-1/3; pi3<-pi/3; s3<-sqrt(3)
u<-(q-((p^2)/3))/3; v<-(r-(p*q/3)+((2*(p^3))/27))/2
d<-(u^3)+(v^2)
if (is.nul(d)) {
	if (is.nul(v)) return(list(-p3))
	else {
		vv<-(-v)^r3
	    return(list((2*vv)-p3,-vv-p3))
		}
	}
if (is.neg(d)) {
	s<-2*sqrt(-u); t<-acos(-v/sqrt(-u^3))
	y1<-(s*cos(t))-p3
	y2<-(-s*cos(t+pi3))-p3
	y3<-(-s*cos(t-pi3))-p3
	return(list(y1,y2,y3))
	}
if (is.pos(d)) {	
	w<-sqrt(d); a<-(w-v)^r3; b<--(w+v)^r3
	return(list(a+b-p3))
	}
}

findQuadraticRoots<-function(a) {
if (is.nul(a[3]))
	return(findLinearRoots(a[1:2]))
p<-a[2]/a[3]; q<-a[1]/a[3]; d<-(p^2)-(4*q)
if (is.neg(d)) 
	return(list())
if (is.pos(d)) 
	return(list(((-p)+sqrt(d))/2,((-p)-sqrt(d))/2))
if (is.nul(d)) 
	return(list(-p/2))
}

findLinearRoots<-function(a) {
if (!is.nul(a[2])) return(list(-a[1]/a[2]))
else return(list())
}

minimizeQuartic<-function(a) {
#
# D.J. Jeffrey: Formulae, algorithms, and quartic extrema. 
# Mathematics Magazine 70, pp349 - 356 (1997).
# http://www.apmaths.uwo.ca/~djeffrey/Offprints/matmag97.ps 
#
a0<-a[1]; a1<-a[2]; a2<-a[3]; a3<-a[4]; a4<-a[5]
if (!is.pos(a4)) stop("no minimum")
b2<-(a2/(3*a4))-((a3*a3)/(8*a4*a4))
b1<-(a1/(2*a4))+((a3*a3*a3)/(16*a4*a4*a4))-((a2*a3)/(4*a4*a4))
if (is.nul(b1)) {
    bm<-min(0,b2)
	return(c(sqrt(-3*bm/2),-9*bm*bm/4))
	}
b0<-((a2*a3*a3)/(16*a4*a4))-((a1*a3)/(4*a4))-((3*a3*a3*a3*a3)/(256*a4*a4*a4))
ss<-(b1*b1)+(b2*b2*b2)+sqrt((b1*b1*b1*b1)+(2*b1*b1*b2*b2*b2))
kf<-(ss^(1/3))+(b2*b2*(ss^(-1/3)))+b2
kk<--(3/4)*(kf-b2)*(kf-(3*b2))
return(c(-b1/kf,(a4*kk)+a0+b0))
}

is.nul<-function(x) {
return(abs(x) < 1e-10)
}

is.pos<-function(x) {
return(x > 1e-10)
}

is.neg<-function(x) {
return(x < -1e-10)
}