Solved

Solve quartic using Ferrari's method

Posted on 2006-11-03
3
1,772 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
[X]
Welcome to Experts Exchange

Add your voice to the tech community where 5M+ people just like you are talking about what matters.

  • Help others & share knowledge
  • Earn cash & points
  • Learn & ask questions
  • 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

Free Tool: Subnet Calculator

The subnet calculator helps you design networks by taking an IP address and network mask and returning information such as network, broadcast address, and host range.

One of a set of tools we're offering as a way of saying thank you for being a part of the community.

Question has a verified solution.

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

Suggested Solutions

Title # Comments Views Activity
C++ to C# code conversion issue 4 134
draw a Christmas tree by using a nested loop? 26 106
Indy 10 not Receiving UDP broadcast 3 57
Coding for the first time 9 104
This article shows you how to optimize memory allocations in C++ using placement new. Applicable especially to usecases dealing with creation of large number of objects. A brief on problem: Lets take example problem for simplicity: - I have a G…
Go is an acronym of golang, is a programming language developed Google in 2007. Go is a new language that is mostly in the C family, with significant input from Pascal/Modula/Oberon family. Hence Go arisen as low-level language with fast compilation…
The goal of the video will be to teach the user the concept of local variables and scope. An example of a locally defined variable will be given as well as an explanation of what scope is in C++. The local variable and concept of scope will be relat…
The viewer will learn how to pass data into a function in C++. This is one step further in using functions. Instead of only printing text onto the console, the function will be able to perform calculations with argumentents given by the user.

752 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