Roots of a quadratic function: Difference between revisions

Content added Content deleted

Inline

Revision as of 23:14, 13 October 2009

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 Moler suggest to try it on a=1, b=-10⁵, c=1. Consider the following implementation in Ada: <lang 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; </lang> 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.

(Note that the above is using single-precision floats; while using double-precision floats will cause the above procedure to produce a nonzero result, a similar failure occurs at about b = -10⁹ = -1e9.)

Task: do it better.

Ada

<lang 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; </lang> 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

ALGOL 68

Translation of: Ada

Works with: ALGOL 68 version Standard - no extensions to language used

Works with: ALGOL 68G version Any - tested with release mk15-0.8b.fc9.i386

<lang algol>quadratic equation: BEGIN

  MODE ROOTS  = UNION([]REAL, []COMPL);
  PROC solve  = (REAL a, b, c)ROOTS:
  BEGIN
     REAL sa = b UP 2 - 4 * a * c;
     IF sa >=0 THEN # handle the +ve case as REAL #
       REAL sd = ( b < 0 | sqrt (sa) | -sqrt(sa));
       REAL x = (- b + sd) / 2 * a;
       []REAL (x, c / (a * x))
    ELSE # handle the -ve case as COMPL conjugate pairs #
       COMPL sd = ( b < 0 | complex sqrt (sa) | -complex sqrt(sa));
       COMPL x = (- b + sd) / 2 * a;
       []COMPL (x, c / (a * x))
    FI
  END # solve #;

  # only a very tiny difference between the 2 examples #
  []ROOTS test = (solve (1, -10e5, 1), solve (1, 0, 1));

  FOR index TO UPB test DO
    ROOTS r = test[index];

Output the two different scenerios #

    CASE r IN
      ([]REAL r): put(stand out,("REAL x1 = ", fixed(r[1], -real width, 8), ", ",
                                      "x2 = ", fixed(r[2], -real width, 8), new line)),
      ([]COMPL c): put(stand out,("COMPL x1,x2 =", fixed(re OF c[1], -real width, 8), " +/-", 
                                                  fixed(-im OF c[1], -real width, 8),"i", new line))
    ESAC

OD

END # quadratic_equation #</lang> Output:

REAL x1 = 999999.99999900, x2 =      0.00000100
COMPL x1,x2 =     0.00000000 +/-     1.00000000i

AutoHotkey

ahk forum: discussion <lang AutoHotkey> MsgBox % quadratic(u,v, 1,-3,2) ", " u ", " v MsgBox % quadratic(u,v, 1,3,2) ", " u ", " v MsgBox % quadratic(u,v, -2,4,-2) ", " u ", " v MsgBox % quadratic(u,v, 1,0,1) ", " u ", " v SetFormat FloatFast, 0.15e MsgBox % quadratic(u,v, 1,-1.0e8,1) ", " u ", " v

quadratic(ByRef x1, ByRef x2, a,b,c) { ; -> #real roots {x1,x2} of ax²+bx+c

  If (a = 0)
     Return -1 ; ERROR: not quadratic
  d := b*b - 4*a*c
  If (d < 0) {
     x1 := x2 := ""
     Return 0
  }
  If (d = 0) {
     x1 := x2 := -b/2/a
     Return 1
  }
  x1 := (-b - (b<0 ? -sqrt(d) : sqrt(d)))/2/a
  x2 := c/a/x1
  Return 2

}</lang>

C

<lang c>#include <stdio.h>

include <math.h>
include <complex.h>

void roots_quadratic_eq(double a, double b, double c, complex double *x) {

 double delta;

 delta = b*b - 4.0*a*c;
 x[0] = (-b + csqrt(delta)) / (2.0*a);
 x[1] = (-b - csqrt(delta)) / (2.0*a);

}</lang>

Translation of: C++

