Roots of a quadratic function: Difference between revisions

Content added Content deleted

Inline

Revision as of 17:46, 28 November 2008

Write a program to find the roots of a quadratic equation. i.e solve the equation ax² + bx + c = 0.

The program must correctly handle complex roots. Error checking on the input (for a = 0) need not be shown. Quadratic equations represents a good example of how dangerous naive programming can be in numeric domain. A naive implementation of the well known formula for explicit roots of a quadratic equation suffers catastrophic loss of accuracy. Consider a naive implementation of. In their classic textbook on numeric methods Computer Methods for Mathematical Computations, George Forsythe, Michael Malcolm and Cleve Mole suggest to try it on a=1, b=-10⁵, c=1. The result of the following implementation <ada> with Ada.Text_IO; use Ada.Text_IO; with Ada.Numerics.Elementary_Functions; use Ada.Numerics.Elementary_Functions;

procedure Quadratic_Equation is

  type Roots is array (1..2) of Float;
  function Solve (A, B, C : Float) return Roots is
     SD : constant Float := sqrt (B**2 - 4.0 * A * C);
     AA : constant Float := 2.0 * A;
  begin
     return ((- B + SD) / AA, (- B - SD) / AA);
  end Solve;

  R : constant Roots := Solve (1.0, -10.0E5, 1.0);

begin

  Put_Line ("X1 =" & Float'Image (R (1)) & " X2 =" & Float'Image (R (2)));

end Quadratic_Equation; </ada> Sample output:

X1 = 1.00000E+06 X2 = 0.00000E+00

As we see the second root has suffered a total precision loss. The right answer is that x₂ is about 10^-6. The naive method is numerically unstable.

Task: do it better.

Ada

<ada> with Ada.Text_IO; use Ada.Text_IO; with Ada.Numerics.Elementary_Functions; use Ada.Numerics.Elementary_Functions;

procedure Quadratic_Equation is

  type Roots is array (1..2) of Float;
  function Solve (A, B, C : Float) return Roots is
     SD : constant Float := sqrt (B**2 - 4.0 * A * C);
     X  : Float;
  begin
     if SD > 0.0 xor B > 0.0 then
        X := (- B + SD) / 2.0 * A;
        return (X, C / (A * X));
     else
        X := (- B - SD) / 2.0 * A;
        return (C / (A * X), X);
     end if;
  end Solve;

  R : constant Roots := Solve (1.0, -10.0E5, 1.0);

begin

  Put_Line ("X1 =" & Float'Image (R (1)) & " X2 =" & Float'Image (R (2)));

end Quadratic_Equation; </ada> Here precision loss is prevented by checking signs of operands. On errors, Constraint_Error is propagated on numeric errors and when roots are complex. Sample output:

X1 = 1.00000E+06 X2 = 1.00000E-06

J

J has a built-in polynomial solver:

p.

Example (using inputs from other solutions):

   coeff =. _3 |.\ 3 4 4r3   3 2 _1   3 2 1
   > {:"1 p. coeff
         _0.666667           _0.666667
                _1            0.333333
_0.333333j0.471405 _0.333333j_0.471405

Of course this generalizes to polynomials of any order.

Fortran

Works with: Fortran version 90 and later

PROGRAM QUADRATIC

  IMPLICIT NONE
  REAL :: a, b, c, e, discriminant, rroot1, rroot2
  COMPLEX :: croot1, croot2
 
  WRITE(*,*) "Enter the coefficients of the equation ax^2 + bx + c"
  WRITE(*, "(A)", ADVANCE="NO") "a = "
  READ *, a
  WRITE(*,"(A)", ADVANCE="NO") "b = "
  READ *, b
  WRITE(*,"(A)", ADVANCE="NO") "c = "
  READ *, c
  WRITE(*,"(3(A,F10.6))") "Coefficients are: a=", a, "   b=", b, "   c=", c
 
  e = 1e-5
  discriminant = b*b - 4*a*c
 
  IF (ABS(discriminant) < e) THEN
     rroot1 = -b / (2*a)
     WRITE(*,*) "The roots are real and equal:"
     WRITE(*,*) "Root=", rroot1
  ELSE IF (discriminant > 0) THEN
     rroot1 = (-b + SQRT(discriminant)) / (2*a)
     rroot2 = (-b - SQRT(discriminant)) / (2*a)
     WRITE(*,*) "The roots are real:"
     WRITE(*,*) "Root1=", rroot1, " Root2=", rroot2
  ELSE
     croot1 = (-b + SQRT(CMPLX(discriminant))) / (2*a)	
     croot2 = CONJG(croot1)
     WRITE(*,*) "The roots are complex:"
     WRITE(*,*) "Root1=", croot1, " Root2=", croot2
  END IF

END PROGRAM QUADRATIC

Sample output

Coefficients are: a=  3.000000   b=  4.000000   c=  1.333333
The roots are real and equal:                    
Root=   -0.666667 

Coefficients are: a=  3.000000   b=  2.000000   c= -1.000000
The roots are real:
Root1=    0.333333    Root2=    -1.00000 

Coefficients are: a=  3.000000   b=  2.000000   c=  1.000000
The roots are complex:
Root1= (-0.33333,0.47140)  Root2= (-0.33333,-0.47140)

Python

No attempt is made to categorize the type of output as this is not mentioned as being part of the task. <python>>>> def quad_roots(a,b,c): a,b,c =complex(a), complex(b), complex(c) discriminant = b*b - 4*a*c root1 = (-b + discriminant**0.5)/2./a root2 = (-b - discriminant**0.5)/2./a return root1, root2

>>> for coeffs in ((3, 4, 4/3.), (3, 2, -1), (3, 2, 1)): print "a, b, c =", coeffs print " (root1, root2) =", quad_roots(*coeffs)

a, b, c = (3, 4, 1.3333333333333333)

 (root1, root2) = ((-0.66666666666666663+0j), (-0.66666666666666663+0j))

a, b, c = (3, 2, -1)

 (root1, root2) = ((0.33333333333333331+0j), (-1+0j))

a, b, c = (3, 2, 1)

 (root1, root2) = ((-0.33333333333333331+0.47140452079103173j), (-0.33333333333333331-0.47140452079103173j))

>>> </python>

Categorized roots version: <python>>>> def quad_discriminating_roots(a,b,c, entier = 1e-5): discriminant = b*b - 4*a*c a,b,c,d =complex(a), complex(b), complex(c), complex(discriminant) root1 = (-b + d**0.5)/2./a root2 = (-b - d**0.5)/2./a if abs(discriminant) < entier: return "real and equal", abs(root1), abs(root1) if discriminant > 0: return "real", root1.real, root2.real return "complex", root1, root2

>>> for coeffs in ((3, 4, 4/3.), (3, 2, -1), (3, 2, 1)): print "Roots of: %gX^2 %+gX %+g are" % coeffs print " %s: %s, %s" % quad_discriminating_roots(*coeffs)

Roots of: 3X^2 +4X +1.33333 are

 real and equal: 0.666666666667, 0.666666666667

Roots of: 3X^2 +2X -1 are

 real: 0.333333333333, -1.0

Roots of: 3X^2 +2X +1 are

 complex: (-0.333333333333+0.471404520791j), (-0.333333333333-0.471404520791j)

>>> </python>