Roots of a quadratic function

From Rosetta Code
Jump to: navigation, search
Task
Roots of a quadratic function
You are encouraged to solve this task according to the task description, using any language you may know.
This task has been clarified. Its programming examples are in need of review to ensure that they still fit the requirements of the task.

Write a program to find the roots of a quadratic equation, i.e., solve the equation ax2 + bx + c = 0. Your program must correctly handle non-real roots, but it need not check that a \neq 0.

The problem of solving a quadratic equation is a good example of how dangerous it can be to ignore the peculiarities of floating-point arithmetic. The obvious way to implement the quadratic formula suffers catastrophic loss of accuracy when one of the roots to be found is much closer to 0 than the other. In their classic textbook on numeric methods Computer Methods for Mathematical Computations, George Forsythe, Michael Malcolm, and Cleve Moler suggest trying the naive algorithm with a = 1, b = − 105, and c = 1. (For double-precision floats, set b = − 109.) Consider the following implementation in 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;

Sample output:

X1 = 1.00000E+06 X2 = 0.00000E+00

As we can see, the second root has lost all significant figures. The right answer is that X2 is about 10 − 6. The naive method is numerically unstable.

Suggested by Middlebrook (D-OA), a better numerical method: to define two parameters  q = \sqrt{a c} / b and  f = 1/2 + \sqrt{1 - 4 q^2} /2

and the two roots of the quardratic are:  \frac{-b}{a} f and  \frac{-c}{b f}


Task: do it better. This means that given a = 1, b = − 109, and c = 1, both of the roots your program returns should be greater than 10 − 11. Or, if your language can't do floating-point arithmetic any more precisely than single precision, your program should be able to handle b = − 106. Either way, show what your program gives as the roots of the quadratic in question. See page 9 of "What Every Scientist Should Know About Floating-Point Arithmetic" for a possible algorithm.

Contents

[edit] 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;

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

[edit] ALGOL 68

Translation of: Ada
Works with: ALGOL 68 version Revision 1 - no extensions to language used
Works with: ALGOL 68G version Any - tested with release 1.18.0-9h.tiny
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 #

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;

[edit] AutoHotkey

ahk forum: discussion

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
}

[edit] BBC BASIC

      FOR test% = 1 TO 7
READ a$, b$, c$
PRINT "For a = " ; a$ ", b = " ; b$ ", c = " ; c$ TAB(32) ;
PROCsolvequadratic(EVAL(a$), EVAL(b$), EVAL(c$))
NEXT
END
 
DATA 1, -1E9, 1
DATA 1, 0, 1
DATA 2, -1, -6
DATA 1, 2, -2
DATA 0.5, SQR(2), 1
DATA 1, 3, 2
DATA 3, 4, 5
 
DEF PROCsolvequadratic(a, b, c)
LOCAL d, f
d = b^2 - 4*a*c
CASE SGN(d) OF
WHEN 0:
PRINT "the single root is " ; -b/2/a
WHEN +1:
f = (1 + SQR(1-4*a*c/b^2))/2
PRINT "the real roots are " ; -f*b/a " and " ; -c/b/f
WHEN -1:
PRINT "the complex roots are " ; -b/2/a " +/- " ; SQR(-d)/2/a "*i"
ENDCASE
ENDPROC

Output:

For a = 1, b = -1E9, c = 1      the real roots are 1E9 and 1E-9
For a = 1, b = 0, c = 1         the complex roots are 0 +/- 1*i
For a = 2, b = -1, c = -6       the real roots are 2 and -1.5
For a = 1, b = 2, c = -2        the real roots are -2.73205081 and 0.732050808
For a = 0.5, b = SQR(2), c = 1  the single root is -1.41421356
For a = 1, b = 3, c = 2         the real roots are -2 and -1
For a = 3, b = 4, c = 5         the complex roots are -0.666666667 +/- 1.1055416*i

[edit] C