<lang c>void roots_quadratic_eq2(double a, double b, double c, complex double *x) {

 b /= a;
 c /= a;
 double delta = b*b - 4*c;
 if ( delta < 0 ) {
   x[0] = -b/2 + I*sqrt(-delta)/2.0;
   x[1] = -b/2 - I*sqrt(-delta)/2.0;
 } else {
   double root = sqrt(delta);
   double sol = (b>0) ? (-b - root)/2.0 : (-b + root)/2.0;
   x[0] = sol;
   x[1] = c/sol;
 }

}</lang>

<lang c>int main() {

 complex double x[2];

 roots_quadratic_eq(1, -1e20, 1, x);
 printf("x1 = (%.20le, %.20le)\nx2 = (%.20le, %.20le)\n\n",

creal(x[0]), cimag(x[0]), creal(x[1]), cimag(x[1]));

 roots_quadratic_eq2(1, -1e20, 1, x);
 printf("x1 = (%.20le, %.20le)\nx2 = (%.20le, %.20le)\n\n",

creal(x[0]), cimag(x[0]), creal(x[1]), cimag(x[1]));

 return 0;

}</lang>

x1 = (1.00000000000000000000e+20, 0.00000000000000000000e+00)
x2 = (0.00000000000000000000e+00, 0.00000000000000000000e+00)

x1 = (1.00000000000000000000e+20, 0.00000000000000000000e+00)
x2 = (9.99999999999999945153e-21, 0.00000000000000000000e+00)

C++

include <iostream>
include <utility>
include <complex>

typedef std::complex<double> complex;

std::pair<complex, complex>

solve_quadratic_equation(double a, double b, double c)

{

 b /= a;
 c /= a;
 double discriminant = b*b-4*c;
 if (discriminant < 0)
   return std::make_pair(complex(-b/2, std::sqrt(-discriminant)/2),
                         complex(-b/2, -std::sqrt(-discriminant)/2));

 double root = std::sqrt(discriminant);
 double solution1 = (b > 0)? (-b - root)/2
                           : (-b + root)/2;

 return std::make_pair(solution1, c/solution1);

}

int main() {

 std::pair<complex, complex> result = solve_quadratic_equation(1, -1e20, 1);
 std::cout << result.first << ", " << result.second << std::endl;

} </lang> Output:

(1e+20,0), (1e-20,0)

D

<lang d> import std.stdio; import std.math; void main() {

       writefln("%A",solve_quadratic(1+0i,0+0i,1+0i));

} creal[]solve_quadratic(creal a,creal b,creal c) {

       creal[]ret;
       ret.length = 2;
       // save the square root of the discriminant to be used later
       ret[1] = sqrt(b*b-4*a*c);
       ret[0] = (-b+ret[1])/(2*a);
       ret[1] = (-b-ret[1])/(2*a);
       return ret;

} </lang>

Forth

Without locals: <lang forth>

quadratic ( fa fb fc -- r1 r2 )

 frot frot
 ( c a b )
 fover 3 fpick f* -4e f*  fover fdup f* f+
 ( c a b det )
 fdup f0< if abort" imaginary roots" then
 fsqrt
 fover f0< if fnegate then
 f+ fnegate
 ( c a b-det )
 2e f/ fover f/  
 ( c a r1 )
 frot frot f/ fover f/ ;

</lang> With locals: <lang forth>

quadratic { F: a F: b F: c -- r1 r2 }

 b b f*  4e a f* c f* f-
 fdup f0< if abort" imaginary roots" then
 fsqrt
 b f0< if fnegate then b f+ fnegate 2e f/ a f/
 c a f/ fover f/ ;

\ test 1 set-precision 1e -1e6 1e quadratic fs. fs. \ 1e-6 1e6 </lang>

Fortran

Works with: Fortran version 90 and later

