We help IT Professionals succeed at work.

We've partnered with Certified Experts, Carl Webster and Richard Faulkner, to bring you two Citrix podcasts. Learn about 2020 trends and get answers to your biggest Citrix questions!Listen Now

x

Solve quartic using Ferrari's method

on
Medium Priority
1,875 Views
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)+complex(2)*y);

cp_roots[0] = complex(-B/(4*A)) + (+W+cp_sqrt(-(complex(3*alpha)+complex(2)*y+complex(2*beta)/W)))/complex(2);
cp_roots[1] = complex(-B/(4*A)) + (+W-cp_sqrt(-(complex(3*alpha)+complex(2)*y+complex(2*beta)/W)))/complex(2);
cp_roots[2] = complex(-B/(4*A)) + (-W+cp_sqrt(-(complex(3*alpha)+complex(2)*y-complex(2*beta)/W)))/complex(2);
cp_roots[3] = complex(-B/(4*A)) + (-W-cp_sqrt(-(complex(3*alpha)+complex(2)*y-complex(2*beta)/W)))/complex(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

-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

-282.871516735
-281.442836722
9.73440081856
8.30571753803

What am I missing here?
Comment
Watch Question

View Solution Only

CERTIFIED EXPERT
Most Valuable Expert 2014
Top Expert 2015
Commented:
-3/8 and -5/6 are integer divisions, which = 0

Not the solution you were looking for? Getting a personalized solution is easy.

CERTIFIED EXPERT
Most Valuable Expert 2014
Top Expert 2015

Commented:
how do you define plus some methods ?
What do you get for alpha,beta,gamma,P,Q,R,U,y,W with SolveQuarticFerrari(1.0,0.0,6.0,-60.0,36.0)?
CERTIFIED EXPERT
Most Valuable Expert 2014
Top Expert 2015

Commented:
-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
Thanks for using Experts Exchange.

• View three pieces of content (articles, solutions, posts, and videos)
• Ask the experts questions (counted toward content limit)
• Customize your dashboard and profile