Code that tries to avoid floating point overflow and other unfortunate loss of precissions: (compiled with gcc -std=c99 for complex, though easily adapted to just real numbers)

#include <stdio.h>
#include <stdlib.h>
#include <complex.h>
#include <math.h>
 
typedef double complex cplx;
 
void quad_root
(double a, double b, double c, cplx * ra, cplx *rb)
{
double d, e;
if (!a) {
*ra = b ? -c / b : 0;
*rb = 0;
return;
}
if (!c) {
*ra = 0;
*rb = -b / a;
return;
}
 
b /= 2;
if (fabs(b) > fabs(c)) {
e = 1 - (a / b) * (c / b);
d = sqrt(fabs(e)) * fabs(b);
} else {
e = (c > 0) ? a : -a;
e = b * (b / fabs(c)) - e;
d = sqrt(fabs(e)) * sqrt(fabs(c));
}
 
if (e < 0) {
e = fabs(d / a);
d = -b / a;
*ra = d + I * e;
*rb = d - I * e;
return;
}
 
d = (b >= 0) ? d : -d;
e = (d - b) / a;
d = e ? (c / e) / a : 0;
*ra = d;
*rb = e;
return;
}
 
int main()
{
cplx ra, rb;
quad_root(1, 1e12 + 1, 1e12, &ra, &rb);
printf("(%g + %g i), (%g + %g i)\n",
creal(ra), cimag(ra), creal(rb), cimag(rb));
 
quad_root(1e300, -1e307 + 1, 1e300, &ra, &rb);
printf("(%g + %g i), (%g + %g i)\n",
creal(ra), cimag(ra), creal(rb), cimag(rb));
 
return 0;
}
Output:
(-1e+12 + 0 i), (-1 + 0 i)
(1.00208e+07 + 0 i), (9.9792e-08 + 0 i)
#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);
}
Translation of: 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;
}
}
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;
}
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)

[edit] C#

using System;
using System.Numerics;
 
class QuadraticRoots
{
static Tuple<Complex, Complex> Solve(double a, double b, double c)
{
var q = -(b + Math.Sign(b) * Complex.Sqrt(b * b - 4 * a * c)) / 2;
return Tuple.Create(q / a, c / q);
}
 
static void Main()
{
Console.WriteLine(Solve(1, -1E20, 1));
}
}

Output:

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

[edit] 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;
}

Output:

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

[edit] Clojure

(defn quadratic 
"Compute the roots of a quadratic in the form ax^2 + bx + c = 1.
Returns any of nil, a float, or a vector."

[a b c]
(let [sq-d (Math/sqrt (- (* b b) (* 4 a c)))
f #(/ (% b sq-d) (* 2 a))]
(cond
(neg? sq-d) nil
(zero? sq-d) (f +)
(pos? sq-d) [(f +) (f -)]
 :else nil))) ; maybe our number ended up as NaN

Output:

user=> (quadratic 1.0 1.0 1.0)
nil
user=> (quadratic 1.0 2.0 1.0)
1.0
user=> (quadratic 1.0 3.0 1.0)
[2.618033988749895 0.3819660112501051]
 

[edit] Common 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))))

[edit] D

import std.math, std.traits;
 
CommonType!(T1, T2, T3)[] naiveQR(T1, T2, T3)
(in T1 a, in T2 b, in T3 c)
pure nothrow if (isFloatingPoint!T1) {
alias ReturnT = typeof(typeof(return).init[0]);
if (a == 0)
return [ReturnT(c / b)]; // It's a linear function.
immutable ReturnT det = b ^^ 2 - 4 * a * c;
if (det < 0)
return []; // No real number root.
immutable SD = sqrt(det);
return [(-b + SD) / 2 * a, (-b - SD) / 2 * a];
}
 
CommonType!(T1, T2, T3)[] cautiQR(T1, T2, T3)
(in T1 a, in T2 b, in T3 c)
pure nothrow if (isFloatingPoint!T1) {
alias ReturnT = typeof(typeof(return).init[0]);
if (a == 0)
return [ReturnT(c / b)]; // It's a linear function.
immutable ReturnT det = b ^^ 2 - 4 * a * c;
if (det < 0)
return []; // No real number root.
immutable SD = sqrt(det);
 
if (b * a < 0) {
immutable x = (-b + SD) / 2 * a;
return [x, c / (a * x)];
} else {
immutable x = (-b - SD) / 2 * a;
return [c / (a * x), x];
}
}
 
void main() {
import std.stdio;
writeln("With 32 bit float type:");
writefln(" Naive: [%(%g, %)]", naiveQR(1.0f, -10e5f, 1.0f));
writefln("Cautious: [%(%g, %)]", cautiQR(1.0f, -10e5f, 1.0f));
writeln("\nWith 64 bit double type:");
writefln(" Naive: [%(%g, %)]", naiveQR(1.0, -10e5, 1.0));
writefln("Cautious: [%(%g, %)]", cautiQR(1.0, -10e5, 1.0));
writeln("\nWith real type:");
writefln(" Naive: [%(%g, %)]", naiveQR(1.0L, -10e5L, 1.0L));
writefln("Cautious: [%(%g, %)]", cautiQR(1.0L, -10e5L, 1.0L));
}
Output:
With 32 bit float type:
   Naive: [1e+06, 0]
Cautious: [1e+06, 1e-06]

With 64 bit double type:
   Naive: [1e+06, 1.00001e-06]
Cautious: [1e+06, 1e-06]

With real type:
   Naive: [1e+06, 1e-06]
Cautious: [1e+06, 1e-06]

[edit] Factor

Translation of: Ada
:: 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 * / ;
( scratchpad ) 1 -1.e20 1 quadratic-equation
--- Data stack:
1.0e+20
9.999999999999999e-21

Middlebrook method

:: quadratic-equation2 ( a b c -- x1 x2 )
a c * sqrt b / :> q
1 4 q sq * - sqrt 0.5 * 0.5 + :> f
b neg a / f * c neg b / f / ;
 


( scratchpad ) 1 -1.e20 1 quadratic-equation
--- Data stack:
1.0e+20
1.0e-20

[edit] Forth

Without locals:

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

With locals:

: 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

[edit] 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
rroot1 = -(b + SIGN(SQRT(discriminant), b)) / (2.0_dp * a)
rroot2 = c / (a * rroot1)
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

[edit] 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 ]
 
# This works also with floating-point numbers
QuadraticRoots(2.0, 2.0, -1.0);
# [ 0.366025, -1.36603 ]

[edit] Go

package main
 
import (
"fmt"
"math"
)
 
func qr(a, b, c float64) ([]float64, []complex128) {
d := b*b-4*a*c
switch {
case d == 0:
// single root
return []float64{-b/(2*a)}, nil
case d > 0:
// two real roots
if b < 0 {
d = math.Sqrt(d)-b
} else {
d = -math.Sqrt(d)-b
}
return []float64{d/(2*a), (2*c)/d}, nil
case d < 0:
// two complex roots
 
den := 1/(2*a)
t1 := complex(-b*den, 0)
t2 := complex(0, math.Sqrt(-d)*den)
return nil, []complex128{t1+t2, t1-t2}
}
// otherwise d overflowed or a coefficient was NAN
return []float64{d}, nil
}
 
func test(a, b, c float64) {
fmt.Print("coefficients: ", a, b, c, " -> ")
r, i := qr(a, b, c)
switch len(r) {
case 1:
fmt.Println("one real root:", r[0])
case 2:
fmt.Println("two real roots:", r[0], r[1])
default:
fmt.Println("two complex roots:", i[0], i[1])
}
}
 
func main() {
for _, c := range [][3]float64{
{1, -2, 1},
{1, 0, 1},
{1, -10, 1},
{1, -1000, 1},
{1, -1e9, 1},
} {
test(c[0], c[1], c[2])
}
}

