Roots of a function: Difference between revisions
Added Easylang
m (oct...) |
(Added Easylang) |
||
| (207 intermediate revisions by 86 users not shown) | |||
Line 1:
{{task|Arithmetic operations}}
;Task:
Create a program that finds and outputs the roots of a given function, range and (if applicable) step width.
The program should identify whether the root is exact or approximate.
For this task, use: <big><big> ƒ(x) = x<sup>3</sup> - 3x<sup>2</sup> + 2x </big></big>
<br><br>
=={{header|11l}}==
{{trans|Python}}
<syntaxhighlight lang="11l">F f(x)
R x^3 - 3 * x^2 + 2 * x
-V step = 0.001
-V start = -1.0
-V stop = 3.0
V sgn = f(start) > 0
V x = start
L x <= stop
V value = f(x)
I value == 0
print(‘Root found at ’x)
E I (value > 0) != sgn
print(‘Root found near ’x)
sgn = value > 0
x += step</syntaxhighlight>
{{out}}
<pre>
Root found near 8.812395258e-16
Root found near 1
Root found near 2.001
</pre>
=={{header|Ada}}==
<
procedure Roots_Of_Function is
Line 42 ⟶ 80:
X := X + Step;
end loop;
end Roots_Of_Function;</
=={{header|ALGOL 68}}==
{{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 to extensive use of FORMATted transput}}
Finding 3 roots using the secant method:
<syntaxhighlight lang="algol68">MODE DBL = LONG REAL;
FORMAT dbl = $g(-long real width, long real width-6, -2)$;
Line 116 ⟶ 157:
)
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|Arturo}}==
<syntaxhighlight lang="arturo">f: function [n]->
((n^3) - 3*n^2) + 2*n
step: 0.01
start: neg 1.0
stop: 3.0
sign: positive? f start
x: start
while [x =< stop][
value: f x
if? value = 0 ->
print ["root found at" to :string .format:".5f" x]
else ->
if sign <> value > 0 -> print ["root found near" to :string .format:".5f" x]
sign: value > 0
'x + step
]</syntaxhighlight>
{{out}}
<pre>root found near 0.00000
root found near 1.00000
root found near 2.00000</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:
* roots() searches for intervals within given limits, shifted by a given “step”, where our function has different signs at the endpoints.
* Having found such an interval, the root() function searches for a value where our function is 0, within a given tolerance.
* It also sets ErrorLevel to info about the root found.
[http://www.autohotkey.com/forum/viewtopic.php?t=44657&postdays=0&postorder=asc&start=139 discussion]
<syntaxhighlight lang="autohotkey">MsgBox % roots("poly", -0.99, 2, 0.1, 1.0e-5)
MsgBox % roots("poly", -1, 3, 0.1, 1.0e-5)
roots(f,x1,x2,step,tol) { ; search for roots in intervals of length "step", within tolerance "tol"
x := x1, y := %f%(x), s := (y>0)-(y<0)
Loop % ceil((x2-x1)/step) {
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
Return res
}
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
}
}
poly(x) {
Return ((x-3)*x+2)*x
}</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"> function$ = "x^3-3*x^2+2*x"
rangemin = -1
rangemax = 3
stepsize = 0.001
accuracy = 1E-8
PROCroots(function$, rangemin, rangemax, stepsize, accuracy)
END
DEF PROCroots(func$, min, max, inc, eps)
LOCAL x, sign%, oldsign%
oldsign% = 0
FOR x = min TO max STEP inc
sign% = SGN(EVAL(func$))
IF sign% = 0 THEN
PRINT "Root found at x = "; x
sign% = -oldsign%
ELSE IF sign% <> oldsign% AND oldsign% <> 0 THEN
IF inc < eps THEN
PRINT "Root found near x = "; x
ELSE
PROCroots(func$, x-inc, x+inc/8, inc/8, eps)
ENDIF
ENDIF
ENDIF
oldsign% = sign%
NEXT x
ENDPROC</syntaxhighlight>
Output:
<pre>Root found near x = 2.29204307E-9
Root found near x = 1
Root found at x = 2</pre>
=={{header|C}}==
=== Secant Method ===
<syntaxhighlight lang="c">#include <math.h>
#include <stdio.h>
double f(double x)
{
return x*x*x-3.0*x*x +2.0*x;
}
double secant( double xA, double xB, double(*f)(double) )
{
double e = 1.0e-12;
double fA, fB;
double d;
int i;
int limit = 50;
fA=(*f)(xA);
for (i=0; i<limit; i++) {
fB=(*f)(xB);
d = (xB - xA) / (fB - fA) * fB;
if (fabs(d) < e)
break;
xA = xB;
fA = fB;
xB -= d;
}
if (i==limit) {
printf("Function is not converging near (%7.4f,%7.4f).\n", xA,xB);
return -99.0;
}
return xB;
}
int main(int argc, char *argv[])
{
double step = 1.0e-2;
double e = 1.0e-12;
double x = -1.032; // just so we use secant method
double xx, value;
int s = (f(x)> 0.0);
while (x < 3.0) {
value = f(x);
if (fabs(value) < e) {
printf("Root found at x= %12.9f\n", x);
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 ===
<syntaxhighlight lang="c">#include <gsl/gsl_poly.h>
#include <stdio.h>
int main(int argc, char *argv[])
{
/* 0 + 2x - 3x^2 + 1x^3 */
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#}}==
{{trans|C++}}
<syntaxhighlight lang="csharp">using System;
class Program
{
public static void Main(string[] args)
{
Func<double, double> f = x => { return x * x * x - 3 * x * x + 2 * x; };
double step = 0.001; // Smaller step values produce more accurate and precise results
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
{
private static int Sign(double x)
{
return x < 0.0 ? -1 : x > 0.0 ? 1 : 0;
}
public static void PrintRoots(Func<double, double> f, double lowerBound,
double upperBound, double step)
{
double x = lowerBound, ox = x;
double y = f(x), oy = y;
int s = Sign(y), os = s;
for (; x <= upperBound; x += step)
{
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;
oy = y;
os = s;
}
}
public static void Main(string[] args)
{
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===
{{trans|C++}}
<syntaxhighlight lang="csharp">using System;
class Program
{
public static void Main(string[] args)
{
Func<double, double> f = x => { return x * x * x - 3 * x * x + 2 * x; };
double root = BrentsFun(f, lower: -1.0, upper: 4, tol: 0.002, maxIter: 100);
}
private static void Swap<T>(ref T a, ref T b)
{
var tmp = a;
a = b;
b = tmp;
}
public static double BrentsFun(Func<double, double> f, double lower, double upper, double tol, uint maxIter)
{
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;
if (!(fa * fb < 0))
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)
{
Swap(ref a, ref b);
Swap(ref fa, ref fb);
}
double c = a; // c now equals the largest magnitude of the lower and upper bounds
double fc = fa; // precompute function evalutation for point c by assigning it the same value as fa
bool mflag = true; // boolean flag used to evaluate if statement later 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)
{
// 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++}}==
<
double f(double x)
Line 160 ⟶ 664:
sign = ( value > 0 );
}
}</
===Brent's Method===
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.
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).
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.
<syntaxhighlight lang="cpp">#include <iostream>
#include <cmath>
#include <algorithm>
#include <functional>
double brents_fun(std::function<double (double)> f, double lower, double upper, double tol, unsigned int max_iter)
{
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))
{
std::cout << "Signs of f(lower_bound) and f(upper_bound) must be opposites" << std::endl; // throws exception if root isn't bracketed
return -11;
}
if (std::abs(fa) < std::abs(b)) // if magnitude of f(lower_bound) is less than magnitude of f(upper_bound)
{
std::swap(a,b);
std::swap(fa,fb);
}
double c = a; // c now equals the largest magnitude of the lower and upper bounds
double fc = fa; // precompute function evalutation for point c by assigning it the same value as fa
bool mflag = true; // boolean flag used to evaluate if statement later on
double s = 0; // Our Root that will be returned
double d = 0; // Only used if mflag is unset (mflag == false)
for (unsigned int iter = 1; iter < max_iter; ++iter)
{
// stop if converged on root or error is less than tolerance
if (std::abs(b-a) < tol)
{
std::cout << "After " << iter << " iterations the root is: " << s << std::endl;
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 && (std::abs(s-b) >= (std::abs(b-c) * 0.5)) ) ||
( !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;
}
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 (std::abs(fa) < std::abs(fb)) // if magnitude of fa is less than magnitude of fb
{
std::swap(a,b); // swap a and b
std::swap(fa,fb); // make sure f(a) and f(b) are correct after swap
}
} // 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>
<pre>
> (findRoots #(+ (* % % %) (* -3 % %) (* 2 %)) -1.0 3.0 0.0001 0.00000001)
(-9.381755897326649E-14 0.9999999999998124 1.9999999999997022)
</pre>
=={{header|CoffeeScript}}==
{{trans|Python}}
<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 ->
# Smaller steps produce more accurate/precise results in general,
# but for many functions we'll never get exact roots, either due
# to imperfect binary representation or irrational roots.
step = 1 / 256
f1 = (x) -> x*x*x - 3*x*x + 2*x
print_roots f1, -1, 5, step
f2 = (x) -> x*x - 4*x + 3
print_roots f2, -1, 5, step
f3 = (x) -> x - 1.5
print_roots f3, 0, 4, step
f4 = (x) -> x*x - 2
print_roots f4, -2, 2, step
</syntaxhighlight>
output
<syntaxhighlight lang="text">
> coffee roots.coffee
-----
Root found at 0
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}}==
{{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.
<syntaxhighlight lang="lisp">(defun find-roots (function start end &optional (step 0.0001))
(let* ((roots '())
(value (funcall function start))
(plusp (plusp value)))
(when (zerop value)
(format t "~&Root found at ~W." start))
(do ((x (+ start step) (+ x step)))
((> x end) (nreverse roots))
(setf value (funcall function x))
(cond
((zerop value)
(format t "~&Root found at ~w." x)
(push x roots))
((not (eql plusp (plusp value)))
(format t "~&Root found near ~w." x)
(push (cons (- x step) x) roots)))
(setf plusp (plusp value)))))</syntaxhighlight>
<pre>> (find-roots #'(lambda (x) (+ (* x x x) (* -3 x x) (* 2 x))) -1 3)
Root found near 5.3588345E-5.
Root found near 1.0000072.
Root found near 2.000073.
((-4.6411653E-5 . 5.3588345E-5)
(0.99990714 . 1.0000072)
(1.9999729 . 2.000073))</pre>
=={{header|D}}==
<syntaxhighlight lang="d">import std.stdio, std.math, std.algorithm;
return abs(a) <= b;
}
T[] findRoot(T)(immutable T function(in T) pure nothrow fi,
in T start, in T end, in T step=T(0.001L),
T tolerance = T(1e-4L)) {
if (step.nearZero)
writefln("WARNING: step size may be too small.");
/// Search root by simple bisection.
T searchRoot(T a, T b) pure nothrow {
T root;
int limit = 49;
T gap = b - a;
if (fi(a).nearZero)
return a;
if (fi(b).nearZero)
return b;
root = (b + a) / 2.0L;
if (fi(root).nearZero)
return root;
((fi(a) * fi(root) < 0) ? b : a) = root;
gap = b - a;
}
}
immutable dir = T(end > start ? 1.0 : -1.0);
immutable step2 = (end > start) ? abs(step) : -abs(step);
T[T] result;
immutable T r = searchRoot(x, x + step2);
return result.keys.sort
}
void report(T)(in T[] r, immutable T function(in T) pure f,
in T tolerance = T(1e-4L)) {
if (r.length) {
writefln("Root found (tolerance = %1.4g):", tolerance);
foreach (const x; r) {
immutable T y = f(x);
if (nearZero(y))
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() {
static real f(in real x) pure nothrow {
return x ^^ 3 - (3 * x ^^ 2) + 2 * x;
}
findRoot(&f, -1.0L, 3.0L, 0.001L).report(&f);
}</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}}==
{{trans|Scala}}
<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* {
for (double x = start; x < stop; x = x + step) {
if (fn(x).abs() < epsilon) yield x;
}
}
main() {
// Vector(-9.381755897326649E-14, 0.9999999999998124, 1.9999999999997022)
print(findRoots(fn, -1.0, 3.0, 0.0001, 0.000000001));
}</syntaxhighlight>
=={{header|Delphi}}==
See [https://rosettacode.org/wiki/Roots_of_a_function#Pascal Pascal].
=={{header|DWScript}}==
{{trans|C}}
<syntaxhighlight lang="delphi">type TFunc = function (x : Float) : Float;
function f(x : Float) : Float;
begin
Result := x*x*x-3.0*x*x +2.0*x;
end;
const e = 1.0e-12;
function Secant(xA, xB : Float; f : TFunc) : Float;
const
limit = 50;
var
fA, fB : Float;
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;
var x := -1.032; // just so we use secant method
var xx, value : Float;
var s := f(x)>0.0;
while (x < 3.0) do begin
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|EasyLang}}==
<syntaxhighlight>
func f x .
return x * x * x - 3 * x * x + 2 * x
.
numfmt 6 0
proc findroot start stop step . .
x = start
while x <= stop
val = f x
if val = 0
print x & " (exact)"
elif sign val <> sign0 and sign0 <> 0
print x & " (err = " & step & ")"
.
sign0 = sign val
x += step
.
.
proc drawfunc start stop . .
linewidth 0.3
drawgrid
x = start
while x <= stop
line x * 10 + 50 f x * 10 + 50
x += 0.1
.
.
drawfunc -1 3
findroot -1 3 pow 2 -20
print ""
findroot -1 3 1e-6
</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>
{{out}}
<pre>
Root found at 8.81239525796218e-16
Root found at 1.0000000000000016
Root found at 1.9999999999998914
</pre>
=={{header|Erlang}}==
<syntaxhighlight lang="erlang">% Implemented by Arjun Sunel
-module(roots).
-export([main/0]).
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) ->
Value = F(X),
if
Value == 0 -> % We hit a root
io:format("Root found at ~p~n",[X]),
while(X+Step, Step, Start, Stop, Value > 0,F);
(Value < 0) == Sign -> % We passed a root
io:format("Root found near ~p~n",[X]),
while(X+Step , Step, Start, Stop, Value > 0,F);
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}}==
<syntaxhighlight lang="erre">
PROGRAM ROOTS_FUNCTION
!VAR E,X,STP,VALUE,S%,I%,LIMIT%,X1,X2,D
FUNCTION F(X)
F=X*X*X-3*X*X+2*X
END FUNCTION
BEGIN
X=-1
STP=1.0E-6
E=1.0E-9
S%=(F(X)>0)
PRINT("VERSION 1: SIMPLY STEPPING X")
WHILE X<3.0 DO
VALUE=F(X)
IF ABS(VALUE)<E THEN
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
PRINT("VERSION 2: SECANT METHOD")
X1=-1.0
X2=3.0
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
EXACT ROOT FOUND AT X = 1
</pre>
=={{header|Fortran}}==
{{works with|Fortran|90 and later}}
<syntaxhighlight 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 <
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
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) = 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.
<syntaxhighlight lang="go">package main
import (
"fmt"
"math"
)
func main() {
example := func(x float64) float64 { return x*x*x - 3*x*x + 2*x }
findroots(example, -.5, 2.6, 1)
}
func findroots(f func(float64) float64, lower, upper, step float64) {
for x0, x1 := lower, lower+step; x0 < upper; x0, x1 = x1, x1+step {
x1 = math.Min(x1, upper)
r, status := secant(f, x0, x1)
if status != "" && r >= x0 && r < x1 {
fmt.Printf(" %6.3f %s\n", r, status)
}
}
}
func secant(f func(float64) float64, x0, x1 float64) (float64, string) {
var f0 float64
f1 := f(x0)
for i := 0; i < 100; i++ {
f0, f1 = f1, f(x1)
switch {
case f1 == 0:
return x1, "exact"
case math.Abs(x1-x0) < 1e-6:
return x1, "approximate"
}
x0, x1 = x1, x1-f1*(x1-x0)/(f1-f0)
}
return 0, ""
}</syntaxhighlight>
Output:
<pre>
0.000 approximate
1.000 exact
2.000 approximate
</pre>
=={{header|Haskell}}==
<syntaxhighlight lang="haskell">f x = x^3-3*x^2+2*x
findRoots start stop step eps =
[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.
<syntaxhighlight lang="haskell">import Numeric.GSL.Polynomials
import Data.Complex
*Main> mapM_ print $ polySolve [0,2,-3,1]
(-5.421010862427522e-20) :+ 0.0
2.000000000000001 :+ 0.0
0.9999999999999996 :+ 0.0</syntaxhighlight>
No complex roots, so:
<syntaxhighlight lang="haskell">*Main> mapM_ (print.realPart) $ polySolve [0,2,-3,1]
-5.421010862427522e-20
2.000000000000001
0.9999999999999996</syntaxhighlight>
It is possible to solve the problem directly and elegantly using robust bisection method and Alternative type class.
<syntaxhighlight lang="haskell">import Control.Applicative
data Root a = Exact a | Approximate a deriving (Show, Eq)
-- looks for roots on an interval
bisection :: (Alternative f, Floating a, Ord a) =>
(a -> a) -> a -> a -> f (Root a)
bisection f a b | f a * f b > 0 = empty
| 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
findRoots :: (Alternative f, Floating a, Ord a) =>
(a -> a) -> [a] -> а (Root a)
findRoots f [] = empty
findRoots f [x] = if f x == 0 then pure (Exact x) else empty
findRoots f (a:b:xs) = bisection f a b <|> findRoots f (b:xs)</syntaxhighlight>
It is possible to use these functions with different Alternative functors: IO, Maybe or List:
<pre>λ> bisection (\x -> x*x-2) 1 2
Approximate 1.414213562373094
λ> bisection (\x -> x-1) 1 2
Exact 1.0
λ> bisection (\x -> x*x-2) 2 3 :: Maybe (Root Double)
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`
<pre>λ> Data.List.nub $ findRoots (\x -> x^3 - 3*x^2 + 2*x) [-3..3]
[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}}==
HicEst's [http://www.HicEst.com/SOLVE.htm SOLVE] function employs the Levenberg-Marquardt method:
<syntaxhighlight lang="hicest">OPEN(FIle='test.txt')
1 DLG(NameEdit=x0, DNum=3)
x = x0
chi2 = SOLVE(NUL=x^3 - 3*x^2 + 2*x, Unknown=x, I=iterations, NumDiff=1E-15)
EDIT(Text='approximate exact ', Word=(chi2 == 0), Parse=solution)
WRITE(FIle='test.txt', LENgth=6, Name) x0, x, solution, chi2, iterations
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|Icon}} and {{header|Unicon}}==
{{trans|Java}}
Works in both languages:
<syntaxhighlight lang="unicon">procedure main()
showRoots(f, -1.0, 4, 0.002)
end
procedure f(x)
return x^3 - 3*x^2 + 2*x
end
procedure showRoots(f, lb, ub, step)
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)
return (x<0, -1) | (x>0, 1) | 0
end</syntaxhighlight>
Output:
<pre>
->roots
~2.616794878713638e-18
~1.0
~2.0
->
</pre>
=={{header|J}}==
J has builtin a root-finding operator, '''<tt>p.</tt>''', whose input is the
<syntaxhighlight lang="j"> 1{::p. 0 2 _3 1
We can determine whether the roots are exact or approximate by evaluating the polynomial at the candidate roots, and testing for zero:
<syntaxhighlight lang="j"> (0=]p.1{::p.) 0 2 _3 1
As you can see, <tt>p.</tt> is also the operator which evaluates polynomials. This is not a coincidence.
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.
<syntaxhighlight lang="j"> blackbox=: 0 2 _3 1&p.
(#~ (=<./)@:|@blackbox) i.&.(1e6&*)&.(1&+) 3
0 1 2
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).
=={{header|Java}}==
<
public interface Function {
public double f(double x);
}
}
public static void printRoots(Function f, double lowerBound,
double upperBound, double step) {
double x = lowerBound, ox = x;
double y = f.f(x), oy = y;
int s = sign(y), os = s;
for (; x <= upperBound ; x += step) {
s = sign(y = f.f(x));
if (s == 0) {
System.out.println(x);
} else if (s != os) {
double dx = x - ox;
double dy = y - oy;
double cx = x - dx * (y / dy);
System.out.println("~" + cx);
}
ox = x; oy = y; os = s;
}
}
Function poly = new Function () {
public double f(double x) {
return x*x*x - 3*x*x + 2*x;
}
};
printRoots(poly, -1.0, 4, 0.002);
}
}</syntaxhighlight>
Produces this output:
<pre>~2.616794878713638E-18
~1.0000000000000002
~2.000000000000001</pre>
=={{header|JavaScript}}==
{{trans|Java}}
{{works with|SpiderMonkey|22}}
{{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) {
return (x < 0.0) ? -1 : (x > 0.0) ? 1 : 0;
}
function printRoots(f, lowerBound, upperBound, step) {
var x = lowerBound, ox = x,
y = f(x), oy = y,
s = sign(y), os = s;
for (; x <= upperBound ; x += step) {
s = sign(y = f(x));
if (s == 0) {
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);
</syntaxhighlight>
=={{header|jq}}==
printRoots(f; lower; upper; step) finds approximations to the roots
of an arbitrary continuous real-valued function, f, in the range
[lower, upper], assuming step is small enough.
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 function, f, may be an expression (as in the example below) or a defined filter.
printRoots/3 emits an array of results, each of which is either a
number (representing an exact root within the limits of machine arithmetic) or a string consisting of "~" followed by an approximation to the root.
<syntaxhighlight lang="jq">def sign:
if . < 0 then -1 elif . > 0 then 1 else 0 end;
def printRoots(f; lowerBound; upperBound; step):
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)
[
0,
1,
2
]
printRoots( .*.*. - 3*.*. + 2*.; -1.0; 4; .001)
[
"~1.320318770141425e-18",
"~1.0000000000000002",
"~1.9999999999999993"
]</syntaxhighlight>
=={{header|Julia}}==
Assuming that one has the Roots package installed:
<syntaxhighlight lang="julia">using Roots
println(find_zero(x -> x^3 - 3x^2 + 2x, (-100, 100)))</syntaxhighlight>
{{out}}
<pre>[0.0,1.0,2.0]</pre>
Without the Roots package, Newton's method may be defined in this manner:
<syntaxhighlight lang="julia">function newton(f, fp, x::Float64,tol=1e-14::Float64,maxsteps=100::Int64)
##f: the function of x
##fp: the derivative of f
local xnew, xold = x, Inf
local fn, fo = f(xnew), Inf
local counter = 1
while (counter < maxsteps) && (abs(xnew - xold) > tol) && ( abs(fn - fo) > tol )
x = xnew - f(xnew)/fp(xnew) ## update x
xnew, xold = x, xnew
fn, fo = f(xnew), fn
counter += 1
end
if counter >= maxsteps
error("Did not converge in ", string(maxsteps), " steps")
else
xnew, counter
end
end
</syntaxhighlight>
Finding the roots of f(x) = x3 - 3x2 + 2x:
<syntaxhighlight lang="julia">
f(x) = x^3 - 3*x^2 + 2*x
fp(x) = 3*x^2-6*x+2
x_s, count = newton(f,fp,1.00)
</syntaxhighlight>
{{out}}
(1.0,2)
=={{header|Kotlin}}==
{{trans|C}}
<syntaxhighlight lang="scala">// version 1.1.2
typealias DoubleToDouble = (Double) -> Double
fun f(x: Double) = x * x * x - 3.0 * x * x + 2.0 * x
fun secant(x1: Double, x2: Double, f: DoubleToDouble): Double {
val e = 1.0e-12
val limit = 50
var xa = x1
var xb = x2
var fa = f(xa)
var i = 0
while (i++ < limit) {
var fb = f(xb)
val d = (xb - xa) / (fb - fa) * fb
if (Math.abs(d) < e) break
xa = 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 e = 1.0e-12
var x = -1.032
var s = f(x) > 0.0
while (x < 3.0) {
val value = f(x)
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>
{{out}}
<pre>
Root found at x = 0.000000000
Root found at x = 1.000000000
Root found at x = 2.000000000
</pre>
=={{header|Lambdatalk}}==
<syntaxhighlight lang="scheme">
1) defining the function:
{def func {lambda {:x} {+ {* 1 :x :x :x} {* -3 :x :x} {* 2 :x}}}}
-> func
2) printing roots:
{S.map {lambda {:x}
{if {< {abs {func :x}} 0.0001}
then {br}- a root found at :x else}}
{S.serie -1 3 0.01}}
->
- a root found at 7.528699885739343e-16
- a root found at 1.0000000000000013
- a root found at 2.000000000000002
3) printing the roots of the "sin" function between -720° to +720°;
{S.map {lambda {:x}
{if {< {abs {sin {* {/ {PI} 180} :x}}} 0.01}
then {br}- a root found at :x° else}}
{S.serie -720 +720 10}}
->
- a root found at -720°
- a root found at -540°
- a root found at -360°
- a root found at -180°
- a root found at 0°
- a root found at 180°
- a root found at 360°
- a root found at 540°
- a root found at 720°
</syntaxhighlight>
=={{header|Liberty BASIC}}==
<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
mainwin 80 12
xMin =-1
xMax = 3
y =f( xMin) ' Since Liberty BASIC has an 'eval(' function the fn
' and limits would be better entered via 'input'.
LastY =y
eps =1E-12 ' closeness acceptable
bigH=0.01
print
print " Checking for roots of x^3 -3 *x^2 +2 *x =0 over range -1 to +3"
print
x=xMin: dx = bigH
do
x=x+dx
y = f(x)
'print x, dx, y
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)
f =x^3 -3 *x^2 +2 *x
end function</syntaxhighlight>
=={{header|Lua}}==
<syntaxhighlight lang="lua">-- Function to have roots found
function f (x) return x^3 - 3*x^2 + 2*x end
-- Find roots of f within x=[start, stop] or approximations thereof
function root (f, start, stop, step)
local roots, x, sign, foundExact, value = {}, start, f(start) > 0
while x <= stop do
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
print("Root (to 12DP)\tMax. Error\n")
for _, r in pairs(root(f, -1, 3, 10^-6)) do
print(string.format("%0.12f", r.val), r.err)
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}}==
<syntaxhighlight
outputs:
<syntaxhighlight
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).
Line 333 ⟶ 1,919:
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.
=={{header|Mathematica}}/{{header|Wolfram Language}}==
There are multiple obvious ways to do this in Mathematica.
===Solve===
This requires a full equation and will perform symbolic operations only:
Output
<pre> {{x->0},{x->1},{x->2}}</pre>
===NSolve===
This requires merely the polynomial and will perform numerical operations if needed:
Output
<pre> {{x->0.},{x->1.},{x->2.}}</pre>
(note that the results here are floats)
===FindRoot===
This will numerically try to find one(!) local root from a given starting point:
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===
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:
Output
<pre>{{x->0}}</pre>
===Reduce===
This will (symbolically) reduce a given expression to the simplest possible form, solving equations and performing substitutions in the process:
(note that this doesn't yield a "solution" but a different expression that expresses the same thing as the original)
=={{header|Maxima}}==
<syntaxhighlight lang="maxima">e: x^3 - 3*x^2 + 2*x$
/* Number of roots in a real interval, using Sturm sequences */
nroots(e, -10, 10);
3
solve(e, x);
[x=1, x=2, x=0]
/* 'solve sets the system variable 'multiplicities */
solve(x^4 - 2*x^3 + 2*x - 1, x);
[x=-1, x=1]
multiplicities;
[1, 3]
/* Rational approximation of roots using Sturm sequences and bisection */
realroots(e);
[x=1, x=2, x=0]
/* 'realroots also sets the system variable 'multiplicities */
multiplicities;
[1, 1, 1]
/* Numerical root using Brent's method (here with another equation) */
find_root(sin(t) - 1/2, t, 0, %pi/2);
0.5235987755983
fpprec: 60$
bf_find_root(sin(t) - 1/2, t, 0, %pi/2);
5.23598775598298873077107230546583814032861566562517636829158b-1
/* Numerical root using Newton's method */
load(newton1)$
newton(e, x, 1.1, 1e-6);
1.000000017531147
/* For polynomials, Jenkins–Traub algorithm */
allroots(x^3 + x + 1);
[x=1.161541399997252*%i+0.34116390191401,
x=0.34116390191401-1.161541399997252*%i,
x=-0.68232780382802]
bfallroots(x^3 + x + 1);
[x=1.16154139999725193608791768724717407484314725802151429063617b0*%i + 3.41163901914009663684741869855524128445594290948999288901864b-1,
x=3.41163901914009663684741869855524128445594290948999288901864b-1 - 1.16154139999725193608791768724717407484314725802151429063617b0*%i,
x=-6.82327803828019327369483739711048256891188581897998577803729b-1]</syntaxhighlight>
=={{header|Nim}}==
<syntaxhighlight lang="nim">import math
import strformat
func f(x: float): float = x ^ 3 - 3 * x ^ 2 + 2 * x
var
step = 0.01
start = -1.0
stop = 3.0
sign = f(start) > 0
x = start
while x <= stop:
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}}
<pre>
Root found near 0.00000
Root found near 1.00000
Root found near 2.00000
</pre>
=={{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}}==
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 =
((u > 0.0) && (v < 0.0)) || ((u < 0.0) && (v > 0.0));;
let xtol a b = (a = b);; (* or use |a-b| < epsilon *)
let rec regula_falsi a b fa fb f =
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 =
let rec next x fx =
if x > hi then [] else
let y = x +. step in
let fy = f y 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) =
Printf.printf "f(%.17f) = %.17f [%s]\n"
x fx (if fx = 0.0 then "exact" else "approx") in
let f x = ((x -. 3.0)*.x +. 2.0)*.x in
List.iter showroot (search (-5.0) 5.0 0.1 f);;</syntaxhighlight>
Output:
<pre>
f(0.00000000000000000) = 0.00000000000000000 [exact]
f(1.00000000000000022) = 0.00000000000000000 [exact]
f(1.99999999999999978) = 0.00000000000000000 [exact]
</pre>
Note these roots are exact solutions with floating-point calculation.
=={{header|Octave}}==
Line 371 ⟶ 2,134:
If the equation is a polynomial, we can put the coefficients in a vector and use ''roots'':
<
r = roots(a);
% let's print it
Line 381 ⟶ 2,144:
endif
printf(" exact)\n");
endfor</
Otherwise we can program our (simple) method:
{{trans|Python}}
<
y = x.^3 -3.*x.^2 + 2.*x;
endfunction
Line 406 ⟶ 2,169:
se = sign(v);
x = x + step;
endwhile</
=={{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>
findRoots(#f, -1, 3, 0.0001)
Root found at 0
Root found at 1
Root found at 2
findRoots(#f, -1.000001, 3, 0.0001)
Root near 9.90000000000713e-005
Root near 1.000099
Root near 2.000099
</pre>
=={{header|ooRexx}}==
<syntaxhighlight lang="oorexx">/* REXX program to solve a cubic polynom equation
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
Parse Arg a
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). -->
<syntaxhighlight lang="parigp">polroots(x^3-3*x^2+2*x)</syntaxhighlight>
===Newton's method===
This uses a modified version of the Newton–Raphson method.
<syntaxhighlight lang="parigp">polroots(x^3-3*x^2+2*x,1)</syntaxhighlight>
===Brent's method===
<syntaxhighlight lang="parigp">solve(x=-.5,.5,x^3-3*x^2+2*x)
solve(x=.5,1.5,x^3-3*x^2+2*x)
solve(x=1.5,2.5,x^3-3*x^2+2*x)</syntaxhighlight>
===Factorization to linear factors===
<syntaxhighlight lang="parigp">findRoots(P)={
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)")
)
)
)
)
};
findRoots(x^3-3*x^2+2*x)</syntaxhighlight>
===Factorization to quadratic factors===
Of course this process could be continued to degrees 3 and 4 with sufficient additional work.
<syntaxhighlight lang="parigp">findRoots(P)={
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}}
<syntaxhighlight lang="pascal">Program RootsFunction;
var
e, x, step, value: double;
s: boolean;
i, limit: integer;
x1, x2, d: double;
function f(const x: double): double;
begin
f := x*x*x - 3*x*x + 2*x;
end;
begin
x := -1;
step := 1.0e-6;
e := 1.0e-9;
s := (f(x) > 0);
writeln('Version 1: simply stepping x:');
while x < 3.0 do
begin
value := f(x);
if abs(value) < e then
begin
writeln ('root found at x = ', x);
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
if i > limit then
begin
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>
Output:
<pre>
Version 1: simply stepping x:
root found at x = 7.91830063542152E-012
root found at x = 1.00000000001584E+000
root found at x = 1.99999999993357E+000
Version 2: secant method:
Exact root found at x = 1.00000000000000E+000
</pre>
=={{header|Perl}}==
<
{
my $x = shift;
Line 446 ⟶ 2,464:
# Update our sign
$sign = ( $value > 0 );
}</
=={{header|Phix}}==
{{trans|CoffeeScript}}
<!--<syntaxhighlight lang="phix">(phixonline)-->
<span style="color: #008080;">procedure</span> <span style="color: #000000;">print_roots</span><span style="color: #0000FF;">(</span><span style="color: #004080;">integer</span> <span style="color: #000000;">f</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;">--
-- Print approximate roots of f between x=start and x=stop, using
-- sign changes as an indicator that a root has been encountered.
--</span>
<span style="color: #004080;">atom</span> <span style="color: #000000;">x</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">start</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">y</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span>
<span style="color: #7060A8;">puts</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"-----\n"</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">while</span> <span style="color: #000000;">x</span><span style="color: #0000FF;"><=</span><span style="color: #000000;">stop</span> <span style="color: #008080;">do</span>
<span style="color: #004080;">atom</span> <span style="color: #000000;">last_y</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">y</span>
<span style="color: #000000;">y</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">f</span><span style="color: #0000FF;">(</span><span style="color: #000000;">x</span><span style="color: #0000FF;">)</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: #008080;">or</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">last_y</span><span style="color: #0000FF;"><</span><span style="color: #000000;">0</span> <span style="color: #008080;">and</span> <span style="color: #000000;">y</span><span style="color: #0000FF;">></span><span style="color: #000000;">0</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">or</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">last_y</span><span style="color: #0000FF;">></span><span style="color: #000000;">0</span> <span style="color: #008080;">and</span> <span style="color: #000000;">y</span><span style="color: #0000FF;"><</span><span style="color: #000000;">0</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">then</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"Root found %s %.10g\n"</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">{</span><span style="color: #008080;">iff</span><span style="color: #0000FF;">(</span><span style="color: #000000;">y</span><span style="color: #0000FF;">=</span><span style="color: #000000;">0</span><span style="color: #0000FF;">?</span><span style="color: #008000;">"at"</span><span style="color: #0000FF;">:</span><span style="color: #008000;">"near"</span><span style="color: #0000FF;">),</span><span style="color: #000000;">x</span><span style="color: #0000FF;">})</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #000000;">x</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">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,
-- but for many functions we'll never get exact roots, either due
-- to imperfect binary representation or irrational roots.</span>
<span style="color: #008080;">constant</span> <span style="color: #000000;">step</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">/</span><span style="color: #000000;">256</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">f1</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;">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>
<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>
-----
Root found at 0
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
</pre>
=={{header|PicoLisp}}==
{{trans|Clojure}}
<syntaxhighlight lang="picolisp">(de findRoots (F Start Stop Step Eps)
(filter
'((N) (> Eps (abs (F N))))
(range Start Stop Step) ) )
(scl 12)
(mapcar round
(findRoots
'((X) (+ (*/ X X X `(* 1.0 1.0)) (*/ -3 X X 1.0) (* 2 X)))
-1.0 3.0 0.0001 0.00000001 ) )</syntaxhighlight>
Output:
<pre>-> ("0.000" "1.000" "2.000")</pre>
=={{header|PL/I}}==
<syntaxhighlight lang="pl/i">
f: procedure (x) returns (float (18));
declare x float (18);
return (x**3 - 3*x**2 + 2*x );
end f;
declare eps float, (x, y) float (18);
declare dx fixed decimal (15,13);
eps = 1e-12;
do dx = -5.03 to 5 by 0.1;
x = dx;
if sign(f(x)) ^= sign(f(dx+0.1)) then
call locate_root;
end;
locate_root: procedure;
declare (left, mid, right) float (18);
put skip list ('Looking for root in [' || x, x+0.1 || ']' );
left = x; right = dx+0.1;
PUT SKIP LIST (F(LEFT), F(RIGHT) );
if abs(f(left) ) < eps then
do; put skip list ('Found a root at x=', left); return; end;
else if abs(f(right) ) < eps then
do; put skip list ('Found a root at x=', right); return; end;
do forever;
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}}==
{{trans|C++}}
<syntaxhighlight lang="purebasic">Procedure.d f(x.d)
ProcedureReturn x*x*x-3*x*x+2*x
EndProcedure
Procedure main()
OpenConsole()
Define.d StepSize= 0.001
Define.d Start=-1, stop=3
Define.d value=f(start), x=start
Define.i oldsign=Sign(value)
If value=0
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>
=={{header|Python}}==
<
step = 0.001 # Smaller step values produce more accurate and precise results
Line 473 ⟶ 2,635:
sign = value > 0
x += step</
=={{header|R}}==
{{trans|Octave}}
<syntaxhighlight lang="r">f <- function(x) x^3 -3*x^2 + 2*x
findroots <- function(f, begin, end, tol = 1e-20, step = 0.001) {
se <- ifelse(sign(f(begin))==0, 1, sign(f(begin)))
x <- 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
}
}
findroots(f, -1, 3)</syntaxhighlight>
=={{header|Racket}}==
<syntaxhighlight lang="racket">
#lang racket
;; Attempts to find all roots of a real-valued function f
;; in a given interval [a b] by dividing the interval into N parts
;; 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
;; in a given interval [a b].
(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
;; with absolute accuracy (tolerance).
(define (root? f x) (almost-equal? 0 (f x)))
;; Returns #t if interval (a b) contains a root
;; (or the odd number of roots) of a real-valued function f.
(define (includes-root? f a b) (< (* (f a) (f b)) 0))
;; Returns #t if a and b are equal with respect to
;; the relative accuracy (tolerance).
(define (almost-equal? a b)
(or (< (abs (+ b a)) (tolerance))
(< (abs (/ (- b a) (+ b a))) (tolerance))))
(define tolerance (make-parameter 5e-16))
</syntaxhighlight>
Different root-finding methods
<syntaxhighlight lang="racket">
(define (secant f a b)
(let next ([x1 a] [y1 (f a)] [x2 b] [y2 (f b)] [n 50])
(define x3 (/ (- (* x1 y2) (* x2 y1)) (- y2 y1)))
(cond
; if the method din't converge within given interval
; switch to more robust bisection method
[(or (not (< a x3 b)) (zero? n)) (bisection f a b)]
[(almost-equal? x3 x2) x3]
[else (next x2 y2 x3 (f x3) (sub1 n))])))
(define (bisection f x1 x2)
(let divide ([a x1] [b x2])
(and (<= (* (f a) (f b)) 0)
(let ([c (* 0.5 (+ a b))])
(if (almost-equal? a b)
c
(or (divide a c) (divide c b)))))))
</syntaxhighlight>
Examples:
<syntaxhighlight lang="racket">
-> (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)))
-> (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.
The functional nature of Racket allows to perform the optimization without changing the main code using memoization.
Simple memoization operator
<syntaxhighlight lang="racket">
(define (memoized f)
(define tbl (make-hash))
(λ 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
<syntaxhighlight lang="racket">
-> (find-roots (memoized f) -3 4 #:divisions 50)
'(2.4932181969624796e-33 1.0 2.0)
</syntaxhighlight>
The profiling shows that memoization reduces the number of function calls
in this example from 184 to 67 (50 calls for primary interval division and about 6 calls for each point refinement).
=={{header|Raku}}==
(formerly Perl 6)
Uses exact arithmetic.
<syntaxhighlight lang="raku" line>sub f(\x) { x³ - 3*x² + 2*x }
my $start = -1;
my $stop = 3;
my $step = 0.001;
for $start, * + $step ... $stop -> $x {
state $sign = 0;
given f($x) {
my $next = .sign;
when 0.0 {
say "Root found at $x";
}
when $sign and $next != $sign {
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}}==
Both of these REXX versions use the '''bisection method'''.
===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= 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 coded in-line===
This version is about '''40%''' faster than the 1<sup>st</sup> REXX version.
<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 /* " " " " " " */
x= bot - inc /*compute 1st value to start compares. */
z= x * (x * (x-3) + 2) /*formula used ──► x^3 - 3x^2 + 2x */
!= sign(z) /*obtain the sign of the initial value.*/
do x=bot to top by inc /*traipse through the specified range. */
z= x * (x * (x-3) + 2); $= sign(z) /*compute new value; obtain the sign. */
if z=0 then say 'found an exact root at' x/1
else if !\==$ then if !\==0 then say 'passed a root at' x/1
!= $ /*use the new sign for the next compare*/
end /*x*/ /*dividing by unity normalizes X [↑] */</syntaxhighlight>
{{out|output|text= is the same as the 1<sup>st</sup> REXX version.}} <br><br>
=={{header|Ring}}==
<syntaxhighlight lang="ring">
load "stdlib.ring"
function = "return pow(x,3)-3*pow(x,2)+2*x"
rangemin = -1
rangemax = 3
stepsize = 0.001
accuracy = 0.1
roots(function, rangemin, rangemax, stepsize, accuracy)
func roots funct, min, max, inc, eps
oldsign = 0
for x = min to max step inc
num = sign(eval(funct))
if num = 0
see "root found at x = " + x + nl
num = -oldsign
else if num != oldsign and oldsign != 0
if inc < eps
see "root found near x = " + x + nl
else roots(funct, x-inc, x+inc/8, inc/8, eps) ok ok ok
oldsign = num
next
</syntaxhighlight>
Output:
<pre>
root found near x = 0.00
root found near x = 1.00
root found near x = 2.00
</pre>
=={{header|RLaB}}==
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.
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''
using the bisection method on the interval -5 to 5 is,
<syntaxhighlight lang="rlab">
f = function(x)
{
rval = x .^ 3 - 3 * x .^ 2 + 2 * x;
return rval;
};
>> findroot(f, , [-5,5])
0
</syntaxhighlight>
For a detailed description of the solver and its parameters interested reader is directed to the ''rlabplus'' manual.
=={{header|Ruby}}==
{{trans|Python}}
<syntaxhighlight lang="ruby">def sign(x)
x <=> 0
end
def find_roots(f, range, step=0.001)
sign = sign(f[range.begin])
range.step(step) do |x|
value = f[x]
if value == 0
puts "Root found at #{x}"
elsif sign(value) == -sign
puts "Root found between #{x-step} and #{x}"
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
Root found at 1.0
Root found at 2.0
</pre>
Or we could use Enumerable#inject, monkey patching and block:
<syntaxhighlight lang="ruby">class Numeric
def sign
self <=> 0
end
end
def find_roots(range, step = 1e-3)
range.step( step ).inject( yield(range.begin).sign ) do |sign, x|
value = yield(x)
if value == 0
puts "Root found at #{x}"
elsif value.sign == -sign
puts "Root found between #{x-step} and #{x}"
end
value.sign
end
end
find_roots(-1..3) { |x| x**3 - 3*x**2 + 2*x }</syntaxhighlight>
=={{header|Rust}}==
<syntaxhighlight lang="rust">// 202100315 Rust programming solution
use roots::find_roots_cubic;
fn main() {
let roots = find_roots_cubic(1f32, -3f32, 2f32, 0f32);
println!("Result : {:?}", roots);
}</syntaxhighlight>
{{out}}
<pre>
Result : Three([0.000000059604645, 0.99999994, 2.0])
</pre>
Another without external crates:
<syntaxhighlight lang="rust">
use num::Float;
/// Note: We cannot use `range_step` here because Floats don't implement
/// 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() {
let roots = find_roots(
|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);
}
</syntaxhighlight>
{{out}}
<pre>
roots of f(x) = x^3 - 3x^2 + 2x are: [-0.00000000000009381755897326649, 0.9999999999998124, 1.9999999999997022]
</pre>
=={{header|Scala}}==
===Imperative version (Ugly, side effects)===
{{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,
lowerBound: Double,
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) {
val y = f(x)
val s = math.signum(y)
if (s == 0) println(x)
else if (s != os) println(s"~${x - (x - ox) * (y / (y - oy))}")
ox = x
oy = y
os = s
}
}
printRoots(poly, -1.0, 4, 0.002)
}</syntaxhighlight>
===Functional version (Recommended)===
<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
println(findRoots(fn, -1.0, 3.0, 0.0001, 0.000000001))
}</syntaxhighlight>
{{out}}
Vector(-9.381755897326649E-14, 0.9999999999998124, 1.9999999999997022)
=={{header|Scheme}}==
For R7RS Scheme.
<syntaxhighlight lang="scheme">
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;
;;; Finding (real) roots of a function.
;;;
;;; I follow the model that breaks the task into two distinct ones:
;;; isolating real roots, and then finding the isolated roots. The
;;; former task I will call "isolating roots", and the latter I will
;;; call "rootfinding".
;;;
;;; Isolating real roots of a polynomial can be done exactly, and the
;;; methods can handle infinite domains. Scheme (because it has exact
;;; rationals) is a relatively easy language in which to write such
;;; code.
;;;
;;; I have also isolated the roots of low-degree polynomials on finite
;;; intervals by the following method, in floating-point arithmetic:
;;; rewrite the polynomial in the Bernstein polynomials basis, then
;;; take derivatives to get critical points, working your way back up
;;; in degree from a straight line. This method goes back and forth
;;; between the two subtasks. (You could use the quadratic formula
;;; once you got down to degree two, but I wouldn’t bother. The cubic
;;; and quartic formulas are numerically very poor and should be
;;; avoided.)
;;;
;;; However, these methods require that the function be a
;;; polynomial. Here I will simply use a step size and the
;;; intermediate value theorem.
;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(cond-expand
(r7rs)
(chicken (import r7rs)))
;;;
;;; Step 1. Isolation of roots.
;;;
;;; I will simply step over the domain and look for intervals that
;;; contain at least one root. There is a risk of getting the error
;;; message "the root is not bracketed" when you run the
;;; rootfinder. (One could have the rootfinder raise a recoverable
;;; exception or return a special value, instead.)
;;;
(define-library (isolate-roots)
(export isolate-roots)
(import (scheme base))
(begin
(define (isolate-roots f x-min x-max x-step)
(define (ith-x i) (+ x-min (* i x-step)))
(let ((x0 x-min)
(y0 (f x-min)))
(let loop ((i1 1)
(x0 x0)
(y0 y0)
(accum (if (zero? y0)
`((,x0 . ,x0))
'())))
(let ((x1 (ith-x i1)))
(if (< x-max x1)
(reverse accum)
(let ((y1 (f x1)))
(cond
((zero? y1)
(loop (+ i1 1) x1 y1 `((,x1 . ,x1) . ,accum)))
((negative? (* y0 y1))
(loop (+ i1 1) x1 y1 `((,x0 . ,x1) . ,accum)))
(else
(loop (+ i1 1) x1 y1 accum)))))))))
)) ;; end library roots-isolator
;;;
;;; Step 2. Rootfinding.
;;;
;;; I will use the ITP method. I wrote this implementation shortly
;;; after the algorithm was published. See
;;; https://sourceforge.net/p/chemoelectric/itp-root-finder
;;;
;;; Reference:
;;;
;;; I.F.D. Oliveira and R.H.C. Takahashi. 2020. An Enhancement of
;;; the Bisection Method Average Performance Preserving Minmax
;;; Optimality. ACM Trans. Math. Softw. 47, 1, Article 5 (December
;;; 2020), 24 pages. https://doi.org/10.1145/3423597
;;;
(define-library (itp-root-finder)
(export itp-root-finder-epsilon
itp-root-finder-extra-steps
itp-root-finder-kappa1
itp-root-finder-kappa2
;; itp-root-bracket-finder returns two values that form a
;; bracket no wider than 2ϵ.
itp-root-bracket-finder
;; itp-root-finder returns the point midway between the ends
;; of the final bracket.
itp-root-finder)
(import (scheme base))
(import (scheme inexact))
(import (scheme case-lambda))
(import (srfi 143)) ; Fixnums.
(import (srfi 144)) ; Flonums.
(begin
(define ϕ
;; The Golden Ratio, (1 + √5)/2, rounded down by about
;; 0.00003398875.
1.618)
(define 1+ϕ (+ 1.0 ϕ))
(define itp-root-finder-epsilon
(make-parameter
(* 1000.0 fl-epsilon)
(lambda (ϵ)
(if (positive? ϵ)
ϵ
(error 'itp-root-finder-epsilon
"a positive value was expected"
ϵ)))))
(define itp-root-finder-extra-steps
;; Increase extra-steps above zero, if you wish to try to speed
;; up convergence, at the expense of that many more steps in the
;; worst case.
(make-parameter
0
(lambda (n₀)
(if (or (negative? n₀) (not (integer? n₀)))
(error 'itp-root-finder-extra-steps
"a non-negative integer was expected"
n₀)
n₀))))
(define itp-root-finder-kappa1
(make-parameter
0.1
(lambda (κ₁)
(if (positive? κ₁)
κ₁
(error 'itp-root-finder-kappa1
"a positive value was expected"
κ₁)))))
(define itp-root-finder-kappa2
(make-parameter
2.0
(lambda (κ₂)
(if (or (< κ₂ 1) (< 1+ϕ κ₂))
;; We allow <= 1+ϕ (instead of ‘< 1+ϕ’) because we
;; already rounded ϕ down.
(error 'itp-root-finder-kappa2
(string-append "a value 1 <= kappa2 <= "
(number->string 1+ϕ)
" was expected")
κ₂)
κ₂))))
(define (sign x)
(cond ((negative? x) -1)
((positive? x) 1)
(else 0)))
(define (apply-sign σ x)
(cond ((fxnegative? σ) (- x))
((fxpositive? σ) x)
(else 0)))
(define (itp-root-bracket-finder%% f a b ϵ n₀ κ₁ κ₂)
(let* ((2ϵ (inexact (* 2 ϵ)))
(n½ (exact (ceiling (log (/ (inexact (- b a)) 2ϵ) 2))))
(n_max (+ n½ n₀))
(ya (f a))
(yb (f b))
(σ_ya (sign ya))
(σ_yb (sign yb)))
(cond
((fxzero? σ_ya) (values a a))
((fxzero? σ_yb) (values b b))
(else
(when (fxpositive? (* σ_ya σ_yb))
(error 'itp-root-bracket-finder
"the root is not bracketed"
a b))
(let loop ((pow2 (expt 2 n_max))
(a (inexact a))
(b (inexact b))
(ya (inexact ya))
(yb (inexact yb)))
(if (or (= pow2 1) (fl<=? (fl- b a) 2ϵ))
(values a b)
(let* ( ;; x½ – the bisection.
(x½ (fl* 0.5 (fl+ a b)))
;; xf – interpolation by regula falsi.
(xf (fl/ (fl- (fl* yb a) (fl* ya b))
(fl- yb ya)))
(b-a (fl- b a))
(δ (fl* κ₁ (flabs (expt b-a κ₂))))
(x½-xf (fl- x½ xf))
(σ (sign x½-xf))
;; xt – the ‘truncation’ of xf.
(xt (if (fl<=? δ (flabs x½-xf))
(fl+ xf (apply-sign σ δ))
x½))
(r (- (* pow2 ϵ) (fl* 0.5 b-a)))
;; xp – the projection of xt onto [x½-r,x½+r].
(xp (if (fl<=? (flabs (fl- xt x½)) r)
xt
(fl- x½ (apply-sign σ r))))
(yp (inexact (f xp))))
(let ((pow2/2 (truncate-quotient pow2 2))
(σ_yp (sign yp)))
(cond ((fx=? σ_ya σ_yp)
;; yp has the same sign as ya. Make it the
;; new ya.
(loop pow2/2 xp b yp yb))
((fx=? σ_yb σ_yp)
;; yp has the same sign as yb. Make it the
;; new yb.
(loop pow2/2 a xp ya yp))
(else
;; yp is zero.
(values xp xp)))))))))))
(define (itp-root-bracket-finder% f a b ϵ n₀ κ₁ κ₂)
(cond
((< b a) (itp-root-bracket-finder% b a f ϵ n₀ κ₁ κ₂))
(else
(let* ((ϵ (or ϵ (itp-root-finder-epsilon)))
(n₀ (or n₀ (itp-root-finder-extra-steps)))
(κ₁ (or κ₁ (itp-root-finder-kappa1)))
(κ₂ (or κ₂ (itp-root-finder-kappa2))))
(when (negative? ϵ)
(error 'itp-root-bracket-finder
"a positive value was expected" ϵ))
(when (negative? κ₁)
(error 'itp-root-bracket-finder
"a positive value was expected" κ₁))
(when (or (< κ₂ 1) (< 1+ϕ κ₂))
;; We allow <= 1+ϕ (instead of ‘< 1+ϕ’) because we already
;; rounded ϕ down.
(error 'itp-root-bracket-finder
(string-append "a value 1 <= kappa2 <= "
(number->string 1+ϕ)
" was expected")
κ₂))
(when (or (negative? n₀) (not (integer? n₀)))
(error 'itp-root-bracket-finder
"a non-negative integer was expected" n₀))
(itp-root-bracket-finder%% f a b ϵ n₀ κ₁ κ₂)))))
(define itp-root-bracket-finder
(case-lambda
((f a b)
(itp-root-bracket-finder% f a b #f #f #f #f))
((f a b ϵ)
(itp-root-bracket-finder% f a b ϵ #f #f #f))
((f a b ϵ n₀)
(itp-root-bracket-finder% f a b ϵ n₀ #f #f))
((f a b ϵ n₀ κ₁)
(itp-root-bracket-finder% f a b ϵ n₀ κ₁ #f))
((f a b ϵ n₀ κ₁ κ₂)
(itp-root-bracket-finder% f a b ϵ n₀ κ₁ κ₂))))
(define (itp-root-finder% f a b ϵ n₀ κ₁ κ₂)
(call-with-values
(lambda ()
(itp-root-bracket-finder f a b ϵ n₀ κ₁ κ₂))
(lambda (a b)
(if (= a b)
a
(* 0.5 (+ a b))))))
(define itp-root-finder
(case-lambda
((f a b)
(itp-root-finder% f a b #f #f #f #f))
((f a b ϵ)
(itp-root-finder% f a b ϵ #f #f #f))
((f a b ϵ n₀)
(itp-root-finder% f a b ϵ n₀ #f #f))
((f a b ϵ n₀ κ₁)
(itp-root-finder% f a b ϵ n₀ κ₁ #f))
((f a b ϵ n₀ κ₁ κ₂)
(itp-root-finder% f a b ϵ n₀ κ₁ κ₂))))
)) ;; end library itp-root-finder
(import (scheme base))
(import (scheme write))
(import (isolate-roots))
(import (itp-root-finder))
(define (f x)
;; x³ - 3x² + 2x, written in Horner form.
(* x (+ 2 (* x (+ -3 x)))))
(define (find-root f interval)
(define (display-exactness root)
(display (if (and (exact? root)
(exact? (f root))
(zero? (f root)))
" (exact) "
" (inexact) ")))
(let ((x0 (car interval))
(x1 (cdr interval)))
(if (= x0 x1)
(begin
(let ((root (if (exact? x0) x0 x1)))
(display-exactness root)
(display "(rootfinder not used) ")
(display root)
(newline)))
(begin
;;
;; I am not careful here to avoid accidentally excluding the
;; root from the bracketing interval [x0,x1]. Floating point
;; is very tricky to work with.
;;
(let ((root (itp-root-finder f x0 x1)))
(display-exactness root)
(display "(rootfinder used) ")
(display root)
(newline))))))
;;; The following two demonstrations find all three roots exactly, as
;;; exact rationals, without the need for a rootfinding step.
(newline)
(display "Stepping by 1/1000 from 0 to 2:")
(newline)
(do ((p (isolate-roots f 0 2 1/1000) (cdr p)))
((not (pair? p)))
(find-root f (car p)))
(newline)
(display "Stepping by 1/1000 from -10 to 10:")
(newline)
(do ((p (isolate-roots f -10 10 1/1000) (cdr p)))
((not (pair? p)))
(find-root f (car p)))
;;; The following demonstration gives inexact results, because the
;;; step size is an inexact number.
(newline)
(display "Stepping by 0.001 from -10.0 to 10.0:")
(newline)
(do ((p (isolate-roots f -10.0 10.0 0.001) (cdr p)))
((not (pair? p)))
(find-root f (car p)))
;;; The following demonstration gives inexact results, because the
;;; rootfinder is needed.
(newline)
(display "Stepping by 13/3333 from -2111/1011 to 33/13:")
(newline)
(do ((p (isolate-roots f -2111/1011 33/13 13/3333) (cdr p)))
((not (pair? p)))
(find-root f (car p)))
(newline)
</syntaxhighlight>
{{out}}
<pre>$ gosh roots-of-a-function.scm
Stepping by 1/1000 from 0 to 2:
(exact) (rootfinder not used) 0
(exact) (rootfinder not used) 1
(exact) (rootfinder not used) 2
Stepping by 1/1000 from -10 to 10:
(exact) (rootfinder not used) 0
(exact) (rootfinder not used) 1
(exact) (rootfinder not used) 2
Stepping by 0.001 from -10.0 to 10.0:
(inexact) (rootfinder not used) 0.0
(inexact) (rootfinder not used) 1.0
(inexact) (rootfinder not used) 2.0
Stepping by 13/3333 from -2111/1011 to 33/13:
(inexact) (rootfinder used) -1.3380129580295458e-15
(inexact) (rootfinder used) 0.9999999999998657
(inexact) (rootfinder used) 1.9999999999999998
</pre>
=={{header|Sidef}}==
<syntaxhighlight lang="ruby">func f(x) {
x*x*x - 3*x*x + 2*x
}
var step = 0.001
var start = -1
var stop = 3
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
Root found at 2</pre>
=={{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).
<syntaxhighlight lang="tcl">proc froots {lambda {start -3} {end 3} {step 0.0001}} {
set res {}
set lastsign [sgn [apply $lambda $start]]
for {set x $start} {$x <= $end} {set x [expr {$x + $step}]} {
set sign [sgn [apply $lambda $x]]
if {$sign != $lastsign} {
lappend res [format ~%.11f $x]
}
set lastsign $sign
}
return $res
}
proc sgn x {expr {($x>0) - ($x<0)}}
puts [froots {x {expr {$x**3 - 3*$x**2 + 2*$x}}}]</syntaxhighlight>
Result and timing:
<pre>/Tcl $ time ./froots.tcl
~0.00000000000 ~1.00000000000 ~2.00000000000
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 \
{x {expr {$x**3 - 3*$x**2 + 2*$x}}} \
{x {expr {3*$x**2 - 6*$x + 2}}}]</syntaxhighlight>
=={{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>.
In this case, the roots are exact; inexact results are marked by decimal points.
=={{header|Wren}}==
{{trans|Go}}
{{libheader|Wren-fmt}}
<syntaxhighlight lang="wren">import "./fmt" for Fmt
var secant = Fn.new { |f, x0, x1|
var f0 = 0
var f1 = f.call(x0)
for (i in 0...100) {
f0 = 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)
}
return [0, ""]
}
var findRoots = Fn.new { |f, lower, upper, step|
var x0 = lower
var x1 = lower + step
while (x0 < upper) {
x1 = (x1 < upper) ? x1 : upper
var res = secant.call(f, x0, x1)
var r = res[0]
var status = res[1]
if (status != "" && r >= x0 && r < x1) {
Fmt.print(" $6.3f $s", r, status)
}
x0 = x1
x1 = x1 + step
}
}
var example = Fn.new { |x| x*x*x - 3*x*x + 2*x }
findRoots.call(example, -0.5, 2.6, 1)</syntaxhighlight>
{{out}}
<pre>
0.000 approximate
1.000 exact
2.000 approximate
</pre>
=={{header|XPL0}}==
{{trans|Wren}}
<syntaxhighlight lang "XPL0">include xpllib; \for Print
func real F(X);
real X;
return X*X*X - 3.*X*X + 2.*X;
char Status;
func real Secant(X0, X1);
real X0, X1, F0, F1, T;
int I;
[F1:= F(X0);
for I:= 0 to 100-1 do
[F0:= F1;
F1:= F(X1);
if F1 = 0. then [Status:= "exact"; return X1];
if abs(X1-X0) < 1e-6 then [Status:= "approximate"; return X1];
T:= X0;
X0:= X1;
X1:= X1 - F1*(X1-T)/(F1-F0);
];
Status:= 0; return 0.;
];
func FindRoots(Lower, Upper, Step);
real Lower, Upper, Step;
real X0, X1, R;
[X0:= Lower;
X1:= Lower + Step;
while X0 < Upper do
[X1:= if X1 < Upper then X1 else Upper;
R:= Secant(X0, X1);
if Status # 0 and R >= X0 and R < X1 then
Print(" %2.3f %s\n", R, Status);
X0:= X1;
X1:= X1 + Step;
];
];
FindRoots(-0.5, 2.6, 1.)</syntaxhighlight>
{{out}}
<pre>
0.000 approximate
1.000 exact
2.000 approximate
</pre>
=={{header|zkl}}==
{{trans|Haskell}}
<syntaxhighlight lang="zkl">fcn findRoots(f,start,stop,step,eps){
[start..stop,step].filter('wrap(x){ f(x).closeTo(0.0,eps) })
}</syntaxhighlight>
<syntaxhighlight lang="zkl">fcn f(x){ x*x*x - 3.0*x*x + 2.0*x }
findRoots(f, -1.0, 3.0, 0.0001, 0.00000001).println();</syntaxhighlight>
{{out}}
<pre>L(-9.38176e-14,1,2)</pre>
{{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 "Function is not converging near (%7.4f,%7.4f).".fmt(xA,xB);
}</syntaxhighlight>
<syntaxhighlight lang="zkl">step:=0.1;
xs:=findRoots(f, -1.032, 3.0, step, 0.1);
xs.println(" --> ",xs.apply('wrap(x){ secant(f,x-step,x+step) }));</syntaxhighlight>
{{out}}
<pre>L(-0.032,0.968,1.068,1.968) --> L(1.87115e-19,1,1,2)</pre>
{{omit from|M4}}
| |||