Copyright 1984 by ABComputing May 15, 1984 ÉÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍ» º Roots of 3rd and 4th Order Polynomials º º º º by º º º º William Squire º ÈÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍͼ ÚÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄ¿ ³ Introduction ³ ÀÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÙ In the 16th century, Italian mathematicians Tartaglia, Ferrari, and Cardan, developed formulas for solving third and fourth order polynomial equations, with real coefficients. Today, very few computer centers have programs implementing these solutions. (They rely on iterative programs that will work for any degree polynomial.) This is not entirely due to a desire to minimize the number of subroutines in a library; the "nominally" exact solutions frequently do not give results as accurate as iterative solutions. This column discusses a FORTRAN implementation of the exact solutions for solving third and fourth order polynomials. These solutions can be found in an "Introduction to the Theory of Equations," by N.B. Conkwright, Ginn Publishers. ÚÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄ¿ ³ Exact Solution of the Cubic ³ ÀÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÙ If P3(Z) is a cubic polynomial, with real coefficients 2 3 P3(Z) = C + C X + C X + X 0 1 2 then one root, "X," is given by X = U + V - (1/3)C (1) 2 where 3 G = ( 2C - 9 C C + 27 C ) / 54 2 1 2 0 2 H = ( 3C - C ) / 9 1 2 2 3 1/2 1/3 U = [ -G + ( G + H ) ] (2) 2 3 1/2 1/3 V = [ -G - ( G + H ) ] (3) The quantity 2 3 D = G + H appearing in (2) and (3) is proportional to the discriminant of the cubic, and indicates the nature of the roots of P3(Z). In particular: if ( D > 0 ) one root is real and two are complex conjugates elseif ( D = 0 ) all roots are real and at least two are equal elseif ( D < 0 ) all roots are real and unequal endif Note in the case D < 0, all three roots are real, and that it is necessary to take the cube root of a complex number in U and V. This may have been an important factor in convincing mathematicians of the "reality" of so called "imaginary" numbers. When D<0, U, and V are complex conjugates, so X is real. The cause of inaccuracy of the exact solutions is, in many cases, the loss of significant digits in U or V when 2 3 G << H This effect may be minimized by calculating: 2 3 1/2 1/3 U = [ -G + ( G + H ) ] if -G is positive, or 2 3 1/2 1/3 V = [ -G - ( G + H ) ] if -G is negative. These forms are chosen so that essentially no subtraction, only addition, occurs. When U or V is determined the mate may be found from UV = - H and a root determined by (1). Dividing the cubic by (X-root) yields a quadratic equation that can be solved easily. ÚÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄ¿ ³ Exact Solution of the Quartic ³ ÀÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÙ The quartic equation 2 3 4 P4(Z) = C + C X + C X + C X + X 0 1 2 3 is solved by finding a real root, "Y" of the resolvent cubic: 2 2 2 3 ( 4C C - C - C C )/8 + 1/4 ( C C - 4 C ) X - 1/2 C X + X = 0 2 0 1 3 0 3 1 0 2 Using this real root, the quartic is factored as the product of two quadratics with real coefficients, and readily solved. The real coefficients mentioned are expressions involving "Y". ÚÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄ¿ ³ Programming Considerations ³ ÀÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÙ The program accompanying this article implements these exact solutions. The user is asked whether a cubic or quartic polynomial is involved, and to provide the coefficients. The values of the roots are determined and printed, and the absolute values of P3 and P4 at each root is printed. ÚÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄ¿ ³ Cubic Polynomial ³ ÀÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÙ If a cubic polynomial is requested, subroutine CUBIC is called. CUBIC calls subroutine CUR to find a real root of the polynomial, and then QDRTC to find the two remaining roots. CUR is a novelty. By treating the various possibilities separately it produces an accurate result. 1. If D < 0, the exponentiation operator (**) is used to find the cube root. CUR then uses the complex conjugate of the cube root in the calculation of the real root. The difficulty is that one must properly match the branch of the cube root used for U and V, as the FORTRAN compiler does not automatically do this. For example, if asked for a cube root of (-8.,0), the FORTRAN processor may not give (-2.,0), but one of the complex roots (1,+1.73205) or (1,-1.73205). Another difficulty is that the "**" operation does not work on negative real arguments. The function CBRT(X) is called from CUR, and it uses SIGN (ABS(X) ** .333333, X) to find the cube root for both negative and positive X. Strictly speaking, CBRT is not needed. The necessary code could be easily inserted in the two places where CBRT is called, but I believe using CBRT makes the program clearer. 2. If D = 0, the root of the cubic becomes: 1/3 2(-G) - C / 3 2 3. If D > 0, and -G is negative, then V is evaluated and U is determined from: U = - H ÄÄÄÄÄ V If D > 0, and -G is positive, then U is evaluated and V is determined from: V = - H ÄÄÄÄÄ U The idea is to match the sign of D to the sign of -G. In this manner, no subtraction, only addition, is performed. This approach minimizes round-off error. ÚÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄ¿ ³ Quartic Polynomial ³ ÀÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÙ Routine QURTC is called to find the roots of P4. It uses CUR to find a real root of the resolvent cubic, and QDRTC is called twice to find the roots of the the two quadratics mentioned in the section on quartic polynomials. QDRTC also treats positive, negative, and vanishing values of the discriminant separately, calculates the first root (R1), and finds the second root by using (R1) (R2) = constant term of quadratic ÚÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄ¿ ³ Conclusion ³ ÀÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÙ The set of subroutines presented gives accurate results for a variety of cubics and quartics with real coefficients. However, three paths were required, depending upon the sign of the discriminant. Perhaps in the future, improvements in the complex arithmetic capabilities of compilers will make it possible to implement the exact solutions for cubics and quartics with complex coefficients directly. ÚÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄ¿ ³ File Name: ÛÛ fort1.txt ÛÛ ³ ÀÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÙