Roots of a quadratic function: Difference between revisions
m
→{{header|Wren}}: Minor tidy
m (→{{header|Ruby}}: Remove phantom category 'complex.rb'. This is part of standard library.) |
m (→{{header|Wren}}: Minor tidy) |
||
| (124 intermediate revisions by 52 users not shown) | |||
Line 1:
{{task|Arithmetic operations}}{{Clarified-review}}Write a program to find the roots of a quadratic equation, i.e., solve the equation <math>ax^2 + bx + c = 0</math>.
Your program must correctly handle non-real roots, but it need not check that <math>a \neq 0</math>.
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 ''[http://www.pdas.com/fmm.htm Computer Methods for Mathematical Computations]'', George Forsythe, Michael Malcolm, and Cleve Moler suggest trying the naive algorithm with <math>a = 1</math>, <math>b = -10^5</math>, and <math>c = 1</math>.
(For double-precision floats, set <math>b = -10^9</math>.)
Consider the following implementation in [[Ada]]:
<syntaxhighlight lang="ada">with Ada.Text_IO; use Ada.Text_IO;
with Ada.Numerics.Elementary_Functions; use Ada.Numerics.Elementary_Functions;
Line 17 ⟶ 22:
begin
Put_Line ("X1 =" & Float'Image (R (1)) & " X2 =" & Float'Image (R (2)));
end Quadratic_Equation;</
{{out}}
<pre>X1 = 1.00000E+06 X2 = 0.00000E+00</pre>
As we can see, the second root has lost all significant figures. The right answer is that <code>X2</code> is about <math>10^{-6}</math>. The naive method is numerically unstable.
Line 27 ⟶ 32:
'''Task''': do it better. This means that given <math>a = 1</math>, <math>b = -10^9</math>, and <math>c = 1</math>, both of the roots your program returns should be greater than <math>10^{-11}</math>. Or, if your language can't do floating-point arithmetic any more precisely than single precision, your program should be able to handle <math>b = -10^6</math>. Either way, show what your program gives as the roots of the quadratic in question. See page 9 of
[https://www.validlab.com/goldberg/paper.pdf "What Every Scientist Should Know About Floating-Point Arithmetic"] for a possible algorithm.
=={{header|11l}}==
<syntaxhighlight lang="11l">F quad_roots(a, b, c)
V sqd = Complex(b^2 - 4*a*c) ^ 0.5
R ((-b + sqd) / (2 * a),
(-b - sqd) / (2 * a))
V testcases = [(3.0, 4.0, 4 / 3),
(3.0, 2.0, -1.0),
(3.0, 2.0, 1.0),
(1.0, -1e9, 1.0),
(1.0, -1e100, 1.0)]
L(a, b, c) testcases
V (r1, r2) = quad_roots(a, b, c)
print(r1, end' ‘ ’)
print(r2)</syntaxhighlight>
{{out}}
<pre>
-0.666667+0i -0.666667+0i
0.333333+0i -1+0i
-0.333333+0.471405i -0.333333-0.471405i
1e+09+0i 0i
1e+100+0i 0i
</pre>
=={{header|Ada}}==
<
with Ada.Numerics.Elementary_Functions; use Ada.Numerics.Elementary_Functions;
Line 40 ⟶ 72:
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;
Line 51 ⟶ 83:
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.
{{out}}
<pre>
X1 = 1.00000E+06 X2 = 1.00000E-06
Line 64 ⟶ 97:
{{works with|ALGOL 68G|Any - tested with release [http://sourceforge.net/projects/algol68/files/algol68g/algol68g-1.18.0/algol68g-1.18.0-9h.tiny.el5.centos.fc11.i386.rpm/download 1.18.0-9h.tiny]}}
{{wont work with|ELLA ALGOL 68|Any (with appropriate job cards) - tested with release [http://sourceforge.net/projects/algol68/files/algol68toc/algol68toc-1.8.8d/algol68toc-1.8-8d.fc9.i386.rpm/download 1.8-8d] - probably need to "USE" the compl ENVIRON}}
<
BEGIN
Line 118 ⟶ 151:
ESAC
OD
END # quadratic_equation #</
{{out}}
<pre>
Quadratic: 1.00000000 x**2 + -1000000.00000000 x + 1.00000000 = 0
Line 142 ⟶ 175:
=={{header|AutoHotkey}}==
ahk forum: [http://www.autohotkey.com/forum/viewtopic.php?p=276617#276617 discussion]
<
MsgBox % quadratic(u,v, 1,3,2) ", " u ", " v
MsgBox % quadratic(u,v, -2,4,-2) ", " u ", " v
Line 164 ⟶ 197:
x2 := c/a/x1
Return 2
}</
=={{header|BBC BASIC}}==
<
READ a$, b$, c$
PRINT "For a = " ; a$ ", b = " ; b$ ", c = " ; c$ TAB(32) ;
Line 194 ⟶ 227:
PRINT "the complex roots are " ; -b/2/a " +/- " ; SQR(-d)/2/a "*i"
ENDCASE
ENDPROC</
{{out}}
<pre>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
Line 206 ⟶ 239:
=={{header|C}}==
Code that tries to avoid floating point overflow and other unfortunate loss of precissions: (compiled with <code>gcc -std=c99</code> for <code>complex</code>, though easily adapted to just real numbers)
<
#include <stdlib.h>
#include <complex.h>
Line 266 ⟶ 299:
return 0;
}</syntaxhighlight>
{{out}}<pre>(
(1.00208e+07 + 0 i), (9.9792e-08 + 0 i)</pre>
<
#include <math.h>
#include <complex.h>
Line 280 ⟶ 314:
x[0] = (-b + csqrt(delta)) / (2.0*a);
x[1] = (-b - csqrt(delta)) / (2.0*a);
}</
{{trans|C++}}
<
{
b /= a;
Line 296 ⟶ 330:
x[1] = c/sol;
}
}</
<
{
complex double x[2];
Line 312 ⟶ 346:
return 0;
}</
<pre>x1 = (1.00000000000000000000e+20, 0.00000000000000000000e+00)
Line 321 ⟶ 355:
=={{header|C sharp|C#}}==
<
using System.Numerics;
Line 336 ⟶ 370:
Console.WriteLine(Solve(1, -1E20, 1));
}
}</
{{out}}
<
=={{header|C++}}==
<
#include <utility>
#include <complex>
Line 367 ⟶ 402:
std::pair<complex, complex> result = solve_quadratic_equation(1, -1e20, 1);
std::cout << result.first << ", " << result.second << std::endl;
}</
{{out}}
(1e+20,0), (1e-20,0)
=={{header|Clojure}}==
<syntaxhighlight lang="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</syntaxhighlight>
{{out}}
<syntaxhighlight lang="clojure">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]
</syntaxhighlight>
=={{header|Common Lisp}}==
<
(list
(
(/ (- (- b) (sqrt (
=={{header|D}}==
<syntaxhighlight lang
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) {
return [x, c / (a * x)]
} else {
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));
}</syntaxhighlight>
{{out}}
<pre>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]</pre>
=={{header|Delphi}}==
See [https://rosettacode.org/wiki/Roots_of_a_quadratic_function#Pascal Pascal].
=={{header|Elixir}}==
<syntaxhighlight lang="elixir">defmodule Quadratic do
def roots(a, b, c) do
IO.puts "Roots of a quadratic function (#{a}, #{b}, #{c})"
d = b * b - 4 * a * c
a2 = a * 2
cond do
d > 0 ->
sd = :math.sqrt(d)
IO.puts " the real roots are #{(- b + sd) / a2} and #{(- b - sd) / a2}"
d == 0 ->
IO.puts " the single root is #{- b / a2}"
true ->
sd = :math.sqrt(-d)
IO.puts " the complex roots are #{- b / a2} +/- #{sd / a2}*i"
end
end
end
Quadratic.roots(1, -2, 1)
Quadratic.roots(1, -3, 2)
Quadratic.roots(1, 0, 1)
Quadratic.roots(1, -1.0e10, 1)
Quadratic.roots(1, 2, 3)
Quadratic.roots(2, -1, -6)</syntaxhighlight>
{{out}}
<pre>
Roots of a quadratic function (1, -2, 1)
the single root is 1.0
Roots of a quadratic function (1, -3, 2)
the real roots are 2.0 and 1.0
Roots of a quadratic function (1, 0, 1)
the complex roots are 0.0 +/- 1.0*i
Roots of a quadratic function (1, -1.0e10, 1)
the real roots are 1.0e10 and 0.0
Roots of a quadratic function (1, 2, 3)
the complex roots are -1.0 +/- 1.4142135623730951*i
Roots of a quadratic function (2, -1, -6)
the real roots are 2.0 and -1.5
</pre>
=={{header|ERRE}}==
<syntaxhighlight lang="text">PROGRAM QUADRATIC
PROCEDURE SOLVE_QUADRATIC
D=B*B-4*A*C
IF ABS(D)<1D-6 THEN D=0 END IF
CASE SGN(D) OF
0->
PRINT("the single root is ";-B/2/A)
END ->
1->
F=(1+SQR(1-4*A*C/(B*B)))/2
PRINT("the real roots are ";-F*B/A;"and ";-C/B/F)
END ->
-1->
PRINT("the complex roots are ";-B/2/A;"+/-";SQR(-D)/2/A;"*i")
END ->
END CASE
END PROCEDURE
BEGIN
PRINT(CHR$(12);) ! CLS
FOR TEST%=1 TO 7 DO
READ(A,B,C)
PRINT("For a=";A;",b=";B;",c=";C;TAB(32);)
SOLVE_QUADRATIC
END FOR
DATA(1,-1E9,1)
DATA(1,0,1)
DATA(2,-1,-6)
DATA(1,2,-2)
DATA(0.5,1.4142135,1)
DATA(1,3,2)
DATA(3,4,5)
END PROGRAM</syntaxhighlight>
{{out}}
<pre>For a= 1 ,b=-1E+09 ,c= 1 the real roots are 1E+09 and 1E-09
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.732051 and .7320508
For a= .5 ,b= 1.414214 ,c= 1 the single root is -1.414214
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 -.6666667 +/- 1.105542 *i
</pre>
=={{header|Factor}}==
{{trans|Ada}}
<
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 * / ;</
<
--- Data stack:
1.0e+20
9.999999999999999e-21</
Middlebrook method
<
a c * sqrt b / :> q
1 4 q sq * - sqrt 0.5 * 0.5 + :> f
b neg a / f * c neg b / f / ;
</syntaxhighlight>
<
--- Data stack:
1.0e+20
1.0e-20</
=={{header|Forth}}==
Without locals:
<
frot frot
( c a b )
Line 458 ⟶ 627:
2e f/ fover f/
( c a r1 )
frot frot f/ fover f/ ;</
With locals:
<
b b f* 4e a f* c f* f-
fdup f0< if abort" imaginary roots" then
Line 469 ⟶ 638:
\ test
1 set-precision
1e -1e6 1e quadratic fs. fs. \ 1e-6 1e6</
=={{header|Fortran}}==
===Fortran 90===
{{works with|Fortran|90 and later}}
<
IMPLICIT NONE
Line 493 ⟶ 663:
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)
WRITE(*,*) "The roots are real:"
WRITE(*,"(2(A,E23.15))") "Root1 = ", rroot1, " Root2 = ", rroot2
Line 511 ⟶ 676:
WRITE(*,*) "The roots are complex:"
WRITE(*,"(2(A,2E23.15,A))") "Root1 = ", croot1, "j ", " Root2 = ", croot2, "j"
END IF</
{{out}}
Coefficients are: a = 0.300000000000000E+01 b = 0.400000000000000E+01 c = 0.133333333330000E+01
The roots are real and equal:
Line 528 ⟶ 693:
The roots are real:
Root1 = 0.999999999999000E+06 Root2 = 0.100000000000100E-05
===Fortran I===
Source code written in FORTRAN I (october 1956) for the IBM 704.
<syntaxhighlight lang="fortran">
COMPUTE ROOTS OF A QUADRATIC FUNCTION - 1956
READ 100,A,B,C
100 FORMAT(3F8.3)
PRINT 100,A,B,C
DISC=B**2-4.*A*C
IF(DISC),1,2,3
1 XR=-B/(2.*A)
XI=SQRT(-DISC)/(2.*A)
XJ=-XI
PRINT 311
PRINT 312,XR,XI,XR,XJ
311 FORMAT(13HCOMPLEX ROOTS)
312 FORMAT(4HX1=(,2E12.4,6H),X2=(,2E12.4,1H))
GO TO 999
2 X1=-B/(2.*A)
X2=X1
PRINT 321
PRINT 332,X1,X2
321 FORMAT(16HEQUAL REAL ROOTS)
GO TO 999
3 X1= (-B+SQRT(DISC)) / (2.*A)
X2= (-B-SQRT(DISC)) / (2.*A)
PRINT 331
PRINT 332,X1,X2
331 FORMAT(10HREAL ROOTS)
332 FORMAT(3HX1=,E12.5,4H,X2=,E12.5)
999 STOP
</syntaxhighlight>
=={{header|FreeBASIC}}==
{{libheader|GMP}}
<syntaxhighlight lang="freebasic">' version 20-12-2020
' compile with: fbc -s console
#Include Once "gmp.bi"
Sub solvequadratic_n(a As Double ,b As Double, c As Double)
Dim As Double f, d = b ^ 2 - 4 * a * c
Select Case Sgn(d)
Case 0
Print "1: the single root is "; -b / 2 / a
Case 1
Print "1: the real roots are "; (-b + Sqr(d)) / 2 * a; " and ";(-b - Sqr(d)) / 2 * a
Case -1
Print "1: the complex roots are "; -b / 2 / a; " +/- "; Sqr(-d) / 2 / a; "*i"
End Select
End Sub
Sub solvequadratic_c(a As Double ,b As Double, c As Double)
Dim As Double f, d = b ^ 2 - 4 * a * c
Select Case Sgn(d)
Case 0
Print "2: the single root is "; -b / 2 / a
Case 1
f = (1 + Sqr(1 - 4 * a *c / b ^ 2)) / 2
Print "2: the real roots are "; -f * b / a; " and "; -c / b / f
Case -1
Print "2: the complex roots are "; -b / 2 / a; " +/- "; Sqr(-d) / 2 / a; "*i"
End Select
End Sub
Sub solvequadratic_gmp(a_ As Double ,b_ As Double, c_ As Double)
#Define PRECISION 1024 ' about 300 digits
#Define MAX 25
Dim As ZString Ptr text
text = Callocate (1000)
Mpf_set_default_prec(PRECISION)
Dim As Mpf_ptr a, b, c, d, t
a = Allocate(Len(__mpf_struct)) : Mpf_init_set_d(a, a_)
b = Allocate(Len(__mpf_struct)) : Mpf_init_set_d(b, b_)
c = Allocate(Len(__mpf_struct)) : Mpf_init_set_d(c, c_)
d = Allocate(Len(__mpf_struct)) : Mpf_init(d)
t = Allocate(Len(__mpf_struct)) : Mpf_init(t)
mpf_mul(d, b, b)
mpf_set_ui(t, 4)
mpf_mul(t, t, a)
mpf_mul(t, t, c)
mpf_sub(d, d, t)
Select Case mpf_sgn(d)
Case 0
mpf_neg(t, b)
mpf_div_ui(t, t, 2)
mpf_div(t, t, a)
Gmp_sprintf(text,"%.*Fe", MAX, t)
Print "3: the single root is "; *text
Case Is > 0
mpf_sqrt(d, d)
mpf_add(a, a, a)
mpf_neg(t, b)
mpf_add(t, t, d)
mpf_div(t, t, a)
Gmp_sprintf(text,"%.*Fe", MAX, t)
Print "3: the real roots are "; *text; " and ";
mpf_neg(t, b)
mpf_sub(t, t, d)
mpf_div(t, t, a)
Gmp_sprintf(text,"%.*Fe", MAX, t)
Print *text
Case Is < 0
mpf_neg(t, b)
mpf_div_ui(t, t, 2)
mpf_div(t, t, a)
Gmp_sprintf(text,"%.*Fe", MAX, t)
Print "3: the complex roots are "; *text; " +/- ";
mpf_neg(t, d)
mpf_sqrt(t, t)
mpf_div_ui(t, t, 2)
mpf_div(t, t, a)
Gmp_sprintf(text,"%.*Fe", MAX, t)
Print *text; "*i"
End Select
End Sub
' ------=< MAIN >=------
Dim As Double a, b, c
Print "1: is the naieve way"
Print "2: is the cautious way"
Print "3: is the naieve way with help of GMP"
Print
For i As Integer = 1 To 10
Read a, b, c
Print "Find root(s) for "; Str(a); "X^2"; IIf(b < 0, "", "+");
Print Str(b); "X"; IIf(c < 0, "", "+"); Str(c)
solvequadratic_n(a, b , c)
solvequadratic_c(a, b , c)
solvequadratic_gmp(a, b , c)
Print
Next
' empty keyboard buffer
While Inkey <> "" : Wend
Print : Print "hit any key to end program"
Sleep
End
Data 1, -1E9, 1
Data 1, 0, 1
Data 2, -1, -6
Data 1, 2, -2
Data 0.5, 1.4142135623731, 1
Data 1, 3, 2
Data 3, 4, 5
Data 1, -1e100, 1
Data 1, -1e200, 1
Data 1, -1e300, 1</syntaxhighlight>
{{out}}
<pre>1: is the naieve way
2: is the cautious way
3: is the naieve way with help of GMP
Find root(s) for 1X^2-1000000000X+1
1: the real roots are 1000000000 and 0
2: the real roots are 1000000000 and 1e-009
3: the real roots are 9.9999999999999999900000000e+08 and 1.0000000000000000010000000e-09
Find root(s) for 1X^2+0X+1
1: the complex roots are -0 +/- 1*i
2: the complex roots are -0 +/- 1*i
3: the complex roots are 0.0000000000000000000000000e+00 +/- 1.0000000000000000000000000e+00*i
Find root(s) for 2X^2-1X-6
1: the real roots are 8 and -6
2: the real roots are 2 and -1.5
3: the real roots are 2.0000000000000000000000000e+00 and -1.5000000000000000000000000e+00
Find root(s) for 1X^2+2X-2
1: the real roots are 0.7320508075688772 and -2.732050807568877
2: the real roots are -2.732050807568877 and 0.7320508075688773
3: the real roots are 7.3205080756887729352744634e-01 and -2.7320508075688772935274463e+00
Find root(s) for 0.5X^2+1.4142135623731X+1
1: the real roots are -0.3535533607909526 and -0.3535534203955974
2: the real roots are -1.414213681582389 and -1.414213443163811
3: the real roots are -1.4142134436707580875788206e+00 and -1.4142136810754419733330398e+00
Find root(s) for 1X^2+3X+2
1: the real roots are -1 and -2
2: the real roots are -2 and -0.9999999999999999
3: the real roots are -1.0000000000000000000000000e+00 and -2.0000000000000000000000000e+00
Find root(s) for 3X^2+4X+5
1: the complex roots are -0.6666666666666666 +/- 1.105541596785133*i
2: the complex roots are -0.6666666666666666 +/- 1.105541596785133*i
3: the complex roots are -6.6666666666666666666666667e-01 +/- 1.1055415967851332830383109e+00*i
Find root(s) for 1X^2-1e+100X+1
1: the real roots are 1e+100 and 0
2: the real roots are 1e+100 and 1e-100
3: the real roots are 1.0000000000000000159028911e+100 and 9.9999999999999998409710889e-101
Find root(s) for 1X^2-1e+200X+1
1: the real roots are 1.#INF and -1.#INF
2: the real roots are 1e+200 and 1e-200
3: the real roots are 9.9999999999999996973312221e+199 and 0.0000000000000000000000000e+00
Find root(s) for 1X^2-1e+300X+1
1: the real roots are 1.#INF and -1.#INF
2: the real roots are 1e+300 and 1e-300
3: the real roots are 1.0000000000000000525047603e+300 and 0.0000000000000000000000000e+00</pre>
=={{header|GAP}}==
<
local d;
d := Sqrt(b*b - 4*a*c);
Line 538 ⟶ 919:
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 ]</syntaxhighlight>
=={{header|Go}}==
<
import (
Line 587 ⟶ 973:
func main() {
for _, c := range [][3]float64{
{1, -2, 1},
Line 597 ⟶ 982:
test(c[0], c[1], c[2])
}
}</
{{out}}
<pre>
coefficients: 1 -2 1 -> one real root: 1
Line 608 ⟶ 993:
=={{header|Haskell}}==
<
type CD = Complex Double
quadraticRoots :: (CD, CD, CD) -> (CD, CD)
quadraticRoots (a, b, c)
| 0
)
| otherwise =
( (- b + d) / (2 * a),
(2 * c) / (- b + d)
)
where
d = sqrt $ b ^ 2 - 4 * a * c
main :: IO ()
main =
mapM_
(print . quadraticRoots)
[ (3, 4, 4 / 3),
(3, 2, -1),
(3, 2, 1),
(1, -10e5, 1),
(1, -10e9, 1)
]</syntaxhighlight>
{{Out}}
<pre>((-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)</pre>
=={{header|Icon}} and {{header|Unicon}}==
{{trans|Ada}}
Works in both languages.
<syntaxhighlight lang="unicon">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</syntaxhighlight>
{{out}}
<pre>
->rqf 1.0 -0.000000001 1.0
1.0*x^2 + -1000000.0*x + 1.0 has roots 999999.999999 and 1.000000000001e-06
->
</pre>
=={{header|IDL}}==
<
print, "input a, press enter, input b, press enter, input c, press enter"
Line 643 ⟶ 1,067:
e= y/z
print, d,e</
=={{header|IS-BASIC}}==
<syntaxhighlight lang="is-basic">
100 PROGRAM "Quadratic.bas"
110 PRINT "Enter coefficients a, b and c:":INPUT PROMPT "a= ,b= ,c= ":A,B,C
120 IF A=0 THEN
130 PRINT "The coefficient of x^2 can not be 0."
140 ELSE
150 LET D=B^2-4*A*C
160 SELECT CASE SGN(D)
170 CASE 0
180 PRINT "The single root is ";-B/2/A
190 CASE 1
200 PRINT "The real roots are ";(-B+SQR(D))/(2*A);"and ";(-B-SQR(D))/(2*A)
210 CASE -1
220 PRINT "The complex roots are ";-B/2/A;"+/- ";STR$(SQR(-D)/2/A);"*i"
230 END SELECT
240 END IF</syntaxhighlight>
=={{header|J}}==
'''Solution''' use J's built-in polynomial solver:
p.
This primitive converts between the coefficient form of a polynomial (with the exponents being the array indices of the coefficients) and the multiplier-and-roots for of a polynomial (with two boxes, the first containing the multiplier and the second containing the roots).
'''Example''' using inputs from other solutions and the unstable example from the task description:
<
> {:"1 p. coeff
_0.666667 _0.666667
_1 0.333333
_0.333333j0.471405 _0.333333j_0.471405
1e6 1e_6
1e9 1e_9</syntaxhighlight>
Of course <code>p.</code> generalizes to polynomials of
<
+-+-----+
|2|4 2 0|
+-+-----+
p. 2 ; 4 2 0 NB. return coefficients
0 16 _12 2</
Exploring the limits of precision:
<
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</
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
<syntaxhighlight lang="j">q_r=: verb define
'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</syntaxhighlight>
=={{header|Java}}==
<
private static class Complex {
double re, im;
Line 759 ⟶ 1,212:
}
}
}</
{{out}}
<pre>
a = 1.00000 b = 22.0000 c = -1323.00
Line 784 ⟶ 1,237:
X1 = -0.500000 + 0.866025i
X2 = -0.500000 - 0.866025i</pre>
=={{header|jq}}==
{{ works with |jq|1.4}}
Currently jq does not include support for complex number operations, so
a small library is included in the first section.
The second section defines <tt>quadratic_roots(a;b;c)</tt>,
which emits a stream of 0 or two solutions, or the value <tt>true</tt> if a==b==c==0.
The third section defines a function for producing a table showing (i, error, solution) for solutions to x^2 - 10^i + 1 = 0 for various values of i.
'''Section 1''': Complex numbers (scrolling window)
<div style="overflow:scroll; height:200px;">
<syntaxhighlight lang="jq"># Complex numbers as points [x,y] in the Cartesian plane
def real(z): if (z|type) == "number" then z else z[0] end;
def imag(z): if (z|type) == "number" then 0 else z[1] end;
def plus(x; y):
if (x|type) == "number" then
if (y|type) == "number" then [ x+y, 0 ]
else [ x + y[0], y[1]]
end
elif (y|type) == "number" then plus(y;x)
else [ x[0] + y[0], x[1] + y[1] ]
end;
def multiply(x; y):
if (x|type) == "number" then
if (y|type) == "number" then [ x*y, 0 ]
else [x * y[0], x * y[1]]
end
elif (y|type) == "number" then multiply(y;x)
else [ x[0] * y[0] - x[1] * y[1], x[0] * y[1] + x[1] * y[0]]
end;
def negate(x): multiply(-1; x);
def minus(x; y): plus(x; multiply(-1; y));
def conjugate(z):
if (z|type) == "number" then [z, 0]
else [z[0], -(z[1]) ]
end;
def invert(z):
if (z|type) == "number" then [1/z, 0]
else
( (z[0] * z[0]) + (z[1] * z[1]) ) as $d
# use "0 + ." to convert -0 back to 0
| [ z[0]/$d, (0 + -(z[1]) / $d)]
end;
def divide(x;y): multiply(x; invert(y));
def magnitude(z):
real( multiply(z; conjugate(z))) | sqrt;
# exp^z
def complex_exp(z):
def expi(x): [ (x|cos), (x|sin) ];
if (z|type) == "number" then z|exp
elif z[0] == 0 then expi(z[1]) # for efficiency
else multiply( (z[0]|exp); expi(z[1]) )
end ;
def complex_sqrt(z):
if imag(z) == 0 and real(z) >= 0 then [(real(z)|sqrt), 0]
else
magnitude(z) as $r
| if $r == 0 then [0,0]
else
(real(z)/$r) as $a
| (imag(z)/$r) as $b
| $r | sqrt as $r
| (($a | acos) / 2)
| [ ($r * cos), ($r * sin)]
end
end ;</syntaxhighlight></div>
'''Section 2:''' quadratic_roots(a;b;c)
<syntaxhighlight lang="jq"># If there are infinitely many solutions, emit true;
# if none, emit empty;
# otherwise always emit two.
# For numerical accuracy, Middlebrook's approach is adopted:
def quadratic_roots(a; b; c):
if a == 0 and b == 0 then
if c == 0 then true # infinitely many
else empty # none
end
elif a == 0 then [-c/b, 0]
elif b == 0 then (complex_sqrt(1/a) | (., negate(.)))
else
divide( plus(1.0; complex_sqrt( minus(1.0; (4 * a * c / (b*b))))); 2) as $f
| negate(divide(multiply(b; $f); a)),
negate(divide(c; multiply(b; $f)))
end
;</syntaxhighlight>
'''Section 3''':
Produce a table showing [i, error, solution] for solutions to x^2 - 10^i + 1 = 0
<syntaxhighlight lang="jq">def example:
def pow(i): . as $in | reduce range(0;i) as $i (1; . * $in);
def poly(a;b;c): plus( plus( multiply(a; multiply(.;.)); multiply(b;.)); c);
def abs: if . < 0 then -. else . end;
def zero(z):
if z == 0 then 0
else (magnitude(z)|abs) as $zero
| if $zero < 1e-10 then "+0" else $zero end
end;
def lpad(n): tostring | (n - length) * " " + .;
range(0; 13) as $i
| (- (10|pow($i))) as $b
| quadratic_roots(1; $b; 1) as $x
| $x | poly(1; $b; 1) as $zero
| "\($i|lpad(4)): error = \(zero($zero)|lpad(18)) x=\($x)"
;
example</syntaxhighlight>
{{Out}} (scrolling window)
<div style="overflow:scroll; height:200px;">
<syntaxhighlight lang="sh">
$ jq -M -r -n -f Roots_of_a_quadratic_function.jq
0: error = +0 x=[0.5,0.8660254037844386]
0: error = +0 x=[0.5000000000000001,-0.8660254037844387]
1: error = +0 x=[9.898979485566356,0]
1: error = +0 x=[0.10102051443364382,-0]
2: error = +0 x=[99.98999899979995,0]
2: error = +0 x=[0.010001000200050014,-0]
3: error = 1.1641532182693481e-10 x=[999.998999999,0]
3: error = +0 x=[0.0010000010000019998,-0]
4: error = +0 x=[9999.999899999999,0]
4: error = +0 x=[0.00010000000100000003,-0]
5: error = +0 x=[99999.99999,0]
5: error = +0 x=[1.0000000001e-05,-0]
6: error = 0.0001220703125 x=[999999.9999989999,0]
6: error = +0 x=[1.000000000001e-06,-0]
7: error = 0.015625 x=[9999999.9999999,0]
7: error = +0 x=[1.0000000000000101e-07,-0]
8: error = 1 x=[99999999.99999999,0]
8: error = +0 x=[1e-08,-0]
9: error = 1 x=[1000000000,0]
9: error = +0 x=[1e-09,-0]
10: error = 1 x=[10000000000,0]
10: error = +0 x=[1e-10,-0]
11: error = 1 x=[100000000000,0]
11: error = +0 x=[1e-11,-0]
12: error = 1 x=[1000000000000,0]
12: error = +0 x=[1e-12,-0]</syntaxhighlight></div>
=={{header|Julia}}==
This solution is an implementation of algorithm from the Goldberg paper cited in the task description. It does check for <tt>a=0</tt> and returns the linear solution in that case. Julia's <tt>sqrt</tt> throws a domain error for negative real inputs, so negative discriminants are converted to complex by adding <tt>0im</tt> prior to taking the square root.
Alternative solutions might make use of Julia's Polynomials or Roots packages.
<syntaxhighlight lang="julia">using Printf
function quadroots(x::Real, y::Real, z::Real)
a, b, c = promote(float(x), y, z)
if a ≈ 0.0 return [-c / b] end
Δ = b ^ 2 - 4a * c
if Δ ≈ 0.0 return [-sqrt(c / a)] end
if Δ < 0.0 Δ = complex(Δ) end
d = sqrt(Δ)
if b < 0.0
d -= b
return [d / 2a, 2c / d]
else
d = -d - b
return [2c / d, d / 2a]
end
end
a = [1, 1, 1.0, 10]
b = [10, 2, -10.0 ^ 9, 1]
c = [1, 1, 1, 1]
for (x, y, z) in zip(a, b, c)
@printf "The roots of %.2fx² + %.2fx + %.2f\n\tx₀ = (%s)\n" x y z join(round.(quadroots(x, y, z), 2), ", ")
end</syntaxhighlight>
{{out}}
<pre>The roots of 1.00x² + 10.00x + 1.00
x₀ = (-0.1, -9.9)
The roots of 1.00x² + 2.00x + 1.00
x₀ = (-1.0)
The roots of 1.00x² + -1000000000.00x + 1.00
x₀ = (1.0e9, 0.0)
The roots of 10.00x² + 1.00x + 1.00
x₀ = (-0.05 + 0.31im, -0.05 - 0.31im)</pre>
=={{header|K}}==
===K6===
{{works with|ngn/k}}
<syntaxhighlight lang="k"> / naive method
/ sqr[x] and sqrt[x] must be provided
quf:{[a;b;c]; s:sqrt[sqr[b]-4*a*c];(-b+s;-b-s)%2*a}
quf[0.5;-2.5;2]
1.0 4.0
quf[1;8;15]
-5.0 -3.0
quf[1;10;1]
-9.898979485566356 -0.10102051443364424
</syntaxhighlight>
=={{header|Kotlin}}==
{{trans|Java}}
<syntaxhighlight lang="scala">import java.lang.Math.*
data class Equation(val a: Double, val b: Double, val c: Double) {
data class Complex(val r: Double, val i: Double) {
override fun toString() = when {
i == 0.0 -> r.toString()
r == 0.0 -> "${i}i"
else -> "$r + ${i}i"
}
}
data class Solution(val x1: Any, val x2: Any) {
override fun toString() = when(x1) {
x2 -> "X1,2 = $x1"
else -> "X1 = $x1, X2 = $x2"
}
}
val quadraticRoots by lazy {
val _2a = a + a
val d = b * b - 4.0 * a * c // discriminant
if (d < 0.0) {
val r = -b / _2a
val i = sqrt(-d) / _2a
Solution(Complex(r, i), Complex(r, -i))
} else {
// avoid calculating -b +/- sqrt(d), to avoid any
// subtractive cancellation when it is near zero.
val r = if (b < 0.0) (-b + sqrt(d)) / _2a else (-b - sqrt(d)) / _2a
Solution(r, c / (a * r))
}
}
}
fun main(args: Array<String>) {
val equations = listOf(Equation(1.0, 22.0, -1323.0), // two distinct real roots
Equation(6.0, -23.0, 20.0), // with a != 1.0
Equation(1.0, -1.0e9, 1.0), // with one root near zero
Equation(1.0, 2.0, 1.0), // one real root (double root)
Equation(1.0, 0.0, 1.0), // two imaginary roots
Equation(1.0, 1.0, 1.0)) // two complex roots
equations.forEach { println("$it\n" + it.quadraticRoots) }
}</syntaxhighlight>
{{out}}
<pre>Equation(a=1.0, b=22.0, c=-1323.0)
X1 = -49.0, X2 = 27.0
Equation(a=6.0, b=-23.0, c=20.0)
X1 = 2.5, X2 = 1.3333333333333333
Equation(a=1.0, b=-1.0E9, c=1.0)
X1 = 1.0E9, X2 = 1.0E-9
Equation(a=1.0, b=2.0, c=1.0)
X1,2 = -1.0
Equation(a=1.0, b=0.0, c=1.0)
X1 = 1.0i, X2 = -1.0i
Equation(a=1.0, b=1.0, c=1.0)
X1 = -0.5 + 0.8660254037844386i, X2 = -0.5 + -0.8660254037844386i</pre>
=={{header|lambdatalk}}==
<syntaxhighlight lang="scheme">
1) using lambdas:
{def equation
{lambda {:a :b :c}
{b equation :a*x{sup 2}+:b*x+:c=0}
{{lambda {:a' :b' :d}
{if {> :d 0}
then {{lambda {:b' :d'}
{equation.disp {+ :b' :d'} {- :b' :d'} 2 real roots}
} :b' {/ {sqrt :d} :a'}}
else {if {< :d 0}
then {{lambda {:b' :d'}
{equation.disp [:b',:d'] [:b',-:d'] 2 complex roots}
} :b' {/ {sqrt {- :d}} :a'} }
else {equation.disp :b' :b' one real double root}
}}
} {* 2 :a} {/ {- :b} {* 2 :a}} {- {* :b :b} {* 4 :a :c}} } }}
2) using let:
{def equation
{lambda {:a :b :c}
{b equation :a*x{sup 2}+:b*x+:c=0}
{let { {:a' {* 2 :a}}
{:b' {/ {- :b} {* 2 :a}}}
{:d {- {* :b :b} {* 4 :a :c}}} }
{if {> :d 0}
then {let { {:b' :b'}
{:d' {/ {sqrt :d} :a'}} }
{equation.disp {+ :b' :d'} {- :b' :d'} 2 real roots} }
else {if {< :d 0}
then {let { {:b' :b'}
{:d' {/ {sqrt {- :d}} :a'}} }
{equation.disp [:b',:d'] [:b',-:d'] 2 complex roots} }
else {equation.disp :b' :b' one real double root} }} }}}
3) a function to display results in an HTML table format
{def equation.disp
{lambda {:x1 :x2 :txt}
{table {@ style="background:#ffa"}
{tr {td :txt: }}
{tr {td x1 = :x1 }}
{tr {td x2 = :x2 }} } }}
4) testing:
equation 1*x2+1*x+-1=0
2 real roots:
x1 = 0.6180339887498949
x2 = -1.618033988749895
equation 1*x2+1*x+1=0
2 complex roots:
x1 = [-0.5,0.8660254037844386]
x2 = [-0.5,-0.8660254037844386]
equation 1*x2+-2*x+1=0
one real double root:
x1 = 1
x2 = 1
</syntaxhighlight>
=={{header|Liberty BASIC}}==
<
'assume a<>0
print quad$(a,b,c)
Line 799 ⟶ 1,581:
quad$=str$(x/(2*a)+sqr(D)/(2*a));" , ";str$(x/(2*a)-sqr(D)/(2*a))
end if
end function</
=={{header|Logo}}==
<
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</
=={{header|Lua}}==
Line 812 ⟶ 1,594:
like that from the Complex Numbers article. However, this should be enough to demonstrate its accuracy:
<
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
Line 820 ⟶ 1,602:
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.
=={{header|
<syntaxhighlight lang="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);</syntaxhighlight>
{{out}}
<pre> (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 </pre>
=={{header|Mathematica}}/{{header|Wolfram Language}}==
Possible ways to do this are (symbolic and numeric examples):
<
Solve[x^2 - 10^5 x + 1 == 0, x]
Root[#1^2 - 10^5 #1 + 1 &, 1]
Line 833 ⟶ 1,634:
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:
Line 872 ⟶ 1,673:
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.
=={{header|MATLAB}} / {{header|Octave}}==
<
=={{header|Maxima}}==
<syntaxhighlight lang="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] */</syntaxhighlight>
=={{header|МК-61/52}}==
<syntaxhighlight lang="text">П2 С/П /-/ <-> / 2 / П3 x^2 С/П
ИП2 / - Вx <-> КвКор НОП x>=0 28 ИП3
x<0 24 <-> /-/ + / Вx С/П /-/ КвКор
ИП3 С/П</syntaxhighlight>
''Input:'' a С/П b С/П c С/П
{{out}} x<sub>1</sub> - РX; x<sub>2</sub> - РY (or error message, if D < 0).
=={{header|Modula-3}}==
{{trans|Ada}}
<
IMPORT IO, Fmt, Math;
Line 890 ⟶ 1,719:
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;
Line 901 ⟶ 1,730:
r := Solve(1.0D0, -10.0D5, 1.0D0);
IO.Put("X1 = " & Fmt.LongReal(r[1]) & " X2 = " & Fmt.LongReal(r[2]) & "\n");
END Quad.</
=={{header|Nim}}==
<syntaxhighlight lang="nim">import math, complex, strformat
const Epsilon = 1e-15
type
SolKind = enum solDouble, solFloat, solComplex
Roots = object
case kind: SolKind
of solDouble:
fvalue: float
of solFloat:
fvalues: (float, float)
of solComplex:
cvalues: (Complex64, Complex64)
func quadRoots(a, b, c: float): Roots =
if a == 0:
raise newException(ValueError, "first coefficient cannot be null.")
let den = a * 2
let Δ = b * b - a * c * 4
if abs(Δ) < Epsilon:
result = Roots(kind: solDouble, fvalue: -b / den)
elif Δ < 0:
let r = -b / den
let i = sqrt(-Δ) / den
result = Roots(kind: solComplex, cvalues: (complex64(r, i), complex64(r, -i)))
else:
let r = (if b < 0: -b + sqrt(Δ) else: -b - sqrt(Δ)) / den
result = Roots(kind: solFloat, fvalues: (r, c / (a * r)))
func `$`(r: Roots): string =
case r.kind
of solDouble:
result = $r.fvalue
of solFloat:
result = &"{r.fvalues[0]}, {r.fvalues[1]}"
of solComplex:
result = &"{r.cvalues[0].re} + {r.cvalues[0].im}i, {r.cvalues[1].re} + {r.cvalues[1].im}i"
when isMainModule:
const Equations = [(1.0, -2.0, 1.0),
(10.0, 1.0, 1.0),
(1.0, -10.0, 1.0),
(1.0, -1000.0, 1.0),
(1.0, -1e9, 1.0)]
for (a, b, c) in Equations:
echo &"Equation: {a=}, {b=}, {c=}"
let roots = quadRoots(a, b, c)
let plural = if roots.kind == solDouble: "" else: "s"
echo &" root{plural}: {roots}"</syntaxhighlight>
{{out}}
<pre>Equation: a=1.0, b=-2.0, c=1.0
root: 1.0
Equation: a=10.0, b=1.0, c=1.0
roots: -0.05 + 0.3122498999199199i, -0.05 + -0.3122498999199199i
Equation: a=1.0, b=-10.0, c=1.0
roots: 9.898979485566356, 0.1010205144336438
Equation: a=1.0, b=-1000.0, c=1.0
roots: 999.998999999, 0.001000001000002
Equation: a=1.0, b=-1000000000.0, c=1.0
roots: 1000000000.0, 1e-09</pre>
=={{header|OCaml}}==
<
| RealRoots of float * float
| ComplexRoots of Complex.t * Complex.t ;;
Line 918 ⟶ 1,818:
{ Complex.re = r; Complex.im = (-.i) })
else
let r =
else ((sqrt d) +. b) /. (-2.0 *. a)
in
RealRoots (r, c /. (r *. a))
;;</syntaxhighlight>
{{out}}
<
- : quadroots = RealRoots (-1.4142135623730951, 1.4142135623730949)
Line 935 ⟶ 1,836:
# quadsolve 1.0 (-1.0e5) 1.0 ;;
- : quadroots = RealRoots (99999.99999, 1.0000000001000001e-005)</
=={{header|Octave}}==
Line 941 ⟶ 1,842:
=={{header|PARI/GP}}==
{{works with|PARI/GP|2.8.0+}}
<syntaxhighlight lang="parigp">roots(a,b,c)=polrootsreal(Pol([a,b,c]))</syntaxhighlight>
{{trans|C}}
Otherwise, coding directly:
<syntaxhighlight lang="parigp">roots(a,b,c)={
b /= a;
c /= a;
Line 952 ⟶ 1,857:
[sol,c/sol]
)
};</
Either way,
<syntaxhighlight lang="parigp">roots(1,-1e9,1)</syntaxhighlight>
gives one root around 0.000000001000000000000000001 and one root around 999999999.999999999.
=={{header|Pascal}}==
some parts translated from Modula2
<syntaxhighlight lang="pascal">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.
</syntaxhighlight>
{{out}}
<pre>
Version 1:
x1: 1.00000000E+010, x2: 1.00000000E-010
Version 2:
x1: 1.00000000E+010, x2: 1.00000000E-010
</pre>
=={{header|Perl}}==
When using [http://perldoc.perl.org/Math/Complex.html Math::Complex] perl automatically convert numbers when necessary.
<
($x1,$x2) = solveQuad(1,2,3);
Line 967 ⟶ 1,915:
my $root = sqrt($b**2 - 4*$a*$c);
return ( -$b + $root )/(2*$a), ( -$b - $root )/(2*$a);
}</
=={{header|
{{trans|ERRE}}
<!--<syntaxhighlight lang="phix">-->
<span style="color: #008080;">procedure</span> <span style="color: #000000;">solve_quadratic</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">t3</span><span style="color: #0000FF;">)</span>
<span style="color: #004080;">atom</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">a</span><span style="color: #0000FF;">,</span><span style="color: #000000;">b</span><span style="color: #0000FF;">,</span><span style="color: #000000;">c</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">t3</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">d</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">b</span><span style="color: #0000FF;">*</span><span style="color: #000000;">b</span><span style="color: #0000FF;">-</span><span style="color: #000000;">4</span><span style="color: #0000FF;">*</span><span style="color: #000000;">a</span><span style="color: #0000FF;">*</span><span style="color: #000000;">c</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">f</span>
<span style="color: #004080;">string</span> <span style="color: #000000;">s</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sprintf</span><span style="color: #0000FF;">(</span><span style="color: #008000;">"for a=%g,b=%g,c=%g"</span><span style="color: #0000FF;">,</span><span style="color: #000000;">t3</span><span style="color: #0000FF;">),</span> <span style="color: #000000;">t</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">u</span>
<span style="color: #008080;">if</span> <span style="color: #7060A8;">abs</span><span style="color: #0000FF;">(</span><span style="color: #000000;">d</span><span style="color: #0000FF;">)<</span><span style="color: #000000;">1e-6</span> <span style="color: #008080;">then</span> <span style="color: #000000;">d</span><span style="color: #0000FF;">=</span><span style="color: #000000;">0</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">switch</span> <span style="color: #7060A8;">sign</span><span style="color: #0000FF;">(</span><span style="color: #000000;">d</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span>
<span style="color: #008080;">case</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">:</span> <span style="color: #000000;">t</span> <span style="color: #0000FF;">=</span> <span style="color: #008000;">"single root is %g"</span>
<span style="color: #000000;">u</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{-</span><span style="color: #000000;">b</span><span style="color: #0000FF;">/</span><span style="color: #000000;">2</span><span style="color: #0000FF;">/</span><span style="color: #000000;">a</span><span style="color: #0000FF;">}</span>
<span style="color: #008080;">case</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">:</span> <span style="color: #000000;">t</span> <span style="color: #0000FF;">=</span> <span style="color: #008000;">"real roots are %g and %g"</span>
<span style="color: #000000;">f</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">+</span><span style="color: #7060A8;">sqrt</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">-</span><span style="color: #000000;">4</span><span style="color: #0000FF;">*</span><span style="color: #000000;">a</span><span style="color: #0000FF;">*</span><span style="color: #000000;">c</span><span style="color: #0000FF;">/(</span><span style="color: #000000;">b</span><span style="color: #0000FF;">*</span><span style="color: #000000;">b</span><span style="color: #0000FF;">)))/</span><span style="color: #000000;">2</span>
<span style="color: #000000;">u</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{-</span><span style="color: #000000;">f</span><span style="color: #0000FF;">*</span><span style="color: #000000;">b</span><span style="color: #0000FF;">/</span><span style="color: #000000;">a</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">c</span><span style="color: #0000FF;">/</span><span style="color: #000000;">b</span><span style="color: #0000FF;">/</span><span style="color: #000000;">f</span><span style="color: #0000FF;">}</span>
<span style="color: #008080;">case</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">:</span> <span style="color: #000000;">t</span> <span style="color: #0000FF;">=</span> <span style="color: #008000;">"complex roots are %g +/- %g*i"</span>
<span style="color: #000000;">u</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{-</span><span style="color: #000000;">b</span><span style="color: #0000FF;">/</span><span style="color: #000000;">2</span><span style="color: #0000FF;">/</span><span style="color: #000000;">a</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">sqrt</span><span style="color: #0000FF;">(-</span><span style="color: #000000;">d</span><span style="color: #0000FF;">)/</span><span style="color: #000000;">2</span><span style="color: #0000FF;">/</span><span style="color: #000000;">a</span><span style="color: #0000FF;">}</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">switch</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"%-25s the %s\n"</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">s</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">sprintf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">t</span><span style="color: #0000FF;">,</span><span style="color: #000000;">u</span><span style="color: #0000FF;">)})</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span>
<span style="color: #008080;">constant</span> <span style="color: #000000;">tests</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{{</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">1E9</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1</span><span style="color: #0000FF;">},</span>
<span style="color: #0000FF;">{</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1</span><span style="color: #0000FF;">},</span>
<span style="color: #0000FF;">{</span><span style="color: #000000;">2</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">6</span><span style="color: #0000FF;">},</span>
<span style="color: #0000FF;">{</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">2</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">2</span><span style="color: #0000FF;">},</span>
<span style="color: #0000FF;">{</span><span style="color: #000000;">0.5</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1.4142135</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1</span><span style="color: #0000FF;">},</span>
<span style="color: #0000FF;">{</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">3</span><span style="color: #0000FF;">,</span><span style="color: #000000;">2</span><span style="color: #0000FF;">},</span>
<span style="color: #0000FF;">{</span><span style="color: #000000;">3</span><span style="color: #0000FF;">,</span><span style="color: #000000;">4</span><span style="color: #0000FF;">,</span><span style="color: #000000;">5</span><span style="color: #0000FF;">}}</span>
<span style="color: #7060A8;">papply</span><span style="color: #0000FF;">(</span><span style="color: #000000;">tests</span><span style="color: #0000FF;">,</span><span style="color: #000000;">solve_quadratic</span><span style="color: #0000FF;">)</span>
<!--</syntaxhighlight>-->
<pre>
for a=1,b=-1e+9,c=1 the real roots are 1e+9 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.73205 and 0.732051
for a=0.5,b=1.41421,c=1 the single root is -1.41421
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.666667 +/- 1.10554*i
</pre>
=={{header|PicoLisp}}==
<syntaxhighlight lang="picolisp">(scl 40)
(de solveQuad (A B C)
(let SD (sqrt (- (* B B) (*
(if (lt0
(*/ 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 .) )</syntaxhighlight>
{{out}}
<pre>-> ("999,999.999999" "0.000001")</pre>
=={{header|PL/I}}==
<syntaxhighlight lang="pl/i">
declare (c1, c2) float complex,
(a, b, c, x1, x2) float;
Line 1,021 ⟶ 1,995:
put data (x1, x2);
end;
</syntaxhighlight>
=={{header|
{{libheader|NumPy}}
This solution compares the naïve method with three "better" methods.
<syntaxhighlight lang="python">#!/usr/bin/env python3
import math
import cmath
import numpy
def quad_discriminating_roots(a,b,c, entier = 1e-5):
"""For reference, the naive algorithm which shows complete loss of
precision on the quadratic in question. (This function also returns a
characterization of the roots.)"""
discriminant = b*b - 4*a*c
a,b,c,d =complex(a), complex(b), complex(c), complex(discriminant)
root1 = (-b + cmath.sqrt(d))/2./a
root2 = (-b - cmath.sqrt(d))/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
def middlebrook(a, b, c):
try:
q = math.sqrt(a*c)/b
f = .5+ math.sqrt(1-4*q*q)/2
except ValueError:
q = cmath.sqrt(a*c)/b
f = .5+ cmath.sqrt(1-4*q*q)/2
return (-b/a)*f, -c/(b*f)
def whatevery(a, b, c):
try:
d = math.sqrt(b*b-4*a*c)
except ValueError:
if b > 0:
return div(2*c, (-b-d)), div((-b-d), 2*a)
else:
return div((-b+d), 2*a), div(2*c, (-b+d))
def div(n, d):
"""Divide, with a useful interpretation of division by zero."""
try:
return n/d
except ZeroDivisionError:
if n:
return n*float('inf')
return float('nan')
testcases = [
(3, 4, 4/3), # real, equal
(3, 2, -1), # real, unequal
(3, 2, 1), # complex
(1, -1e9, 1), # ill-conditioned "quadratic in question" required by task.
(1, -1e100, 1),
(1, -1e200, 1),
(1, -1e300, 1),
]
print('Naive:')
for c in testcases:
print("{} {:.5} {:.5}".format(*quad_discriminating_roots(*c)))
print('\nMiddlebrook:')
for c in testcases:
print(("{:.5} "*2).format(*middlebrook(*c)))
print('\nWhat Every...')
for c in testcases:
print(("{:.5} "*2).format(*whatevery(*c)))
print('\nNumpy:')
for c in testcases:
print(("{:.5} "*2).format(*numpy.roots(c)))</syntaxhighlight>
{{out}}
<pre>
Naive:
real and equal 0.66667 0.66667
real 0.33333 -1.0
complex (-0.33333+0.4714j) (-0.33333-0.4714j)
real 1e+09 0.0
real 1e+100 0.0
real nan nan
real nan nan
Middlebrook:
-0.66667 -0.66667
(-1+0j) (0.33333+0j)
(-0.33333-0.4714j) (-0.33333+0.4714j)
1e+09 1e-09
1e+100 1e-100
1e+200 1e-200
1e+300 1e-300
What Every...
-0.66667 -0.66667
0.33333 -1.0
(-0.33333+0.4714j) (-0.33333-0.4714j)
1e+09 1e-09
1e+100 1e-100
inf 0.0
inf 0.0
Numpy:
-0.66667 -0.66667
-1.0 0.33333
(-0.33333+0.4714j) (-0.33333-0.4714j)
1e+09 1e-09
1e+100 1e-100
1e+200 1e-200
1e+300 0.0
</pre>
=={{header|R}}==
<syntaxhighlight lang="r">qroots <- function(a, b, c) {
r <- sqrt(b * b - 4 * a * c + 0i)
} else {
z <- (-b - r) / (2 * a)
}
c(z, c / (z * a))
}
qroots(1, 0, 2i)
[1] -1+1i 1-1i
qroots(1, -1e9, 1)
[1] 1e+09+0i 1e-09+0i</syntaxhighlight>
Using the builtin '''polyroot''' function (note the order of coefficients is reversed):
<syntaxhighlight lang="r">polyroot(c(2i, 0, 1))
[1] -1+1i 1-1i
polyroot(c(1, -1e9, 1))
[1] 1e-09+0i 1e+09+0i</syntaxhighlight>
=={{header|Racket}}==
<syntaxhighlight lang="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)</syntaxhighlight>
=={{header|Raku}}==
(formerly Perl 6)
{{Works with|Rakudo|2022.12}}
''Works with previous versions also but will return slightly less precise results.''
Raku has complex number handling built in.
<syntaxhighlight lang="raku" line>for
[1, 2, 1],
[1, 2, 3],
[1, -2, 1],
[1, 0, -4],
[1, -10⁶, 1]
-> @coefficients {
printf "Roots for %d, %d, %d\t=> (%s, %s)\n",
|@coefficients, |quadroots(@coefficients);
}
sub quadroots (*[$a, $b, $c]) {
( -$b + $_ ) / (2 × $a),
( -$b - $_ ) / (2 × $a)
given
($b² - 4 × $a × $c ).Complex.sqrt.narrow
}</syntaxhighlight>
{{out}}
<pre>Roots for 1, 2, 1 => (-1, -1)
Roots for 1, 2, 3 => (-1+1.4142135623730951i, -1-1.4142135623730951i)
Roots for 1, -2, 1 => (1, 1)
Roots for 1, 0, -4 => (2, -2)
Roots for 1, -1000000, 1 => (999999.999999, 1.00000761449337e-06)</pre>
=={{header|REXX}}==
===version 1===
The REXX language doesn't have a '''sqrt''' function, nor does it support complex numbers natively.
Since "unlimited" decimal precision is part of the REXX language, the '''numeric digits''' was increased
<br>(from a default of '''9''') to '''200''' to accommodate when a root is closer to zero than the other root.
Note that only nine decimal digits (precision) are shown in the ''displaying'' of the output.
This REXX version supports ''complex numbers'' for the result.
<syntaxhighlight lang="rexx">/*REXX program finds the roots (which may be complex) of a quadratic function. */
parse arg a b c . /*obtain the specified arguments: A B C*/
call quad a,b,c /*solve quadratic function via the sub.*/
r1= r1/1; r2= r2/1; a= a/1; b= b/1; c= c/1 /*normalize numbers to a new precision.*/
if r1j\=0 then r1=r1||left('+',r1j>0)(r1j/1)"i" /*Imaginary part? Handle complex number*/
if r2j\=0 then r2=r2||left('+',r2j>0)(r2j/1)"i" /* " " " " " */
say ' a =' a /*display the normalized value of A. */
say ' b =' b /* " " " " " B. */
say ' c =' c /* " " " " " C. */
say; say 'root1 =' r1 /* " " " " 1st root*/
say 'root2 =' r2 /* " " " " 2nd root*/
exit 0 /*stick a fork in it, we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
quad: parse arg aa,bb,cc; numeric digits 200 /*obtain 3 args; use enough dec. digits*/
$= sqrt(bb**2-4*aa*cc); L= length($) /*compute SQRT (which may be complex).*/
r= 1 /(aa+aa); ?= right($, 1)=='i' /*compute reciprocal of 2*aa; Complex?*/
if ? then do; r1= -bb *r; r2=r1; r1j= left($,L-1)*r; r2j=-r1j; end
else do; r1=(-bb+$)*r; r2=(-bb-$)*r; r1j= 0; r2j= 0; end
return
/*──────────────────────────────────────────────────────────────────────────────────────*/
sqrt: procedure; parse arg x 1 ox; if x=0 then return 0; d= digits(); m.= 9; numeric form
numeric digits 9; h= d+6; x=abs(x); parse value format(x,2,1,,0) 'E0' with g 'E' _ .
g=g*.5'e'_%2; do j=0 while h>9; m.j=h; h=h%2+1; end /*j*/
do k=j+5 to 0 by -1; numeric digits m.k; g=(g+x/g)*.5; end /*k*/
numeric digits d; return (g/1)left('i', ox<0) /*make complex if OX<0. */</syntaxhighlight>
{{out|output|text= when using the input of: <tt> 1 -10e5 1 </tt>}}
<pre>
a = 1
b = -1000000
c = 1
root1 = 1000000
root2 = 0.000001
</pre>
The following output is when Regina 3.9.3 REXX is used.
{{out|output|text= when using the input of: <tt> 1 -10e9 1 </tt>}}
<pre>
a = 1
b = -1.0E+10
c = 1
root1 = 1.000000000E+10
root2 = 1E-10
</pre>
The following output is when R4 REXX is used.
{{out|output|text= when using the input of: <tt> 1 -10e9 1 </tt>}}
<pre>
a = 1
b = -1E+10
c = 1
root1 = 1E+10
root2 = 0.0000000001
</pre>
{{out|output|text= when using the input of: <tt> 3 2 1 </tt>}}
<pre>
a = 3
b = 2
c = 1
root1 = -0.333333333+0.471404521i
root2 = -0.333333333-0.471404521i
</pre>
{{out|output|text= when using the input of: <tt> 1 0 1 </tt>
<pre>
a = 1
b = 0
c = 1
root1 = 0+1i
root2 = 0-1i
</pre>
=== Version 2 ===
<syntaxhighlight lang="rexx">/* 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)
Exit</syntaxhighlight>
{{out}}
<pre>
Version 1:
a = 1
b = -1
c = 0
root1 = 1
root2 = 0
Version 2:
1*x**2-1.0000000001*x+1.e-9 = 0
x1=0.9999999991000000000025
x2=0.0000000009999999999975
</pre>
=={{header|Ring}}==
<syntaxhighlight lang="text">
x1 = 0
x2 = 0
quadratic(3, 4, 4/3.0) # [-2/3]
see "x1 = " + x1 + " x2 = " + x2 + nl
quadratic(3, 2, -1) # [1/3, -1]
see "x1 = " + x1 + " x2 = " + x2 + nl
quadratic(-2, 7, 15) # [-3/2, 5]
see "x1 = " + x1 + " x2 = " + x2 + nl
quadratic(1, -2, 1) # [1]
see "x1 = " + x1 + " x2 = " + x2 + nl
func quadratic a, b, c
sqrtDiscriminant = sqrt(pow(b,2) - 4*a*c)
x1 = (-b + sqrtDiscriminant) / (2.0*a)
x2 = (-b - sqrtDiscriminant) / (2.0*a)
return [x1, x2]
</syntaxhighlight>
=={{header|RPL}}==
RPL can solve quadratic functions directly :
'x^2-1E9*x+1' 'x' QUAD
returns
1: '(1000000000+s1*1000000000)/2'
which can then be turned into roots by storing 1 or -1 in the <code>s1</code> variable and evaluating the formula:
DUP 1 's1' STO EVAL SWAP -1 's1' STO EVAL
hence returning
2: 1000000000
1: 0
So let's implement the algorithm proposed by the task:
{| class="wikitable"
! RPL code
! Comment
|-
|
≪
'''IF''' DUP TYPE 1 == '''THEN'''
'''IF''' DUP IM NOT '''THEN''' RE '''END END'''
≫ '<span style="color:blue">REALZ</span>' STO
.
≪ → a b c
≪ '''IF''' b NOT '''THEN''' c a / NEG √ DUP NEG '''ELSE'''
a c * √ b /
1 SWAP SQ 4 * - √ 2 / 0.5 +
b * NEG
DUP a / <span style="color:blue">REALZ</span>
c ROT / <span style="color:blue">REALZ</span> '''END'''
≫ ≫ '<span style="color:blue">QROOT</span>' STO
|
<span style="color:blue">REALZ</span> ''( number → number )''
if number is a complex
with no imaginary part, then turn it into a real
<span style="color:blue">QROOT</span> ''( a b c → r1 r2 ) ''
if b=0 then roots are obvious, else
q = sqrt(a*c)/b
f = 1/2+sqrt(1-4*q^2)/2
get -b*f
root1 = -b/a*f
root2 = -c/(b*f)
|}
1 -1E9 1 <span style="color:blue">QROOT</span>
actually returns a more correct answer:
2: 1000000000
1: .000000001
=={{header|Ruby}}==
{{works with|Ruby|1.9.3+}}
The CMath#sqrt method will return a Complex instance if necessary.
<syntaxhighlight lang="ruby">require 'cmath'
def quadratic(a, b, c)
sqrt_discriminant =
[(-b + sqrt_discriminant) / (2.0*a), (-b - sqrt_discriminant) / (2.0*a)]
end
p quadratic(3, 4, 4/3.0) #
p quadratic(3, 2, -1) #
p quadratic(3, 2, 1) #
p quadratic(1, 0, 1) #
p quadratic(1, -1e6, 1) #
p quadratic(-2, 7, 15) #
p quadratic(1, -2, 1) #
p quadratic(1, 3, 3) #
{{out}}
<pre>
[-0.6666666666666666, -0.6666666666666666]
[
[(-0.3333333333333333+0.47140452079103173i), (-0.3333333333333333-0.47140452079103173i)]
[(0.0+1.0i), (0.0-1.0i)]
[999999.999999, 1.00000761449337e-06]
[-1.5, 5.0]
[1.0, 1.0]
[
</pre>
=={{header|Run BASIC}}==
<syntaxhighlight lang="runbasic">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</syntaxhighlight><pre>FOR 1,2,3 => -1 +i1.41421356 , -1 -i1.41421356
FOR 4,5,6 => -0.625 +i1.05326872 , -0.625 -i1.05326872</pre>
=={{header|Scala}}==
Using [[Arithmetic/Complex#Scala|Complex]] class from task Arithmetic/Complex.
<
object QuadraticRoots {
def solve(a:Double, b:Double, c:Double)={
Line 1,222 ⟶ 2,483:
}
}
}</
Usage:
<
(1.0, 22.0, -1323.0), // two distinct real roots
(6.0, -23.0, 20.0), // with a != 1.0
Line 1,240 ⟶ 2,501:
if(roots._1 != roots._2) println("x2="+roots._2)
println
}</
{{out}}
<pre>a=1.00000 b=22.0000 c=-1323.00
x1=-49.0
Line 1,266 ⟶ 2,527:
=={{header|Scheme}}==
<
(if (= a 0)
(if (= b 0) 'fail (- (/ c b)))
Line 1,305 ⟶ 2,566:
(quadratic 1 -1e5 1)
; (99999.99999 1.0000000001000001e-05)</
=={{header|Seed7}}==
{{trans|Ada}}
<syntaxhighlight lang="seed7">$ 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;</syntaxhighlight>
{{out}}
<pre>
X1 = 1000000.000000 X2 = 0.000001
</pre>
=={{header|Sidef}}==
<syntaxhighlight lang="ruby">var sets = [
[1, 2, 1],
[1, 2, 3],
[1, -2, 1],
[1, 0, -4],
[1, -1e6, 1],
]
func quadroots(a, b, c) {
var root = sqrt(b**2 - 4*a*c)
[(-b + root) / (2 * a),
(-b - root) / (2 * a)]
}
sets.each { |coefficients|
say ("Roots for #{coefficients}",
"=> (#{quadroots(coefficients...).join(', ')})")
}</syntaxhighlight>
{{out}}
<pre>
Roots for [1, 2, 1]=> (-1, -1)
Roots for [1, 2, 3]=> (-1+1.41421356237309504880168872420969807856967187538i, -1-1.41421356237309504880168872420969807856967187538i)
Roots for [1, -2, 1]=> (1, 1)
Roots for [1, 0, -4]=> (2, -2)
Roots for [1, -1000000, 1]=> (999999.999998999999999998999999999997999999999995, 0.00000100000000000100000000000200000000000500000000002)
</pre>
=={{header|Stata}}==
<syntaxhighlight lang="stata">mata
: polyroots((-2,0,1))
1 2
+-----------------------------+
1 | 1.41421356 -1.41421356 |
+-----------------------------+
: polyroots((2,0,1))
1 2
+-------------------------------+
1 | -1.41421356i 1.41421356i |
+-------------------------------+</syntaxhighlight>
=={{header|Tcl}}==
{{tcllib|math::complexnumbers}}
<
namespace import math::complexnumbers::complex math::complexnumbers::tostring
Line 1,345 ⟶ 2,692:
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)}</
{{out}}
<pre>3x**2 + 4x + 1.3333333333333333 = 0
x = -0.6666666666666666
Line 1,372 ⟶ 2,719:
=={{header|TI-89 BASIC}}==
TI-89 BASIC has built-in numeric and algebraic solvers.
<syntaxhighlight lang="text">solve(x^2-1E9x+1.0)</syntaxhighlight>
returns
<pre>x=1.E-9 or x=1.E9</pre>
=={{header|Wren}}==
{{trans|Go}}
{{libheader|Wren-complex}}
<syntaxhighlight lang="wren">import "./complex" for Complex
var quadratic = Fn.new { |a, b, c|
var d = b*b - 4*a*c
if (d == 0) {
// single root
return [[-b/(2*a)], null]
}
if (d > 0) {
// two real roots
var sr = d.sqrt
d = (b < 0) ? sr - b : -sr - b
return [[d/(2*a), 2*c/d], null]
}
// two complex roots
var den = 1 / (2*a)
var t1 = Complex.new(-b*den, 0)
var t2 = Complex.new(0, (-d).sqrt * den)
return [[], [t1+t2, t1-t2]]
}
var test = Fn.new { |a, b, c|
System.write("coefficients: %(a), %(b), %(c) -> ")
var roots = quadratic.call(a, b, c)
var r = roots[0]
if (r.count == 1) {
System.print("one real root: %(r[0])")
} else if (r.count == 2) {
System.print("two real roots: %(r[0]) and %(r[1])")
} else {
var i = roots[1]
System.print("two complex roots: %(i[0]) and %(i[1])")
}
}
var coeffs = [
[1, -2, 1],
[1, 0, 1],
[1, -10, 1],
[1, -1000, 1],
[1, -1e9, 1]
]
for (c in coeffs) test.call(c[0], c[1], c[2])</syntaxhighlight>
{{out}}
<pre>
coefficients: 1, -2, 1 -> one real root: 1
coefficients: 1, 0, 1 -> two complex roots: 0 + i and 0 - i
coefficients: 1, -10, 1 -> two real roots: 9.8989794855664 and 0.10102051443364
coefficients: 1, -1000, 1 -> two real roots: 999.998999999 and 0.001000001000002
coefficients: 1, -1000000000, 1 -> two real roots: 1000000000 and 1e-09
</pre>
=={{header|XPL0}}==
{{trans|Go}}
<syntaxhighlight lang "XPL0">include xpllib; \for Print
func real QuadRoots(A, B, C); \Return roots of quadratic equation
real A, B, C;
real D, E, R;
[R:= [0., 0., 0.];
R(0):= 0.; R(1):= 0.; R(2):= 0.;
D:= B*B - 4.*A*C;
case of
D = 0.: [R(0):= -B / (2.*A); \single root
R(1):= R(0);
];
D > 0.: [if B < 0. then \two real roots
E:= sqrt(D) - B
else E:= -sqrt(D) - B;
R(0):= E / (2.*A);
R(1):= 2. * C / E;
];
D < 0.: [R(0):= -B / (2.*A); \real
R(2):= sqrt(-D) /(2.*A); \imaginary
]
other []; \D overflowed or a coefficient was NaN
return R;
];
func Test(A, B, C);
real A, B, C;
real R;
[Print("coefficients: %g, %g, %g -> ", A, B, C);
R:= QuadRoots(A, B, C);
if R(2) # 0. then
Print("two complex roots: %g+%gi, %g-%gi\n", R(0), R(2), R(0), R(2))
else [if R(0) = R(1) then
Print("one real root: %g\n", R(0))
else Print("two real roots: %15.15g, %15.15g\n", R(0), R(1));
];
];
real C; int I;
[C:= [ [1., -2., 1.],
[1., 0., 1.],
[1., -10., 1.],
[1., -1000., 1.],
[1., -1e9, 1.],
[1., -4., 6.] ];
for I:= 0 to 5 do
Test(C(I,0), C(I,1), C(I,2));
]</syntaxhighlight>
{{out}}
<pre>
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.89897948556636, 0.101020514433644
coefficients: 1, -1000, 1 -> two real roots: 999.998999999, 0.001000001000002
coefficients: 1, -1e9, 1 -> two real roots: 1000000000, 0.000000001
coefficients: 1, -4, 6 -> two complex roots: 2+1.41421i, 2-1.41421i
</pre>
=={{header|zkl}}==
zkl doesn't have a complex number package.
{{trans|Elixir}}
<syntaxhighlight lang="zkl">fcn quadratic(a,b,c){ b=b.toFloat();
println("Roots of a quadratic function %s, %s, %s".fmt(a,b,c));
d,a2:=(b*b - 4*a*c), a+a;
if(d>0){
sd:=d.sqrt();
println(" the real roots are %s and %s".fmt((-b + sd)/a2,(-b - sd)/a2));
}
else if(d==0) println(" the single root is ",-b/a2);
else{
sd:=(-d).sqrt();
println(" the complex roots are %s and \U00B1;%si".fmt(-b/a2,sd/a2));
}
}</syntaxhighlight>
<syntaxhighlight lang="zkl">foreach a,b,c in (T( T(1,-2,1), T(1,-3,2), T(1,0,1), T(1,-1.0e10,1), T(1,2,3), T(2,-1,-6)) ){
quadratic(a,b,c)
}</syntaxhighlight>
{{out}}
<pre>
Roots of a quadratic function 1, -2, 1
the single root is 1
Roots of a quadratic function 1, -3, 2
the real roots are 2 and 1
Roots of a quadratic function 1, 0, 1
the complex roots are 0 and ±1i
Roots of a quadratic function 1, -1e+10, 1
the real roots are 1e+10 and 0
Roots of a quadratic function 1, 2, 3
the complex roots are -1 and ±1.41421i
Roots of a quadratic function 2, -1, -6
the real roots are 2 and -1.5
</pre>
{{omit from|M4}}
| |||