Rosetta Code

Roots of a quadratic function: Difference between revisions

Page Discussion

Roots of a quadratic function (view source)

Revision as of 12:55, 28 August 2022

7,377 bytes added ,  1 year ago
m
syntax highlighting fixup automation
m (syntax highlighting fixup automation)
Line 1:
{{task|Arithmetic operations}}
{{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>.
 
;Task:
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.
Create a program that finds and outputs the roots of a given function, range and (if applicable) step width.
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]]:
<lang ada>with Ada.Text_IO; use Ada.Text_IO;
with Ada.Numerics.Elementary_Functions; use Ada.Numerics.Elementary_Functions;
 
The program should identify whether the root is exact or approximate.
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>
{{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.
 
For this task, use: &nbsp; &nbsp; <big><big> ƒ(x) &nbsp; = &nbsp; x<sup>3</sup> - 3x<sup>2</sup> + 2x </big></big>
Suggested by Middlebrook (D-OA), a better numerical method: to define two parameters <math> q = \sqrt{a c} / b </math> and <math> f = 1/2 + \sqrt{1 - 4 q^2} /2 </math>
<br><br>
 
=={{header|11l}}==
and the two roots of the quardratic are: <math> \frac{-b}{a} f </math> and <math> \frac{-c}{b f} </math>
{{trans|Python}}
 
<syntaxhighlight lang="11l">F f(x)
R x^3 - 3 * x^2 + 2 * x
 
-V step = 0.001
'''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
-V start = -1.0
[https://www.validlab.com/goldberg/paper.pdf "What Every Scientist Should Know About Floating-Point Arithmetic"] for a possible algorithm.
-V stop = 3.0
 
V sgn = f(start) > 0
=={{header|11l}}==
V x = start
<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))
 
L x <= stop
V testcases = [(3.0, 4.0, 4 / 3),
V value = f(3.0, 2.0, -1.0x),
(3.0, 2.0, 1.0),
(1.0, -1e9, 1.0),
(1.0, -1e100, 1.0)]
 
I value == 0
L(a, b, c) testcases
print(‘Root found at ’x)
V (r1, r2) = quad_roots(a, b, c)
E I (value > 0) != sgn
print(r1, end' ‘ ’)
print(r2‘Root found near ’x)</lang>
 
sgn = value > 0
x += step</syntaxhighlight>
 
{{out}}
<pre>
Root found near 8.812395258e-16
-0.666667+0i -0.666667+0i
Root found near 1
0.333333+0i -1+0i
Root found near 2.001
-0.333333+0.471405i -0.333333-0.471405i
1e+09+0i 0i
1e+100+0i 0i
</pre>
 
=={{header|Ada}}==
<syntaxhighlight lang ="ada">with Ada.Text_IOText_Io; use Ada.Text_IOText_Io;
with Ada.Numerics.Elementary_Functions; use Ada.Numerics.Elementary_Functions;
procedure Roots_Of_Function is
 
package Real_Io is new Ada.Text_Io.Float_Io(Long_Float);
procedure Quadratic_Equation is
use Real_Io;
type Roots is array (1..2) of Float;
function Solve (A, B, C : Float) return Roots is
function F(X : Long_Float) return Long_Float is
SD : constant Float := sqrt (B**2 - 4.0 * A * C);
X : Float;
begin
ifreturn B(X**3 <- 03.0*X*X then+ 2.0*X);
end F;
X := (- B + SD) / (2.0 * A);
return (X, C / (A * X));
Step : constant Long_Float := 1.0E-6;
else
Start : constant XLong_Float := (- B - SD) / (21.0 * A);
Stop : constant Long_Float := 3.0;
return (C / (A * X), X);
Value : Long_Float := F(Start);
Sign : Boolean := Value > 0.0;
X : Long_Float := Start + Step;
begin
if Value = 0.0 then
Put("Root found at ");
Put(Item => Start, Fore => 1, Aft => 6, Exp => 0);
New_Line;
end if;
while X <= Stop loop
Value := F(X);
if (Value > 0.0) /= Sign then
Put("Root found near ");
Put(Item => X, Fore => 1, Aft => 6, Exp => 0);
New_Line;
elsif Value = 0.0 then
Put("Root found at ");
Put(Item => X, Fore => 1, Aft => 6, Exp => 0);
New_Line;
end if;
Sign := Value > 0.0;
end Solve;
X := X + Step;
 
end loop;
R : constant Roots := Solve (1.0, -10.0E5, 1.0);
end Roots_Of_Function;</syntaxhighlight>
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.
{{out}}
<pre>
X1 = 1.00000E+06 X2 = 1.00000E-06
</pre>
 
=={{header|ALGOL 68}}==
{{trans|Ada}}
 
{{works with|ALGOL 68|Revision 1 - no extensions to language used}}
 
{{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] - due probablyto need toextensive "USE"use theof complFORMATted ENVIRONtransput}}
Finding 3 roots using the secant method:
<lang algol68>quadratic equation:
<syntaxhighlight lang="algol68">MODE DBL = LONG REAL;
BEGIN
FORMAT dbl = $g(-long real width, long real width-6, -2)$;
 
MODE ROOTS XY = UNIONSTRUCT([]REALDBL x, []COMPLy);
FORMAT xy root = $f(dbl)" ("b("Exactly", "Approximately")")"$;
MODE QUADRATIC = STRUCT(REAL a,b,c);
 
MODE DBLOPT = UNION(DBL, VOID);
PROC solve = (QUADRATIC q)ROOTS:
MODE XYRES = UNION(XY, VOID);
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;
 
PROC find root = (PROC (DBL)DBL f, DBLOPT in x1, in x2, in x error, in y error)XYRES:(
# only a very tiny difference between the 2 examples #
INT limit = ENTIER (long real width / log(2)); # worst case of a binary search) #
[]QUADRATIC test = ((1, -10e5, 1), (1, 0, 1), (1,-3,2), (1,3,2), (4,0,4), (3,4,5));
DBL x1 := (in x1|(DBL x1):x1|-5.0), # if x1 is EMPTY then -5.0 #
x2 := (in x2|(DBL x2):x2|+5.0),
FORMAT real fmt = $g(-0,8)$;
x error := (in x error|(DBL x error):x error|small real),
FORMAT compl fmt = $f(real fmt)"+"f(real fmt)"i"$;
y error := (in y error|(DBL y error):y error|small real);
FORMAT quadratic fmt = $f(real fmt)" x**2 + "f(real fmt)" x + "f(real fmt)" = 0"$;
DBL y1 := f(x1), y2;
DBL dx := x1 - x2, dy;
 
FORIF indexy1 TO= UPB0 test DOTHEN
XY(x1, y1) # we already have a solution! #
QUADRATIC quadratic = test[index];
ELSE
ROOTS r = solve(quadratic);
FOR i WHILE
y2 := f(x2);
# Output the two different scenerios #
IF y2 = 0 THEN stop iteration FI;
printf(($"Quadratic: "$, quadratic fmt, quadratic, $l$));
IF i = limit THEN value error FI;
CASE r IN
IF y1 = y2 THEN value error FI;
([]REAL r):
dy printf(($"REAL x1 := "$, realy1 fmt,- r[1],y2;
dx $", x2 := "$,dx real/ fmt,dy r[2],* $"y2; "$,
x1 := x2; y1 := y2; # retain for next iteration #
$"REAL y1 = "$, real fmt, real evaluate(quadratic,r[1]),
x2 -:= dx;
$", y2 = "$, real fmt, real evaluate(quadratic,r[2]), $";"ll$
# WHILE # ABS dx > x error AND ABS dy > y error DO
)),
([]COMPL c):SKIP
OD;
printf(($"COMPL x1,x2 = "$, real fmt, re OF c[1], $"+/-"$,
stop iteration:
real fmt, ABS im OF c[1], $"; "$,
XY(x2, y2) EXIT
$"COMPL y1 = "$, compl fmt, compl evaluate(quadratic,c[1]),
value error:
$", y2 = "$, compl fmt, compl evaluate(quadratic,c[2]), $";"ll$
))EMPTY
ESACFI
);
OD
END # quadratic_equation #</lang>
{{out}}
<pre>
Quadratic: 1.00000000 x**2 + -1000000.00000000 x + 1.00000000 = 0
REAL x1 = 999999.99999900, x2 = .00000100; REAL y1 = -.00000761, y2 = -.00000761;
 
PROC f = (DBL x)DBL: x UP 3 - LONG 3.1 * x UP 2 + LONG 2.0 * x;
Quadratic: 1.00000000 x**2 + .00000000 x + 1.00000000 = 0
COMPL x1,x2 = .00000000+/-1.00000000; COMPL y1 = .00000000+.00000000i, y2 = .00000000+.00000000i;
 
DBL first root, second root, third root;
Quadratic: 1.00000000 x**2 + -3.00000000 x + 2.00000000 = 0
REAL x1 = 2.00000000, x2 = 1.00000000; REAL y1 = .00000000, y2 = .00000000;
 
XYRES first result = find root(f, LENG -1.0, LENG 3.0, EMPTY, EMPTY);
Quadratic: 1.00000000 x**2 + 3.00000000 x + 2.00000000 = 0
CASE first result IN
REAL x1 = -2.00000000, x2 = -1.00000000; REAL y1 = .00000000, y2 = .00000000;
(XY first result): (
printf(($"1st root found at x = "f(xy root)l$, x OF first result, y OF first result=0));
first root := x OF first result
)
OUT printf($"No first root found"l$); stop
ESAC;
 
XYRES second result = find root( (DBL x)DBL: f(x) / (x - first root), EMPTY, EMPTY, EMPTY, EMPTY);
Quadratic: 4.00000000 x**2 + .00000000 x + 4.00000000 = 0
CASE second result IN
COMPL x1,x2 = .00000000+/-1.00000000; COMPL y1 = .00000000+.00000000i, y2 = .00000000+.00000000i;
(XY second result): (
printf(($"2nd root found at x = "f(xy root)l$, x OF second result, y OF second result=0));
second root := x OF second result
)
OUT printf($"No second root found"l$); stop
ESAC;
 
XYRES third result = find root( (DBL x)DBL: f(x) / (x - first root) / ( x - second root ), EMPTY, EMPTY, EMPTY, EMPTY);
Quadratic: 3.00000000 x**2 + 4.00000000 x + 5.00000000 = 0
CASE third result IN
COMPL x1,x2 = -.66666667+/-1.10554160; COMPL y1 = .00000000+.00000000i, y2 = .00000000+-.00000000i;
(XY third result): (
printf(($"3rd root found at x = "f(xy root)l$, x OF third result, y OF third result=0));
third root := x OF third result
)
OUT printf($"No third root found"l$); stop
ESAC</syntaxhighlight>
Output:
<pre>1st root found at x = 9.1557112297752398099031e-1 (Approximately)
2nd root found at x = 2.1844288770224760190097e 0 (Approximately)
3rd root found at x = 0.0000000000000000000000e 0 (Exactly)
</pre>
 
=={{header|ATS}}==
<syntaxhighlight lang="ats">
#include
"share/atspre_staload.hats"
 
typedef d = double
 
fun
findRoots
(
start: d, stop: d, step: d, f: (d) -> d, nrts: int, A: d
) : void = (
//
if
start < stop
then let
val A2 = f(start)
var nrts: int = nrts
val () =
if A2 = 0.0
then (
nrts := nrts + 1;
$extfcall(void, "printf", "An exact root is found at %12.9f\n", start)
) (* end of [then] *)
// end of [if]
val () =
if A * A2 < 0.0
then (
nrts := nrts + 1;
$extfcall(void, "printf", "An approximate root is found at %12.9f\n", start)
) (* end of [then] *)
// end of [if]
in
findRoots(start+step, stop, step, f, nrts, A2)
end // end of [then]
else (
if nrts = 0
then $extfcall(void, "printf", "There are no roots found!\n")
// end of [if]
) (* end of [else] *)
//
) (* end of [findRoots] *)
 
(* ****** ****** *)
 
implement
main0 () =
findRoots (~1.0, 3.0, 0.001, lam (x) => x*x*x - 3.0*x*x + 2.0*x, 0, 0.0)
</syntaxhighlight>
 
=={{header|AutoHotkey}}==
Poly(x) is a test function of one variable, here we are searching for its roots:
ahk forum: [http://www.autohotkey.com/forum/viewtopic.php?p=276617#276617 discussion]
* roots() searches for intervals within given limits, shifted by a given “step”, where our function has different signs at the endpoints.
<lang AutoHotkey>MsgBox % quadratic(u,v, 1,-3,2) ", " u ", " v
* Having found such an interval, the root() function searches for a value where our function is 0, within a given tolerance.
MsgBox % quadratic(u,v, 1,3,2) ", " u ", " v
* It also sets ErrorLevel to info about the root found.
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
 
[http://www.autohotkey.com/forum/viewtopic.php?t=44657&postdays=0&postorder=asc&start=139 discussion]
quadratic(ByRef x1, ByRef x2, a,b,c) { ; -> #real roots {x1,x2} of ax²+bx+c
<syntaxhighlight lang="autohotkey">MsgBox % roots("poly", -0.99, 2, 0.1, 1.0e-5)
If (a = 0)
MsgBox % roots("poly", -1, 3, 0.1, 1.0e-5)
Return -1 ; ERROR: not quadratic
 
d := b*b - 4*a*c
roots(f,x1,x2,step,tol) { ; search for roots in intervals of length "step", within tolerance "tol"
If (d < 0) {
x := x1, y := x2%f%(x), s := ""(y>0)-(y<0)
Loop % ceil((x2-x1)/step) {
Return 0
x += step, y := %f%(x), t := (y>0)-(y<0)
If (s=0 || s!=t)
res .= root(f, x-step, x, tol) " [" ErrorLevel "]`n"
s := t
}
Sort res, UN ; remove duplicate endpoints
If (d = 0) {
Return res
x1 := x2 := -b/2/a
}
Return 1
 
root(f,x1,x2,d) { ; find x in [x1,x2]: f(x)=0 within tolerance d, by bisection
If (!y1 := %f%(x1))
Return x1, ErrorLevel := "Exact"
If (!y2 := %f%(x2))
Return x2, ErrorLevel := "Exact"
If (y1*y2>0)
Return "", ErrorLevel := "Need different sign ends!"
Loop {
x := (x2+x1)/2, y := %f%(x)
If (y = 0 || x2-x1 < d)
Return x, ErrorLevel := y ? "Approximate" : "Exact"
If ((y>0) = (y1>0))
x1 := x, y1 := y
Else
x2 := x, y2 := y
}
}
x1 := (-b - (b<0 ? -sqrt(d) : sqrt(d)))/2/a
 
x2 := c/a/x1
poly(x) {
Return 2
Return ((x-3)*x+2)*x
}</lang>
}</syntaxhighlight>
 
=={{header|Axiom}}==
Using a polynomial solver:
<syntaxhighlight lang="axiom">expr := x^3-3*x^2+2*x
solve(expr,x)</syntaxhighlight>
Output:
<syntaxhighlight lang="axiom"> (1) [x= 2,x= 1,x= 0]
Type: List(Equation(Fraction(Polynomial(Integer))))</syntaxhighlight>
Using the secant method in the interpreter:
<syntaxhighlight lang="axiom">digits(30)
secant(eq: Equation Expression Float, binding: SegmentBinding(Float)):Float ==
eps := 1.0e-30
expr := lhs eq - rhs eq
x := variable binding
seg := segment binding
x1 := lo seg
x2 := hi seg
fx1 := eval(expr, x=x1)::Float
abs(fx1)<eps => return x1
for i in 1..100 repeat
fx2 := eval(expr, x=x2)::Float
abs(fx2)<eps => return x2
(x1, fx1, x2) := (x2, fx2, x2 - fx2 * (x2 - x1) / (fx2 - fx1))
error "Function not converging."</syntaxhighlight>
The example can now be called using:
<syntaxhighlight lang="axiom">secant(expr=0,x=-0.5..0.5)</syntaxhighlight>
 
=={{header|BBC BASIC}}==
<syntaxhighlight lang ="bbcbasic"> FOR test%function$ = 1 TO 7"x^3-3*x^2+2*x"
rangemin = READ a$, b$, c$-1
rangemax = 3
PRINT "For a = " ; a$ ", b = " ; b$ ", c = " ; c$ TAB(32) ;
stepsize = 0.001
PROCsolvequadratic(EVAL(a$), EVAL(b$), EVAL(c$))
NEXTaccuracy = 1E-8
PROCroots(function$, rangemin, rangemax, stepsize, accuracy)
END
DATADEF 1PROCroots(func$, -1E9min, 1max, inc, eps)
DATALOCAL 1x, 0sign%, 1oldsign%
DATAoldsign% 2,= -1, -60
DATAFOR 1,x 2,= -2min TO max STEP inc
DATA 0.5, SQRsign% = SGN(2EVAL(func$)), 1
DATA 1, 3,IF 2sign% = 0 THEN
PRINT "Root found at x = "; x
DATA 3, 4, 5
sign% = -oldsign%
ELSE IF sign% <> oldsign% AND oldsign% <> 0 THEN
DEF PROCsolvequadratic(a, b, c)
LOCAL d, f IF inc < eps THEN
PRINT "Root found near x = "; x
d = b^2 - 4*a*c
CASE SGN(d) OF ELSE
PROCroots(func$, x-inc, x+inc/8, inc/8, eps)
WHEN 0:
PRINT "the single root is " ; -b/2/aENDIF
WHEN +1: ENDIF
f = (1 + SQR(1-4*a*c/b^2))/2ENDIF
oldsign% = sign%
PRINT "the real roots are " ; -f*b/a " and " ; -c/b/f
NEXT WHEN -1:x
ENDPROC</syntaxhighlight>
PRINT "the complex roots are " ; -b/2/a " +/- " ; SQR(-d)/2/a "*i"
Output:
ENDCASE
<pre>Root found near x = 2.29204307E-9
ENDPROC</lang>
Root found near x = 1
{{out}}
Root found at x = 2</pre>
<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
For a = 2, b = -1, c = -6 the real roots are 2 and -1.5
For a = 1, b = 2, c = -2 the real roots are -2.73205081 and 0.732050808
For a = 0.5, b = SQR(2), c = 1 the single root is -1.41421356
For a = 1, b = 3, c = 2 the real roots are -2 and -1
For a = 3, b = 4, c = 5 the complex roots are -0.666666667 +/- 1.1055416*i</pre>
 
=={{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)
<lang C>#include <stdio.h>
#include <stdlib.h>
#include <complex.h>
#include <math.h>
 
=== Secant Method ===
typedef double complex cplx;
 
<syntaxhighlight lang="c">#include <math.h>
void quad_root
#include <stdio.h>
(double a, double b, double c, cplx * ra, cplx *rb)
 
double f(double x)
{
return x*x*x-3.0*x*x +2.0*x;
double d, e;
}
if (!a) {
*ra = b ? -c / b : 0;
*rb = 0;
return;
}
if (!c) {
*ra = 0;
*rb = -b / a;
return;
}
 
double secant( double xA, double xB, double(*f)(double) )
b /= 2;
{
if (fabs(b) > fabs(c)) {
e = 1 - (adouble /e b)= * (c / b)1.0e-12;
double fA, fB;
d = sqrt(fabs(e)) * fabs(b);
double d;
} else {
e = (c > 0)int ? a : -ai;
e = b * (bint /limit fabs(c)) -= e50;
d = sqrt(fabs(e)) * sqrt(fabs(c));
}
 
fA=(*f)(xA);
if (e < 0) {
for (i=0; i<limit; i++) {
e = fabs(d / a);
fB=(*f)(xB);
d = -b / a;
*ra = d += (xB - xA) / (fB - IfA) * efB;
*rb = d - I * if (fabs(d) < e;)
break;
return;
xA = xB;
}
fA = fB;
 
d = (b >= 0) ? d : xB -= d;
e = (d - b) / a;}
d = e ? if (c / ei==limit) / a : 0;{
printf("Function is not converging near (%7.4f,%7.4f).\n", xA,xB);
*ra = d;
return -99.0;
*rb = e;
}
return;
return xB;
}
 
int main(int argc, char *argv[])
{
double step = 1.0e-2;
cplx ra, rb;
double e = 1.0e-12;
quad_root(1, 1e12 + 1, 1e12, &ra, &rb);
double x = -1.032; // just so we use secant method
printf("(%g + %g i), (%g + %g i)\n",
double xx, value;
creal(ra), cimag(ra), creal(rb), cimag(rb));
 
int s = (f(x)> 0.0);
quad_root(1e300, -1e307 + 1, 1e300, &ra, &rb);
printf("(%g + %g i), (%g + %g i)\n",
creal(ra), cimag(ra), creal(rb), cimag(rb));
 
while (x < 3.0) {
return 0;
value = f(x);
}</lang>
if (fabs(value) < e) {
{{out}}<pre>(-1e+12 + 0 i), (-1 + 0 i)
printf("Root found at x= %12.9f\n", x);
(1.00208e+07 + 0 i), (9.9792e-08 + 0 i)</pre>
s = (f(x+.0001)>0.0);
}
else if ((value > 0.0) != s) {
xx = secant(x-step, x,&f);
if (xx != -99.0) // -99 meaning secand method failed
printf("Root found at x= %12.9f\n", xx);
else
printf("Root found near x= %7.4f\n", x);
s = (f(x+.0001)>0.0);
}
x += step;
}
return 0;
}</syntaxhighlight>
 
=== GNU Scientific Library ===
<lang c>#include <stdio.h>
#include <math.h>
#include <complex.h>
 
<syntaxhighlight lang="c">#include <gsl/gsl_poly.h>
void roots_quadratic_eq(double a, double b, double c, complex double *x)
#include <stdio.h>
 
int main(int argc, char *argv[])
{
/* 0 + 2x - 3x^2 + 1x^3 */
double delta;
double p[] = {0, 2, -3, 1};
double z[6];
gsl_poly_complex_workspace *w = gsl_poly_complex_workspace_alloc(4);
gsl_poly_complex_solve(p, 4, w, z);
gsl_poly_complex_workspace_free(w);
 
for(int i = 0; i < 3; ++i)
printf("%.12f\n", z[2 * i]);
 
return 0;
}</syntaxhighlight>
 
One can also use the GNU Scientific Library to find roots of functions. Compile with <pre>gcc roots.c -lgsl -lcblas -o roots</pre>
 
=={{header|C sharp|C#}}==
 
delta = b*b - 4.0*a*c;
x[0] = (-b + csqrt(delta)) / (2.0*a);
x[1] = (-b - csqrt(delta)) / (2.0*a);
}</lang>
{{trans|C++}}
 
<lang c>void roots_quadratic_eq2(double a, double b, double c, complex double *x)
<syntaxhighlight lang="csharp">using System;
 
class Program
{
public static void Main(string[] args)
b /= a;
c /= a;{
Func<double, deltadouble> f = bx => { return x *b x * x - 43 *c x * x + 2 * x; };
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>
 
double step = 0.001; // Smaller step values produce more accurate and precise results
<lang c>int main()
double start = -1;
double stop = 3;
double value = f(start);
int sign = (value > 0) ? 1 : 0;
// Check for root at start
if (value == 0)
Console.WriteLine("Root found at {0}", start);
 
for (var x = start + step; x <= stop; x += step)
{
value = f(x);
if (((value > 0) ? 1 : 0) != sign)
// We passed a root
Console.WriteLine("Root found near {0}", x);
else if (value == 0)
// We hit a root
Console.WriteLine("Root found at {0}", x);
// Update our sign
sign = (value > 0) ? 1 : 0;
}
}
}</syntaxhighlight>
 
{{trans|Java}}
<syntaxhighlight lang="csharp">using System;
 
class Program
{
complex private static int Sign(double x[2];)
{
return x < 0.0 ? -1 : x > 0.0 ? 1 : 0;
}
 
public static void PrintRoots(Func<double, double> f, double lowerBound,
roots_quadratic_eq(1, -1e20, 1, x);
double upperBound, double step)
printf("x1 = (%.20le, %.20le)\nx2 = (%.20le, %.20le)\n\n",
{
creal(x[0]), cimag(x[0]),
double x = lowerBound, ox = x;
creal(x[1]), cimag(x[1]));
double y = f(x), oy = y;
roots_quadratic_eq2(1, -1e20, 1, x);
int s = Sign(y), os = s;
printf("x1 = (%.20le, %.20le)\nx2 = (%.20le, %.20le)\n\n",
creal(x[0]), cimag(x[0]),
creal(x[1]), cimag(x[1]));
 
for (; x <= upperBound; x += step)
return 0;
{
}</lang>
s = Sign(y = f(x));
if (s == 0)
{
Console.WriteLine(x);
}
else if (s != os)
{
var dx = x - ox;
var dy = y - oy;
var cx = x - dx * (y / dy);
Console.WriteLine("~{0}", cx);
}
 
ox = x;
<pre>x1 = (1.00000000000000000000e+20, 0.00000000000000000000e+00)
oy = y;
x2 = (0.00000000000000000000e+00, 0.00000000000000000000e+00)
os = s;
}
}
 
public static void Main(string[] args)
x1 = (1.00000000000000000000e+20, 0.00000000000000000000e+00)
{
x2 = (9.99999999999999945153e-21, 0.00000000000000000000e+00)</pre>
Func<double, double> f = x => { return x * x * x - 3 * x * x + 2 * x; };
PrintRoots(f, -1.0, 4, 0.002);
}
}</syntaxhighlight>
 
===Brent's Method===
=={{header|C sharp|C#}}==
 
<lang csharp>using System;
{{trans|C++}}
using System.Numerics;
<syntaxhighlight lang="csharp">using System;
 
class QuadraticRootsProgram
{
public static void Main(string[] args)
static Tuple<Complex, Complex> Solve(double a, double b, double c)
{
varFunc<double, qdouble> f = -(bx +=> Math.Sign(b){ return x * Complex.Sqrt(bx * bx - 43 * ax * c))x /+ 2 * x; };
returndouble Tuple.Createroot = BrentsFun(qf, /lower: a-1.0, cupper: 4, tol: 0.002, /maxIter: q100);
}
 
private static void MainSwap<T>(ref T a, ref T b)
{
Console.WriteLine(Solve(1,var -1E20,tmp 1))= a;
a = b;
b = tmp;
}
}</lang>
{{out}}
<pre>((1E+20, 0), (1E-20, 0))</pre>
 
public static double BrentsFun(Func<double, double> f, double lower, double upper, double tol, uint maxIter)
=={{header|C++}}==
{
<lang cpp>#include <iostream>
double a = lower;
#include <utility>
double b = upper;
#include <complex>
double fa = f(a); // calculated now to save function calls
double fb = f(b); // calculated now to save function calls
double fs;
 
if (!(fa * fb < 0))
typedef std::complex<double> complex;
throw new ArgumentException("Signs of f(lower_bound) and f(upper_bound) must be opposites");
 
if (Math.Abs(fa) < Math.Abs(b)) // if magnitude of f(lower_bound) is less than magnitude of f(upper_bound)
std::pair<complex, complex>
{
solve_quadratic_equation(double a, double b, double c)
Swap(ref a, ref b);
{
Swap(ref fa, ref fb);
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 c = a; // c now equals the largest magnitude of the lower and upper bounds
double root = std::sqrt(discriminant);
double fc = fa; // precompute function evalutation for point c by assigning it the same value as fa
double solution1 = (b > 0)? (-b - root)/2
bool mflag = true; // boolean flag used to evaluate if statement later : (-b + root)/2;on
double s = 0; // Our Root that will be returned
double d = 0; // Only used if mflag is unset (mflag == false)
 
for (uint iter = 1; iter < maxIter; ++iter)
return std::make_pair(solution1, c/solution1);
{
// stop if converged on root or error is less than tolerance
if (Math.Abs(b - a) < tol)
{
Console.WriteLine("After {0} iterations the root is: {1}", iter, s);
return s;
} // end if
 
if (fa != fc && fb != fc)
{
// use inverse quadratic interopolation
s = (a * fb * fc / ((fa - fb) * (fa - fc)))
+ (b * fa * fc / ((fb - fa) * (fb - fc)))
+ (c * fa * fb / ((fc - fa) * (fc - fb)));
}
else
{
// secant method
s = b - fb * (b - a) / (fb - fa);
}
 
// checks to see whether we can use the faster converging quadratic && secant methods or if we need to use bisection
if ( ( (s < (3 * a + b) * 0.25) || (s > b)) ||
( mflag && (Math.Abs(s - b) >= (Math.Abs(b - c) * 0.5)) ) ||
( !mflag && (Math.Abs(s - b) >= (Math.Abs(c - d) * 0.5)) ) ||
( mflag && (Math.Abs(b - c) < tol) ) ||
( !mflag && (Math.Abs(c - d) < tol)) )
{
// bisection method
s = (a + b) * 0.5;
 
mflag = true;
}
else
{
mflag = false;
}
 
fs = f(s);// calculate fs
d = c; // first time d is being used (wasnt used on first iteration because mflag was set)
c = b; // set c equal to upper bound
fc = fb; // set f(c) = f(b)
 
if (fa * fs < 0) // fa and fs have opposite signs
{
b = s;
fb = fs; // set f(b) = f(s)
}
else
{
a = s;
fa = fs; // set f(a) = f(s)
}
 
if (Math.Abs(fa) < Math.Abs(fb)) // if magnitude of fa is less than magnitude of fb
{
Swap(ref a, ref b); // swap a and b
Swap(ref fa, ref fb); // make sure f(a) and f(b) are correct after swap
}
} // end for
 
throw new AggregateException("The solution does not converge or iterations are not sufficient");
}
// end brents_fun
}
</syntaxhighlight>
 
=={{header|C++}}==
<syntaxhighlight lang="cpp">#include <iostream>
 
double f(double x)
{
return (x*x*x - 3*x*x + 2*x);
}
 
int main()
{
double step = 0.001; // Smaller step values produce more accurate and precise results
std::pair<complex, complex> result = solve_quadratic_equation(1, -1e20, 1);
double start = -1;
std::cout << result.first << ", " << result.second << std::endl;
double stop = 3;
}</lang>
double value = f(start);
{{out}}
double sign = (value > 0);
(1e+20,0), (1e-20,0)
// Check for root at start
if ( 0 == value )
std::cout << "Root found at " << start << std::endl;
 
for( double x = start + step;
=={{header|Clojure}}==
x <= stop;
x += step )
{
value = f(x);
if ( ( value > 0 ) != sign )
// We passed a root
std::cout << "Root found near " << x << std::endl;
else if ( 0 == value )
// We hit a root
std::cout << "Root found at " << x << std::endl;
// Update our sign
sign = ( value > 0 );
}
}</syntaxhighlight>
 
===Brent's Method===
<lang clojure>(defn quadratic
Brent's Method uses a combination of the bisection method, inverse quadratic interpolation, and the secant method to find roots. It has a guaranteed run time equal to that of the bisection method (which always converges in a known number of steps (log2[(upper_bound-lower_bound)/tolerance] steps to be precise ) unlike the other methods), but the algorithm uses the much faster inverse quadratic interpolation and secant method whenever possible. The algorithm is robust and commonly used in libraries with a roots() function built in.
"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</lang>
 
The algorithm is coded as a function that returns a double value for the root. The function takes an input that requires the function being evaluated, the lower and upper bounds, the tolerance one is looking for before converging (i recommend 0.0001) and the maximum number of iterations before giving up on finding the root (the root will always be found if the root is bracketed and a sufficient number of iterations is allowed).
{{out}}
<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]
</lang>
 
The implementation is taken from the pseudo code on the wikipedia page for Brent's Method found here: https://en.wikipedia.org/wiki/Brent%27s_method.
=={{header|Common Lisp}}==
<syntaxhighlight lang="cpp">#include <iostream>
<lang lisp>(defun quadratic (a b c)
#include <cmath>
(list
#include <algorithm>
(/ (+ (- b) (sqrt (- (expt b 2) (* 4 a c)))) (* 2 a))
#include <functional>
(/ (- (- b) (sqrt (- (expt b 2) (* 4 a c)))) (* 2 a))))</lang>
 
double brents_fun(std::function<double (double)> f, double lower, double upper, double tol, unsigned int max_iter)
=={{header|D}}==
{
<lang d>import std.math, std.traits;
double a = lower;
double b = upper;
double fa = f(a); // calculated now to save function calls
double fb = f(b); // calculated now to save function calls
double fs = 0; // initialize
 
if (!(fa * fb < 0))
CommonType!(T1, T2, T3)[] naiveQR(T1, T2, T3)
{
(in T1 a, in T2 b, in T3 c)
std::cout << "Signs of f(lower_bound) and f(upper_bound) must be opposites" << std::endl; // throws exception if root isn't bracketed
pure nothrow if (isFloatingPoint!T1) {
return -11;
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];
}
 
if (std::abs(fa) < std::abs(b)) // if magnitude of f(lower_bound) is less than magnitude of f(upper_bound)
CommonType!(T1, T2, T3)[] cautiQR(T1, T2, T3)
{
(in T1 a, in T2 b, in T3 c)
std::swap(a,b);
pure nothrow if (isFloatingPoint!T1) {
std::swap(fa,fb);
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);
 
double c = a; // c now equals the largest magnitude of the lower and upper bounds
if (b * a < 0) {
double fc = fa; // precompute function evalutation for point c by assigning it the same value as fa
immutable x = (-b + SD) / 2 * a;
bool mflag = true; // boolean flag used to evaluate if statement later on
return [x, c / (a * x)];
double s = 0; // Our Root that will be returned
} else {
double d = 0; // Only used if mflag is unset (mflag == false)
immutable x = (-b - SD) / 2 * a;
return [c / (a * x), x];
}
}
 
for (unsigned int iter = 1; iter < max_iter; ++iter)
void main() {
{
import std.stdio;
// stop if converged on root or error is less than tolerance
writeln("With 32 bit float type:");
if (std::abs(b-a) < tol)
writefln(" Naive: [%(%g, %)]", naiveQR(1.0f, -10e5f, 1.0f));
{
writefln("Cautious: [%(%g, %)]", cautiQR(1.0f, -10e5f, 1.0f));
std::cout << "After " << iter << " iterations the root is: " << s << std::endl;
writeln("\nWith 64 bit double type:");
return s;
writefln(" Naive: [%(%g, %)]", naiveQR(1.0, -10e5, 1.0));
} // end if
writefln("Cautious: [%(%g, %)]", cautiQR(1.0, -10e5, 1.0));
writeln("\nWith real type:");
if (fa != fc && fb != fc)
writefln(" Naive: [%(%g, %)]", naiveQR(1.0L, -10e5L, 1.0L));
{
writefln("Cautious: [%(%g, %)]", cautiQR(1.0L, -10e5L, 1.0L));
// use inverse quadratic interopolation
}</lang>
s = ( a * fb * fc / ((fa - fb) * (fa - fc)) )
{{out}}
+ ( b * fa * fc / ((fb - fa) * (fb - fc)) )
<pre>With 32 bit float type:
+ ( c * fa * fb / ((fc - fa) * (fc - fb)) );
Naive: [1e+06, 0]
}
Cautious: [1e+06, 1e-06]
else
{
// secant method
s = b - fb * (b - a) / (fb - fa);
}
 
// checks to see whether we can use the faster converging quadratic && secant methods or if we need to use bisection
With 64 bit double type:
if ( ( (s < (3 * a + b) * 0.25) || (s > b) ) ||
Naive: [1e+06, 1.00001e-06]
( mflag && (std::abs(s-b) >= (std::abs(b-c) * 0.5)) ) ||
Cautious: [1e+06, 1e-06]
( !mflag && (std::abs(s-b) >= (std::abs(c-d) * 0.5)) ) ||
( mflag && (std::abs(b-c) < tol) ) ||
( !mflag && (std::abs(c-d) < tol)) )
{
// bisection method
s = (a+b)*0.5;
 
mflag = true;
With real type:
}
Naive: [1e+06, 1e-06]
else
Cautious: [1e+06, 1e-06]</pre>
{
mflag = false;
}
 
fs = f(s); // calculate fs
=={{header|Delphi}}==
d = c; // first time d is being used (wasnt used on first iteration because mflag was set)
See [https://rosettacode.org/wiki/Roots_of_a_quadratic_function#Pascal Pascal].
c = b; // set c equal to upper bound
fc = fb; // set f(c) = f(b)
 
if ( fa * fs < 0) // fa and fs have opposite signs
=={{header|Elixir}}==
{
<lang elixir>defmodule Quadratic do
b = s;
def roots(a, b, c) do
fb = fs; // set f(b) = f(s)
IO.puts "Roots of a quadratic function (#{a}, #{b}, #{c})"
}
d = b * b - 4 * a * c
else
a2 = a * 2
{
cond do
a = s;
d > 0 ->
fa = fs; // set sdf(a) = :math.sqrtf(ds)
}
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
 
if (std::abs(fa) < std::abs(fb)) // if magnitude of fa is less than magnitude of fb
Quadratic.roots(1, -2, 1)
{
Quadratic.roots(1, -3, 2)
std::swap(a,b); // swap a and b
Quadratic.roots(1, 0, 1)
std::swap(fa,fb); // make sure f(a) and f(b) are correct after swap
Quadratic.roots(1, -1.0e10, 1)
}
Quadratic.roots(1, 2, 3)
 
Quadratic.roots(2, -1, -6)</lang>
} // end for
 
std::cout<< "The solution does not converge or iterations are not sufficient" << std::endl;
 
} // end brents_fun
 
</syntaxhighlight>
 
=={{header|Clojure}}==
 
{{trans|Haskell}}
<syntaxhighlight lang="clojure">
 
(defn findRoots [f start stop step eps]
(filter #(-> (f %) Math/abs (< eps)) (range start stop step)))
</syntaxhighlight>
 
{{out}}
<pre>
> (findRoots #(+ (* % % %) (* -3 % %) (* 2 %)) -1.0 3.0 0.0001 0.00000001)
Roots of a quadratic function (1, -2, 1)
(-9.381755897326649E-14 0.9999999999998124 1.9999999999997022)
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|ERRECoffeeScript}}==
{{trans|Python}}
<lang>PROGRAM QUADRATIC
<syntaxhighlight lang="coffeescript">
print_roots = (f, begin, end, step) ->
# Print approximate roots of f between x=begin and x=end,
# using sign changes as an indicator that a root has been
# encountered.
x = begin
y = f(x)
last_y = y
cross_x_axis = ->
(last_y < 0 and y > 0) or (last_y > 0 and y < 0)
console.log '-----'
while x <= end
y = f(x)
if y == 0
console.log "Root found at", x
else if cross_x_axis()
console.log "Root found near", x
x += step
last_y = y
 
do ->
PROCEDURE SOLVE_QUADRATIC
# Smaller steps produce more accurate/precise results in general,
D=B*B-4*A*C
# but for many functions we'll never get exact roots, either due
IF ABS(D)<1D-6 THEN D=0 END IF
# to imperfect binary representation or irrational roots.
CASE SGN(D) OF
step = 0->1 / 256
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
 
f1 = (x) -> x*x*x - 3*x*x + 2*x
BEGIN
print_roots f1, -1, 5, step
PRINT(CHR$(12);) ! CLS
FORf2 TEST%=1 TO(x) 7-> DOx*x - 4*x + 3
print_roots f2, -1, 5, step
READ(A,B,C)
f3 = (x) -> x - 1.5
PRINT("For a=";A;",b=";B;",c=";C;TAB(32);)
print_roots f3, 0, 4, step
SOLVE_QUADRATIC
f4 = (x) -> x*x - 2
END FOR
print_roots f4, -2, 2, step
DATA(1,-1E9,1)
</syntaxhighlight>
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</lang>
{{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>
 
output
=={{header|Factor}}==
{{trans|Ada}}
<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>
 
<syntaxhighlight lang="text">
<lang factor>( scratchpad ) 1 -1.e20 1 quadratic-equation
> coffee roots.coffee
--- Data stack:
-----
1.0e+20
Root found at 0
9.999999999999999e-21</lang>
Root found at 1
Root found at 2
-----
Root found at 1
Root found at 3
-----
Root found at 1.5
-----
Root found near -1.4140625
Root found near 1.41796875
</syntaxhighlight>
 
=={{header|Common Lisp}}==
Middlebrook method
<lang factor>:: quadratic-equation2 ( a b c -- x1 x2 )
a c * sqrt b / :> q
1 4 q sq * - sqrt 0.5 * 0.5 + :> f
b neg a / f * c neg b / f / ;
</lang>
 
{{trans|Perl}}
 
<code>find-roots</code> prints roots (and values near roots) and returns a list of root designators, each of which is either a number <code><var>n</var></code>, in which case <code>(zerop (funcall function <var>n</var>))</code> is true, or a <code>cons</code> whose <code>car</code> and <code>cdr</code> are such that the sign of function at car and cdr changes.
<lang factor>( scratchpad ) 1 -1.e20 1 quadratic-equation
--- Data stack:
1.0e+20
1.0e-20</lang>
 
<syntaxhighlight lang="lisp">(defun find-roots (function start end &optional (step 0.0001))
=={{header|Forth}}==
(let* ((roots '())
Without locals:
(value (funcall function start))
<lang forth>: quadratic ( fa fb fc -- r1 r2 )
(plusp (plusp value)))
frot frot
( c a(when b(zerop value)
(format t "~&Root found at ~W." start))
fover 3 fpick f* -4e f* fover fdup f* f+
(do c((x a(+ bstart step) (+ detx step)))
((> x end) (nreverse roots))
fdup f0< if abort" imaginary roots" then
(setf value (funcall function x))
fsqrt
fover f0< if fnegate then(cond
((zerop value)
f+ fnegate
(format t "~&Root found at ~w." x)
( c a b-det )
2e f/ fover f/ (push x roots))
((not c(eql aplusp r1(plusp value)))
(format t "~&Root found near ~w." x)
frot frot f/ fover f/ ;</lang>
(push (cons (- x step) x) roots)))
With locals:
(setf plusp (plusp value)))))</syntaxhighlight>
<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/ ;
 
<pre>> (find-roots #'(lambda (x) (+ (* x x x) (* -3 x x) (* 2 x))) -1 3)
\ test
Root found near 5.3588345E-5.
1 set-precision
Root found near 1.0000072.
1e -1e6 1e quadratic fs. fs. \ 1e-6 1e6</lang>
Root found near 2.000073.
((-4.6411653E-5 . 5.3588345E-5)
(0.99990714 . 1.0000072)
(1.9999729 . 2.000073))</pre>
 
=={{header|FortranD}}==
<syntaxhighlight lang="d">import std.stdio, std.math, std.algorithm;
===Fortran 90===
{{works with|Fortran|90 and later}}
<lang fortran>PROGRAM QUADRATIC
 
bool nearZero(T)(in T a, in T b = T.epsilon * 4) pure nothrow {
IMPLICIT NONE
return abs(a) <= b;
INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(15)
}
REAL(dp) :: a, b, c, e, discriminant, rroot1, rroot2
COMPLEX(dp) :: croot1, croot2
 
T[] findRoot(T)(immutable T function(in T) pure nothrow fi,
WRITE(*,*) "Enter the coefficients of the equation ax^2 + bx + c"
in T start, in T end, in T step=T(0.001L),
WRITE(*, "(A)", ADVANCE="NO") "a = "
T tolerance = T(1e-4L)) {
READ *, a
if (step.nearZero)
WRITE(*,"(A)", ADVANCE="NO") "b = "
writefln("WARNING: step size may be too small.");
READ *, b
WRITE(*,"(A)", ADVANCE="NO") "c = "
READ *, c
WRITE(*,"(3(A,E23.15))") "Coefficients are: a = ", a, " b = ", b, " c = ", c
e = 1.0e-9_dp
discriminant = b*b - 4.0_dp*a*c
IF (ABS(discriminant) < e) THEN
rroot1 = -b / (2.0_dp * a)
WRITE(*,*) "The roots are real and equal:"
WRITE(*,"(A,E23.15)") "Root = ", rroot1
ELSE IF (discriminant > 0) THEN
rroot1 = -(b + SIGN(SQRT(discriminant), b)) / (2.0_dp * a)
rroot2 = c / (a * rroot1)
WRITE(*,*) "The roots are real:"
WRITE(*,"(2(A,E23.15))") "Root1 = ", rroot1, " Root2 = ", rroot2
ELSE
croot1 = (-b + SQRT(CMPLX(discriminant))) / (2.0_dp*a)
croot2 = CONJG(croot1)
WRITE(*,*) "The roots are complex:"
WRITE(*,"(2(A,2E23.15,A))") "Root1 = ", croot1, "j ", " Root2 = ", croot2, "j"
END IF</lang>
{{out}}
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
 
/// Search root by simple bisection.
===Fortran I===
T searchRoot(T a, T b) pure nothrow {
Source code written in FORTRAN I (october 1956) for the IBM 704.
T root;
<lang fortran>
int limit = 49;
COMPUTE ROOTS OF A QUADRATIC FUNCTION - 1956
READ 100,A,B,C T gap = b - a;
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
</lang>
 
while (!nearZero(gap) && limit--) {
=={{header|FreeBASIC}}==
if (fi(a).nearZero)
{{libheader|GMP}}
return a;
<lang freebasic>' version 20-12-2020
if (fi(b).nearZero)
' compile with: fbc -s console
return b;
root = (b + a) / 2.0L;
if (fi(root).nearZero)
return root;
((fi(a) * fi(root) < 0) ? b : a) = root;
gap = b - a;
}
 
return root;
#Include Once "gmp.bi"
}
 
immutable dir = T(end > start ? 1.0 : -1.0);
Sub solvequadratic_n(a As Double ,b As Double, c As Double)
immutable step2 = (end > start) ? abs(step) : -abs(step);
T[T] result;
for (T x = start; (x * dir) <= (end * dir); x += step2)
if (fi(x) * fi(x + step2) <= 0) {
immutable T r = searchRoot(x, x + step2);
result[r] = fi(r);
}
 
return result.keys.sort().release;
Dim As Double f, d = b ^ 2 - 4 * a * c
}
 
void report(T)(in T[] r, immutable T function(in T) pure f,
Select Case Sgn(d)
in T tolerance = T(1e-4L)) {
Case 0
if (r.length) {
Print "1: the single root is "; -b / 2 / a
writefln("Root found (tolerance = %1.4g):", tolerance);
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
 
foreach (const x; r) {
End Sub
immutable T y = f(x);
 
if (nearZero(y))
Sub solvequadratic_c(a As Double ,b As Double, c As Double)
writefln("... EXACTLY at %+1.20f, f(x) = %+1.4g",x,y);
else if (nearZero(y, tolerance))
writefln(".... MAY-BE at %+1.20f, f(x) = %+1.4g",x,y);
else
writefln("Verify needed, f(%1.4g) = " ~
"%1.4g > tolerance in magnitude", x, y);
}
} else
writefln("No root found.");
}
 
void main() {
Dim As Double f, d = b ^ 2 - 4 * a * c
static real f(in real x) pure nothrow {
Select Case Sgn(d)
return x ^^ 3 - (3 * x ^^ 2) + 2 * x;
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
 
findRoot(&f, -1.0L, 3.0L, 0.001L).report(&f);
Sub solvequadratic_gmp(a_ As Double ,b_ As Double, c_ As Double)
}</syntaxhighlight>
{{out}}
<pre>Root found (tolerance = 0.0001):
.... MAY-BE at -0.00000000000000000080, f(x) = -1.603e-18
... EXACTLY at +1.00000000000000000020, f(x) = -2.168e-19
.... MAY-BE at +1.99999999999999999950, f(x) = -8.674e-19</pre>
NB: smallest increment for real type in D is real.epsilon = 1.0842e-19.
 
=={{header|Dart}}==
#Define PRECISION 1024 ' about 300 digits
{{trans|Scala}}
#Define MAX 25
<syntaxhighlight lang="dart">double fn(double x) => x * x * x - 3 * x * x + 2 * x;
 
findRoots(Function(double) f, double start, double stop, double step, double epsilon) sync* {
Dim As ZString Ptr text
for (double x = start; x < stop; x = x + step) {
text = Callocate (1000)
if (fn(x).abs() < epsilon) yield x;
Mpf_set_default_prec(PRECISION)
}
}
 
main() {
Dim As Mpf_ptr a, b, c, d, t
// Vector(-9.381755897326649E-14, 0.9999999999998124, 1.9999999999997022)
a = Allocate(Len(__mpf_struct)) : Mpf_init_set_d(a, a_)
print(findRoots(fn, -1.0, 3.0, 0.0001, 0.000000001));
b = Allocate(Len(__mpf_struct)) : Mpf_init_set_d(b, b_)
}</syntaxhighlight>
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)
 
=={{header|Delphi}}==
mpf_mul(d, b, b)
See [https://rosettacode.org/wiki/Roots_of_a_function#Pascal Pascal].
mpf_set_ui(t, 4)
mpf_mul(t, t, a)
mpf_mul(t, t, c)
mpf_sub(d, d, t)
 
=={{header|DWScript}}==
Select Case mpf_sgn(d)
{{trans|C}}
Case 0
<syntaxhighlight lang="delphi">type TFunc = function (x : Float) : Float;
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
 
function f(x : Float) : Float;
End Sub
begin
Result := x*x*x-3.0*x*x +2.0*x;
end;
 
const e = 1.0e-12;
' ------=< MAIN >=------
 
function Secant(xA, xB : Float; f : TFunc) : Float;
Dim As Double a, b, c
const
Print "1: is the naieve way"
limit = 50;
Print "2: is the cautious way"
var
Print "3: is the naieve way with help of GMP"
fA, fB : Float;
Print
d : Float;
i : Integer;
begin
fA := f(xA);
for i := 0 to limit do begin
fB := f(xB);
d := (xB-xA)/(fB-fA)*fB;
if Abs(d) < e then
Exit(xB);
xA := xB;
fA := fB;
xB -= d;
end;
PrintLn(Format('Function is not converging near (%7.4f,%7.4f).', [xA, xB]));
Result := -99.0;
end;
 
const fstep = 1.0e-2;
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
 
var x := -1.032; // just so we use secant method
' empty keyboard buffer
var xx, value : Float;
While Inkey <> "" : Wend
var s := f(x)>0.0;
Print : Print "hit any key to end program"
 
Sleep
while (x < 3.0) do begin
End
value := f(x);
if Abs(value)<e then begin
PrintLn(Format("Root found at x= %12.9f", [x]));
s := (f(x+0.0001)>0.0);
end else if (value>0.0) <> s then begin
xx := Secant(x-fstep, x, f);
if xx <> -99.0 then // -99 meaning secand method failed
PrintLn(Format('Root found at x = %12.9f', [xx]))
else PrintLn(Format('Root found near x= %7.4f', [xx]));
s := (f(x+0.0001)>0.0);
end;
x += fstep;
end;</syntaxhighlight>
 
=={{header|EchoLisp}}==
We use the 'math' library, and define f(x) as the polynomial : x<sup>3</sup> -3x<sup>2</sup> +2x
 
<syntaxhighlight lang="lisp">
(lib 'math.lib)
Lib: math.lib loaded.
(define fp ' ( 0 2 -3 1))
(poly->string 'x fp) → x^3 -3x^2 +2x
(poly->html 'x fp) → x<sup>3</sup> -3x<sup>2</sup> +2x
(define (f x) (poly x fp))
(math-precision 1.e-6) → 0.000001
 
(root f -1000 1000) → 2.0000000133245677 ;; 2
(root f -1000 (- 2 epsilon)) → 1.385559938161431e-7 ;; 0
(root f epsilon (- 2 epsilon)) → 1.0000000002190812 ;; 1
</syntaxhighlight>
 
=={{header|Elixir}}==
{{trans|Ruby}}
<syntaxhighlight lang="elixir">defmodule RC do
def find_roots(f, range, step \\ 0.001) do
first .. last = range
max = last + step / 2
Stream.iterate(first, &(&1 + step))
|> Stream.take_while(&(&1 < max))
|> Enum.reduce(sign(first), fn x,sn ->
value = f.(x)
cond do
abs(value) < step / 100 ->
IO.puts "Root found at #{x}"
0
sign(value) == -sn ->
IO.puts "Root found between #{x-step} and #{x}"
-sn
true -> sign(value)
end
end)
end
defp sign(x) when x>0, do: 1
defp sign(x) when x<0, do: -1
defp sign(0) , do: 0
end
 
f = fn x -> x*x*x - 3*x*x + 2*x end
RC.find_roots(f, -1..3)</syntaxhighlight>
 
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</lang>
{{out}}
<pre>
<pre>1: is the naieve way
Root found at 8.81239525796218e-16
2: is the cautious way
Root found at 1.0000000000000016
3: is the naieve way with help of GMP
Root found at 1.9999999999998914
</pre>
 
=={{header|Erlang}}==
Find root(s) for 1X^2-1000000000X+1
<syntaxhighlight lang="erlang">% Implemented by Arjun Sunel
1: the real roots are 1000000000 and 0
-module(roots).
2: the real roots are 1000000000 and 1e-009
-export([main/0]).
3: the real roots are 9.9999999999999999900000000e+08 and 1.0000000000000000010000000e-09
main() ->
F = fun(X)->X*X*X - 3*X*X + 2*X end,
Step = 0.001, % Using smaller steps will provide more accurate results
Start = -1,
Stop = 3,
Sign = F(Start) > 0,
X = Start,
while(X, Step, Start, Stop, Sign,F).
 
while(X, Step, Start, Stop, Sign,F) ->
Find root(s) for 1X^2+0X+1
Value = F(X),
1: the complex roots are -0 +/- 1*i
if
2: the complex roots are -0 +/- 1*i
Value == 0 -> % We hit a root
3: the complex roots are 0.0000000000000000000000000e+00 +/- 1.0000000000000000000000000e+00*i
io:format("Root found at ~p~n",[X]),
while(X+Step, Step, Start, Stop, Value > 0,F);
 
(Value < 0) == Sign -> % We passed a root
Find root(s) for 2X^2-1X-6
io:format("Root found near ~p~n",[X]),
1: the real roots are 8 and -6
while(X+Step , Step, Start, Stop, Value > 0,F);
2: the real roots are 2 and -1.5
3: the real roots are 2.0000000000000000000000000e+00 and -1.5000000000000000000000000e+00
X > Stop ->
io:format("") ;
true ->
while(X+Step, Step, Start, Stop, Value > 0,F)
end.
</syntaxhighlight>
{{out}}
<pre>Root found near 8.81239525796218e-16
Root found near 1.0000000000000016
Root found near 2.0009999999998915
ok</pre>
 
=={{header|ERRE}}==
Find root(s) for 1X^2+2X-2
<syntaxhighlight lang="erre">
1: the real roots are 0.7320508075688772 and -2.732050807568877
PROGRAM ROOTS_FUNCTION
2: the real roots are -2.732050807568877 and 0.7320508075688773
3: the real roots are 7.3205080756887729352744634e-01 and -2.7320508075688772935274463e+00
 
!VAR E,X,STP,VALUE,S%,I%,LIMIT%,X1,X2,D
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
 
FUNCTION F(X)
Find root(s) for 1X^2+3X+2
F=X*X*X-3*X*X+2*X
1: the real roots are -1 and -2
END FUNCTION
2: the real roots are -2 and -0.9999999999999999
3: the real roots are -1.0000000000000000000000000e+00 and -2.0000000000000000000000000e+00
 
BEGIN
Find root(s) for 3X^2+4X+5
X=-1
1: the complex roots are -0.6666666666666666 +/- 1.105541596785133*i
STP=1.0E-6
2: the complex roots are -0.6666666666666666 +/- 1.105541596785133*i
E=1.0E-9
3: the complex roots are -6.6666666666666666666666667e-01 +/- 1.1055415967851332830383109e+00*i
S%=(F(X)>0)
 
PRINT("VERSION 1: SIMPLY STEPPING X")
Find root(s) for 1X^2-1e+100X+1
WHILE X<3.0 DO
1: the real roots are 1e+100 and 0
VALUE=F(X)
2: the real roots are 1e+100 and 1e-100
IF ABS(VALUE)<E THEN
3: the real roots are 1.0000000000000000159028911e+100 and 9.9999999999999998409710889e-101
PRINT("ROOT FOUND AT X =";X)
S%=NOT S%
ELSE
IF ((VALUE>0)<>S%) THEN
PRINT("ROOT FOUND AT X =";X)
S%=NOT S%
END IF
END IF
X=X+STP
END WHILE
 
PRINT
Find root(s) for 1X^2-1e+200X+1
PRINT("VERSION 2: SECANT METHOD")
1: the real roots are 1.#INF and -1.#INF
X1=-1.0
2: the real roots are 1e+200 and 1e-200
X2=3.0
3: the real roots are 9.9999999999999996973312221e+199 and 0.0000000000000000000000000e+00
E=1.0E-15
I%=1
LIMIT%=300
LOOP
IF I%>LIMIT% THEN
PRINT("ERROR: FUNCTION NOT CONVERGING")
EXIT
END IF
D=(X2-X1)/(F(X2)-F(X1))*F(X2)
IF ABS(D)<E THEN
IF D=0 THEN
PRINT("EXACT ";)
ELSE
PRINT("APPROXIMATE ";)
END IF
PRINT("ROOT FOUND AT X =";X2)
EXIT
END IF
X1=X2
X2=X2-D
I%=I%+1
END LOOP
END PROGRAM
</syntaxhighlight>
Note: Outputs are calculated in single precision.
{{out}}
<pre>
VERSION 1: SIMPLY STEPPING X
ROOT FOUND AT X = 8.866517E-07
ROOT FOUND AT X = 1.000001
ROOT FOUND AT X = 2
 
VERSION 2: SECANT METHOD
Find root(s) for 1X^2-1e+300X+1
EXACT ROOT FOUND AT X = 1
1: the real roots are 1.#INF and -1.#INF
</pre>
2: the real roots are 1e+300 and 1e-300
3: the real roots are 1.0000000000000000525047603e+300 and 0.0000000000000000000000000e+00</pre>
 
=={{header|GAPFortran}}==
{{works with|Fortran|90 and later}}
<lang gap>QuadraticRoots := function(a, b, c)
<syntaxhighlight lang="fortran">PROGRAM ROOTS_OF_A_FUNCTION
local d;
d := Sqrt(b*b - 4*a*c);
return [ (-b+d)/(2*a), (-b-d)/(2*a) ];
end;
 
IMPLICIT NONE
# 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 ]
 
INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(15)
# This works also with floating-point numbers
REAL(dp) :: f, e, x, step, value
QuadraticRoots(2.0, 2.0, -1.0);
LOGICAL :: s
# [ 0.366025, -1.36603 ]</lang>
f(x) = x*x*x - 3.0_dp*x*x + 2.0_dp*x
x = -1.0_dp ; step = 1.0e-6_dp ; e = 1.0e-9_dp
s = (f(x) > 0)
DO WHILE (x < 3.0)
value = f(x)
IF(ABS(value) < e) THEN
WRITE(*,"(A,F12.9)") "Root found at x =", x
s = .NOT. s
ELSE IF ((value > 0) .NEQV. s) THEN
WRITE(*,"(A,F12.9)") "Root found near x = ", x
s = .NOT. s
END IF
x = x + step
END DO
END PROGRAM ROOTS_OF_A_FUNCTION</syntaxhighlight>
The following approach uses the [[wp:Secant_method|Secant Method]] to numerically find one root. Which root is found will depend on the start values x1 and x2 and if these are far from a root this method may not converge.
<syntaxhighlight lang="fortran">INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(15)
INTEGER :: i=1, limit=100
REAL(dp) :: d, e, f, x, x1, x2
f(x) = x*x*x - 3.0_dp*x*x + 2.0_dp*x
x1 = -1.0_dp ; x2 = 3.0_dp ; e = 1.0e-15_dp
DO
IF (i > limit) THEN
WRITE(*,*) "Function not converging"
EXIT
END IF
d = (x2 - x1) / (f(x2) - f(x1)) * f(x2)
IF (ABS(d) < e) THEN
WRITE(*,"(A,F18.15)") "Root found at x = ", x2
EXIT
END IF
x1 = x2
x2 = x2 - d
i = i + 1
END DO</syntaxhighlight>
 
=={{header|FreeBASIC}}==
Simple bisection method.
<syntaxhighlight lang="freebasic">#Include "crt.bi"
const iterations=20000000
 
sub bisect( f1 as function(as double) as double,min as double,max as double,byref O as double,a() as double)
dim as double last,st=(max-min)/iterations,v
for n as double=min to max step st
v=f1(n)
if sgn(v)<>sgn(last) then
redim preserve a(1 to ubound(a)+1)
a(ubound(a))=n
O=n+st:exit sub
end if
last=v
next
end sub
 
function roots(f1 as function(as double) as double,min as double,max as double, a() as double) as long
redim a(0)
dim as double last,O,st=(max-min)/iterations,v
for n as double=min to max step st
v=f1(n)
if sgn(v)<>sgn(last) and n>min then bisect(f1,n-st,n,O,a()):n=O
last=v
next
return ubound(a)
end function
 
Function CRound(Byval x As Double,Byval precision As Integer=30) As String
If precision>30 Then precision=30
Dim As zstring * 40 z:Var s="%." &str(Abs(precision)) &"f"
sprintf(z,s,x)
If Val(z) Then Return Rtrim(Rtrim(z,"0"),".")Else Return "0"
End Function
 
function defn(x as double) as double
return x^3-3*x^2+2*x
end function
 
redim as double r()
 
print
if roots(@defn,-20,20,r()) then
print "in range -20 to 20"
print "All roots approximate"
print "number","root to 6 dec places","function value at root"
for n as long=1 to ubound(r)
print n,CRound(r(n),6),,defn(r(n))
next n
end if
sleep</syntaxhighlight>
{{out}}
<pre>in range -20 to 20
All roots approximate
number root to 6 dec places function value at root
1 0 -2.929925652002424e-009
2 1 1.477781779325033e-009
3 2 -2.897852187377925e-009</pre>
 
=={{header|Go}}==
Secant method. No error checking.
<lang go>package main
<syntaxhighlight lang="go">package main
 
import (
"fmt"
"math"
)
 
func main() {
func qr(a, b, c float64) ([]float64, []complex128) {
example := func(x float64) float64 { return x*x*x - 3*x*x + 2*x }
d := b*b-4*a*c
findroots(example, -.5, 2.6, 1)
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 testfindroots(af func(float64) float64, lower, bupper, cstep float64) {
for x0, x1 := lower, lower+step; x0 < upper; x0, x1 = x1, x1+step {
fmt.Print("coefficients: ", a, b, c, " -> ")
x1 = math.Min(x1, upper)
r, i := qr(a, b, c)
r, status := secant(f, x0, x1)
switch len(r) {
if status != "" && r >= x0 && r < x1 {
case 1:
fmt.PrintlnPrintf("one real root:%6.3f %s\n", r[0], status)
}
case 2:
}
fmt.Println("two real roots:", r[0], r[1])
default:
fmt.Println("two complex roots:", i[0], i[1])
}
}
 
func secant(f func(float64) float64, x0, x1 float64) (float64, string) {
func main() {
var f0 for _, c := range [][3]float64{
f1 := f(x0)
{1, -2, 1},
for i := 0; i < 100; i++ {1, 0, 1},
f0, f1 = f1, f(x1)
{1, -10, 1},
switch {
{1, -1000, 1},
case f1 == 0:
{1, -1e9, 1},
return x1, "exact"
} {
case math.Abs(x1-x0) < 1e-6:
test(c[0], c[1], c[2])
return x1, "approximate"
}
}
}</lang>
x0, x1 = x1, x1-f1*(x1-x0)/(f1-f0)
{{out}}
}
return 0, ""
}</syntaxhighlight>
Output:
<pre>
0.000 approximate
coefficients: 1 -2 1 -> one real root: 1
1.000 exact
coefficients: 1 0 1 -> two complex roots: (0+1i) (-0-1i)
2.000 approximate
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
</pre>
 
=={{header|Haskell}}==
<syntaxhighlight lang="haskell">f x = x^3-3*x^2+2*x
<lang haskell>import Data.Complex (Complex, realPart)
 
findRoots start stop step eps =
type CD = Complex Double
[x | x <- [start, start+step .. stop], abs (f x) < eps]</syntaxhighlight>
Executed in GHCi:
<syntaxhighlight lang="haskell">*Main> findRoots (-1.0) 3.0 0.0001 0.000000001
[-9.381755897326649e-14,0.9999999999998124,1.9999999999997022]</syntaxhighlight>
 
Or using package [http://hackage.haskell.org/package/hmatrix hmatrix] from HackageDB.
quadraticRoots :: (CD, CD, CD) -> (CD, CD)
<syntaxhighlight lang="haskell">import Numeric.GSL.Polynomials
quadraticRoots (a, b, c)
import Data.Complex
| 0 < realPart b =
( (2 * c) / (- b - d),
(- b - d) / (2 * a)
)
| otherwise =
( (- b + d) / (2 * a),
(2 * c) / (- b + d)
)
where
d = sqrt $ b ^ 2 - 4 * a * c
 
*Main> mapM_ print $ polySolve [0,2,-3,1]
main :: IO ()
(-5.421010862427522e-20) :+ 0.0
main =
2.000000000000001 :+ 0.0
mapM_
0.9999999999999996 :+ 0.0</syntaxhighlight>
(print . quadraticRoots)
No complex roots, so:
[ (3, 4, 4 / 3),
<syntaxhighlight lang="haskell">*Main> mapM_ (print.realPart) $ polySolve [0,2,-3,1]
(3, 2, -1),
-5.421010862427522e-20
(3, 2, 1),
2.000000000000001
(1, -10e5, 1),
0.9999999999999996</syntaxhighlight>
(1, -10e9, 1)
]</lang>
{{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>
 
It is possible to solve the problem directly and elegantly using robust bisection method and Alternative type class.
=={{header|Icon}} and {{header|Unicon}}==
<syntaxhighlight lang="haskell">import Control.Applicative
 
data Root a = Exact a | Approximate a deriving (Show, Eq)
{{trans|Ada}}
 
-- looks for roots on an interval
Works in both languages.
bisection :: (Alternative f, Floating a, Ord a) =>
<lang unicon>procedure main()
(a -> a) -> a -> a -> f (Root a)
solve(1.0, -10.0e5, 1.0)
bisection f a b | f a * f b > 0 = empty
end
| f a == 0 = pure (Exact a)
| f b == 0 = pure (Exact b)
| smallInterval = pure (Approximate c)
| otherwise = bisection f a c <|> bisection f c b
where c = (a + b) / 2
smallInterval = abs (a-b) < 1e-15 || abs ((a-b)/c) < 1e-15
 
-- looks for roots on a grid
procedure solve(a,b,c)
findRoots :: (Alternative f, Floating a, Ord a) =>
d := sqrt(b*b - 4.0*a*c)
roots := if b < 0 then [r1 := (a -b+d)/(2.0*> a), c/-> [a] -> а (Root a*r1)]
findRoots f [] = empty
else [r1 := (-b-d)/(2.0*a), c/(a*r1)]
findRoots f [x] = if f x == 0 then pure (Exact x) else empty
write(a,"*x^2 + ",b,"*x + ",c," has roots ",roots[1]," and ",roots[2])
findRoots f (a:b:xs) = bisection f a b <|> findRoots f (b:xs)</syntaxhighlight>
end</lang>
 
It is possible to use these functions with different Alternative functors: IO, Maybe or List:
{{out}}
<pre>λ> bisection (\x -> x*x-2) 1 2
<pre>
Approximate 1.414213562373094
->rqf 1.0 -0.000000001 1.0
λ> bisection (\x -> x-1) 1 2
1.0*x^2 + -1000000.0*x + 1.0 has roots 999999.999999 and 1.000000000001e-06
Exact 1.0
->
λ> bisection (\x -> x*x-2) 2 3 :: Maybe (Root Double)
</pre>
Nothing
λ> findRoots (\x -> x^3 - 3*x^2 + 2*x) [-3..3] :: Maybe (Root Double)
Just (Exact 0.0)
λ> findRoots (\x -> x^3 - 3*x^2 + 2*x) [-3..3] :: [Root Double]
[Exact 0.0,Exact 0.0,Exact 1.0,Exact 2.0]</pre>
 
To get rid of repeated roots use `Data.List.nub`
=={{header|IDL}}==
<pre>λ> Data.List.nub $ findRoots (\x -> x^3 - 3*x^2 + 2*x) [-3..3]
<lang idl>compile_OPT IDL2
[Exact 0.0,Exact 1.0,Exact 2.0]
λ> Data.List.nub $ findRoots (\x -> x^3 - 3*x^2 + x) [-3..3]
[Exact 0.0,Approximate 2.6180339887498967]</pre>
 
=={{header|HicEst}}==
print, "input a, press enter, input b, press enter, input c, press enter"
HicEst's [http://www.HicEst.com/SOLVE.htm SOLVE] function employs the Levenberg-Marquardt method:
read,a,b,c
<syntaxhighlight lang="hicest">OPEN(FIle='test.txt')
Promt='Enter values of a,b,c and hit enter'
 
1 DLG(NameEdit=x0, DNum=3)
a0=0.0
b0=0.0
c0=0.0 ;make them floating point variables
 
x = x0
x=-b+sqrt((b^2)-4*a*c)
chi2 = SOLVE(NUL=x^3 - 3*x^2 + 2*x, Unknown=x, I=iterations, NumDiff=1E-15)
y=-b-sqrt((b^2)-4*a*c)
EDIT(Text='approximate exact ', Word=(chi2 == 0), Parse=solution)
z=2*a
d= x/z
e= y/z
 
WRITE(FIle='test.txt', LENgth=6, Name) x0, x, solution, chi2, iterations
print, d,e</lang>
GOTO 1</syntaxhighlight>
<syntaxhighlight lang="hicest">x0=0.5; x=1; solution=exact; chi2=79E-32 iterations=65;
x0=0.4; x=2E-162 solution=exact; chi2=0; iterations=1E4;
x0=0.45; x=1; solution=exact; chi2=79E-32 iterations=67;
x0=0.42; x=2E-162 solution=exact; chi2=0; iterations=1E4;
x0=1.5; x=1.5; solution=approximate; chi2=0.1406; iterations=14:
x0=1.54; x=1; solution=exact; chi2=44E-32 iterations=63;
x0=1.55; x=2; solution=exact; chi2=79E-32 iterations=55;
x0=1E10; x=2; solution=exact; chi2=18E-31 iterations=511;
x0=-1E10; x=0; solution=exact; chi2=0; iterations=1E4;</syntaxhighlight>
 
=={{header|IS-BASICIcon}} and {{header|Unicon}}==
{{trans|Java}}
<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</lang>
 
Works in both languages:
=={{header|J}}==
<syntaxhighlight lang="unicon">procedure main()
'''Solution''' use J's built-in polynomial solver:
showRoots(f, -1.0, 4, 0.002)
p.
end
 
procedure f(x)
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).
return x^3 - 3*x^2 + 2*x
end
 
procedure showRoots(f, lb, ub, step)
'''Example''' using inputs from other solutions and the unstable example from the task description:
ox := x := lb
oy := f(x)
os := sign(oy)
while x <= ub do {
if (s := sign(y := f(x))) = 0 then write(x)
else if s ~= os then {
dx := x-ox
dy := y-oy
cx := x-dx*(y/dy)
write("~",cx)
}
(ox := x, oy := y, os := s)
x +:= step
}
end
 
procedure sign(x)
<lang j> coeff =. _3 |.\ 3 4 4r3 3 2 _1 3 2 1 1 _1e6 1 1 _1e9 1
return (x<0, -1) | (x>0, 1) | 0
> {:"1 p. coeff
end</syntaxhighlight>
_0.666667 _0.666667
_1 0.333333
_0.333333j0.471405 _0.333333j_0.471405
1e6 1e_6
1e9 1e_9</lang>
 
Output:
Of course <code>p.</code> generalizes to polynomials of arbitrary order (which isn't as great as that might sound, because of practical limitations). Given the coefficients <code>p.</code> returns the multiplier and roots of the polynomial. Given the multiplier and roots it returns the coefficients. For example using the cubic <math>0 + 16x - 12x^2 + 2x^3</math>:
<pre>
<lang j> p. 0 16 _12 2 NB. return multiplier ; roots
->roots
+-+-----+
~2.616794878713638e-18
|2|4 2 0|
~1.0
+-+-----+
~2.0
p. 2 ; 4 2 0 NB. return coefficients
->
0 16 _12 2</lang>
</pre>
 
=={{header|J}}==
Exploring the limits of precision:
 
J has builtin a root-finding operator, '''<tt>p.</tt>''', whose input is the coeffiecients of the polynomial (where the exponent of the indeterminate variable matches the index of the coefficient: 0 1 2 would be 0 + x + (2 times x squared)). Hence:
<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>
 
<syntaxhighlight lang="j"> 1{::p. 0 2 _3 1
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.
2 1 0</syntaxhighlight>
 
We can determine whether the roots are exact or approximate by evaluating the polynomial at the candidate roots, and testing for zero:
Middlebrook formula implemented in J
 
<syntaxhighlight lang="j"> (0=]p.1{::p.) 0 2 _3 1
<lang j>q_r=: verb define
1 1 1</syntaxhighlight>
'a b c' =. y
q=. b %~ %: a * c
f=. 0.5 + 0.5 * %:(1-4*q*q)
(-b*f%a),(-c%b*f)
)
 
As you can see, <tt>p.</tt> is also the operator which evaluates polynomials. This is not a coincidence.
q_r 1 _1e6 1
1e6 1e_6</lang>
 
That said, we could also implement the technique used by most others here. Specifically: we can implement the function as a black box and check every 1 millionth of a unit between minus one and three, and we can test that result for exactness.
=={{header|Java}}==
<lang java>public class QuadraticRoots {
private static class Complex {
double re, im;
 
<syntaxhighlight lang="j"> blackbox=: 0 2 _3 1&p.
public Complex(double re, double im) {
(#~ (=<./)@:|@blackbox) i.&.(1e6&*)&.(1&+) 3
this.re = re;
0 1 2
this.im = im;
0=blackbox 0 1 }2
1 1 1</syntaxhighlight>
 
Here, we see that each of the results (0, 1 and 2) are as accurate as we expect our computer arithmetic to be. (The = returns 1 where paired values are equal and 0 where they are not equal).
@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);
}
 
=={{header|Java}}==
@Override
<syntaxhighlight lang="java">public class Roots {
public String toString() {
public interface Function {
if (im == 0.0) {return String.format("%g", re);}
public double f(double x);
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[]int quadraticRootssign(double a, double b, double cx) {
return (x < 0.0) ? -1 : (x Complex[]> roots0.0) =? new1 Complex[2]: 0;
}
double d = b * b - 4.0 * a * c; // discriminant
double aa = a + a;
 
public static void printRoots(Function f, double lowerBound,
if (d < 0.0) {
double upperBound, double re = -b /step) aa;{
double x = lowerBound, ox = x;
double im = Math.sqrt(-d) / aa;
double y = f.f(x), oy = y;
roots[0] = new Complex(re, im);
int s = sign(y), os = s;
roots[1] = new Complex(re, -im);
 
} else if (b < 0.0) {
for (; x <= upperBound ; x += step) {
// Avoid calculating -b - Math.sqrt(d), to avoid any
s = sign(y = f.f(x));
// subtractive cancellation when it is near zero.
if (s == 0) {
double re = (-b + Math.sqrt(d)) / aa;
System.out.println(x);
roots[0] = new Complex(re, 0.0);
} else if (s != os) {
roots[1] = new Complex(c / (a * re), 0.0);
double dx = x - ox;
} else {
double dy = y - oy;
// Avoid calculating -b + Math.sqrt(d).
double recx = (-bx - Math.sqrtdx * (d))y / aady);
System.out.println("~" + cx);
roots[1] = new Complex(re, 0.0);
}
roots[0] = new Complex(c / (a * re), 0.0);
ox = x; oy }= y; os = s;
}
return roots;
}
 
public static void main(String[] args) {
Function poly = new Function () {
double[][] equations = {
public double f(double x) {
{1.0, 22.0, -1323.0}, // two distinct real roots
return x*x*x - 3*x*x + 2*x;
{6.0, -23.0, 20.0}, // with a != 1.0
}
{1.0, -1.0e9, 1.0}, // with one root near zero
};
{1.0, 2.0, 1.0}, // one real root (double root)
printRoots(poly, -1.0, 4, 0.002);
{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]);
}
}
}
}</langsyntaxhighlight>
Produces this output:
{{out}}
<pre>~2.616794878713638E-18
<pre>
~1.0000000000000002
a = 1.00000 b = 22.0000 c = -1323.00
~2.000000000000001</pre>
X1 = 27.0000
X2 = -49.0000
 
=={{header|JavaScript}}==
a = 6.00000 b = -23.0000 c = 20.0000
{{trans|Java}}
X1 = 2.50000
{{works with|SpiderMonkey|22}}
X2 = 1.33333
{{works with|Firefox|22}}
<syntaxhighlight lang="javascript">
// This function notation is sorta new, but useful here
// Part of the EcmaScript 6 Draft
// developer.mozilla.org/en-US/docs/Web/JavaScript/Reference/Functions_and_function_scope
var poly = (x => x*x*x - 3*x*x + 2*x);
 
function sign(x) {
a = 1.00000 b = -1.00000e+09 c = 1.00000
return (x < 0.0) ? -1 : (x > 0.0) ? 1 : 0;
X1 = 1.00000e+09
}
X2 = 1.00000e-09
 
function printRoots(f, lowerBound, upperBound, step) {
a = 1.00000 b = 2.00000 c = 1.00000
var x = lowerBound, ox = x,
X1,2 = -1.00000
y = f(x), oy = y,
s = sign(y), os = s;
 
for (; x <= upperBound ; x += step) {
a = 1.00000 b = 0.00000 c = 1.00000
s = sign(y = f(x));
X1 = 1.00000i
if (s == 0) {
X2 = -1.00000i
console.log(x);
}
else if (s != os) {
var dx = x - ox;
var dy = y - oy;
var cx = x - dx * (y / dy);
console.log("~" + cx);
}
ox = x; oy = y; os = s;
}
}
 
printRoots(poly, -1.0, 4, 0.002);
a = 1.00000 b = 1.00000 c = 1.00000
</syntaxhighlight>
X1 = -0.500000 + 0.866025i
X2 = -0.500000 - 0.866025i</pre>
 
=={{header|jq}}==
printRoots(f; lower; upper; step) finds approximations to the roots
{{ works with |jq|1.4}}
of an arbitrary continuous real-valued function, f, in the range
Currently jq does not include support for complex number operations, so
[lower, upper], assuming step is small enough.
a small library is included in the first section.
 
The algorithm is similar to that used for example in the Javascript section on this page, except that a bug has been removed at the point when the previous and current signs are compared.
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 function, f, may be an expression (as in the example below) or a defined filter.
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.
 
printRoots/3 emits an array of results, each of which is either a
'''Section 1''': Complex numbers (scrolling window)
number (representing an exact root within the limits of machine arithmetic) or a string consisting of "~" followed by an approximation to the root.
<div style="overflow:scroll; height:200px;">
<syntaxhighlight lang="jq">def sign:
<lang jq># Complex numbers as points [x,y] in the Cartesian plane
def real(z): if (z|type). ==< "number"0 then -1 elif . > 0 then z1 else z[0] end;
 
def printRoots(f; lowerBound; upperBound; step):
def imag(z): if (z|type) == "number" then 0 else z[1] end;
lowerBound as $x
| ($x|f) as $y
| ($y|sign) as $s
| reduce range($x; upperBound+step; step) as $x
# state: [ox, oy, os, roots]
( [$x, $y, $s, [] ];
.[0] as $ox | .[1] as $oy | .[2] as $os
| ($x|f) as $y
| ($y | sign) as $s
| if $s == 0 then [$x, $y, $s, (.[3] + [$x] )]
elif $s != $os and $os != 0 then
($x - $ox) as $dx
| ($y - $oy) as $dy
| ($x - ($dx * $y / $dy)) as $cx # by geometry
| [$x, $y, $s, (.[3] + [ "~\($cx)" ])] # an approximation
else [$x, $y, $s, .[3] ]
end )
| .[3] ;
</syntaxhighlight>
We present two examples, one where step is a power of 1/2, and one where it is not:
{{Out}}
<syntaxhighlight lang="jq">printRoots( .*.*. - 3*.*. + 2*.; -1.0; 4; 1/256)
 
[
def plus(x; y):
0,
if (x|type) == "number" then
1,
if (y|type) == "number" then [ x+y, 0 ]
2
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;
 
printRoots( .*.*. - 3*.*. + 2*.; -1.0; 4; .001)
def negate(x): multiply(-1; x);
[
"~1.320318770141425e-18",
"~1.0000000000000002",
"~1.9999999999999993"
]</syntaxhighlight>
 
=={{header|Julia}}==
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;
 
Assuming that one has the Roots package installed:
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;
 
<syntaxhighlight lang="julia">using Roots
def divide(x;y): multiply(x; invert(y));
 
println(find_zero(x -> x^3 - 3x^2 + 2x, (-100, 100)))</syntaxhighlight>
def magnitude(z):
real( multiply(z; conjugate(z))) | sqrt;
 
{{out}}
# 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 ;</lang></div>
'''Section 2:''' quadratic_roots(a;b;c)
<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
;</lang>
'''Section 3''':
Produce a table showing [i, error, solution] for solutions to x^2 - 10^i + 1 = 0
<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) * " " + .;
 
<pre>[0.0,1.0,2.0]</pre>
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</lang>
{{Out}} (scrolling window)
<div style="overflow:scroll; height:200px;">
<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]</lang></div>
 
Without the Roots package, Newton's method may be defined in this manner:
=={{header|Julia}}==
<syntaxhighlight lang="julia">function newton(f, fp, x::Float64,tol=1e-14::Float64,maxsteps=100::Int64)
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.
##f: the function of x
 
##fp: the derivative of f
Alternative solutions might make use of Julia's Polynomials or Roots packages.
 
local xnew, xold = x, Inf
<lang julia>using Printf
local fn, fo = f(xnew), Inf
 
local counter = 1
function quadroots(x::Real, y::Real, z::Real)
a, b, c = promote(float(x), y, z)
while (counter < maxsteps) && (abs(xnew - xold) > tol) && ( abs(fn - fo) > tol )
if a ≈ 0.0 return [-c / b] end
Δx = bxnew ^- 2 -f(xnew)/fp(xnew) 4a## *update cx
xnew, xold = x, xnew
if Δ ≈ 0.0 return [-sqrt(c / a)] end
if Δ < 0.0 Δ fn, fo = complexf(Δxnew), endfn
counter d += sqrt(Δ)1
end
if b < 0.0
d -= b
if counter >= maxsteps
return [d / 2a, 2c / d]
error("Did not converge in ", string(maxsteps), " steps")
else
d = -d - belse
xnew, counter
return [2c / d, d / 2a]
end
end
</syntaxhighlight>
 
Finding the roots of f(x) = x3 - 3x2 + 2x:
a = [1, 1, 1.0, 10]
b = [10, 2, -10.0 ^ 9, 1]
c = [1, 1, 1, 1]
 
<syntaxhighlight lang="julia">
for (x, y, z) in zip(a, b, c)
f(x) = x^3 - 3*x^2 + 2*x
@printf "The roots of %.2fx² + %.2fx + %.2f\n\tx₀ = (%s)\n" x y z join(round.(quadroots(x, y, z), 2), ", ")
fp(x) = 3*x^2-6*x+2
end</lang>
 
x_s, count = newton(f,fp,1.00)
</syntaxhighlight>
{{out}}
 
<pre>The roots of 1.00x² + 10.00x + 1.00
(1.0,2)
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|Kotlin}}==
{{trans|JavaC}}
<syntaxhighlight lang ="scala">import// java.langversion 1.Math1.*2
 
datatypealias classDoubleToDouble = 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"
}
}
 
fun f(x: Double) = x * x * x - 3.0 * x * x + 2.0 * x
data class Solution(val x1: Any, val x2: Any) {
override fun toString() = when(x1) {
x2 -> "X1,2 = $x1"
else -> "X1 = $x1, X2 = $x2"
}
}
 
fun secant(x1: Double, x2: Double, f: DoubleToDouble): Double {
val quadraticRoots by lazy {
val _2ae = a + a1.0e-12
val limit = 50
val d = b * b - 4.0 * a * c // discriminant
var xa = if (d < 0.0) {x1
var val rxb = -b / _2ax2
var val ifa = sqrtf(-dxa) / _2a
var i = 0
Solution(Complex(r, i), Complex(r, -i))
while (i++ < } elselimit) {
var fb = f(xb)
// avoid calculating -b +/- sqrt(d), to avoid any
val d = (xb //- subtractivexa) cancellation/ when(fb it- isfa) near* zero.fb
val r = if (b < 0Math.0) (-b + sqrtabs(d)) /< _2a else (-b - sqrt(de)) / _2abreak
xa = Solution(r, c / (a * r))xb
}fa = fb
xb -= d
}
if (i == limit) {
println("Function is not converging near (${"%7.4f".format(xa)}, ${"%7.4f".format(xb)}).")
return -99.0
}
return xb
}
 
fun main(args: Array<String>) {
val step = 1.0e-2
val equations = listOf(Equation(1.0, 22.0, -1323.0), // two distinct real roots
val e = 1.0e-12
Equation(6.0, -23.0, 20.0), // with a != 1.0
var x = -1.032
Equation(1.0, -1.0e9, 1.0), // with one root near zero
var s = f(x) > 0.0
Equation(1.0, 2.0, 1.0), // one real root (double root)
while (x < 3.0) {
Equation(1.0, 0.0, 1.0), // two imaginary roots
val value = f(x)
Equation(1.0, 1.0, 1.0)) // two complex roots
if (Math.abs(value) < e) {
println("Root found at x = ${"%12.9f".format(x)}")
s = f(x + 0.0001) > 0.0
}
else if ((value > 0.0) != s) {
val xx = secant(x - step, x, ::f)
if (xx != -99.0)
println("Root found at x = ${"%12.9f".format(xx)}")
else
println("Root found near x = ${"%7.4f".format(x)}")
s = f(x + 0.0001) > 0.0
}
x += step
}
}</syntaxhighlight>
 
equations.forEach { println("$it\n" + it.quadraticRoots) }
}</lang>
{{out}}
<pre>
<pre>Equation(a=1.0, b=22.0, c=-1323.0)
Root found at x = 0.000000000
X1 = -49.0, X2 = 27.0
Root found at x = 1.000000000
Equation(a=6.0, b=-23.0, c=20.0)
Root found at x = 2.000000000
X1 = 2.5, X2 = 1.3333333333333333
</pre>
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|lambdatalkLambdatalk}}==
<langsyntaxhighlight lang="scheme">
1) defining the function:
1) using lambdas:
{def func {lambda {:x} {+ {* 1 :x :x :x} {* -3 :x :x} {* 2 :x}}}}
-> func
 
2) printing roots:
{def equation
{S.map {lambda {:a :b :cx}
{if {< {abs {func :x}} 0.0001}
{b equation :a*x{sup 2}+:b*x+:c=0}
then {br}- a root found at :x else}}
{{lambda {:a' :b' :d}
{if {>S.serie -1 :d3 0.01}}
->
then {{lambda {:b' :d'}
- a root found at 7.528699885739343e-16
{equation.disp {+ :b' :d'} {- :b' :d'} 2 real roots}
- a root found at 1.0000000000000013
} :b' {/ {sqrt :d} :a'}}
- a root found at 2.000000000000002
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}} } }}
 
3) printing the roots of the "sin" function between -720° to +720°;
2) using let:
 
{S.map {lambda {:x}
{def equation
{if {< {abs {sin {* {/ {PI} 180} :x}}} 0.01}
{lambda {:a :b :c}
then {br}- a root found at :x° else}}
{b equation :a*x{sup 2}+:b*x+:c=0}
{let { {:a' {*S.serie -720 2+720 :a10}}
->
{:b' {/ {- :b} {* 2 :a}}}
- a root found at -720° {:d {- {* :b :b} {* 4 :a :c}}} }
- a root found at -540°
{if {> :d 0}
- a root found at -360°
then {let { {:b' :b'}
- a root found at -180° {:d' {/ {sqrt :d} :a'}} }
- a root found at 0° {equation.disp {+ :b' :d'} {- :b' :d'} 2 real roots} }
- a root found at 180°
else {if {< :d 0}
- a root found at 360°
then {let { {:b' :b'}
- a root found at 540° {:d' {/ {sqrt {- :d}} :a'}} }
- a root found at 720°
{equation.disp [:b',:d'] [:b',-:d'] 2 complex roots} }
</syntaxhighlight>
else {equation.disp :b' :b' one real double root} }} }}}
 
=={{header|Liberty BASIC}}==
3) a function to display results in an HTML table format
<syntaxhighlight lang="lb">' Finds and output the roots of a given function f(x),
' within a range of x values.
 
' [RC]Roots of an function
{def equation.disp
{lambda {:x1 :x2 :txt}
{table {@ style="background:#ffa"}
{tr {td :txt: }}
{tr {td x1 = :x1 }}
{tr {td x2 = :x2 }} } }}
 
mainwin 80 12
4) testing:
 
xMin =-1
equation 1*x2+1*x+-1=0
xMax = 3
2 real roots:
y =f( xMin) ' Since Liberty BASIC has an 'eval(' function the fn
x1 = 0.6180339887498949
' and limits would be better entered via 'input'.
x2 = -1.618033988749895
LastY =y
equation 1*x2+1*x+1=0
2 complex roots:
x1 = [-0.5,0.8660254037844386]
x2 = [-0.5,-0.8660254037844386]
 
eps =1E-12 ' closeness acceptable
equation 1*x2+-2*x+1=0
one real double root:
x1 = 1
x2 = 1
</lang>
 
bigH=0.01
=={{header|Liberty BASIC}}==
<lang lb>a=1:b=2:c=3
'assume a<>0
print quad$(a,b,c)
end
 
print
function quad$(a,b,c)
print " Checking for roots of x^3 -3 *x^2 +2 *x =0 over range -1 to +3"
D=b^2-4*a*c
x=-1*bprint
 
if D<0 then
x=xMin: dx = bigH
quad$=str$(x/(2*a));" +i";str$(sqr(abs(D))/(2*a));" , ";str$(x/(2*a));" -i";str$(sqr(abs(D))/abs(2*a))
elsedo
x=x+dx
quad$=str$(x/(2*a)+sqr(D)/(2*a));" , ";str$(x/(2*a)-sqr(D)/(2*a))
end if y = f(x)
'print x, dx, y
end function</lang>
if y*LastY <0 then 'there is a root, should drill deeper
if dx < eps then 'we are close enough
print " Just crossed axis, solution f( x) ="; y; " at x ="; using( "#.#####", x)
LastY = y
dx = bigH 'after closing on root, continue with big step
else
x=x-dx 'step back
dx = dx/10 'repeat with smaller step
end if
end if
loop while x<xMax
 
print
print " Finished checking in range specified."
 
end
 
function f( x)
=={{header|Logo}}==
f =x^3 -3 *x^2 +2 *x
<lang logo>to quadratic :a :b :c
end function</syntaxhighlight>
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>
 
=={{header|Lua}}==
<syntaxhighlight lang="lua">-- Function to have roots found
In order to correctly handle complex roots, qsolve must be given objects from a suitable complex number library,
function f (x) return x^3 - 3*x^2 + 2*x end
like that from the Complex Numbers article. However, this should be enough to demonstrate its accuracy:
 
-- Find roots of f within x=[start, stop] or approximations thereof
<lang lua>function qsolve(a, b, c)
function root (f, start, stop, step)
if b < 0 then return qsolve(-a, -b, -c) end
local roots, x, sign, foundExact, value = {}, start, f(start) > 0
val = b + (b^2 - 4*a*c)^(1/2) --this never exhibits instability if b > 0
while x <= stop do
return -val / (2 * a), -2 * c / val --2c / val is the same as the "unstable" second root
value = f(x)
if value == 0 then
table.insert(roots, {val = x, err = 0})
foundExact = true
end
if value > 0 ~= sign then
if foundExact then
foundExact = false
else
table.insert(roots, {val = x, err = step})
end
end
sign = value > 0
x = x + step
end
return roots
end
 
-- Main procedure
for i = 1, 12 do
print(qsolve"Root (1,to 0-10^i,12DP)\tMax. 1)Error\n")
for _, r in pairs(root(f, -1, 3, 10^-6)) do
end</lang>
print(string.format("%0.12f", r.val), r.err)
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.
end</syntaxhighlight>
{{out}}
<pre>Root (to 12DP) Max. Error
 
0.000000000008 1e-06
1.000000000016 1e-06
2.000000999934 1e-06</pre>
Note that the roots found are all near misses because fractional numbers that seem nice and 'round' in decimal (such as 10^-6) often have some rounding error when represented in binary. To increase the chances of finding exact integer roots, try using an integer start value with a step value that is a power of two.
<syntaxhighlight lang="lua">-- Main procedure
print("Root (to 12DP)\tMax. Error\n")
for _, r in pairs(root(f, -1, 3, 2^-10)) do
print(string.format("%0.12f", r.val), r.err)
end</syntaxhighlight>
{{out}}
<pre>Root (to 12DP) Max. Error
 
0.000000000000 0
1.000000000000 0
2.000000000000 0</pre>
 
=={{header|Maple}}==
<lang Maple>solve(a*x^2+b*x+c,x);
 
<syntaxhighlight lang="maple">f := x^3-3*x^2+2*x;
solve(1.0*x^2-10.0^9*x+1.0,x,explicit,allsolutions);
roots(f,x);</syntaxhighlight>
 
outputs:
fsolve(x^2-10^9*x+1,x,complex);</lang>
{{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
 
<syntaxhighlight lang="maple">[[0, 1], [1, 1], [2, 1]]</syntaxhighlight>
-9 9
1.000000000 10 , 1.000000000 10 </pre>
 
which means there are three roots. Each root is named as a pair where the first element is the value (0, 1, and 2), the second one the multiplicity (=1 for each means none of the three are degenerate).
=={{header|Mathematica}}/{{header|Wolfram Language}}==
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:
 
By itself (i.e. unless specifically asked to do so), Maple will only perform exact (symbolic) operations and not attempt to do any kind of numerical approximation.
<math>\left\{\left\{x\to \frac{-b-\sqrt{b^2-4 a c}}{2 a}\right\},\left\{x\to \frac{-b+\sqrt{b^2-4 a c}}{2 a}\right\}\right\}</math>
 
=={{header|Mathematica}}/{{header|Wolfram Language}}==
<math>\left\{\left\{x\to \frac{1}{50000+\sqrt{2499999999}}\right\},\left\{x\to 50000+\sqrt{2499999999}\right\}\right\}</math>
There are multiple obvious ways to do this in Mathematica.
===Solve===
This requires a full equation and will perform symbolic operations only:
<syntaxhighlight lang="mathematica">Solve[x^3-3*x^2+2*x==0,x]</syntaxhighlight>
Output
<pre> {{x->0},{x->1},{x->2}}</pre>
 
===NSolve===
<math>50000-\sqrt{2499999999}</math>
This requires merely the polynomial and will perform numerical operations if needed:
<syntaxhighlight lang="mathematica"> NSolve[x^3 - 3*x^2 + 2*x , x]</syntaxhighlight>
Output
<pre> {{x->0.},{x->1.},{x->2.}}</pre>
(note that the results here are floats)
 
===FindRoot===
<math>50000+\sqrt{2499999999}</math>
This will numerically try to find one(!) local root from a given starting point:
<syntaxhighlight lang="mathematica">FindRoot[x^3 - 3*x^2 + 2*x , {x, 1.5}]</syntaxhighlight>
Output
<pre> {x->0.}</pre>
From a different start point:
<syntaxhighlight lang="mathematica">FindRoot[x^3 - 3*x^2 + 2*x , {x, 1.1}]</syntaxhighlight>
Output
<pre>{x->1.}</pre>
(note that there is no guarantee which one is found).
 
===FindInstance===
<math>\begin{align}
This finds a value (optionally out of a given domain) for the given variable (or a set of values for a set of given variables) that satisfy a given equality or inequality:
\Biggl(
<syntaxhighlight lang="mathematica"> FindInstance[x^3 - 3*x^2 + 2*x == 0, x]</syntaxhighlight>
a & \neq 0 \And \And
Output
\left(
<pre>{{x->0}}</pre>
x==\frac{-b-\sqrt{b^2-4 a c}}{2 a}
\|
x==\frac{-b+\sqrt{b^2-4 a c}}{2 a}
\right)
\Biggr)\\
&\biggl\|
\left(
a==0 \And\And b\neq 0 \And\And x==-\frac{c}{b}
\right)\\
&\biggr\|
(c==0 \And \And b==0 \And \And a==0)
\end{align}
</math>
 
===Reduce===
<math>x==\frac{1}{50000+\sqrt{2499999999}}\|x==50000+\sqrt{2499999999}</math>
This will (symbolically) reduce a given expression to the simplest possible form, solving equations and performing substitutions in the process:
<syntaxhighlight lang="mathematica">Reduce[x^3 - 3*x^2 + 2*x == 0, x]</syntaxhighlight>
<pre> x==0||x==1||x==2</pre>
(note that this doesn't yield a "solution" but a different expression that expresses the same thing as the original)
 
=={{header|Maxima}}==
<math>\left\{\left\{x\to \frac{1}{50000+\sqrt{2499999999}}\right\},\left\{x\to 50000+\sqrt{2499999999}\right\}\right\}</math>
<syntaxhighlight lang="maxima">e: x^3 - 3*x^2 + 2*x$
 
/* Number of roots in a real interval, using Sturm sequences */
<math>\{x\to 0.00001\}</math>
nroots(e, -10, 10);
3
 
solve(e, x);
<math>\{x\to 100000.\}</math>
[x=1, x=2, x=0]
 
/* 'solve sets the system variable 'multiplicities */
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.
 
solve(x^4 - 2*x^3 + 2*x - 1, x);
=={{header|MATLAB}} / {{header|Octave}}==
[x=-1, x=1]
<lang Matlab>roots([1 -3 2]) % coefficients in decreasing order of power e.g. [x^n ... x^2 x^1 x^0]</lang>
 
multiplicities;
=={{header|Maxima}}==
[1, 3]
<lang maxima>solve(a*x^2 + b*x + c = 0, x);
 
/* Rational approximation of roots using Sturm sequences and bisection */
/* 2 2
sqrt(b - 4 a c) + b sqrt(b - 4 a c) - b
[x = - --------------------, x = --------------------]
2 a 2 a */
 
realroots(e);
fpprec: 40$
[x=1, x=2, x=0]
 
/* 'realroots also sets the system variable 'multiplicities */
solve(x^2 - 10^9*x + 1 = 0, x);
/* [x = 500000000 - sqrt(249999999999999999),
x = sqrt(249999999999999999) + 500000000] */
 
multiplicities;
bfloat(%);
[1, 1, 1]
/* [x = 1.0000000000000000009999920675269450501b-9,
x = 9.99999999999999998999999999999999999b8] */</lang>
 
/* Numerical root using Brent's method (here with another equation) */
=={{header|МК-61/52}}==
<lang>П2 С/П /-/ <-> / 2 / П3 x^2 С/П
ИП2 / - Вx <-> КвКор НОП x>=0 28 ИП3
x<0 24 <-> /-/ + / Вx С/П /-/ КвКор
ИП3 С/П</lang>
 
find_root(sin(t) - 1/2, t, 0, %pi/2);
''Input:'' a С/П b С/П c С/П
0.5235987755983
 
fpprec: 60$
{{out}} x<sub>1</sub> - РX; x<sub>2</sub> - РY (or error message, if D < 0).
 
bf_find_root(sin(t) - 1/2, t, 0, %pi/2);
=={{header|Modula-3}}==
5.23598775598298873077107230546583814032861566562517636829158b-1
{{trans|Ada}}
<lang modula3>MODULE Quad EXPORTS Main;
 
/* Numerical root using Newton's method */
IMPORT IO, Fmt, Math;
 
load(newton1)$
TYPE Roots = ARRAY [1..2] OF LONGREAL;
newton(e, x, 1.1, 1e-6);
1.000000017531147
 
/* For polynomials, Jenkins–Traub algorithm */
VAR r: Roots;
 
allroots(x^3 + x + 1);
PROCEDURE Solve(a, b, c: LONGREAL): Roots =
[x=1.161541399997252*%i+0.34116390191401,
VAR sd: LONGREAL := Math.sqrt(b * b - 4.0D0 * a * c);
x=0.34116390191401-1.161541399997252*%i,
x: LONGREAL;
x=-0.68232780382802]
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;
 
bfallroots(x^3 + x + 1);
BEGIN
[x=1.16154139999725193608791768724717407484314725802151429063617b0*%i + 3.41163901914009663684741869855524128445594290948999288901864b-1,
r := Solve(1.0D0, -10.0D5, 1.0D0);
x=3.41163901914009663684741869855524128445594290948999288901864b-1 - 1.16154139999725193608791768724717407484314725802151429063617b0*%i,
IO.Put("X1 = " & Fmt.LongReal(r[1]) & " X2 = " & Fmt.LongReal(r[2]) & "\n");
x=-6.82327803828019327369483739711048256891188581897998577803729b-1]</syntaxhighlight>
END Quad.</lang>
 
=={{header|Nim}}==
<syntaxhighlight lang Nim="nim">import math, complex, strformat
import strformat
 
func f(x: float): float = x ^ 3 - 3 * x ^ 2 + 2 * x
const Epsilon = 1e-15
 
var
type
step = 0.01
start = -1.0
stop = 3.0
sign = f(start) > 0
x = start
 
while x <= stop:
SolKind = enum solDouble, solFloat, solComplex
var value = f(x)
if value == 0:
echo fmt"Root found at {x:.5f}"
elif (value > 0) != sign:
echo fmt"Root found near {x:.5f}"
sign = value > 0
x += step</syntaxhighlight>
 
{{out}}
Roots = object
<pre>
case kind: SolKind
Root found near 0.00000
of solDouble:
Root found near 1.00000
fvalue: float
Root found near 2.00000
of solFloat:
</pre>
fvalues: (float, float)
of solComplex:
cvalues: (Complex64, Complex64)
 
=={{header|Objeck}}==
{{trans|C++}}
<syntaxhighlight lang="objeck">
bundle Default {
class Roots {
function : f(x : Float) ~ Float
{
return (x*x*x - 3.0*x*x + 2.0*x);
}
function : Main(args : String[]) ~ Nil
{
step := 0.001;
start := -1.0;
stop := 3.0;
value := f(start);
sign := (value > 0);
if(0.0 = value) {
start->PrintLine();
};
for(x := start + step; x <= stop; x += step;) {
value := f(x);
if((value > 0) <> sign) {
IO.Console->Instance()->Print("~")->PrintLine(x);
}
else if(0 = value) {
IO.Console->Instance()->Print("~")->PrintLine(x);
};
sign := (value > 0);
};
}
}
}
</syntaxhighlight>
 
=={{header|OCaml}}==
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)))
 
A general root finder using the False Position (Regula Falsi) method, which will find all simple roots given a small step size.
 
<syntaxhighlight lang="ocaml">let bracket u v =
func `$`(r: Roots): string =
((u > 0.0) && (v < 0.0)) || ((u < 0.0) && (v > 0.0));;
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"
 
let xtol a b = (a = b);; (* or use |a-b| < epsilon *)
 
let rec regula_falsi a b fa fb f =
when isMainModule:
if xtol a b then (a, fa) else
let c = (fb*.a -. fa*.b) /. (fb -. fa) in
let fc = f c in
if fc = 0.0 then (c, fc) else
if bracket fa fc then
regula_falsi a c fa fc f
else
regula_falsi c b fc fb f;;
 
let search lo hi step f =
const Equations = [(1.0, -2.0, 1.0),
let rec next x fx =
(10.0, 1.0, 1.0),
if x > hi then [] else
(1.0, -10.0, 1.0),
let y = x (1+.0, -1000.0,step 1.0),in
let fy = f y (1.0, -1e9, 1.0)]in
if fx = 0.0 then
(x,fx) :: next y fy
else if bracket fx fy then
(regula_falsi x y fx fy f) :: next y fy
else
next y fy in
next lo (f lo);;
 
let showroot (x,fx) =
for (a, b, c) in Equations:
Printf.printf "f(%.17f) = %.17f [%s]\n"
echo &"Equation: {a=}, {b=}, {c=}"
x fx (if fx = 0.0 then "exact" else "approx") in
let roots = quadRoots(a, b, c)
let pluralf x = if((x roots-.kind ==3.0)*.x solDouble:+. ""2.0)*.x else: "s"in
List.iter showroot (search (-5.0) 5.0 0.1 f);;</syntaxhighlight>
echo &" root{plural}: {roots}"</lang>
 
Output:
{{out}}
<pre>
<pre>Equation: a=1.0, b=-2.0, c=1.0
f(0.00000000000000000) = 0.00000000000000000 [exact]
root: 1.0
f(1.00000000000000022) = 0.00000000000000000 [exact]
Equation: a=10.0, b=1.0, c=1.0
f(1.99999999999999978) = 0.00000000000000000 [exact]
roots: -0.05 + 0.3122498999199199i, -0.05 + -0.3122498999199199i
</pre>
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>
 
Note these roots are exact solutions with floating-point calculation.
=={{header|OCaml}}==
 
=={{header|Octave}}==
<lang ocaml>type quadroots =
| RealRoots of float * float
| ComplexRoots of Complex.t * Complex.t ;;
 
If the equation is a polynomial, we can put the coefficients in a vector and use ''roots'':
let quadsolve a b c =
 
let d = (b *. b) -. (4.0 *. a *. c) in
<syntaxhighlight lang="octave">a = [ 1, -3, 2, 0 ];
if d < 0.0
r = roots(a);
then
% let's print it
let r = -. b /. (2.0 *. a)
for i = 1:3
and i = sqrt(-. d) /. (2.0 *. a) in
n = polyval(a, r(i));
ComplexRoots ({ Complex.re = r; Complex.im = i },
printf("x%d = %f (%f", i, r(i), n);
{ Complex.re = r; Complex.im = (-.i) })
if (n != 0.0)
else
letprintf(" r =not");
endif
if b < 0.0
printf(" exact)\n");
then ((sqrt d) -. b) /. (2.0 *. a)
endfor</syntaxhighlight>
else ((sqrt d) +. b) /. (-2.0 *. a)
 
in
Otherwise we can program our (simple) method:
RealRoots (r, c /. (r *. a))
 
;;</lang>
{{trans|Python}}
<syntaxhighlight lang="octave">function y = f(x)
y = x.^3 -3.*x.^2 + 2.*x;
endfunction
 
step = 0.001;
tol = 10 .* eps;
start = -1;
stop = 3;
se = sign(f(start));
 
x = start;
while (x <= stop)
v = f(x);
if ( (v < tol) && (v > -tol) )
printf("root at %f\n", x);
elseif ( sign(v) != se )
printf("root near %f\n", x);
endif
se = sign(v);
x = x + step;
endwhile</syntaxhighlight>
 
=={{header|Oforth}}==
 
<syntaxhighlight lang="oforth">: findRoots(f, a, b, st)
| x y lasty |
a f perform dup ->y ->lasty
 
a b st step: x [
x f perform -> y
y ==0 ifTrue: [ System.Out "Root found at " << x << cr ]
else: [ y lasty * sgn -1 == ifTrue: [ System.Out "Root near " << x << cr ] ]
y ->lasty
] ;
 
: f(x) x 3 pow x sq 3 * - x 2 * + ; </syntaxhighlight>
 
{{out}}
<pre>
<lang ocaml># quadsolve 1.0 0.0 (-2.0) ;;
findRoots(#f, -1, 3, 0.0001)
- : quadroots = RealRoots (-1.4142135623730951, 1.4142135623730949)
Root found at 0
Root found at 1
Root found at 2
 
findRoots(#f, -1.000001, 3, 0.0001)
# quadsolve 1.0 0.0 2.0 ;;
Root near 9.90000000000713e-005
- : quadroots =
Root near 1.000099
ComplexRoots ({Complex.re = 0.; Complex.im = 1.4142135623730951},
Root near 2.000099
{Complex.re = 0.; Complex.im = -1.4142135623730951})
</pre>
 
=={{header|ooRexx}}==
# quadsolve 1.0 (-1.0e5) 1.0 ;;
<syntaxhighlight lang="oorexx">/* REXX program to solve a cubic polynom equation
- : quadroots = RealRoots (99999.99999, 1.0000000001000001e-005)</lang>
a*x**3+b*x**2+c*x+d =(x-x1)*(x-x2)*(x-x3)
*/
Numeric Digits 16
pi3=Rxcalcpi()/3
Parse Value '1 -3 2 0' with a b c d
p=3*a*c-b**2
q=2*b**3-9*a*b*c+27*a**2*d
det=q**2+4*p**3
say 'p='p
say 'q='q
Say 'det='det
If det<0 Then Do
phi=Rxcalcarccos(-q/(2*rxCalcsqrt(-p**3)),16,'R')
Say 'phi='phi
phi3=phi/3
y1=rxCalcsqrt(-p)*2*Rxcalccos(phi3,16,'R')
y2=rxCalcsqrt(-p)*2*Rxcalccos(phi3+2*pi3,16,'R')
y3=rxCalcsqrt(-p)*2*Rxcalccos(phi3+4*pi3,16,'R')
End
Else Do
t=q**2+4*p**3
tu=-4*q+4*rxCalcsqrt(t)
tv=-4*q-4*rxCalcsqrt(t)
u=qroot(tu)/2
v=qroot(tv)/2
y1=u+v
y2=-(u+v)/2 (u+v)/2*rxCalcsqrt(3)
y3=-(u+v)/2 (-(u+v)/2*rxCalcsqrt(3))
End
say 'y1='y1
say 'y2='y2
say 'y3='y3
x1=y2x(y1)
x2=y2x(y2)
x3=y2x(y3)
Say 'x1='x1
Say 'x2='x2
Say 'x3='x3
Exit
 
qroot: Procedure
=={{header|Octave}}==
Parse Arg a
See [[Quadratic Equation#MATLAB|MATLAB]].
return sign(a)*rxcalcpower(abs(a),1/3,16)
 
y2x: Procedure Expose a b
Parse Arg real imag
xr=(real-b)/(3*a)
If imag<>'' Then Do
xi=(imag-b)/(3*a)
Return xr xi'i'
End
Else
Return xr
::requires 'rxmath' LIBRARY</syntaxhighlight>
{{out}}
<pre>p=-3
q=0
det=-108
phi=1.570796326794897
y1=2.999999999999999
y2=-3.000000000000000
y3=0.000000000000002440395154978758
x1=2
x2=0
x3=1.000000000000001</pre>
 
=={{header|PARI/GP}}==
===Gourdon–Schönhage algorithm===<!-- X. Gourdon, "Algorithmique du théorème fondamental de l'algèbre" (1993). -->
{{works with|PARI/GP|2.8.0+}}
<syntaxhighlight lang="parigp">polroots(x^3-3*x^2+2*x)</syntaxhighlight>
<lang parigp>roots(a,b,c)=polrootsreal(Pol([a,b,c]))</lang>
 
===Newton's method===
{{trans|C}}
This uses a modified version of the Newton–Raphson method.
Otherwise, coding directly:
<syntaxhighlight lang="parigp">polroots(x^3-3*x^2+2*x,1)</syntaxhighlight>
<lang parigp>roots(a,b,c)={
 
b /= a;
===Brent's method===
c /= a;
<syntaxhighlight lang="parigp">solve(x=-.5,.5,x^3-3*x^2+2*x)
my (delta = b^2 - 4*c, root=sqrt(delta));
solve(x=.5,1.5,x^3-3*x^2+2*x)
if (delta < 0,
solve(x=1.5,2.5,x^3-3*x^2+2*x)</syntaxhighlight>
[root-b,-root-b]/2
 
,
===Factorization to linear factors===
my(sol=if(b>0, -b - root,-b + root)/2);
<syntaxhighlight lang="parigp">findRoots(P)={
[sol,c/sol]
my(f=factor(P),t);
for(i=1,#f[,1],
if(poldegree(f[i,1]) == 1,
for(j=1,f[i,2],
print(-polcoeff(f[i,1], 0), " (exact)")
)
);
if(poldegree(f[i,1]) > 1,
t=polroots(f[i,1]);
for(j=1,#t,
for(k=1,f[i,2],
print(if(imag(t[j]) == 0.,real(t[j]),t[j]), " (approximate)")
)
)
)
)
};</lang>
findRoots(x^3-3*x^2+2*x)</syntaxhighlight>
 
===Factorization to quadratic factors===
Either way,
Of course this process could be continued to degrees 3 and 4 with sufficient additional work.
<lang parigp>roots(1,-1e9,1)</lang>
<syntaxhighlight lang="parigp">findRoots(P)={
gives one root around 0.000000001000000000000000001 and one root around 999999999.999999999.
my(f=factor(P),t);
for(i=1,#f[,1],
if(poldegree(f[i,1]) == 1,
for(j=1,f[i,2],
print(-polcoeff(f[i,1], 0), " (exact)")
)
);
if(poldegree(f[i,1]) == 2,
t=solveQuadratic(polcoeff(f[i,1],2),polcoeff(f[i,1],1),polcoeff(f[i,1],0));
for(j=1,f[i,2],
print(t[1]" (exact)\n"t[2]" (exact)")
)
);
if(poldegree(f[i,1]) > 2,
t=polroots(f[i,1]);
for(j=1,#t,
for(k=1,f[i,2],
print(if(imag(t[j]) == 0.,real(t[j]),t[j]), " (approximate)")
)
)
)
)
};
solveQuadratic(a,b,c)={
my(t=-b/2/a,s=b^2/4/a^2-c/a,inner=core(numerator(s))/core(denominator(s)),outer=sqrtint(s/inner));
if(inner < 0,
outer *= I;
inner *= -1
);
s=if(inner == 1,
outer
,
if(outer == 1,
Str("sqrt(", inner, ")")
,
Str(outer, " * sqrt(", inner, ")")
)
);
if (t,
[Str(t, " + ", s), Str(t, " - ", s)]
,
[s, Str("-", s)]
)
};
findRoots(x^3-3*x^2+2*x)</syntaxhighlight>
 
=={{header|Pascal}}==
{{trans|Fortran}}
some parts translated from Modula2
<langsyntaxhighlight lang="pascal">Program QuadraticRootsRootsFunction;
 
var
ae, bx, cstep, q, fvalue: double;
s: boolean;
i, limit: integer;
begin
ax1, x2, d:= 1double;
b := -10e9;
c := 1;
q := sqrt(a * c) / b;
f := (1 + sqrt(1 - 4 * q * q)) / 2;
 
function f(const x: double): double;
writeln ('Version 1:');
begin
writeln ('x1: ', (-b * f / a):16, ', x2: ', (-c / (b * f)):16);
f := x*x*x - 3*x*x + 2*x;
end;
 
begin
writeln ('Version 2:');
qx := sqrt(b * b:= - 4 * a * c)1;
ifstep b:= < 0 then1.0e-6;
e := 1.0e-9;
s := (f(x) > 0);
 
writeln('Version 1: simply stepping x:');
while x < 3.0 do
begin
fvalue := f(-b + qx) / 2 * a;
if abs(value) < e then
writeln ('x1: ', f:16, ', x2: ', (c / (a * f)):16);
end begin
writeln ('root found at x = ', x);
else
s := not s;
end
else if ((value > 0) <> s) then
begin
writeln ('root found at x = ', x);
s := not s;
end;
x := x + step;
end;
writeln('Version 2: secant method:');
x1 := -1.0;
x2 := 3.0;
e := 1.0e-15;
i := 1;
limit := 300;
while true do
begin
fif :=i (-b> -limit q) / 2 * a;then
begin
writeln ('x1: ', (c / (a * f)):16, ', x2: ', f:16);
writeln('Error: function not converging');
exit;
end;
d := (x2 - x1) / (f(x2) - f(x1)) * f(x2);
if abs(d) < e then
begin
if d = 0 then
write('Exact ')
else
write('Approximate ');
writeln('root found at x = ', x2);
exit;
end;
x1 := x2;
x2 := x2 - d;
i := i + 1;
end;
end.
</syntaxhighlight>
</lang>
Output:
{{out}}
<pre>
Version 1: simply stepping x:
root found at x = 7.91830063542152E-012
x1: 1.00000000E+010, x2: 1.00000000E-010
root found at x = 1.00000000001584E+000
Version 2:
root found at x = 1.99999999993357E+000
x1: 1.00000000E+010, x2: 1.00000000E-010
Version 2: secant method:
Exact root found at x = 1.00000000000000E+000
</pre>
 
=={{header|Perl}}==
<syntaxhighlight lang="perl">sub f
When using [http://perldoc.perl.org/Math/Complex.html Math::Complex] perl automatically convert numbers when necessary.
{
<lang perl>use Math::Complex;
my $x = shift;
 
return ($x * $x * $x - 3*$x*$x + 2*$x);
}
 
my $step = 0.001; # Smaller step values produce more accurate and precise results
my $start = -1;
my $stop = 3;
my $value = &f($start);
my $sign = $value > 0;
 
# Check for root at start
($x1,$x2) = solveQuad(1,2,3);
 
print "x1Root =found at $x1,start\n" x2if ( 0 == $x2\n"value );
 
for( my $x = $start + $step;
sub solveQuad
$x <= $stop;
$x += $step )
{
$value = &f($x);
my ($a,$b,$c) = @_;
 
my $root = sqrt($b**2 - 4*$a*$c);
return ( -$b + $root )/(2*$a), if ( -$b0 -== $rootvalue )/(2*$a);
{
}</lang>
# We hit a root
print "Root found at $x\n";
}
elsif ( ( $value > 0 ) != $sign )
{
# We passed a root
print "Root found near $x\n";
}
 
# Update our sign
$sign = ( $value > 0 );
}</syntaxhighlight>
 
=={{header|Phix}}==
{{trans|ERRECoffeeScript}}
<!--<langsyntaxhighlight Phixlang="phix">(phixonline)-->
<span style="color: #008080;">procedure</span> <span style="color: #000000;">solve_quadraticprint_roots</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequenceinteger</span> <span style="color: #000000;">t3f</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">atom</span> <span style="color: #000000;">start</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">stop</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">step</span><span style="color: #0000FF;">)</span>
<span style="color: #000080;font-style:italic;">--
<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>
-- Print approximate roots of f between x=start and x=stop, using
<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>
-- sign changes as an indicator that a root has been encountered.
<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>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">u</span>
<span style="color: #008080004080;">ifatom</span> <span style="color: #7060A8000000;">absx</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">d</span><span style="color: #0000FF;">)<</span><span style="color: #000000;">1e-6start</span> <span style="color: #0080800000FF;">then,</span> <span style="color: #000000;">dy</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: #0080807060A8;">switchputs</span> <span style="color: #7060A80000FF;">sign(</span><span style="color: #0000FF000000;">(1</span><span style="color: #0000000000FF;">d,</span><span style="color: #0000FF008000;">)"-----\n"</span> <span style="color: #0080800000FF;">do)</span>
<span style="color: #008080;">casewhile</span> <span style="color: #000000;">0x</span><span style="color: #0000FF;">:</span> <span style="color: #000000;">t</span> <span style="color: #0000FF000000;">=stop</span> <span style="color: #008000008080;">"single root is %g"do</span>
<span style="color: #000000004080;">uatom</span> <span style="color: #0000FF000000;">=last_y</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;">}y</span>
<span style="color: #008080000000;">casey</span> <span style="color: #0000000000FF;">1=</span> <span style="color: #0000FF000000;">:f</span> <span style="color: #0000000000FF;">t(</span> <span style="color: #0000FF000000;">=x</span> <span style="color: #0080000000FF;">"real roots are %g and %g")</span>
<span style="color: #008080;">if</span> <span style="color: #000000;">y</span><span style="color: #0000FF;">=</span><span style="color: #000000;">0</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: #000000008080;">uor</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{-(</span><span style="color: #000000;">flast_y</span><span style="color: #0000FF;">*<</span><span style="color: #000000;">b0</span><span style="color: #0000FF;">/</span><span style="color: #000000008080;">aand</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">cy</span><span style="color: #0000FF;">/</span><span style="color: #000000;">b</span><span style="color: #0000FF;">/</span><span style="color: #000000;">f0</span><span style="color: #0000FF;">})</span>
<span style="color: #008080;">caseor</span> <span style="color: #0000FF;">-(</span><span style="color: #000000;">1last_y</span><span style="color: #0000FF;">:></span> <span style="color: #000000;">t0</span> <span style="color: #0000FF008080;">=and</span> <span style="color: #008000000000;">y</span><span style="complexcolor: roots#0000FF;"><</span><span arestyle="color: %g#000000;">0</span><span +style="color: #0000FF;">)</-span> %g*i<span style="color: #008080;">then</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">u1</span> <span style="color: #0000FF;">=,</span> <span style="color: #0000FF008000;">{-"Root found %s %.10g\n"</span><span style="color: #0000000000FF;">b,</span> <span style="color: #0000FF;">/{</span><span style="color: #000000008080;">2iff</span><span style="color: #0000FF;">/(</span><span style="color: #000000;">ay</span><span style="color: #0000FF;">,=</span><span style="color: #7060A8000000;">sqrt0</span><span style="color: #0000FF;">(-?</span><span style="color: #000000008000;">d"at"</span><span style="color: #0000FF;">)/:</span><span style="color: #000000008000;">2"near"</span><span style="color: #0000FF;">/),</span><span style="color: #000000;">ax</span><span style="color: #0000FF;">})</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">switchif</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1x</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;">)})step</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">while</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span>
<span style="color: #000080;font-style:italic;">-- Smaller steps produce more accurate/precise results in general,
<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>
-- but for many functions we'll never get exact roots, either due
<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>
-- to imperfect binary representation or irrational roots.</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: #0000FF008080;">{constant</span> <span style="color: #000000;">1step</span> <span style="color: #0000FF;">,=</span> <span style="color: #000000;">21</span><span style="color: #0000FF;">,-/</span><span style="color: #000000;">2</span><span style="color: #0000FF;">},256</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: #7060A8008080;">papplyfunction</span> <span style="color: #000000;">f1</span><span style="color: #0000FF;">(</span><span style="color: #004080;">atom</span> <span style="color: #000000;">testsx</span><span style="color: #0000FF;">,)</span> <span style="color: #008080;">return</span> <span style="color: #000000;">solve_quadraticx</span><span style="color: #0000FF;">)*</span><span style="color: #000000;">x</span><span style="color: #0000FF;">*</span><span style="color: #000000;">x</span><span style="color: #0000FF;">-</span><span style="color: #000000;">3</span><span style="color: #0000FF;">*</span><span style="color: #000000;">x</span><span style="color: #0000FF;">*</span><span style="color: #000000;">x</span><span style="color: #0000FF;">+</span><span style="color: #000000;">2</span><span style="color: #0000FF;">*</span><span style="color: #000000;">x</span> <span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">f2</span><span style="color: #0000FF;">(</span><span style="color: #004080;">atom</span> <span style="color: #000000;">x</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">return</span> <span style="color: #000000;">x</span><span style="color: #0000FF;">*</span><span style="color: #000000;">x</span><span style="color: #0000FF;">-</span><span style="color: #000000;">4</span><span style="color: #0000FF;">*</span><span style="color: #000000;">x</span><span style="color: #0000FF;">+</span><span style="color: #000000;">3</span> <span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<!--</lang>-->
<span style="color: #008080;">function</span> <span style="color: #000000;">f3</span><span style="color: #0000FF;">(</span><span style="color: #004080;">atom</span> <span style="color: #000000;">x</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">return</span> <span style="color: #000000;">x</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1.5</span> <span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">f4</span><span style="color: #0000FF;">(</span><span style="color: #004080;">atom</span> <span style="color: #000000;">x</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">return</span> <span style="color: #000000;">x</span><span style="color: #0000FF;">*</span><span style="color: #000000;">x</span><span style="color: #0000FF;">-</span><span style="color: #000000;">2</span> <span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #000000;">print_roots</span><span style="color: #0000FF;">(</span><span style="color: #000000;">f1</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;">5</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">step</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">print_roots</span><span style="color: #0000FF;">(</span><span style="color: #000000;">f2</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;">5</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">step</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">print_roots</span><span style="color: #0000FF;">(</span><span style="color: #000000;">f3</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">4</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">step</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">print_roots</span><span style="color: #0000FF;">(</span><span style="color: #000000;">f4</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;">2</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">step</span><span style="color: #0000FF;">)</span>
<!--</syntaxhighlight>-->
{{out}}
<pre>
-----
for a=1,b=-1e+9,c=1 the real roots are 1e+9 and 1e-9
Root found at 0
for a=1,b=0,c=1 the complex roots are 0 +/- 1*i
Root found at 1
for a=2,b=-1,c=-6 the real roots are 2 and -1.5
Root found at 2
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
Root found at 1
for a=1,b=3,c=2 the real roots are -2 and -1
Root found at 3
for a=3,b=4,c=5 the complex roots are -0.666667 +/- 1.10554*i
-----
Root found at 1.5
-----
Root found near -1.4140625
Root found near 1.41796875
</pre>
 
=={{header|PicoLisp}}==
{{trans|Clojure}}
<lang PicoLisp>(scl 40)
<syntaxhighlight lang="picolisp">(de findRoots (F Start Stop Step Eps)
(filter
'((N) (> Eps (abs (F N))))
(range Start Stop Step) ) )
 
(scl 12)
(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
(findRoots
(solveQuad 1.0 -1000000.0 1.0)
'((X) (+ (*/ X X X `(* 1.0 1.0)) (*/ -3 X X 1.0) (* 2 X)))
(6 .) )</lang>
-1.0 3.0 0.0001 0.00000001 ) )</syntaxhighlight>
{{out}}
Output:
<pre>-> ("999,999.999999" "0.000001")</pre>
<pre>-> ("0.000" "1.000" "2.000")</pre>
 
=={{header|PL/I}}==
<syntaxhighlight lang="pl/i">
<lang PL/I>
f: declareprocedure (c1, c2x) returns (float complex,(18));
declare x float (a, b, c, x1, x218) float;
return (x**3 - 3*x**2 + 2*x );
end f;
 
declare eps get listfloat, (ax, b,y) cfloat (18);
declare dx fixed decimal (15,13);
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>
 
eps = 1e-12;
=={{header|Python}}==
{{libheader|NumPy}}
This solution compares the naïve method with three "better" methods.
<lang python>#!/usr/bin/env python3
 
do dx = -5.03 to 5 by 0.1;
import math
x = dx;
import cmath
if sign(f(x)) ^= sign(f(dx+0.1)) then
import numpy
call locate_root;
end;
 
locate_root: procedure;
def quad_discriminating_roots(a,b,c, entier = 1e-5):
declare (left, mid, right) float (18);
"""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
 
put skip list ('Looking for root in [' || x, x+0.1 || ']' );
def middlebrook(a, b, c):
left = x; right = dx+0.1;
try:
PUT SKIP LIST (F(LEFT), F(RIGHT) );
q = math.sqrt(a*c)/b
if abs(f(left) ) < eps then
f = .5+ math.sqrt(1-4*q*q)/2
do; put skip list ('Found a root at x=', left); return; end;
except ValueError:
else if abs(f(right) ) < eps then
q = cmath.sqrt(a*c)/b
do; put skip list ('Found a root at x=', right); return; end;
f = .5+ cmath.sqrt(1-4*q*q)/2
do forever;
return (-b/a)*f, -c/(b*f)
mid = (left+right)/2;
if sign(f(mid)) = 0 then
do; put skip list ('Root found at x=', mid); return; end;
else if sign(f(left)) ^= sign(f(mid)) then
right = mid;
else
left = mid;
/* put skip list (left || right); */
if abs(right-left) < eps then
do; put skip list ('There is a root near ' ||
(left+right)/2); return;
end;
end;
end locate_root;
</syntaxhighlight>
 
=={{header|PureBasic}}==
def whatevery(a, b, c):
{{trans|C++}}
try:
<syntaxhighlight lang="purebasic">Procedure.d f(x.d)
d = math.sqrt(b*b-4*a*c)
ProcedureReturn x*x*x-3*x*x+2*x
except ValueError:
EndProcedure
d = cmath.sqrt(b*b-4*a*c)
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))
 
Procedure main()
def div(n, d):
OpenConsole()
"""Divide, with a useful interpretation of division by zero."""
Define.d StepSize= 0.001
try:
Define.d Start=-1, stop=3
return n/d
Define.d value=f(start), x=start
except ZeroDivisionError:
Define.i oldsign=Sign(value)
if n:
return n*float('inf')
If value=0
return float('nan')
PrintN("Root found at "+StrF(start))
EndIf
While x<=stop
value=f(x)
If Sign(value) <> oldsign
PrintN("Root found near "+StrF(x))
ElseIf value = 0
PrintN("Root found at "+StrF(x))
EndIf
oldsign=Sign(value)
x+StepSize
Wend
EndProcedure
 
main()</syntaxhighlight>
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),
]
 
=={{header|Python}}==
print('Naive:')
{{trans|Perl}}
for c in testcases:
<syntaxhighlight lang="python">f = lambda x: x * x * x - 3 * x * x + 2 * x
print("{} {:.5} {:.5}".format(*quad_discriminating_roots(*c)))
 
step = 0.001 # Smaller step values produce more accurate and precise results
print('\nMiddlebrook:')
start = -1
for c in testcases:
stop = 3
print(("{:.5} "*2).format(*middlebrook(*c)))
 
sign = f(start) > 0
print('\nWhat Every...')
for c in testcases:
print(("{:.5} "*2).format(*whatevery(*c)))
 
x = start
print('\nNumpy:')
while x <= stop:
for c in testcases:
value = f(x)
print(("{:.5} "*2).format(*numpy.roots(c)))</lang>
{{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
 
if value == 0:
Middlebrook:
# We hit a root
-0.66667 -0.66667
print "Root found at", x
(-1+0j) (0.33333+0j)
elif (value > 0) != sign:
(-0.33333-0.4714j) (-0.33333+0.4714j)
# We passed a root
1e+09 1e-09
print "Root found near", x
1e+100 1e-100
1e+200 1e-200
1e+300 1e-300
 
# Update our sign
What Every...
sign = value > 0
-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
 
x += step</syntaxhighlight>
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}}==
{{trans|Octave}}
<lang R>qroots <- function(a, b, c) {
<syntaxhighlight lang="r">f <- sqrtfunction(bx) * bx^3 - 4 3*x^2 a+ 2* c + 0i)x
 
if (abs(b - r) > abs(b + r)) {
findroots <- function(f, begin, end, tol = 1e-20, step = 0.001) {
z <- (-b + r) / (2 * a)
se <- ifelse(sign(f(begin))==0, 1, sign(f(begin)))
} else {
zx <- (-b - r) / (2 * a)begin
while ( x <= end ) {
v <- f(x)
if ( abs(v) < tol ) {
print(sprintf("root at %f", x))
} else if ( ifelse(sign(v)==0, 1, sign(v)) != se ) {
print(sprintf("root near %f", x))
}
se <- ifelse( sign(v) == 0 , 1, sign(v))
x <- x + step
}
c(z, c / (z * a))
}
 
findroots(f, -1, 3)</syntaxhighlight>
qroots(1, 0, 2i)
[1] -1+1i 1-1i
 
=={{header|Racket}}==
qroots(1, -1e9, 1)
[1] 1e+09+0i 1e-09+0i</lang>
 
<syntaxhighlight lang="racket">
Using the builtin '''polyroot''' function (note the order of coefficients is reversed):
#lang racket
 
;; Attempts to find all roots of a real-valued function f
<lang R>polyroot(c(2i, 0, 1))
;; in a given interval [a b] by dividing the interval into N parts
[1] -1+1i 1-1i
;; and using the root-finding method on each subinterval
;; which proves to contain a root.
(define (find-roots f a b
#:divisions [N 10]
#:method [method secant])
(define h (/ (- b a) N))
(for*/list ([x1 (in-range a b h)]
[x2 (in-value (+ x1 h))]
#:when (or (root? f x1)
(includes-root? f x1 x2)))
(find-root f x1 x2 #:method method)))
 
;; Finds a root of a real-valued function f
polyroot(c(1, -1e9, 1))
;; in a given interval [a b].
[1] 1e-09+0i 1e+09+0i</lang>
(define (find-root f a b #:method [method secant])
(cond
[(root? f a) a]
[(root? f b) b]
[else (and (includes-root? f a b) (method f a b))]))
 
;; Returns #t if x is a root of a real-valued function f
=={{header|Racket}}==
;; with absolute accuracy (tolerance).
<lang Racket>#lang racket
(define (quadraticroot? af bx) c(almost-equal? 0 (f x)))
(let* ((-b (- b))
(delta (- (expt b 2) (* 4 a c)))
(denominator (* 2 a)))
(list
(/ (+ -b (sqrt delta)) denominator)
(/ (- -b (sqrt delta)) denominator))))
 
;; Returns #t if interval (a b) contains a root
;(quadratic 1 0.0000000000001 -1)
;; (or the odd number of roots) of a real-valued function f.
;'(0.99999999999995 -1.00000000000005)
(define (includes-root? f a b) (< (* (f a) (f b)) 0))
;(quadratic 1 0.0000000000001 1)
;'(-5e-014+1.0i -5e-014-1.0i)</lang>
 
;; Returns #t if a and b are equal with respect to
=={{header|Raku}}==
;; the relative accuracy (tolerance).
(formerly Perl 6)
(define (almost-equal? a b)
(or (< (abs (+ b a)) (tolerance))
(< (abs (/ (- b a) (+ b a))) (tolerance))))
 
(define tolerance (make-parameter 5e-16))
Raku has complex number handling built in.
</syntaxhighlight>
 
Different root-finding methods
<lang perl6>for
[1, 2, 1],
[1, 2, 3],
[1, -2, 1],
[1, 0, -4],
[1, -10**6, 1]
-> @coefficients {
printf "Roots for %d, %d, %d\t=> (%s, %s)\n",
|@coefficients, |quadroots(@coefficients);
}
 
<syntaxhighlight lang="racket">
sub quadroots (*[$a, $b, $c]) {
(define (secant f a b)
( -$b + $_ ) / (2 * $a),
(let next ([x1 -$ba] -[y1 $_(f a)] /[x2 (2b] *[y2 $a(f b)] [n 50])
(define x3 (/ (- (* x1 y2) (* x2 y1)) (- y2 y1)))
given
(cond
($b ** 2 - 4 * $a * $c ).Complex.sqrt.narrow
; if the method din't converge within given interval
}</lang>
; switch to more robust bisection method
{{out}}
[(or (not (< a x3 b)) (zero? n)) (bisection f a b)]
<pre>Roots for 1, 2, 1 => (-1, -1)
[(almost-equal? x3 x2) x3]
Roots for 1, 2, 3 => (-1+1.4142135623731i, -1-1.4142135623731i)
Roots for 1, -2, 1 [else (next x2 y2 x3 (f =>x3) (1,sub1 1n))])))
Roots for 1, 0, -4 => (2, -2)
Roots for 1, -1000000, 1 => (999999.999999, 1.00000761449337e-06)</pre>
 
(define (bisection f x1 x2)
=={{header|REXX}}==
(let divide ([a x1] [b x2])
===version 1===
(and (<= (* (f a) (f b)) 0)
The REXX language doesn't have a &nbsp; '''sqrt''' &nbsp; function, &nbsp; nor does it support complex numbers natively.
(let ([c (* 0.5 (+ a b))])
(if (almost-equal? a b)
c
(or (divide a c) (divide c b)))))))
</syntaxhighlight>
 
Examples:
Since "unlimited" decimal precision is part of the REXX language, the &nbsp; '''numeric digits''' &nbsp; was increased
<syntaxhighlight lang="racket">
<br>(from a default of &nbsp; '''9''') &nbsp; to &nbsp; '''200''' &nbsp; to accommodate when a root is closer to zero than the other root.
-> (find-root (λ (x) (- 2. (* x x))) 1 2)
1.414213562373095
-> (sqrt 2)
1.4142135623730951
 
-> (define (f x) (+ (* x x x) (* -3.0 x x) (* 2.0 x)))
Note that only nine decimal digits (precision) are shown in the &nbsp; ''displaying'' &nbsp; of the output.
-> (find-roots f -3 4 #:divisions 50)
'(2.4932181969624796e-33 1.0 2.0)
</syntaxhighlight>
 
In order to provide a comprehensive code the given solution does not optimize the number of function calls.
This REXX version supports &nbsp; ''complex numbers'' &nbsp; for the result.
The functional nature of Racket allows to perform the optimization without changing the main code using memoization.
<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. */</lang>
{{out|output|text=&nbsp; when using the input of: &nbsp; &nbsp; <tt> 1 &nbsp; -10e5 &nbsp; 1 </tt>}}
<pre>
a = 1
b = -1000000
c = 1
 
Simple memoization operator
root1 = 1000000
<syntaxhighlight lang="racket">
root2 = 0.000001
(define (memoized f)
</pre>
(define tbl (make-hash))
The following output is when Regina 3.9.3 REXX is used.
(λ x
(cond [(hash-ref tbl x #f) => values]
[else (define res (apply f x))
(hash-set! tbl x res)
res])))
</syntaxhighlight>
 
To use memoization just call
{{out|output|text=&nbsp; when using the input of: &nbsp; &nbsp; <tt> 1 &nbsp; -10e9 &nbsp; 1 </tt>}}
<syntaxhighlight lang="racket">
<pre>
-> (find-roots (memoized f) -3 4 #:divisions 50)
a = 1
'(2.4932181969624796e-33 1.0 2.0)
b = -1.0E+10
</syntaxhighlight>
c = 1
 
The profiling shows that memoization reduces the number of function calls
root1 = 1.000000000E+10
in this example from 184 to 67 (50 calls for primary interval division and about 6 calls for each point refinement).
root2 = 1E-10
</pre>
The following output is when R4 REXX is used.
 
=={{header|Raku}}==
{{out|output|text=&nbsp; when using the input of: &nbsp; &nbsp; <tt> 1 &nbsp; -10e9 &nbsp; 1 </tt>}}
(formerly Perl 6)
<pre>
Uses exact arithmetic.
a = 1
<syntaxhighlight lang="raku" line>sub f(\x) { x³ - 3*x² + 2*x }
b = -1E+10
c = 1
 
my $start = -1;
root1 = 1E+10
my $stop = 3;
root2 = 0.0000000001
my $step = 0.001;
</pre>
{{out|output|text=&nbsp; when using the input of: &nbsp; &nbsp; <tt> 3 &nbsp; 2 &nbsp; 1 </tt>}}
<pre>
a = 3
b = 2
c = 1
 
for $start, * + $step ... $stop -> $x {
root1 = -0.333333333+0.471404521i
state $sign = 0;
root2 = -0.333333333-0.471404521i
given f($x) {
</pre>
my $next = .sign;
{{out|output|text=&nbsp; when using the input of: &nbsp; &nbsp; <tt> 1 &nbsp; 0 &nbsp; 1 </tt>
when 0.0 {
<pre>
say "Root found at $x";
a = 1
b = 0 }
when $sign and $next != $sign {
c = 1
say "Root found near $x";
}
NEXT $sign = $next;
}
}</syntaxhighlight>
{{out}}
<pre>Root found at 0
Root found at 1
Root found at 2</pre>
 
=={{header|REXX}}==
root1 = 0+1i
Both of these REXX versions use the &nbsp; '''bisection method'''.
root2 = 0-1i
===function coded as a REXX function===
<syntaxhighlight lang="rexx">/*REXX program finds the roots of a specific function: x^3 - 3*x^2 + 2*x via bisection*/
parse arg bot top inc . /*obtain optional arguments from the CL*/
if bot=='' | bot=="," then bot= -5 /*Not specified? Then use the default.*/
if top=='' | top=="," then top= +5 /* " " " " " " */
if inc=='' | inc=="," then inc= .0001 /* " " " " " " */
z= f(bot - inc) /*compute 1st value to start compares. */
!= sign(z) /*obtain the sign of the initial value.*/
do j=bot to top by inc /*traipse through the specified range. */
z= f(j); $= sign(z) /*compute new value; obtain the sign. */
if z=0 then say 'found an exact root at' j/1
else if !\==$ then if !\==0 then say 'passed a root at' j/1
!= $ /*use the new sign for the next compare*/
end /*j*/ /*dividing by unity normalizes J [↑] */
exit /*stick a fork in it, we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
f: parse arg x; return x * (x * (x-3) +2) /*formula used ──► x^3 - 3x^2 + 2x */
/*with factoring ──► x{ x^2 -3x + 2 } */
/*more " ──► x{ x( x-3 ) + 2 } */</syntaxhighlight>
{{out|output|text=&nbsp; when using the defaults for input:}}
<pre>
found an exact root at 0
found an exact root at 1
found an exact root at 2
</pre>
 
===function Version 2coded in-line===
This version is about &nbsp; '''40%''' &nbsp; faster than the 1<sup>st</sup> REXX version.
<lang rexx>/* REXX ***************************************************************
<syntaxhighlight lang="rexx">/*REXX program finds the roots of a specific function: x^3 - 3*x^2 + 2*x via bisection*/
* 26.07.2913 Walter Pachl
parse arg bot top inc . /*obtain optional arguments from the CL*/
**********************************************************************/
if bot=='' | bot=="," then bot= -5 /*Not specified? Then use the default.*/
Numeric Digits 30
if top=='' | top=="," then top= +5 /* " " " " " " */
Parse Arg a b c 1 alist
if inc=='' | inc=="," then inc= .0001 /* " " " " " " */
Select
x= bot - inc /*compute 1st value to start compares. */
When a='' | a='?' Then
z= x * (x * (x-3) + 2) /*formula used ──► x^3 - 3x^2 + 2x */
Call exit 'rexx qgl a b c solves a*x**2+b*x+c'
!= sign(z) /*obtain the sign of the initial value.*/
When words(alist)<>3 Then
do x=bot to top by inc /*traipse through the specified range. */
Call exit 'three numbers are required'
z= x * (x * (x-3) + 2); $= sign(z) /*compute new value; obtain the sign. */
Otherwise
if z=0 then say 'found an exact root at' x/1
Nop
else if !\==$ then if !\==0 then say 'passed a root at' x/1
End
!= $ /*use the new sign for the next compare*/
gl=a'*x**2'
end /*x*/ /*dividing by unity normalizes X [↑] */</syntaxhighlight>
Select
{{out|output|text=&nbsp; is the same as the 1<sup>st</sup> REXX version.}} <br><br>
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'
 
=={{header|Ring}}==
d=b**2-4*a*c
<syntaxhighlight lang="ring">
If d<0 Then Do
load "stdlib.ring"
dd=sqrt(-d)
function = "return pow(x,3)-3*pow(x,2)+2*x"
r=-b/(2*a)
rangemin = -1
i=dd/(2*a)
rangemax = 3
x1=r'+'i'i'
stepsize = 0.001
x2=r'-'i'i'
accuracy = 0.1
End
roots(function, rangemin, rangemax, stepsize, accuracy)
Else Do
dd=sqrt(d)
func roots funct, min, max, inc, eps
x1=(-b+dd)/(2*a)
oldsign = 0
x2=(-b-dd)/(2*a)
for x = min to max step inc
End
num = sign(eval(funct))
Say 'x1='||x1
if num = 0
Say 'x2='||x2
see "root found at x = " + x + nl
Exit
num = -oldsign
sqrt:
else if num != oldsign and oldsign != 0
/* REXX ***************************************************************
if inc < eps
* EXEC to calculate the square root of x with high precision
see "root found near x = " + x + nl
**********************************************************************/
else roots(funct, x-inc, x+inc/8, inc/8, eps) ok ok ok
Parse Arg x
oldsign = num
prec=digits()
next
prec1=2*prec
</syntaxhighlight>
eps=10**(-prec1)
Output:
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</lang>
{{out}}
<pre>
root found near x = 0.00
Version 1:
root found near ax = 1.00
root found near bx = -12.00
</pre>
c = 0
 
=={{header|RLaB}}==
root1 = 1
RLaB implements a number of solvers from the GSL and the netlib that find the roots of a real or vector function of a real or vector variable.
root2 = 0
The solvers are grouped with respect whether the variable is a scalar, ''findroot'', or a vector, ''findroots''. Furthermore, for each group there are two types of solvers, one that does not require the derivative of the objective function (which root(s) are being sought), and one that does.
 
The script that finds a root of a scalar function <math>f(x) = x^3-3\,x^2 + 2\,x</math> of a scalar variable ''x''
Version 2:
using the bisection method on the interval -5 to 5 is,
1*x**2-1.0000000001*x+1.e-9 = 0
<syntaxhighlight lang="rlab">
x1=0.9999999991000000000025
f = function(x)
x2=0.0000000009999999999975
{
</pre>
rval = x .^ 3 - 3 * x .^ 2 + 2 * x;
return rval;
};
 
>> findroot(f, , [-5,5])
=={{header|Ring}}==
0
<lang>
</syntaxhighlight>
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
 
For a detailed description of the solver and its parameters interested reader is directed to the ''rlabplus'' manual.
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]
</lang>
 
=={{header|Ruby}}==
{{trans|Python}}
{{works with|Ruby|1.9.3+}}
The CMath#sqrt method will return a Complex instance if necessary.
<lang ruby>require 'cmath'
 
<syntaxhighlight lang="ruby">def sign(x)
def quadratic(a, b, c)
x <=> 0
sqrt_discriminant = CMath.sqrt(b**2 - 4*a*c)
[(-b + sqrt_discriminant) / (2.0*a), (-b - sqrt_discriminant) / (2.0*a)]
end
 
def find_roots(f, range, step=0.001)
p quadratic(3, 4, 4/3.0) # [-2/3]
sign = sign(f[range.begin])
p quadratic(3, 2, -1) # [1/3, -1]
range.step(step) do |x|
p quadratic(3, 2, 1) # [(-1/3 + sqrt(2/9)i), (-1/3 - sqrt(2/9)i)]
value = f[x]
p quadratic(1, 0, 1) # [(0+i), (0-i)]
if value == 0
p quadratic(1, -1e6, 1) # [1e6, 1e-6]
puts "Root found at #{x}"
p quadratic(-2, 7, 15) # [-3/2, 5]
elsif sign(value) == -sign
p quadratic(1, -2, 1) # [1]
puts "Root found between #{x-step} and #{x}"
p quadratic(1, 3, 3) # [(-3 + sqrt(3)i)/2), (-3 - sqrt(3)i)/2)]</lang>
end
sign = sign(value)
end
end
 
f = lambda { |x| x**3 - 3*x**2 + 2*x }
find_roots(f, -1..3)</syntaxhighlight>
 
{{out}}
<pre>
Root found at 0.0
[-0.6666666666666666, -0.6666666666666666]
Root found at 1.0
[0.3333333333333333, -1.0]
Root found at 2.0
[(-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]
[(-1.5+0.8660254037844386i), (-1.5-0.8660254037844386i)]
</pre>
 
Or we could use Enumerable#inject, monkey patching and block:
=={{header|Run BASIC}}==
<lang runbasic>print "FOR 1,2,3 => ";quad$(1,2,3)
print "FOR 4,5,6 => ";quad$(4,5,6)
 
<syntaxhighlight lang="ruby">class Numeric
FUNCTION quad$(a,b,c)
def sign
d = b^2-4 * a*c
x self <=> -1*b0
end
if d<0 then
end
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><pre>FOR 1,2,3 => -1 +i1.41421356 , -1 -i1.41421356
FOR 4,5,6 => -0.625 +i1.05326872 , -0.625 -i1.05326872</pre>
 
def find_roots(range, step = 1e-3)
=={{header|Scala}}==
range.step( step ).inject( yield(range.begin).sign ) do |sign, x|
Using [[Arithmetic/Complex#Scala|Complex]] class from task Arithmetic/Complex.
value = yield(x)
<lang scala>import ArithmeticComplex._
if value == 0
object QuadraticRoots {
puts "Root found at #{x}"
def solve(a:Double, b:Double, c:Double)={
valelsif dvalue.sign == b*b-4.0*a*csign
puts "Root found between #{x-step} and #{x}"
val aa = a+a
end
value.sign
if (d < 0.0) { // complex roots
end
val re= -b/aa;
end
val im = math.sqrt(-d)/aa;
(Complex(re, im), Complex(re, -im))
}
else { // real roots
val re=if (b < 0.0) (-b+math.sqrt(d))/aa else (-b -math.sqrt(d))/aa
(re, (c/(a*re)))
}
}
}</lang>
Usage:
<lang scala>val equations=Array(
(1.0, 22.0, -1323.0), // two distinct real roots
(6.0, -23.0, 20.0), // with a != 1.0
(1.0, -1.0e9, 1.0), // with one root near zero
(1.0, 2.0, 1.0), // one real root (double root)
(1.0, 0.0, 1.0), // two imaginary roots
(1.0, 1.0, 1.0) // two complex roots
);
equations.foreach{v =>
val (a,b,c)=v
println("a=%g b=%g c=%g".format(a,b,c))
val roots=solve(a, b, c)
println("x1="+roots._1)
if(roots._1 != roots._2) println("x2="+roots._2)
println
}</lang>
{{out}}
<pre>a=1.00000 b=22.0000 c=-1323.00
x1=-49.0
x2=27.0
 
find_roots(-1..3) { |x| x**3 - 3*x**2 + 2*x }</syntaxhighlight>
a=6.00000 b=-23.0000 c=20.0000
x1=2.5
x2=1.3333333333333333
 
=={{header|Rust}}==
a=1.00000 b=-1.00000e+09 c=1.00000
<syntaxhighlight lang="rust">// 202100315 Rust programming solution
x1=1.0E9
x2=1.0E-9
 
use roots::find_roots_cubic;
a=1.00000 b=2.00000 c=1.00000
x1=-1.0
 
fn main() {
a=1.00000 b=0.00000 c=1.00000
x1=-0.0 + 1.0i
x2=-0.0 + -1.0i
 
let roots = find_roots_cubic(1f32, -3f32, 2f32, 0f32);
a=1.00000 b=1.00000 c=1.00000
x1=-0.5 + 0.8660254037844386i
x2=-0.5 + -0.8660254037844386i</pre>
 
println!("Result : {:?}", roots);
=={{header|Scheme}}==
}</syntaxhighlight>
<lang scheme>(define (quadratic a b c)
{{out}}
(if (= a 0)
<pre>
(if (= b 0) 'fail (- (/ c b)))
Result : Three([0.000000059604645, 0.99999994, 2.0])
(let ((delta (- (* b b) (* 4 a c))))
</pre>
(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))))))
 
Another without external crates:
<syntaxhighlight lang="rust">
use num::Float;
 
/// Note: We cannot use `range_step` here because Floats don't implement
; examples
/// the `CheckedAdd` trait.
fn find_roots<T, F>(f: F, start: T, stop: T, step: T, epsilon: T) -> Vec<T>
where
T: Copy + PartialOrd + Float,
F: Fn(T) -> T,
{
let mut ret = vec![];
let mut current = start;
while current < stop {
if f(current).abs() < epsilon {
ret.push(current);
}
current = current + step;
}
ret
}
 
fn main() {
(quadratic 1 -1 -1)
let roots = find_roots(
; (1.618033988749895 -0.6180339887498948)
|x: f64| x * x * x - 3.0 * x * x + 2.0 * x,
-1.0,
3.0,
0.0001,
0.00000001,
);
 
println!("roots of f(x) = x^3 - 3x^2 + 2x are: {:?}", roots);
(quadratic 1 0 -2)
}
; (-1.4142135623730951 1.414213562373095)
 
</syntaxhighlight>
(quadratic 1 0 2)
{{out}}
; (0+1.4142135623730951i 0-1.4142135623730951i)
<pre>
roots of f(x) = x^3 - 3x^2 + 2x are: [-0.00000000000009381755897326649, 0.9999999999998124, 1.9999999999997022]
</pre>
 
=={{header|Scala}}==
(quadratic 1+1i 2 5)
===Imperative version (Ugly, side effects)===
; (-1.0922677260818898-1.1884256155834088i 0.09226772608188982+2.1884256155834088i)
{{trans|Java}}
{{Out}}Best seen running in your browser either by [https://scalafiddle.io/sf/T63KUsH/0 (ES aka JavaScript, non JVM)] or [https://scastie.scala-lang.org/bh8von94Q1y0tInvEZ3cBQ Scastie (remote JVM)].
<syntaxhighlight lang="scala">object Roots extends App {
val poly = (x: Double) => x * x * x - 3 * x * x + 2 * x
 
private def printRoots(f: Double => Double,
(quadratic 0 4 3)
lowerBound: Double,
; -3/4
upperBound: Double,
step: Double): Unit = {
val y = f(lowerBound)
var (ox, oy, os) = (lowerBound, y, math.signum(y))
 
for (x <- lowerBound to upperBound by step) {
(quadratic 0 0 1)
val y = f(x)
; fail
val s = math.signum(y)
if (s == 0) println(x)
else if (s != os) println(s"~${x - (x - ox) * (y / (y - oy))}")
 
ox = x
(quadratic 1 2 0)
oy = y
; (-2 0)
os = s
}
}
 
printRoots(poly, -1.0, 4, 0.002)
(quadratic 1 2 1)
; (-1 -1)
 
}</syntaxhighlight>
(quadratic 1 -1e5 1)
===Functional version (Recommended)===
; (99999.99999 1.0000000001000001e-05)</lang>
<syntaxhighlight lang="scala">object RootsOfAFunction extends App {
def findRoots(fn: Double => Double, start: Double, stop: Double, step: Double, epsilon: Double) = {
for {
x <- start to stop by step
if fn(x).abs < epsilon
} yield x
}
 
def fn(x: Double) = x * x * x - 3 * x * x + 2 * x
=={{header|Seed7}}==
{{trans|Ada}}
<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;</lang>
 
println(findRoots(fn, -1.0, 3.0, 0.0001, 0.000000001))
}</syntaxhighlight>
{{out}}
Vector(-9.381755897326649E-14, 0.9999999999998124, 1.9999999999997022)
<pre>
X1 = 1000000.000000 X2 = 0.000001
</pre>
 
=={{header|Sidef}}==
<syntaxhighlight lang ="ruby">varfunc sets =f(x) [{
x*x*x - 3*x*x [1, + 2, 1],*x
[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)]
}
 
 
var step = 0.001
sets.each { |coefficients|
var start = -1
say ("Roots for #{coefficients}",
var stop = 3
"=> (#{quadroots(coefficients...).join(', ')})")
 
}</lang>
for x in range(start+step, stop, step) {
static sign = false
given (var value = f(x)) {
when (0) {
say "Root found at #{x}"
}
case (sign && ((value > 0) != sign)) {
say "Root found near #{x}"
}
}
sign = value>0
}</syntaxhighlight>
{{out}}
<pre>Root found at 0
Root found at 1
Roots for [1, 2, 1]=> (-1, -1)
Root found at 2</pre>
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}}==
<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 |
+-------------------------------+</lang>
 
=={{header|Tcl}}==
This simple brute force iteration marks all results, with a leading "~", as approximate. This version always reports its results as approximate because of the general limits of computation using fixed-width floating-point numbers (i.e., IEEE double-precision floats).
{{tcllib|math::complexnumbers}}
<syntaxhighlight lang="tcl">proc froots {lambda {start -3} {end 3} {step 0.0001}} {
<lang tcl>package require math::complexnumbers
set res {}
namespace import math::complexnumbers::complex math::complexnumbers::tostring
set lastsign [sgn [apply $lambda $start]]
 
for {set x $start} {$x <= $end} {set x [expr {$x + $step}]} {
proc quadratic {a b c} {
set discrimsign [exprsgn {$b**2[apply - 4*$a*lambda $c}x]]
if {$sign != $lastsign} {
set roots [list]
lappend res [format ~%.11f $x]
if {$discrim < 0} {
set term1 [expr {(-1.0*$b)/(2*$a)}]
set term2lastsign [expr {sqrt(abs($discrim))/(2*$a)}]sign
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 $rootsres
}
proc sgn x {expr {($x>0) - ($x<0)}}
 
puts [froots {x {expr {$x**3 - 3*$x**2 + 2*$x}}}]</syntaxhighlight>
proc report_quad {a b c} {
Result and timing:
puts [format "%sx**2 + %sx + %s = 0" $a $b $c]
<pre>/Tcl $ time ./froots.tcl
foreach root [quadratic $a $b $c] {
~0.00000000000 ~1.00000000000 ~2.00000000000
puts " x = $root"
 
real 0m0.368s
user 0m0.062s
sys 0m0.030s</pre>
A more elegant solution (and faster, because you can usually make the initial search coarser) is to use brute-force iteration and then refine with [[wp:Newton's method|Newton-Raphson]], but that requires the differential of the function with respect to the search variable.
<syntaxhighlight lang="tcl">proc frootsNR {f df {start -3} {end 3} {step 0.001}} {
set res {}
set lastsign [sgn [apply $f $start]]
for {set x $start} {$x <= $end} {set x [expr {$x + $step}]} {
set sign [sgn [apply $f $x]]
if {$sign != $lastsign} {
lappend res [format ~%.15f [nr $x $f $df]]
}
set lastsign $sign
}
return $res
}
proc sgn x {expr {($x>0) - ($x<0)}}
proc nr {x1 f df} {
# Newton's method converges very rapidly indeed
for {set iters 0} {$iters < 10} {incr iters} {
set x1 [expr {
[set x0 $x1] - [apply $f $x0]/[apply $df $x0]
}]
if {$x0 == $x1} {
break
}
}
return $x1
}
 
puts [frootsNR \
# examples on this page
{x {expr {$x**3 - 3*$x**2 + 2*$x}}} \
report_quad 3 4 [expr {4/3.0}] ;# {-2/3}
{x {expr {3*$x**2 - 6*$x + 2}}}]</syntaxhighlight>
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>
{{out}}
<pre>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</pre>
 
=={{header|TI-89 BASIC}}==
 
Finding roots is a built-in function: <code>zeros(x^3-3x^2+2x, x)</code> returns <code>{0,1,2}</code>.
TI-89 BASIC has built-in numeric and algebraic solvers.
 
<lang>solve(x^2-1E9x+1.0)</lang>
In this case, the roots are exact; inexact results are marked by decimal points.
returns
<pre>x=1.E-9 or x=1.E9</pre>
 
=={{header|Wren}}==
{{trans|Go}}
{{libheader|Wren-complexfmt}}
<langsyntaxhighlight lang="ecmascript">import "/complexfmt" for ComplexFmt
 
var quadraticsecant = Fn.new { |af, bx0, cx1|
var df0 = b*b - 4*a*c0
ifvar (df1 == 0f.call(x0) {
for (i in 0...100) // single root{
returnf0 [[-b/(2*a)],= null]f1
f1 = f.call(x1)
if (f1 == 0) return [x1, "exact"]
if ((x1-x0).abs < 1e-6) return [x1, "approximate"]
var t = x0
x0 = x1
x1 = x1-f1*(x1-t)/(f1-f0)
}
if (d >return [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 testfindRoots = Fn.new { |af, blower, cupper, step|
var x0 = lower
System.write("coefficients: %(a), %(b), %(c) -> ")
var rootsx1 = quadratic.call(a,lower b,+ c)step
varwhile r(x0 =< roots[0]upper) {
if (r.count = x1 = 1(x1 < upper) {? x1 : upper
System.print("onevar realres root:= %secant.call(r[0])"f, x0, x1)
} else if ( var r.count == 2) {res[0]
System.print("twovar realstatus roots:= %(r[0]) and %(rres[1])")
if (status != "" && r >= x0 && r < x1) {
} else {
var i = roots[1] Fmt.print(" $6.3f $s", r, status)
}
System.print("two complex roots: %(i[0]) and %(i[1])")
x0 = x1
x1 = x1 + step
}
}
 
var example = Fn.new { |x| x*x*x - 3*x*x + 2*x }
var coeffs = [
findRoots.call(example, -0.5, 2.6, 1)</syntaxhighlight>
[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])</lang>
 
{{out}}
<pre>
0.000 approximate
coefficients: 1, -2, 1 -> one real root: 1
1.000 exact
coefficients: 1, 0, 1 -> two complex roots: 0 + i and 0 - i
2.000 approximate
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|zkl}}==
{{trans|Haskell}}
zkl doesn't have a complex number package.
<syntaxhighlight lang="zkl">fcn findRoots(f,start,stop,step,eps){
{{trans|Elixir}}
[start..stop,step].filter('wrap(x){ f(x).closeTo(0.0,eps) })
<lang zkl>fcn quadratic(a,b,c){ b=b.toFloat();
}</syntaxhighlight>
println("Roots of a quadratic function %s, %s, %s".fmt(a,b,c));
<syntaxhighlight lang="zkl">fcn f(x){ x*x*x - 3.0*x*x + 2.0*x }
d,a2:=(b*b - 4*a*c), a+a;
findRoots(f, -1.0, 3.0, 0.0001, 0.00000001).println();</syntaxhighlight>
if(d>0){
{{out}}
sd:=d.sqrt();
<pre>L(-9.38176e-14,1,2)</pre>
println(" the real roots are %s and %s".fmt((-b + sd)/a2,(-b - sd)/a2));
{{trans|C}}
<syntaxhighlight lang="zkl">fcn secant(f,xA,xB){
reg e=1.0e-12;
 
fA:=f(xA); if(fA.closeTo(0.0,e)) return(xA);
 
do(50){
fB:=f(xB);
d:=(xB - xA) / (fB - fA) * fB;
if(d.closeTo(0,e)) break;
xA = xB; fA = fB; xB -= d;
}
if(f(xB).closeTo(0.0,e)) xB
else if(d==0) println(" the single root is ",-b/a2);
else "Function is not converging near (%7.4f,%7.4f).".fmt(xA,xB);
else{
}</syntaxhighlight>
sd:=(-d).sqrt();
<syntaxhighlight lang="zkl">step:=0.1;
println(" the complex roots are %s and \U00B1;%si".fmt(-b/a2,sd/a2));
xs:=findRoots(f, -1.032, 3.0, step, 0.1);
}
xs.println(" --> ",xs.apply('wrap(x){ secant(f,x-step,x+step) }));</syntaxhighlight>
}</lang>
<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)
}</lang>
{{out}}
<pre>L(-0.032,0.968,1.068,1.968) --> L(1.87115e-19,1,1,2)</pre>
<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}}
Thundergnat
10,327

edits

Retrieved from "https://rosettacode.org/wiki/Special:MobileDiff/328301"