Roots of a quadratic function: Difference between revisions
(adding scheme) |
(→{{header|Scheme}}: task example) |
||
| Line 870: | Line 870: | ||
(quadratic 1 2 1) |
(quadratic 1 2 1) |
||
; (-1 -1) |
; (-1 -1) |
||
(quadratic 1 -1e5 1) |
|||
; (99999.99999 1.0000000001000001e-05)</lang> |
|||
=={{header|Tcl}}== |
=={{header|Tcl}}== |
||
Revision as of 22:12, 1 January 2011
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 represent 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 when one of the roots is much closer to zero than the other. 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=-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.
(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 = -109 = -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
<lang algol68>quadratic equation: BEGIN
MODE ROOTS = UNION([]REAL, []COMPL); MODE QUADRATIC = STRUCT(REAL a,b,c);
PROC solve = (QUADRATIC q)ROOTS:
BEGIN
REAL a = a OF q, b = b OF q, c = c OF q;
REAL sa = b**2 - 4*a*c;
IF sa >=0 THEN # handle the +ve case as REAL #
REAL sqrt sa = ( b<0 | sqrt(sa) | -sqrt(sa));
REAL r1 = (-b + sqrt sa)/(2*a),
r2 = (-b - sqrt sa)/(2*a);
[]REAL((r1,r2))
ELSE # handle the -ve case as COMPL conjugate pairs #
COMPL compl sqrt sa = ( b<0 | complex sqrt(sa) | -complex sqrt(sa));
COMPL r1 = (-b + compl sqrt sa)/(2*a),
r2 = (-b - compl sqrt sa)/(2*a);
[]COMPL (r1, r2)
FI
END # solve #;
PROC real evaluate = (QUADRATIC q, REAL x )REAL: (a OF q*x + b OF q)*x + c OF q;
PROC compl evaluate = (QUADRATIC q, COMPL x)COMPL: (a OF q*x + b OF q)*x + c OF q;
# only a very tiny difference between the 2 examples # []QUADRATIC test = ((1, -10e5, 1), (1, 0, 1), (1,-3,2), (1,3,2), (4,0,4), (3,4,5)); FORMAT real fmt = $g(-0,8)$; FORMAT compl fmt = $f(real fmt)"+"f(real fmt)"i"$; FORMAT quadratic fmt = $f(real fmt)" x**2 + "f(real fmt)" x + "f(real fmt)" = 0"$;
FOR index TO UPB test DO QUADRATIC quadratic = test[index]; ROOTS r = solve(quadratic);
- Output the two different scenerios #
printf(($"Quadratic: "$, quadratic fmt, quadratic, $l$));
CASE r IN
([]REAL r):
printf(($"REAL x1 = "$, real fmt, r[1],
$", x2 = "$, real fmt, r[2], $"; "$,
$"REAL y1 = "$, real fmt, real evaluate(quadratic,r[1]),
$", y2 = "$, real fmt, real evaluate(quadratic,r[2]), $";"ll$
)),
([]COMPL c):
printf(($"COMPL x1,x2 = "$, real fmt, re OF c[1], $"+/-"$,
real fmt, ABS im OF c[1], $"; "$,
$"COMPL y1 = "$, compl fmt, compl evaluate(quadratic,c[1]),
$", y2 = "$, compl fmt, compl evaluate(quadratic,c[2]), $";"ll$
))
ESAC
OD
END # quadratic_equation #</lang> Output:
Quadratic: 1.00000000 x**2 + -1000000.00000000 x + 1.00000000 = 0 REAL x1 = 999999.99999900, x2 = .00000100; REAL y1 = -.00000761, y2 = -.00000761; Quadratic: 1.00000000 x**2 + .00000000 x + 1.00000000 = 0 COMPL x1,x2 = .00000000+/-1.00000000; COMPL y1 = .00000000+.00000000i, y2 = .00000000+.00000000i; Quadratic: 1.00000000 x**2 + -3.00000000 x + 2.00000000 = 0 REAL x1 = 2.00000000, x2 = 1.00000000; REAL y1 = .00000000, y2 = .00000000; Quadratic: 1.00000000 x**2 + 3.00000000 x + 2.00000000 = 0 REAL x1 = -2.00000000, x2 = -1.00000000; REAL y1 = .00000000, y2 = .00000000; Quadratic: 4.00000000 x**2 + .00000000 x + 4.00000000 = 0 COMPL x1,x2 = .00000000+/-1.00000000; COMPL y1 = .00000000+.00000000i, y2 = .00000000+.00000000i; Quadratic: 3.00000000 x**2 + 4.00000000 x + 5.00000000 = 0 COMPL x1,x2 = -.66666667+/-1.10554160; COMPL y1 = .00000000+.00000000i, y2 = .00000000+-.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>
<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++
<lang cpp>#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)
Common Lisp
<lang lisp>(defun quadratic (a b c)
"Compute the roots of a quadratic in the form ax^2 + bx + c = 1. Evaluates to a list of the two roots."
(let ((discriminant (- (expt b 2) (* 4 a c)))
(denominator (* 2 a))
(neg-b (* b -1)))
(list (/ (+ neg-b (sqrt discriminant)) denominator)
(/ (- neg-b (sqrt discriminant)) denominator))))</lang>
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>
Factor
<lang factor>:: quadratic-equation ( a b c -- x1 x2 )
b sq a c * 4 * - sqrt :> sd b 0 < [ b neg sd + a 2 * / ] [ b neg sd - a 2 * / ] if :> x x c a x * / ;</lang>
<lang factor>( scratchpad ) 1 -1.e20 1 quadratic-equation --- Data stack: 1.0e+20 9.999999999999999e-21</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
<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
GAP
<lang gap>QuadraticRoots := function(a, b, c)
local d; d := Sqrt(b*b - 4*a*c); return [ (-b+d)/(2*a), (-b-d)/(2*a) ];
end;
- Hint : E(12) is a 12th primitive root of 1
QuadraticRoots(2, 2, -1);
- [ 1/2*E(12)^4-1/2*E(12)^7+1/2*E(12)^8+1/2*E(12)^11,
- 1/2*E(12)^4+1/2*E(12)^7+1/2*E(12)^8-1/2*E(12)^11 ]</lang>
Haskell
<lang haskell>import Data.Complex
solveQuadrEquation ::
[Complex Double] -> (Complex Double, Complex Double)
solveQuadrEquation [a, b, c] = (x1,x2) where
discr = sqrt $ b^2 - 4*a*c x1 = ((-b) + discr) / (2*a) x2 = ((-b) - discr) / (2*a)</lang>
Examples: <lang haskell>*Main> mapM_ print $ map solveQuadrEquation [[3 ,4, 4/3],[3,2,-1],[3,2,1],[1,-10^5,1]] ((-0.6666666666666666) :+ 0.0,(-0.6666666666666666) :+ 0.0) (0.3333333333333333 :+ 0.0,(-1.0) :+ 0.0) ((-0.3333333333333333) :+ 0.47140452079103173,(-0.3333333333333333) :+ (-0.47140452079103173)) (99999.99999 :+ 0.0,1.0000003385357559e-5 :+ 0.0)</lang>
IDL
<lang idl>compile_OPT IDL2
print, "input a, press enter, input b, press enter, input c, press enter" read,a,b,c Promt='Enter values of a,b,c and hit enter'
a0=0.0 b0=0.0 c0=0.0 ;make them floating point variables
x=-b+sqrt((b^2)-4*a*c) y=-b-sqrt((b^2)-4*a*c) z=2*a d= x/z e= y/z
print, d,e</lang>
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 p. generalizes to polynomials of any order. Given the coefficients p. returns the multiplier and roots of the polynomial. Given the multiplier and roots it returns the coefficients. For example using the cubic :
<lang j> p. 0 16 _12 2 NB. return multiplier ; roots
+-+-----+
|2|4 2 0|
+-+-----+
p. 2 ; 4 2 0 NB. return coefficients
0 16 _12 2</lang>
Java
<lang java>public class Main {
static class Complex {
double re, im;
public Complex(double re, double im) {
this.re = re;
this.im = im;
}
@Override
public boolean equals(Object obj) {
if (obj == this) {return true;}
if (!(obj instanceof Complex)) {return false;}
Complex other = (Complex) obj;
return (re == other.re) && (im == other.im);
}
@Override
public String toString() {
if (im == 0.0) {return String.format("%g", re);}
if (re == 0.0) {return String.format("%gi", im);}
return String.format("%g %c %gi", re,
(im < 0.0 ? '-' : '+'), Math.abs(im));
}
}
private static Complex[] quadraticRoots(double a, double b, double c) {
Complex[] roots = new Complex[2];
double d = b * b - 4.0 * a * c;
double aa = a + a;
if (d < 0.0) {
double re = -b / aa;
double im = Math.sqrt(-d) / aa;
roots[0] = new Complex(re, im);
roots[1] = new Complex(re, -im);
} else {
roots[0] = new Complex((-b + Math.sqrt(d)) / aa, 0.0);
roots[1] = new Complex((-b - Math.sqrt(d)) / aa, 0.0);
}
return roots;
}
public static void main(String[] args) {
double[][] equations = {
{1.0, -10.0e5, 1.0}, // two distinct real roots
{1.0, 2.0, 1.0}, // one real root (double root)
{1.0, 0.0, 1.0}, // two imaginary roots
{1.0, 1.0, 1.0} // two complex roots
};
for (int i = 0; i < equations.length; i++) {
Complex[] roots = quadraticRoots(
equations[i][0], equations[i][1], equations[i][2]);
System.out.format("%na = %g b = %g c = %g%n",
equations[i][0], equations[i][1], equations[i][2]);
if (roots[0].equals(roots[1])) {
System.out.format("X1,2 = %s%n", roots[0]);
} else {
System.out.format("X1 = %s%n", roots[0]);
System.out.format("X2 = %s%n", roots[1]);
}
}
}
}</lang> Output:
a = 1.00000 b = -1.00000e+06 c = 1.00000 X1 = 1.00000e+06 X2 = 1.00001e-06 a = 1.00000 b = 2.00000 c = 1.00000 X1,2 = -1.00000 a = 1.00000 b = 0.00000 c = 1.00000 X1 = 1.00000i X2 = -1.00000i a = 1.00000 b = 1.00000 c = 1.00000 X1 = -0.500000 + 0.866025i X2 = -0.500000 - 0.866025i
Liberty BASIC
<lang lb>a=1:b=2:c=3
'assume a<>0 print quad$(a,b,c) end
function quad$(a,b,c)
D=b^2-4*a*c
x=-1*b
if D<0 then
quad$=str$(x/(2*a));" +i";str$(sqr(abs(D))/(2*a));" , ";str$(x/(2*a));" -i";str$(sqr(abs(D))/abs(2*a))
else
quad$=str$(x/(2*a)+sqr(D)/(2*a));" , ";str$(x/(2*a)-sqr(D)/(2*a))
end if
end function</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>
Lua
In order to correctly handle complex roots, qsolve must be given objects from a suitable complex number library, like that from the Complex Numbers article. However, this should be enough to demonstrate its accuracy:
<lang lua>function qsolve(a, b, c)
if b < 0 then return qsolve(-a, -b, -c) end val = b + (b^2 - 4*a*c)^(1/2) --this never exhibits instability if b > 0 return -val / (2 * a), -2 * c / val --2c / val is the same as the "unstable" second root
end
for i = 1, 12 do
print(qsolve(1, 0-10^i, 1))
end</lang> The "trick" lies in avoiding subtracting large values that differ by a small amount, which is the source of instability in the "normal" formula. It is trivial to prove that 2c/(b + sqrt(b^2-4ac)) = (b - sqrt(b^2-4ac))/2a.
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:
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 of power e.g. [x^n ... x^2 x^1 x^0]</lang>
Modula-3
<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>
OCaml
<lang ocaml>type quadroots =
| RealRoots of float * float | ComplexRoots of Complex.t * Complex.t ;;
let quadsolve a b c =
let d = (b *. b) -. (4.0 *. a *. c) in
if d < 0.0
then
let r = -. b /. (2.0 *. a)
and i = sqrt(-. d) /. (2.0 *. a) in
ComplexRoots ({ Complex.re = r; Complex.im = i },
{ Complex.re = r; Complex.im = (-.i) })
else
let r = if b < 0.0
then ((sqrt d) -. b) /. (2.0 *. a)
else ((sqrt d) +. b) /. (-2.0 *. a)
in
RealRoots (r, c /. (r *. a))
- </lang>
Sample output: <lang ocaml># quadsolve 1.0 0.0 (-2.0) ;; - : quadroots = RealRoots (-1.4142135623730951, 1.4142135623730949)
- quadsolve 1.0 0.0 2.0 ;;
- : quadroots = ComplexRoots ({Complex.re = 0.; Complex.im = 1.4142135623730951},
{Complex.re = 0.; Complex.im = -1.4142135623730951})
- quadsolve 1.0 (-1.0e5) 1.0 ;;
- : quadroots = RealRoots (99999.99999, 1.0000000001000001e-005)</lang>
Octave
See MATLAB.
PARI/GP
<lang>roots(a,b,c)={
b /= a; c /= a; my (delta = b^2 - 4*c, root=sqrt(delta)); if (delta < 0, [root-b,-root-b]/2 , my(sol=if(b>0, -b - root,-b + root)/2); [sol,c/sol] )
};</lang>
Perl
When using Math::Complex perl automatically convert numbers when necessary. <lang perl>use Math::Complex;
($x1,$x2) = solveQuad(1,2,3);
print "x1 = $x1, x2 = $x2\n";
sub solveQuad { my ($a,$b,$c) = @_; my $root = sqrt($b**2 - 4*$a*$c); return ( -$b + $root )/(2*$a), ( -$b - $root )/(2*$a); }</lang>
Perl 6
Perl 6 has complex number handling built in.
<lang perl6>my @sets = [1, 2, 1],
[1, 2, 3],
[1, -2, 1],
[1, 0, -4],
[1, -10**6, 1];
for @sets -> @coefficients {
say "Roots for @coefficients.join(', ').fmt("%-16s")",
"=> (&quadroots( @coefficients ).join(', '))";
}
multi sub quadroots ($a, $b, $c) {
my $root = (my $t = $b ** 2 - 4 * $a * $c ) < 0
?? $t.Complex.sqrt
!! $t.sqrt;
return ( -$b + $root ) / (2 * $a),
( -$b - $root ) / (2 * $a);
}
multi sub quadroots (@a) {
@a == 3 or die "Expected three elements, got {+@a}";
quadroots |@a;
}</lang> Output:
Roots for 1, 2, 1 => (-1, -1) Roots for 1, 2, 3 => (-1 + 1.4142135623731i, -1 + -1.4142135623731i) Roots for 1, -2, 1 => (1, 1) Roots for 1, 0, -4 => (2, -2) Roots for 1, -1000000, 1 => (999999.999999, 1.00000761449337e-06)
PL/I
<lang PL/I>
declare (c1, c2) float complex,
(a, b, c, x1, x2) float;
get list (a, b, c);
if b**2 < 4*a*c then
do;
c1 = (-b + sqrt(b**2 - 4+0i*a*c)) / (2*a);
c2 = (-b - sqrt(b**2 - 4+0i*a*c)) / (2*a);
put data (c1, c2);
end;
else
do;
x1 = (-b + sqrt(b**2 - 4*a*c)) / (2*a);
x2 = (-b - sqrt(b**2 - 4*a*c)) / (2*a);
put data (x1, x2);
end;
</lang>
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
<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
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)]
Scheme
<lang scheme>(define (quadratic a b c) (if (= a 0) (if (= b 0) 'fail (- (/ c b))) (let ((delta (- (* b b) (* 4 a c)))) (if (and (real? delta) (> delta 0)) (let ((u (+ b (* (if (>= b 0) 1 -1) (sqrt delta))))) (list (/ u -2 a) (/ (* -2 c) u))) (list (/ (- (sqrt delta) b) 2 a) (/ (- (+ (sqrt delta) b)) 2 a))))))
- examples
(quadratic 1 -1 -1)
- (1.618033988749895 -0.6180339887498948)
(quadratic 1 0 -2)
- (-1.4142135623730951 1.414213562373095)
(quadratic 1 0 2)
- (0+1.4142135623730951i 0-1.4142135623730951i)
(quadratic 1+1i 2 5)
- (-1.0922677260818898-1.1884256155834088i 0.09226772608188982+2.1884256155834088i)
(quadratic 0 4 3)
- -3/4
(quadratic 0 0 1)
- fail
(quadratic 1 2 0)
- (-2 0)
(quadratic 1 2 1)
- (-1 -1)
(quadratic 1 -1e5 1)
- (99999.99999 1.0000000001000001e-05)</lang>
Tcl
<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.