Solved

Solve quartic using Ferrari's method

Posted on 2006-11-03
3
1,742 Views
Last Modified: 2012-06-21
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

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?
0
Comment
Question by:jhshukla
  • 3
3 Comments
 
LVL 84

Accepted Solution

by:
ozo earned 125 total points
Comment Utility
-3/8 and -5/6 are integer divisions, which = 0
0
 
LVL 84

Expert Comment

by:ozo
Comment Utility
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)?
0
 
LVL 84

Expert Comment

by:ozo
Comment Utility
-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
0

Featured Post

How to run any project with ease

Manage projects of all sizes how you want. Great for personal to-do lists, project milestones, team priorities and launch plans.
- Combine task lists, docs, spreadsheets, and chat in one
- View and edit from mobile/offline
- Cut down on emails

Join & Write a Comment

Often, when implementing a feature, you won't know how certain events should be handled at the point where they occur and you'd rather defer to the user of your function or class. For example, a XML parser will extract a tag from the source code, wh…
Many modern programming languages support the concept of a property -- a class member that combines characteristics of both a data member and a method.  These are sometimes called "smart fields" because you can add logic that is applied automaticall…
The goal of the tutorial is to teach the user how to use functions in C++. The video will cover how to define functions, how to call functions and how to create functions prototypes. Microsoft Visual C++ 2010 Express will be used as a text editor an…
The viewer will learn additional member functions of the vector class. Specifically, the capacity and swap member functions will be introduced.

744 members asked questions and received personalized solutions in the past 7 days.

Join the community of 500,000 technology professionals and ask your questions.

Join & Ask a Question

Need Help in Real-Time?

Connect with top rated Experts

18 Experts available now in Live!

Get 1:1 Help Now