jhshukla
asked on
Solve quartic using Ferrari's method
I wrote the C++ code for Ferrari's method described at http://en.wikipedia.org/wiki/Quartic_equation (skip to "Summary of Ferrari's method")
void SolveQuarticFerrari(double A, double B, double C, double D, double E, double roots[4])
{
complex cp_roots[4];
#ifdef _DEBUG // test for correctness
cp_roots[0] = complex(0,1)*complex(1,1);
cp_roots[1] = complex(0,1)/complex(1,1);
cp_roots[2] = complex(3,1)-complex(7,2);
cp_roots[3] = complex(2,1)+complex(9,1);
cp_roots[0] = cp_sqrt(cp_roots[0]);
cp_roots[1] = cp_cuberoot(cp_roots[1]);
cp_roots[2] = -cp_roots[2];
cp_roots[3] = cube(cp_roots[3]);
#endif
double alpha = -3/8*sqr(B/A) + C/A;
double beta = (B*B*B)/(8*A*A*A) - (B*C)/(2*A*A) + D/A;
double gamma = (-3*B*B*B*B)/(256*A*A*A*A) + (C*B*B)/(16*A*A*A) - (B*D)/(4*A*A) + E/A;
if ( SimdEqual(beta, 0) )
{
complex u1 = (complex(-alpha) + cp_sqrt(sqr(alpha)-4*gamma ))/complex (2);
complex u2 = (complex(-alpha) - cp_sqrt(sqr(alpha)-4*gamma ))/complex (2);
cp_roots[0] = +cp_sqrt(u1) - complex(B/(4*A));
cp_roots[1] = -cp_sqrt(u1) - complex(B/(4*A));
cp_roots[2] = +cp_sqrt(u2) - complex(B/(4*A));
cp_roots[3] = -cp_sqrt(u2) - complex(B/(4*A));
}
else
{
double P = -sqr(alpha)/12 - gamma;
double Q = -cube(alpha)/108 + alpha*gamma/3 - sqr(beta)/8;
complex R = complex(Q/2)+cp_sqrt(Q*Q/4 +P*P*P/27) ;
complex U = isreal(R) ? cubic_rt(R.real) : cp_cuberoot(R);
complex y = complex(-5/6*alpha) - U;
if ( U )
y = y + complex(P/3)/U;
complex W = cp_sqrt(complex(alpha)+com plex(2)*y) ;
cp_roots[0] = complex(-B/(4*A)) + (+W+cp_sqrt(-(complex(3*al pha)+compl ex(2)*y+co mplex(2*be ta)/W)))/c omplex(2);
cp_roots[1] = complex(-B/(4*A)) + (+W-cp_sqrt(-(complex(3*al pha)+compl ex(2)*y+co mplex(2*be ta)/W)))/c omplex(2);
cp_roots[2] = complex(-B/(4*A)) + (-W+cp_sqrt(-(complex(3*al pha)+compl ex(2)*y-co mplex(2*be ta)/W)))/c omplex(2);
cp_roots[3] = complex(-B/(4*A)) + (-W-cp_sqrt(-(complex(3*al pha)+compl ex(2)*y-co mplex(2*be ta)/W)))/c omplex(2);
}
if ( isreal(cp_roots[0]) ) roots[0] = cp_roots[0].real;
if ( isreal(cp_roots[1]) ) roots[1] = cp_roots[1].real;
if ( isreal(cp_roots[2]) ) roots[2] = cp_roots[2].real;
if ( isreal(cp_roots[3]) ) roots[3] = cp_roots[3].real;
}
complex is defined as
struct complex { double real, imag; /* plus some methods */ };
Here are my two sample inputs, outputs and real answers:
SAMPLE#1 (does not pass beta=0 test)
input:
A 1.0000000000000000 double
B 0.00000000000000000 double
C 6.0000000000000000 double
D -60.000000000000000 double
E 36.000000000000000 double
output:
cp_roots complex [4]
+ [0x0] {real=2.4504888521002162 imag=1.6981366294306341 } complex
+ [0x1] {real=2.4504888521002171 imag=-1.6981366294306341 } complex
+ [0x2] {real=-2.4504888521002162 imag=3.8892317000046215 } complex
+ [0x3] {real=-2.4504888521002171 imag=-3.8892317000046215 } complex
correct answer:
-1.87213664412, +3.8101353368
-1.87213664412, -3.8101353368
3.09987442402, 0
.644398864227, 0
=======================
SAMPLE#2: (passes beta=0 test)
input:
A 0.015624996274709951 double
B 8.5355328884178014 double
C 1086.1359215487923 double
D -21727.922061357851 double
E 100574.00000000000 double
output:
cp_roots complex [4]
+ [0x0] {real=-136.56855877516438 imag=84.916864289141827 } complex
+ [0x1] {real=-136.56855877516438 imag=-84.916864289141827 } complex
+ [0x2] {real=-136.56855877516435 imag=249.60336879045511 } complex
+ [0x3] {real=-136.56855877516441 imag=-249.60336879045511 } complex
correct answer:
-282.871516735
-281.442836722
9.73440081856
8.30571753803
What am I missing here?
void SolveQuarticFerrari(double
{
complex cp_roots[4];
#ifdef _DEBUG // test for correctness
cp_roots[0] = complex(0,1)*complex(1,1);
cp_roots[1] = complex(0,1)/complex(1,1);
cp_roots[2] = complex(3,1)-complex(7,2);
cp_roots[3] = complex(2,1)+complex(9,1);
cp_roots[0] = cp_sqrt(cp_roots[0]);
cp_roots[1] = cp_cuberoot(cp_roots[1]);
cp_roots[2] = -cp_roots[2];
cp_roots[3] = cube(cp_roots[3]);
#endif
double alpha = -3/8*sqr(B/A) + C/A;
double beta = (B*B*B)/(8*A*A*A) - (B*C)/(2*A*A) + D/A;
double gamma = (-3*B*B*B*B)/(256*A*A*A*A)
if ( SimdEqual(beta, 0) )
{
complex u1 = (complex(-alpha) + cp_sqrt(sqr(alpha)-4*gamma
complex u2 = (complex(-alpha) - cp_sqrt(sqr(alpha)-4*gamma
cp_roots[0] = +cp_sqrt(u1) - complex(B/(4*A));
cp_roots[1] = -cp_sqrt(u1) - complex(B/(4*A));
cp_roots[2] = +cp_sqrt(u2) - complex(B/(4*A));
cp_roots[3] = -cp_sqrt(u2) - complex(B/(4*A));
}
else
{
double P = -sqr(alpha)/12 - gamma;
double Q = -cube(alpha)/108 + alpha*gamma/3 - sqr(beta)/8;
complex R = complex(Q/2)+cp_sqrt(Q*Q/4
complex U = isreal(R) ? cubic_rt(R.real) : cp_cuberoot(R);
complex y = complex(-5/6*alpha) - U;
if ( U )
y = y + complex(P/3)/U;
complex W = cp_sqrt(complex(alpha)+com
cp_roots[0] = complex(-B/(4*A)) + (+W+cp_sqrt(-(complex(3*al
cp_roots[1] = complex(-B/(4*A)) + (+W-cp_sqrt(-(complex(3*al
cp_roots[2] = complex(-B/(4*A)) + (-W+cp_sqrt(-(complex(3*al
cp_roots[3] = complex(-B/(4*A)) + (-W-cp_sqrt(-(complex(3*al
}
if ( isreal(cp_roots[0]) ) roots[0] = cp_roots[0].real;
if ( isreal(cp_roots[1]) ) roots[1] = cp_roots[1].real;
if ( isreal(cp_roots[2]) ) roots[2] = cp_roots[2].real;
if ( isreal(cp_roots[3]) ) roots[3] = cp_roots[3].real;
}
complex is defined as
struct complex { double real, imag; /* plus some methods */ };
Here are my two sample inputs, outputs and real answers:
SAMPLE#1 (does not pass beta=0 test)
input:
A 1.0000000000000000 double
B 0.00000000000000000 double
C 6.0000000000000000 double
D -60.000000000000000 double
E 36.000000000000000 double
output:
cp_roots complex [4]
+ [0x0] {real=2.4504888521002162 imag=1.6981366294306341 } complex
+ [0x1] {real=2.4504888521002171 imag=-1.6981366294306341 } complex
+ [0x2] {real=-2.4504888521002162 imag=3.8892317000046215 } complex
+ [0x3] {real=-2.4504888521002171 imag=-3.8892317000046215 } complex
correct answer:
-1.87213664412, +3.8101353368
-1.87213664412, -3.8101353368
3.09987442402, 0
.644398864227, 0
=======================
SAMPLE#2: (passes beta=0 test)
input:
A 0.015624996274709951 double
B 8.5355328884178014 double
C 1086.1359215487923 double
D -21727.922061357851 double
E 100574.00000000000 double
output:
cp_roots complex [4]
+ [0x0] {real=-136.56855877516438 imag=84.916864289141827 } complex
+ [0x1] {real=-136.56855877516438 imag=-84.916864289141827 } complex
+ [0x2] {real=-136.56855877516435 imag=249.60336879045511 } complex
+ [0x3] {real=-136.56855877516441 imag=-249.60336879045511 } complex
correct answer:
-282.871516735
-281.442836722
9.73440081856
8.30571753803
What am I missing here?
ASKER CERTIFIED SOLUTION
membership
This solution is only available to members.
To access this solution, you must be a member of Experts Exchange.
-3/8==0 and -5/6==0 seem to be the problem, I was able to duplicate your error when I choose the right cube root.
SimdEqual(beta, 0) may be problematic too. floating point rounding errors may mean beta is not exactly 0
SimdEqual(beta, 0) may be problematic too. floating point rounding errors may mean beta is not exactly 0
What do you get for alpha,beta,gamma,P,Q,R,U,y