Roots of a quadratic function

Roots of a quadratic function
You are encouraged to solve this task according to the task description, using any language you may know.

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. 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. Consider the following implementation in 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);
     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 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

Solution use J's built-in polynomial solver:

   p.

Example using inputs from other solutions and the unstable example from the task description:

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

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), (1.0, -10.0E5, 1.0)): 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)

Roots of: 1X^2 -1e+06X +1 are

 real: 999999.999999, 1.00000761449e-06

>>> </python>