Roots of a quadratic function: Difference between revisions
m (Corrected order of examples) |
No edit summary |
||
Line 2: | Line 2: | ||
The program must correctly handle complex roots. Error checking on the input (for a = 0) need not be shown. |
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 [http://www.pdas.com/fmm.htm Computer Methods for Mathematical Computations], George Forsythe, Michael Malcolm and Cleve Mole suggest to try it on ''a''=1, ''b''=-10<sup>5</sup>, ''c''=1. Consider the following implementation in [[Ada]]: |
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 [http://www.pdas.com/fmm.htm Computer Methods for Mathematical Computations], George Forsythe, Michael Malcolm and Cleve Mole suggest to try it on ''a''=1, ''b''=-10<sup>5</sup>, ''c''=1. Consider the following implementation in [[Ada]]: |
||
<ada> |
<lang ada> |
||
with Ada.Text_IO; use Ada.Text_IO; |
with Ada.Text_IO; use Ada.Text_IO; |
||
with Ada.Numerics.Elementary_Functions; use Ada.Numerics.Elementary_Functions; |
with Ada.Numerics.Elementary_Functions; use Ada.Numerics.Elementary_Functions; |
||
Line 19: | Line 19: | ||
Put_Line ("X1 =" & Float'Image (R (1)) & " X2 =" & Float'Image (R (2))); |
Put_Line ("X1 =" & Float'Image (R (1)) & " X2 =" & Float'Image (R (2))); |
||
end Quadratic_Equation; |
end Quadratic_Equation; |
||
</ |
</lang> |
||
Sample output: |
Sample output: |
||
<pre> |
<pre> |
||
Line 28: | Line 28: | ||
'''Task''': do it better. |
'''Task''': do it better. |
||
=={{header|Ada}}== |
=={{header|Ada}}== |
||
<ada> |
<lang ada> |
||
with Ada.Text_IO; use Ada.Text_IO; |
with Ada.Text_IO; use Ada.Text_IO; |
||
with Ada.Numerics.Elementary_Functions; use Ada.Numerics.Elementary_Functions; |
with Ada.Numerics.Elementary_Functions; use Ada.Numerics.Elementary_Functions; |
||
Line 51: | Line 51: | ||
Put_Line ("X1 =" & Float'Image (R (1)) & " X2 =" & Float'Image (R (2))); |
Put_Line ("X1 =" & Float'Image (R (1)) & " X2 =" & Float'Image (R (2))); |
||
end Quadratic_Equation; |
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: |
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: |
||
<pre> |
<pre> |
||
Line 131: | Line 131: | ||
=={{header|Python}}== |
=={{header|Python}}== |
||
No attempt is made to categorize the type of output as this is not mentioned as being part of the task. |
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): |
<lang python>>>> def quad_roots(a,b,c): |
||
a,b,c =complex(a), complex(b), complex(c) |
a,b,c =complex(a), complex(b), complex(c) |
||
discriminant = b*b - 4*a*c |
discriminant = b*b - 4*a*c |
||
Line 149: | Line 149: | ||
a, b, c = (3, 2, 1) |
a, b, c = (3, 2, 1) |
||
(root1, root2) = ((-0.33333333333333331+0.47140452079103173j), (-0.33333333333333331-0.47140452079103173j)) |
(root1, root2) = ((-0.33333333333333331+0.47140452079103173j), (-0.33333333333333331-0.47140452079103173j)) |
||
>>> </ |
>>> </lang> |
||
Categorized roots version: |
Categorized roots version: |
||
<python>>>> def quad_discriminating_roots(a,b,c, entier = 1e-5): |
<lang python>>>> def quad_discriminating_roots(a,b,c, entier = 1e-5): |
||
discriminant = b*b - 4*a*c |
discriminant = b*b - 4*a*c |
||
a,b,c,d =complex(a), complex(b), complex(c), complex(discriminant) |
a,b,c,d =complex(a), complex(b), complex(c), complex(discriminant) |
||
Line 176: | Line 176: | ||
Roots of: 1X^2 -1e+06X +1 are |
Roots of: 1X^2 -1e+06X +1 are |
||
real: 999999.999999, 1.00000761449e-06 |
real: 999999.999999, 1.00000761449e-06 |
||
>>> </ |
>>> </lang> |
Revision as of 15:46, 3 February 2009
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
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
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>