Output:

coefficients: 1 -2 1 -> one real root: 1
coefficients: 1 0 1 -> two complex roots: (0+1i) (-0-1i)
coefficients: 1 -10 1 -> two real roots: 9.898979485566356 0.10102051443364381
coefficients: 1 -1000 1 -> two real roots: 999.998999999 0.001000001000002
coefficients: 1 -1e+09 1 -> two real roots: 1e+09 1e-09

[edit] Haskell

import Data.Complex
 
type CD = Complex Double
 
quadraticRoots :: (CD, CD, CD) -> (CD, CD)
quadraticRoots (a, b, c) =
if realPart b > 0
then ((2*c) / (-b - d), (-b - d) / (2*a))
else ((-b + d) / (2*a), (2*c) / (-b + d))
where d = sqrt $ b^2 - 4*a*c
*Main> mapM_ print $ map quadraticRoots [(3, 4, 4/3), (3, 2, -1), (3, 2, 1), (1, -10e5, 1), (1, -10e9, 1)]
((-0.6666666666666666) :+ (-0.0),(-0.6666666666666666) :+ 0.0)
(0.3333333333333333 :+ 0.0,(-1.0) :+ 0.0)
((-0.33333333333333326) :+ 0.4714045207910316,(-0.3333333333333333) :+ (-0.47140452079103173))
(999999.999999 :+ 0.0,1.000000000001e-6 :+ 0.0)
(1.0e10 :+ 0.0,1.0e-10 :+ 0.0)

[edit] 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

[edit] Icon and Unicon

Translation of: Ada

Works in both languages.

procedure main()
solve(1.0, -10.0e5, 1.0)
end
 
procedure solve(a,b,c)
d := sqrt(b*b - 4.0*a*c)
roots := if b < 0 then [r1 := (-b+d)/2.0*a, c/(a*r1)]
else [r1 := (-b-d)/2.0*a, c/(a*r1)]
write(a,"*x^2 + ",b,"*x + ",c," has roots ",roots[1]," and ",roots[2])
end

Output:

->rqf 1.0 -0.000000001 1.0
1.0*x^2 + -1000000.0*x + 1.0 has roots 999999.999999 and 1.000000000001e-06
->