<lang fortran> PROGRAM QUADRATIC

 IMPLICIT NONE
 INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(15)
 REAL(dp) :: a, b, c, e, discriminant, rroot1, rroot2
 COMPLEX(dp) :: 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,E23.15))") "Coefficients are: a = ", a, "   b = ", b, "   c = ", c
 e = 1.0e-9_dp
 discriminant = b*b - 4.0_dp*a*c
 
 IF (ABS(discriminant) < e) THEN
    rroot1 = -b / (2.0_dp*a)
    WRITE(*,*) "The roots are real and equal:"
    WRITE(*,"(A,E23.15)") "Root = ", rroot1
 ELSE IF (discriminant > 0) THEN
    IF (b < 0) THEN
       rroot1 = (-b + SQRT(discriminant)) / (2.0_dp*a)
       rroot2 = c / (a * rroot1)
    ELSE
       rroot1 = (-b - SQRT(discriminant)) / (2.0_dp*a)
       rroot2 = c / (a * rroot1)
    END IF
    WRITE(*,*) "The roots are real:"
    WRITE(*,"(2(A,E23.15))") "Root1 = ", rroot1, "  Root2 = ", rroot2
 ELSE
    croot1 = (-b + SQRT(CMPLX(discriminant))) / (2.0_dp*a) 
    croot2 = CONJG(croot1)
    WRITE(*,*) "The roots are complex:" 
    WRITE(*,"(2(A,2E23.15,A))") "Root1 = ", croot1, "j ", "  Root2 = ", croot2, "j"
 END IF</lang>

Sample output

Coefficients are: a =   0.300000000000000E+01   b =   0.400000000000000E+01   c =   0.133333333330000E+01
The roots are real and equal:
Root =  -0.666666666666667E+00

Coefficients are: a =   0.300000000000000E+01   b =   0.200000000000000E+01   c =  -0.100000000000000E+01
The roots are real:
Root1 =  -0.100000000000000E+01  Root2 =   0.333333333333333E+00

Coefficients are: a =   0.300000000000000E+01   b =   0.200000000000000E+01   c =   0.100000000000000E+01
The roots are complex:
Root1 =  -0.333333333333333E+00  0.471404512723287E+00j   Root2 =  -0.333333333333333E+00 -0.471404512723287E+00j

Coefficients are: a =   0.100000000000000E+01   b =  -0.100000000000000E+07   c =   0.100000000000000E+01
The roots are real:
Root1 =   0.999999999999000E+06  Root2 =   0.100000000000100E-05

J

Solution use J's built-in polynomial solver:

p.

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

