Roots of a quadratic function: Difference between revisions

From Rosetta Code
Content added Content deleted
(→‎[[Quadratic Equation#ALGOL 68]]: return either an array of REAL, or an array of COMPL)
Line 62: Line 62:
{{works with|ALGOL 68|Standard - no extensions to language used}}
{{works with|ALGOL 68|Standard - no extensions to language used}}
{{works with|ALGOL 68G|Any - tested with release mk15-0.8b.fc9.i386}}
{{works with|ALGOL 68G|Any - tested with release mk15-0.8b.fc9.i386}}
{{works with|ELLA ALGOL 68|Any (with appropriate job cards) - tested with release 1.8.8d.fc9.i386}}
<!-- {{works with|ELLA ALGOL 68|Any (with appropriate job cards) - tested with release 1.8.8d.fc9.i386}} probably need to "USE" the compl ENVIRON -->
<lang algol>quadratic equation:
<lang algol>quadratic equation:
BEGIN
BEGIN
MODE ROOTS = [2]REAL;
MODE ROOTS = UNION([]REAL, []COMPL);
PROC solve = (REAL a, b, c)ROOTS:
PROC solve = (REAL a, b, c)ROOTS:
BEGIN
BEGIN
REAL sd := sqrt (b UP 2 - 4 * a * c);
REAL sa = b UP 2 - 4 * a * c;
REAL x;
IF sa >=0 THEN # handle the +ve case as REAL #
IF b < 0 THEN
REAL sd = ( b < 0 | sqrt (sa) | -sqrt(sa));
x := (- b + sd) / 2 * a;
REAL x = (- b + sd) / 2 * a;
(x, c / (a * x))
[]REAL (x, c / (a * x))
ELSE # handle the -ve case as COMPL conjugate pairs #
⚫
ELSE
x := (- b - sd) / 2 * a;
COMPL sd = ( b < 0 | complex sqrt (sa) | -complex sqrt(sa));
(c / (a * x), x)
COMPL x = (- b + sd) / 2 * a;
[]COMPL (x, c / (a * x))
⚫
FI
⚫
FI
END # solve #;
END # solve #;
# only a very tiny difference between the 2 examples #
ROOTS r := solve (1.0, -10e5, 1);
[]ROOTS test = (solve (1, -10e5, 1), solve (1, 0, 1));


FOR index TO UPB test DO
⚫
put(stand out,("x1 =", fixed(r[1], -real width, 8), " x2 =", fixed(r[2], -real width, 8), new line))
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>
END # quadratic_equation #</lang>
Output:
Output:
<pre>
<pre>
x1 =999999.99999900 x2 = 0.00000100
REAL x1 = 999999.99999900, x2 = 0.00000100
COMPL x1,x2 = 0.00000000 +/- 1.00000000i
</pre>
</pre>

=={{header|Fortran}}==
=={{header|Fortran}}==
{{works with|Fortran|90 and later}}
{{works with|Fortran|90 and later}}

Revision as of 11:39, 16 February 2009

Task
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 ax2 + 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=-105, 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 x2 is about 10-6. The naive method is numerically unstable.

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];
  1. 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

Fortran

Works with: Fortran version 90 and later
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

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: 

   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.

Python

No attempt is made to categorize the type of output as this is not mentioned as being part of the task. <lang 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))

>>> </lang>

Categorized roots version: <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>

Retrieved from "https://rosettacode.org/wiki/Roots_of_a_quadratic_function?oldid=89219"