Roots of a quadratic function: Difference between revisions
(→{{header|Java}}: Fix for a = 1, b = -10e9, c = 1, as task requires.) |
|||
Line 556: | Line 556: | ||
1 _1e5 1 p. 1e5 1e_5-1e_5 _1e_15 NB. test displayed roots with adjustment |
1 _1e5 1 p. 1e5 1e_5-1e_5 _1e_15 NB. test displayed roots with adjustment |
||
_3.38436e_7 0</lang> |
_3.38436e_7 0</lang> |
||
Updated example: |
|||
<lang j> 1 {:: p.1 _1e9 1 |
|||
1e9 1e_9</lang> |
|||
=={{header|Java}}== |
=={{header|Java}}== |
Revision as of 19:31, 25 April 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 . Your program must correctly handle non-real roots, but it need not check that .
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 , , and . (For double-precision floats, set .) 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 can see, the second root has lost all significant figures. The right answer is that X2
is about . The naive method is numerically unstable.
Task: do it better. This means that given , , and , both of the roots your program returns should be greater than . Or, if your language can't do floating-point arithmetic any more precisely than single precision, your program should be able to handle . 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.
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 csharp>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)); }
}</lang> Output: <lang>((1E+20, 0), (1E-20, 0))</lang>
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, std.math, std.traits ;
auto naiveQR(T,U,V)(T a, U b, V c) if(isFloatingPoint!T) {
alias CommonType!(T,U,V) RT ; if(a == 0) return [cast(RT)c/b] ; // linear RT det = b*b - 4 * a * c ; if(det < 0) return cast(RT[])[] ; // no real number root immutable SD = sqrt(det) ;
return [(-b + SD) / 2 * a, (-b - SD) / 2 * a] ;
}
auto cautiQR(T,U,V)(T a, U b, V c) if(isFloatingPoint!T) {
alias CommonType!(T,U,V) RT ; if(a == 0) return [cast(RT)c/b] ; // linear RT det = b*b - 4 * a * c ; if(det < 0) return cast(RT[])[] ; // no real number root immutable SD = sqrt(det) ;
if(b*a < 0) { auto x = (-b + SD) / 2 * a ; return [x, c/(a * x)] ; } else { auto x = (-b - SD) / 2 * a ; return [c/(a * x), x] ; }
}
void main() {
writefln(" naive : %6g", naiveQR(1F,-10e5F,1)) ; // force use 32bit float writefln("cautious : %6g", cautiQR(1F,-10e5F,1)) ;
}</lang> output:
naive : [ 1e+06, 0] cautious : [ 1e+06, 1e-06]
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>
Go
<lang 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() {
var zero float64 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]) }
}</lang> 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
Haskell
<lang 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</lang>
*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)
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>
Exploring the limits of precision:
<lang j> 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</lang>
Updated example:
<lang j> 1 {:: p.1 _1e9 1 1e9 1e_9</lang>
Java
<lang 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, -10.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]); } } }
}</lang> 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+10 c = 1.00000 X1 = 1.00000e+10 X2 = 1.00000e-10 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>
PicoLisp
<lang 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 .) )</lang>
Output:
-> ("999,999.999999" "0.000001")
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>
REXX
The REXX language has no SQRT function nor does it have complex numbers.
The SQRT function is included below.
Since unlimited decimal precision is part of the language, the NUMERIC DIGITS
was increased (from a default of 9) to 120 to accomidate roots closer to zero than the other.
<lang rexx>
/*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. */
say 'a=' a /*show value of A. */
say 'b=' b /* " " " B. */
say 'c=' c /* " " " C. */
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
say 'root1=' r1 /*show 1st root (may be complex).*/
say 'root2=' r2 /* " 2nd " " " " */
exit
/*─────────────────────────────────────QUADRATIC subroutine─────────────*/ quadratic: parse arg aa bb cc . ?=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=sqrtguess();
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
sqrtguess: i=left('i',x<0); numeric form scientific; m.=11; p=d+d%4+2;
parse value format(x,2,1,,0) 'E0' with g 'E' _ .; return g*.5'E'_%2
</lang>
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
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
.