<lang j> coeff =. _3 |.\ 3 4 4r3 3 2 _1 3 2 1 1 _1e6 1

  > {:"1 p. coeff
        _0.666667           _0.666667
               _1            0.333333

_0.333333j0.471405 _0.333333j_0.471405

              1e6                1e_6

</lang>

Of course this generalizes to polynomials of any order, for example the roots of $7+6x-2x^{2}-x^{3}$ are: <lang j>

  > {:"1 p. 7 6 _2 _1

_3.1925824 2.1925824 _1 </lang>

Logo

<lang logo> to quadratic :a :b :c

 localmake "d sqrt (:b*:b - 4*:a*:c)
 if :b < 0 [make "d minus :d]
 output list (:d-:b)/(2*:a) (2*:c)/(:d-:b)

end </lang>

Mathematica

Possible ways to do this are (symbolic and numeric examples): <lang Mathematica>

Solve[a x^2 + b x + c == 0, x]
Solve[x^2 - 10^5 x + 1 == 0, x]
Root[#1^2 - 10^5 #1 + 1 &, 1]
Root[#1^2 - 10^5 #1 + 1 &, 2]
Reduce[a x^2 + b x + c == 0, x]
Reduce[x^2 - 10^5 x + 1 == 0, x]
FindInstance[x^2 - 10^5 x + 1 == 0, x, Reals, 2]
FindRoot[x^2 - 10^5 x + 1 == 0, {x, 0}]
FindRoot[x^2 - 10^5 x + 1 == 0, {x, 10^6}]

</lang> gives back:

$\left\{\left\{x\to {\frac {-b-{\sqrt {b^{2}-4ac}}}{2a}}\right\},\left\{x\to {\frac {-b+{\sqrt {b^{2}-4ac}}}{2a}}\right\}\right\}$

$\left\{\left\{x\to {\frac {1}{50000+{\sqrt {2499999999}}}}\right\},\left\{x\to 50000+{\sqrt {2499999999}}\right\}\right\}$

$50000-{\sqrt {2499999999}}$

$50000+{\sqrt {2499999999}}$

$\left(a\neq 0\And \And \left(x=={\frac {-b-{\sqrt {b^{2}-4ac}}}{2a}}\|x=={\frac {-b+{\sqrt {b^{2}-4ac}}}{2a}}\right)\right)\left\|\left(a==0\And \And b\neq 0\And \And x==-{\frac {c}{b}}\right)\right\|(c==0\And \And b==0\And \And a==0)$

$x=={\frac {1}{50000+{\sqrt {2499999999}}}}\|x==50000+{\sqrt {2499999999}}$

$\left\{\left\{x\to {\frac {1}{50000+{\sqrt {2499999999}}}}\right\},\left\{x\to 50000+{\sqrt {2499999999}}\right\}\right\}$

$\{x\to 0.00001\}$

$\{x\to 100000.\}$

Note that some functions do not really give the answer (like reduce) rather it gives another way of writing it (boolean expression). However note that reduce gives the explicit cases for a zero and nonzero, b zero and nonzero, et cetera. Some functions are numeric by nature, other can handle both symbolic and numeric. In generals the solution will be exact if the input is exact. Any exact result can be approximated to arbitrary precision using the function N[expression,number of digits]. Further notice that some functions give back exact answers in a different form then others, however the answers are both correct, the answers are just written differently.

MATLAB

<lang Matlab> roots([1 -3 2]) % coefficients in decreasing order </lang>

Modula-3

Translation of: Ada

<lang modula3>MODULE Quad EXPORTS Main;

IMPORT IO, Fmt, Math;

TYPE Roots = ARRAY [1..2] OF LONGREAL;

VAR r: Roots;

PROCEDURE Solve(a, b, c: LONGREAL): Roots =

 VAR sd: LONGREAL := Math.sqrt(b * b - 4.0D0 * a * c);
     x: LONGREAL;
 BEGIN
   IF b < 0.0D0 THEN
     x := (-b + sd) / 2.0D0 * a;
     RETURN Roots{x, c / (a * x)};
   ELSE
     x := (-b - sd) / 2.0D0 * a;
     RETURN Roots{c / (a * x), x};
   END;
 END Solve;

BEGIN

 r := Solve(1.0D0, -10.0D5, 1.0D0);
 IO.Put("X1 = " & Fmt.LongReal(r[1]) & " X2 = " & Fmt.LongReal(r[2]) & "\n");

END Quad.</lang>

Octave

See MATLAB.

Python

<lang 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

>>> </lang>

R

Translation of: Python

<lang R>quaddiscrroots <- function(a,b,c, tol=1e-5) {

 d <- b*b - 4*a*c + 0i
 root1 <- (-b + sqrt(d))/(2*a)
 root2 <- (-b - sqrt(d))/(2*a)
 if ( abs(Re(d)) < tol ) {
   list("real and equal", abs(root1), abs(root1))
 } else if ( Re(d) > 0 ) {
   list("real", Re(root1), Re(root2))
 } else {
   list("complex", root1, root2)
 }

}

for(coeffs in list(c(3,4,4/3), c(3,2,-1), c(3,2,1), c(1, -1e6, 1)) ) {

 cat(sprintf("roots of %gx^2 %+gx^1 %+g are\n", coeffs[1], coeffs[2], coeffs[3]))
 r <- quaddiscrroots(coeffs[1], coeffs[2], coeffs[3])
 cat(sprintf("  %s: %s, %s\n", r1, r2, r3))

}</lang>

Ruby

With the

Library: complex.rb

package, the Math#sqrt method will return a Complex instance if necessary.

<lang ruby>require 'complex'

def quadratic(a, b, c)

 sqrt_discriminant = Math.sqrt(b**2 - 4*a*c)
 [(-b + sqrt_discriminant) / (2.0*a), (-b - sqrt_discriminant) / (2.0*a)]

end

p quadratic(3, 4, 4/3.0) # {-2/3} p quadratic(3, 2, -1) # {1/3, -1} p quadratic(3, 2, 1) # {(-1/3 + sqrt(2/9)i), (-1/3 - sqrt(2/9)i)} p quadratic(1, 0, 1) # {(0+i), (0-i)} p quadratic(1, -1e6, 1) # {1e6, 1e-6} p quadratic(-2, 7, 15) # {5, -3/2} p quadratic(1, -2, 1) # {1} p quadratic(1, 3, 3) # {(-3 - sqrt(3)i)/2), (-3 + sqrt(3)i)/2)}</lang>

[-0.666666666666667, -0.666666666666667]
[0.333333333333333, -1.0]
[Complex(-0.333333333333333, 0.471404520791032), Complex(-0.333333333333333, -0.471404520791032)]
[Complex(0.0, 1.0), Complex(0.0, -1.0)]
[999999.999999, 1.00000761449337e-06]
[-1.5, 5.0]
[1.0, 1.0]
[Complex(-1.5, 0.866025403784439), Complex(-1.5, -0.866025403784439)]

Tcl

Uses package math::complexnumbers from

Library: tcllib

<lang tcl>package require math::complexnumbers namespace import math::complexnumbers::complex math::complexnumbers::tostring

proc quadratic {a b c} {

   set discrim [expr {$b**2 - 4*$a*$c}]
   set roots [list]
   if {$discrim < 0} {
       set term1 [expr {(-1.0*$b)/(2*$a)}]
       set term2 [expr {sqrt(abs($discrim))/(2*$a)}]
       lappend roots [tostring [complex $term1 $term2]] \
               [tostring [complex $term1 [expr {-1 * $term2}]]]
   } elseif {$discrim == 0} {
       lappend roots [expr {-1.0*$b / (2*$a)}]
   } else {
       lappend roots [expr {(-1*$b + sqrt($discrim)) / (2 * $a)}] \
               [expr {(-1*$b - sqrt($discrim)) / (2 * $a)}]
   }
   return $roots

}

proc report_quad {a b c} {

   puts [format "%sx**2 + %sx + %s = 0" $a $b $c]
   foreach root [quadratic $a $b $c] {
       puts "    x = $root"
   }

}

examples on this page

report_quad 3 4 [expr {4/3.0}] ;# {-2/3} report_quad 3 2 -1 ;# {1/3, -1} report_quad 3 2 1 ;# {(-1/3 + sqrt(2/9)i), (-1/3 - sqrt(2/9)i)} report_quad 1 0 1 ;# {(0+i), (0-i)} report_quad 1 -1e6 1 ;# {1e6, 1e-6}

examples from http://en.wikipedia.org/wiki/Quadratic_equation

report_quad -2 7 15 ;# {5, -3/2} report_quad 1 -2 1 ;# {1} report_quad 1 3 3 ;# {(-3 - sqrt(3)i)/2), (-3 + sqrt(3)i)/2)}</lang> outputs

3x**2 + 4x + 1.3333333333333333 = 0
    x = -0.6666666666666666
3x**2 + 2x + -1 = 0
    x = 0.3333333333333333
    x = -1.0
3x**2 + 2x + 1 = 0
    x = -0.3333333333333333+0.47140452079103173i
    x = -0.3333333333333333-0.47140452079103173i
1x**2 + 0x + 1 = 0
    x = i
    x = -i
1x**2 + -1e6x + 1 = 0
    x = 999999.999999
    x = 1.00000761449337e-6
-2x**2 + 7x + 15 = 0
    x = -1.5
    x = 5.0
1x**2 + -2x + 1 = 0
    x = 1.0
1x**2 + 3x + 3 = 0
    x = -1.5+0.8660254037844386i
    x = -1.5-0.8660254037844386i

TI-89 BASIC

TI-89 BASIC has built-in numeric and algebraic solvers. solve(x^2-1E9x+1.0) returns x=1.E-9 or x=1.E9.