Solved

Solve quartic using Ferrari's method

Posted on 2006-11-03
3
1,750 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
ID: 17871538
-3/8 and -5/6 are integer divisions, which = 0
0
 
LVL 84

Expert Comment

by:ozo
ID: 17871593
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
ID: 17871730
-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

Is Your Active Directory as Secure as You Think?

More than 75% of all records are compromised because of the loss or theft of a privileged credential. Experts have been exploring Active Directory infrastructure to identify key threats and establish best practices for keeping data safe. Attend this month’s webinar to learn more.

Question has a verified solution.

If you are experiencing a similar issue, please ask a related question

Suggested Solutions

Title # Comments Views Activity
object modeling - load and save functions 5 60
Unable to start eclipse ? 17 133
Issues with C++ Class 19 81
Copy output image from TWindowsMediaPlayer 6 23
Written by John Humphreys C++ Threading and the POSIX Library This article will cover the basic information that you need to know in order to make use of the POSIX threading library available for C and C++ on UNIX and most Linux systems.   [s…
This article will show you some of the more useful Standard Template Library (STL) algorithms through the use of working examples.  You will learn about how these algorithms fit into the STL architecture, how they work with STL containers, and why t…
The viewer will learn additional member functions of the vector class. Specifically, the capacity and swap member functions will be introduced.
The viewer will be introduced to the member functions push_back and pop_back of the vector class. The video will teach the difference between the two as well as how to use each one along with its functionality.

867 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

19 Experts available now in Live!

Get 1:1 Help Now