[edit] 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 _1e6 1
> {:"1 p. coeff
_0.666667 _0.666667
_1 0.333333
_0.333333j0.471405 _0.333333j_0.471405
1e6 1e_6

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 0 + 16x − 12x2 + 2x3:

   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

Exploring the limits of precision:

   1{::p. 1 _1e5 1  NB. display roots
100000 1e_5
1 _1e5 1 p. 1{::p. 1 _1e5 1 NB. test roots
_3.38436e_7 0
1 _1e5 1 p. 1e5 1e_5 NB. test displayed roots
1 9.99999e_11
1e5 1e_5 - 1{::p. 1 _1e5 1 NB. find difference
1e_5 _1e_15
1 _1e5 1 p. 1e5 1e_5-1e_5 _1e_15 NB. test displayed roots with adjustment
_3.38436e_7 0

Updated example:

   1 {:: p.1 _1e9 1
1e9 1e_9

As before, when these "roots" are plugged back into the original polynomial, the results are nowhere near zero. However, double precision floating point does not have enough bits to represent the (extremely close) answers that would give a zero result.

Middlebrook formula implemented in J

 
q_r =: 3 : 0
'a b c' =. y
q=. b %~ %: a * c
f=. 0.5 + 0.5*%:(1-4*q*q)
(-b*f%a),(-c%b*f)
)
 
q_r 1 _1e6 1
1e6 1e_6
 

[edit] Java

public class QuadraticRoots {
private 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; // discriminant
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 if (b < 0.0) {
// Avoid calculating -b - Math.sqrt(d), to avoid any
// subtractive cancellation when it is near zero.
double re = (-b + Math.sqrt(d)) / aa;
roots[0] = new Complex(re, 0.0);
roots[1] = new Complex(c / (a * re), 0.0);
} else {
// Avoid calculating -b + Math.sqrt(d).
double re = (-b - Math.sqrt(d)) / aa;
roots[1] = new Complex(re, 0.0);
roots[0] = new Complex(c / (a * re), 0.0);
}
return roots;
}
 
public static void main(String[] args) {
double[][] equations = {
{1.0, 22.0, -1323.0}, // two distinct real roots
{6.0, -23.0, 20.0}, // with a != 1.0
{1.0, -1.0e9, 1.0}, // with one root near zero
{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]);
}
}
}
}

Output:

a = 1.00000   b = 22.0000   c = -1323.00
X1 = 27.0000
X2 = -49.0000

a = 6.00000   b = -23.0000   c = 20.0000
X1 = 2.50000
X2 = 1.33333

a = 1.00000   b = -1.00000e+09   c = 1.00000
X1 = 1.00000e+09
X2 = 1.00000e-09

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

[edit] Liberty BASIC

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

[edit]

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

[edit] 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:

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

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.


[edit] Maple

solve(a*x^2+b*x+c,x);
 
solve(1.0*x^2-10.0^9*x+1.0,x,explicit,allsolutions);
 
fsolve(x^2-10^9*x+1,x,complex);

Output:

                                (1/2)                     (1/2)
                   /          2\             /          2\     
              -b + \-4 a c + b /         b + \-4 a c + b /     
              -----------------------, - ----------------------
                        2 a                       2 a    
      
                                    9                -9
                      1.000000000 10 , 1.000000000 10  

                                    -9                9
                      1.000000000 10  , 1.000000000 10 

[edit] Mathematica

Possible ways to do this are (symbolic and numeric examples):

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}]

gives back:

\left\{\left\{x\to \frac{-b-\sqrt{b^2-4 a c}}{2 a}\right\},\left\{x\to \frac{-b+\sqrt{b^2-4 a c}}{2 a}\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}

\begin{align}
  \Biggl(
  a & \neq 0  \And \And
     \left(
        x==\frac{-b-\sqrt{b^2-4 a c}}{2 a}
        \|
        x==\frac{-b+\sqrt{b^2-4 a c}}{2 a}
      \right)
  \Biggr)\\
  &\biggl\|
    \left(
      a==0 \And\And b\neq 0 \And\And x==-\frac{c}{b}
    \right)\\
  &\biggr\|
  (c==0 \And \And b==0 \And \And a==0)
\end{align}

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.

[edit] MATLAB / Octave

roots([1 -3 2])    % coefficients in decreasing order of power e.g. [x^n ... x^2 x^1 x^0]

[edit] Maxima

solve(a*x^2 + b*x + c = 0, x);
 
/* 2 2
sqrt(b - 4 a c) + b sqrt(b - 4 a c) - b
[x = - --------------------, x = --------------------]
2 a 2 a */
 
fpprec: 40$
 
solve(x^2 - 10^9*x + 1 = 0, x);
/* [x = 500000000 - sqrt(249999999999999999),
x = sqrt(249999999999999999) + 500000000] */
 
bfloat(%);
/* [x = 1.0000000000000000009999920675269450501b-9,
x = 9.99999999999999998999999999999999999b8] */

[edit] МК-61/52

П2	С/П	/-/	<->	/	2	/	П3	x^2	С/П
ИП2 / - Вx <-> КвКор НОП x>=0 28 ИП3
x<0 24 <-> /-/ + / Вx С/П /-/ КвКор
ИП3 С/П

Input: a С/П b С/П c С/П

Output: x1 - РX; x2 - РY (or error message, if D < 0).

[edit] Modula-3

Translation of: Ada
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.

[edit] 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))
;;

Sample output:

# 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)

[edit] Octave

See MATLAB.

[edit] PARI/GP

Translation of: C
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]
)
};

[edit] Pascal

some parts translated from Modula2

Program QuadraticRoots;
 
