2012-10-20 3 views
2

내 프로젝트의 일부로 마감 된 형태 (C++의 경우)에서 4 차 다항식을 풀어야합니다.4 차 루트 (4 차 다항식)에 대한 C++ 해결

A * X 4 + B * X 3 + C * X 2 + D * X + E = 0

I는이를 여러 링크를 발견했다. 그 중 하나는 here입니다. 그러나 그것은 모든 뿌리를 계산하는 반면, 나는 진정한 뿌리를 원합니다. 알고리즘은 주로 페라리의 방법을 사용하여 순서를 줄입니다.

다음
bool solveCubic(double &a, double &b, double &c, double &d, double &root) 
{ 
if(a == 0.0 || abs(a/b) < 1.0e-6) 
    return solveQuadratic(b, c, d, root); 

double B = b/a, C = c/a, D = d/a; 

double Q = (B*B - C*3.0)/9.0, QQQ = Q*Q*Q; 
double R = (2.0*B*B*B - 9.0*B*C + 27.0*D)/54.0, RR = R*R; 

// 3 real roots 
if(RR<QQQ) 
{ 
    /* This sqrt and division is safe, since RR >= 0, so QQQ > RR, */ 
    /* so QQQ > 0. The acos is also safe, since RR/QQQ < 1, and  */ 
    /* thus R/sqrt(QQQ) < 1.          */ 
    double theta = acos(R/sqrt(QQQ)); 
    /* This sqrt is safe, since QQQ >= 0, and thus Q >= 0    */ 
    double r1, r2, r3; 
    r1 = r2 = r3 = -2.0*sqrt(Q); 
    r1 *= cos(theta/3.0); 
    r2 *= cos((theta+2*PI)/3.0); 
    r3 *= cos((theta-2*PI)/3.0); 

    r1 -= B/3.0; 
    r2 -= B/3.0; 
    r3 -= B/3.0; 

    root = 1000000.0; 

    if(r1 >= 0.0) root = r1; 
    if(r2 >= 0.0 && r2 < root) root = r2; 
    if(r3 >= 0.0 && r3 < root) root = r3; 

    return true; 
} 
// 1 real root 
else 
{ 
    double A2 = -pow(fabs®+sqrt(RR-QQQ),1.0/3.0); 
    if (A2!=0.0) { 
     if (R<0.0) A2 = -A2; 
     root = A2 + Q/A2; 
    } 
    root -= B/3.0; 
    return true; 
} 
} 

코드를 설명하는 몇 가지 링크입니다 :

bool solveQuartic(double a, double b, double c, double d, double e, double &root) 
{ 
// I switched to this method, and it seems to be more numerically stable. 
// http://www.gamedev.n...topic_id=451048 

// When a or (a and b) are magnitudes of order smaller than C,D,E 
// just ignore them entirely. This seems to happen because of numerical 
// inaccuracies of the line-circle algorithm. I wanted a robust solver, 
// so I put the fix here instead of there. 
if(a == 0.0 || abs(a/b) < 1.0e-5 || abs(a/c) < 1.0e-5 || abs(a/d) < 1.0e-5) 
    return solveCubic(b, c, d, e, root); 

double B = b/a, C = c/a, D = d/a, E = e/a; 
double BB = B*B; 
double I = -3.0*BB*0.125 + C; 
double J = BB*B*0.125 - B*C*0.5 + D; 
double K = -3*BB*BB/256.0 + C*BB/16.0 - B*D*0.25 + E; 

double z; 
bool foundRoot2 = false, foundRoot3 = false, foundRoot4 = false, foundRoot5 = false; 
if(solveCubic(1.0, I+I, I*I - 4*K, -(J*J), z)) 
{ 
    double value = z*z*z + z*z*(I+I) + z*(I*I - 4*K) - J*J; 

    double p = sqrt(z); 
    double r = -p; 
    double q = (I + z - J/p)*0.5; 
    double s = (I + z + J/p)*0.5; 

    bool foundRoot = false, foundARoot; 
    double aRoot; 
    foundRoot = solveQuadratic(1.0, p, q, root); 
    root -= B/4.0; 

    foundARoot = solveQuadratic(1.0, r, s, aRoot); 
    aRoot -= B/4.0; 
    if((foundRoot && foundARoot && ((aRoot < root && aRoot >= 0.0) 
     || root < 0.0)) || (!foundRoot && foundARoot)) 
    { 
     root = aRoot; 
     foundRoot = true; 
    } 

    foundARoot = solveQuadraticOther(1.0, p, q, aRoot); 
    aRoot -= B/4.0; 
    if((foundRoot && foundARoot && ((aRoot < root && aRoot >= 0.0) 
     || root < 0.0)) || (!foundRoot && foundARoot)) 
    { 
     root = aRoot; 
     foundRoot = true; 
    } 

    foundARoot = solveQuadraticOther(1.0, r, s, aRoot); 
    aRoot -= B/4.0; 
    if((foundRoot && foundARoot && ((aRoot < root && aRoot >= 0.0) 
     || root < 0.0)) || (!foundRoot && foundARoot)) 
    { 
     root = aRoot; 
     foundRoot = true; 
    } 
    return foundRoot; 
} 
return false; 
} 

는 현실과 상상의 두 솔루션을 제공 solveCubic()를 사용합니다. solveCubicsolveQuartic

실제 루트에 대한 4 차 다항식을 해결하기 위해 코드를 수정할 수있는 사람이 있습니까?

최대한 효율적으로 구현하고 싶습니다. BTW 누군가가 LAPACK 같은이 목적을 위해 유용한 라이브러리를 도입하면 감사하겠습니다 (그것은 4 차 다항식의 근원을 직접 계산할 수없는 것 같습니다).

+6

모든 뿌리를 계산 한 다음 삭제 하시겠습니까? –

+0

숫자 솔루션이나 기호 솔루션을 찾고 계십니까? – RBarryYoung

+4

우리는 "코드를 제공하지 않습니다". * 귀하의 * 코드의 오류를 이해하도록 도와드립니다. – alestanis

답변

1

아마도이 방정식을 실제 뿌리에 대해 닫힌 형식으로 풀 수있는 가장 효율적인 방법은 모든 뿌리에 대해 닫힌 형식으로 풀어서 가상의 뿌리를 버리는 것입니다.

try/catch 쌍을 사용하여 허수가 자르기를 결정할 수 있다고 생각할 수도 있지만 실제 루트를 계산할 때 생성하는 중간 값 중 일부는 허수가 될 수 있기 때문에 매우 좋은 전략은 아닙니다.

따라서 C++ 복소수 라이브러리 (here 또는 here 참조)를 사용하여 계산할 수 있습니다.

그런 다음 숫자의 허수 부 분이 0이 아닌지 확인하고 삭제하는 경우 확인하십시오. 그러나 부동 소수점 연산은 정확하지 않으므로 "0"은 0에 매우 가까운 수의 범위를 포함합니다.