四次方程　quartic equation 需修正！

function cubicRoot(x1)
{
if (x1 > 0.0) {
return x1^(1/3);
} else {
return -(-x1)^(1/3);
}
}

b3 = -a3^3/(a4^3) + 4*a2*a3/a4^2 - 8*a1/a4;
b5 = a2^2 - 3*a1*a3 + 12*a0*a4;
b6 = 2*a2^3 - 9*a1*a2*a3 + 27*a0*a3^2 + 27*a1^2*a4 - 72*a0*a2*a4;
b7 = -4*b5^3 + b6^2;
b8 = b6 + sqrt(b7);

crb8 = cubicRoot(b8);

b4 = a3^2/(4*a4^2) + 2^(1/3)*b5/(3*crb8*a4) + crb8/(3*2^(1/3)*a4) - 2*a2/(3*a4);
b2 = -b3/(b4)^(1/4) + a3^2/(2*a4^2) - 2^(1/3)*b5/(3*crb8*a4) - crb8/(3*2^(1/3)*a4) - 4*a2/(3*a4);
b1 = b3/b4^(1/4) + a3^2/(2*a4^2) - 2^(1/3)*b5/(3*crb8*a4) - crb8/(3*2^(1/3)*a4) - 4*a2/(3*a4);

x1 = -a3/(4*a4) - sqrt(b2)/2 - sqrt(b4)/2;
x2 = -a3/(4*a4) + sqrt(b2)/2 - sqrt(b4)/2;
x3 = -a3/(4*a4) - sqrt(b1)/2 + sqrt(b4)/2;
x4 = -a3/(4*a4) + sqrt(b1)/2 + sqrt(b4)/2;

print(b1);
print(b2);
print(b3);
print(b4);
print(b5);
print(b6);
print(b7);
print(b8);
print(crb8);
print(x1);
print(x2);
print(x3);
print(x4);