var
a, b, c, q, f: double;
 
begin
a := 1;
b := -10e9;
c := 1;
q := sqrt(a * c) / b;
f := (1 + sqrt(1 - 4 * q * q)) / 2;
 
writeln ('Version 1:');
writeln ('x1: ', (-b * f / a):16, ', x2: ', (-c / (b * f)):16);
 
writeln ('Version 2:');
q := sqrt(b * b - 4 * a * c);
if b < 0 then
begin
f := (-b + q) / 2 * a;
writeln ('x1: ', f:16, ', x2: ', (c / (a * f)):16);
end
else
begin
f := (-b - q) / 2 * a;
writeln ('x1: ', (c / (a * f)):16, ', x2: ', f:16);
end;
end.
 

Output:

Version 1:
x1:  1.00000000E+010, x2:  1.00000000E-010
Version 2:
x1:  1.00000000E+010, x2:  1.00000000E-010

[edit] Perl

When using Math::Complex perl automatically convert numbers when necessary.

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);
}

[edit] Perl 6

Perl 6 has complex number handling built in.

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;
}

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)

[edit] PicoLisp

(scl 40)
 
(de solveQuad (A B C)
(let SD (sqrt (- (* B B) (* 4 A C)))
(if (lt0 B)
(list
(*/ (- SD B) A 2.0)
(*/ C 2.0 (*/ A A (- SD B) `(* 1.0 1.0))) )
(list
(*/ C 2.0 (*/ A A (- 0 B SD) `(* 1.0 1.0)))
(*/ (- 0 B SD) A 2.0) ) ) ) )
 
(mapcar round
(solveQuad 1.0 -1000000.0 1.0)
(6 .) )

Output:

-> ("999,999.999999" "0.000001")

[edit] 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;
 

[edit] 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
>>>

[edit] R

Translation of: Python
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", r[[1]], r[[2]], r[[3]]))
}

[edit] Racket

#lang racket
(define (quadratic a b c)
(let* ((-b (- b))
(delta (- (expt b 2) (* 4 a c)))
(denominator (* 2 a)))
(list
(/ (+ -b (sqrt delta)) denominator)
(/ (- -b (sqrt delta)) denominator))))
 
;(quadratic 1 0.0000000000001 -1)
;'(0.99999999999995 -1.00000000000005)
;(quadratic 1 0.0000000000001 1)
;'(-5e-014+1.0i -5e-014-1.0i)

[edit] REXX

[edit] version 1

The REXX language has no   sqrt   function nor does it have complex numbers.
Since "unlimited" decimal precision is part of the language, the NUMERIC DIGITS
was increased (from a default of 9) to 120 to accomodate roots closer to zero than the other roots.
Note that only 10 digits are shown in the output (precision).

/*REXX program finds the roots (may be complex) of a quadratic function.*/
numeric digits 120 /*use enough digits for extremes.*/
parse arg a b c . /*get specified arguments: A B C*/
a=a/1; b=b/1; c=c/1 /*normalize the three numbers. */
call quadratic a b c /*solve the quadratic function. */
numeric digits sqrt(digits())%1 /*reduce digits for human beans. */
r1=r1/1 /*normalize to the new digits. */
r2=r2/1 /* " " " " " */
if r1j\=0 then r1=r1 || left('+',r1j>0)(r1j/1)"i" /*handle complex num.*/
if r2j\=0 then r2=r2 || left('+',r2j>0)(r2j/1)"i" /* " " " */
say ' a =' a /*show value of A. */
say ' b =' b /* " " " B. */
say ' c =' c /* " " " C. */
say
say 'root1 =' r1 /*show 1st root (may be complex).*/
say 'root2 =' r2 /* " 2nd " " " " */
exit /*stick a fork in it, we're done.*/
/*──────────────────────────────────QUADRATIC subroutine────────────────*/
quadratic: parse arg aa bb cc . /*obtain the specified arguments.*/
?=sqrt(bb**2-4*aa*cc) /*compute sqrt (might be complex)*/
aa2=1 / (aa+aa) /*compute reciprocal of 2*aa */
if right(?,1)=='i' then do /*are the roots complex? */
 ?i=left(?,length(?)-1)
r1=-bb*aa2; r2=r1; r1j=?i*aa2; r2j=-?i*aa2
end
else do
r1=(-bb+?)*aa2; r2=(-bb-?)*aa2; r1j=0; r2j=0
end
return
/*──────────────────────────────────SQRT subroutine─────────────────────*/
sqrt: procedure; parse arg x,f; if x=0 then return 0; d=digits()
numeric digits 11; g=.sqrtG(); do j=0 while p>9; m.j=p; p=p%2+1; end
do k=j+5 to 0 by -1; if m.k>11 then numeric digits m.k; g=.5*(g+x/g); end
numeric digits d;return (g/1)i
.sqrtG: i=left('i',x<0); numeric form; m.=11; p=d+d%4+2; x=abs(x)
parse value format(x,2,1,,0) 'E0' with g 'E' _ .; return g*.5'E'_%2

output when using the input of: 1 -10e5 1

    a = 1
    b = -1000000
    c = 1

root1 = 1000000
root2 = 0.000001

output when using the input of: 1 -10e9 1

    a = 1
    b = -10000000000
    c = 1

root1 = 1.000000000E+10
root2 = 1E-10

output when using the input of: 3 2 1

    a = 3
    b = 2
    c = 1

root1 = -0.3333333333-0.5333870616i
root2 = -0.3333333333+0.5333870616i

output when using the input of: 1 0 1

    a = 1
    b = 0
    c = 1

root1 = 0+1i
root2 = 0-1i

[edit] Version 2

/* REXX ***************************************************************
* 26.07.2913 Walter Pachl
**********************************************************************/

Numeric Digits 30
Parse Arg a b c 1 alist
Select
When a='' | a='?' Then
Call exit 'rexx qgl a b c solves a*x**2+b*x+c'
When words(alist)<>3 Then
Call exit 'three numbers are required'
Otherwise
Nop
End
gl=a'*x**2'
Select
When b<0 Then gl=gl||b'*x'
When b>0 Then gl=gl||'+'||b'*x'
Otherwise Nop
End
Select
When c<0 Then gl=gl||c
When c>0 Then gl=gl||'+'||c
Otherwise Nop
End
Say gl '= 0'
 
d=b**2-4*a*c
If d<0 Then Do
dd=sqrt(-d)
r=-b/(2*a)
i=dd/(2*a)
x1=r'+'i'i'
x2=r'-'i'i'
End
Else Do
dd=sqrt(d)
x1=(-b+dd)/(2*a)
x2=(-b-dd)/(2*a)
End
Say 'x1='||x1
Say 'x2='||x2
Exit
sqrt:
/* REXX ***************************************************************
* EXEC to calculate the square root of x with high precision
**********************************************************************/

Parse Arg x
prec=digits()
prec1=2*prec
eps=10**(-prec1)
k = 1
Numeric Digits prec1
r0= x
r = 1
Do i=1 By 1 Until r=r0 | (abs(r*r-x)<eps)
r0 = r
r = (r + x/r) / 2
k = min(prec1,2*k)
Numeric Digits (k + 5)
End
Numeric Digits prec
Return (r+0)
exit: Say arg(1)

Output:

Version 1:
a= 1
b= -1.0000000001
c= 0.000000001

root1= 0.9999999991
root2= 0.000000001000000001

Version 2:
1*x**2-1.0000000001*x+1.e-9 = 0
x1=0.9999999991000000000025
x2=0.0000000009999999999975

[edit] Ruby

With the 'complex' package from the standard library, the Math#sqrt method will return a Complex instance if necessary.

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) # [-3/2, 5]
p quadratic(1, -2, 1) # [1]
p quadratic(1, 3, 3) # [(-3 + sqrt(3)i)/2), (-3 - sqrt(3)i)/2)]
Output:
[-0.6666666666666666, -0.6666666666666666]
[0.3333333333333333, -1.0]
[((-2/6)+0.47140452079103173i), ((-2/6)-0.47140452079103173i)]
[((0/2)+1.0i), ((0/2)-1.0i)]
[999999.999999, 1.00000761449337e-06]
[-1.5, 5.0]
[1.0, 1.0]
[((-3/2)+0.8660254037844386i), ((-3/2)-0.8660254037844386i)]

[edit] Run BASIC

print "FOR 1,2,3 => ";quad$(1,2,3)
print "FOR 4,5,6 => ";quad$(4,5,6)
 
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
FOR 1,2,3 => -1 +i1.41421356 , -1 -i1.41421356
FOR 4,5,6 => -0.625 +i1.05326872 , -0.625 -i1.05326872

[edit] Scala

Using Complex class from task Arithmetic/Complex.

import ArithmeticComplex._
object QuadraticRoots {
def solve(a:Double, b:Double, c:Double)={
val d = b*b-4.0*a*c
val aa = a+a
 
if (d < 0.0) { // complex roots
val re= -b/aa;
val im = math.sqrt(-d)/aa;
(Complex(re, im), Complex(re, -im))
}
else { // real roots
val re=if (b < 0.0) (-b+math.sqrt(d))/aa else (-b -math.sqrt(d))/aa
(re, (c/(a*re)))
}
}
}

Usage:

val equations=Array(
(1.0, 22.0, -1323.0), // two distinct real roots
(6.0, -23.0, 20.0), // with a != 1.0
(1.0, -1.0e9, 1.0), // with one root near zero
(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
);
 
equations.foreach{v =>
val (a,b,c)=v
println("a=%g b=%g c=%g".format(a,b,c))
val roots=solve(a, b, c)
println("x1="+roots._1)
if(roots._1 != roots._2) println("x2="+roots._2)
println
}

Output:

a=1.00000   b=22.0000   c=-1323.00
x1=-49.0
x2=27.0

a=6.00000   b=-23.0000   c=20.0000
x1=2.5
x2=1.3333333333333333

a=1.00000   b=-1.00000e+09   c=1.00000
x1=1.0E9
x2=1.0E-9

a=1.00000   b=2.00000   c=1.00000
x1=-1.0

a=1.00000   b=0.00000   c=1.00000
x1=-0.0 + 1.0i
x2=-0.0 + -1.0i

a=1.00000   b=1.00000   c=1.00000
x1=-0.5 + 0.8660254037844386i
x2=-0.5 + -0.8660254037844386i

[edit] 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)

[edit] Seed7

Translation of: Ada
$ include "seed7_05.s7i";
include "float.s7i";
include "math.s7i";
 
const type: roots is new struct
var float: x1 is 0.0;
var float: x2 is 0.0;
end struct;
 
const func roots: solve (in float: a, in float: b, in float: c) is func
result
var roots: solution is roots.value;
local
var float: sd is 0.0;
var float: x is 0.0;
begin
sd := sqrt(b**2 - 4.0 * a * c);
if b < 0.0 then
x := (-b + sd) / 2.0 * a;
solution.x1 := x;
solution.x2 := c / (a * x);
else
x := (-b - sd) / 2.0 * a;
solution.x1 := c / (a * x);
solution.x2 := x;
end if;
end func;
 
const proc: main is func
local
var roots: r is roots.value;
begin
r := solve(1.0, -10.0E5, 1.0);
writeln("X1 = " <& r.x1 digits 6 <& " X2 = " <& r.x2 digits 6);
end func;

Output:

X1 = 1000000.000000 X2 = 0.000001

[edit] Tcl

Library: Tcllib (Package: math::complexnumbers)
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)}

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

[edit] 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.

Personal tools
Namespaces

Variants
Actions
Community
Explore
Misc
Toolbox