Roots of a function: Difference between revisions
Content added Content deleted
(Added Perl version) |
(added D code) |
||
Line 41: | Line 41: | ||
} |
} |
||
</pre> |
</pre> |
||
=={{header|D}}== |
|||
<pre>module findroot ; |
|||
import std.stdio ; |
|||
import std.math ; |
|||
void report(T)(T[] r, T function(T) f, T tolerance = cast(T) 1e-4L) { |
|||
if (r.length) { |
|||
writefln("Root found (tolerance = %1.4g) :", tolerance) ; |
|||
foreach(x ; r) { |
|||
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.") ; |
|||
} |
|||
T nearZero(T)(T a, T b = T.epsilon * 4) { return abs(a) <= b ; } |
|||
T[] findroot(T)(T function(T) f, T start, T end, T step = cast(T) 0.001L, |
|||
T tolerance = cast(T) 1e-4L) { |
|||
T[T] result ; |
|||
if (nearZero(step)) |
|||
writefln("WARNING: step size may be too small.") ; |
|||
T searchRoot(T a, T b) { // search root by simple bisection |
|||
T root ; |
|||
int limit = 49 ; |
|||
T gap = b - a ; |
|||
while (!nearZero(gap) && limit--) { |
|||
if (nearZero(f(a))) return a ; |
|||
if (nearZero(f(b))) return b ; |
|||
root = (b + a)/2.0L ; |
|||
if (nearZero(f(root))) return root ; |
|||
if (f(a) * f(root) < 0) |
|||
b = root ; |
|||
else |
|||
a = root ; |
|||
gap = b - a ; |
|||
} |
|||
return root ; |
|||
} |
|||
T dir = cast(T) (end > start ? 1.0 : -1.0) ; |
|||
step = (end > start) ? abs(step) : - abs(step) ; |
|||
for(T x = start ; x*dir <= end*dir ; x = x + step) |
|||
if (f(x)*f(x + step) <= 0) { |
|||
T r = searchRoot(x, x+ step) ; |
|||
result[r] = f(r) ; |
|||
} |
|||
return result.keys.sort ; // reduce duplacated root, if any |
|||
} |
|||
real f(real x){ |
|||
return x*x*x - 3*x*x + 2*x ; |
|||
} |
|||
void main(){ |
|||
findroot(&f, -1.0L, 3.0L, 0.001L).report(&f) ; |
|||
}</pre> |
|||
Output ( NB:smallest increment for real type in D is real.epsilon = 1.0842e-19 ): |
|||
<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> |
|||
=={{header|Perl}}== |
=={{header|Perl}}== |
||
<pre>sub f |
<pre>sub f |
Revision as of 09:14, 22 February 2008
![Task](http://static.miraheze.org/rosettacodewiki/thumb/b/ba/Rcode-button-task-crushed.png/64px-Rcode-button-task-crushed.png)
You are encouraged to solve this task according to the task description, using any language you may know.
Create a program that finds an 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 example, use f(x)=x^3-3x^2+2x.
C++
#include <iostream> double f(double x) { return (x*x*x - 3*x*x + 2*x); } int main() { double step = 0.001; // Smaller step values produce more accurate and precise results double start = -1; double stop = 3; double value = f(start); double sign = (value > 0); // Check for root at start if ( 0 == value ) std::cout << "Root found at " << start << std::endl; for( double x = start + step; x <= stop; x += step ) { value = f(x); if ( ( value > 0 ) != sign ) // We passed a root std::cout << "Root found near " << x << std::endl; else if ( 0 == value ) // We hit a root std::cout << "Root found at " << x << std::endl; // Update our sign sign = ( value > 0 ); } }
D
module findroot ; import std.stdio ; import std.math ; void report(T)(T[] r, T function(T) f, T tolerance = cast(T) 1e-4L) { if (r.length) { writefln("Root found (tolerance = %1.4g) :", tolerance) ; foreach(x ; r) { 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.") ; } T nearZero(T)(T a, T b = T.epsilon * 4) { return abs(a) <= b ; } T[] findroot(T)(T function(T) f, T start, T end, T step = cast(T) 0.001L, T tolerance = cast(T) 1e-4L) { T[T] result ; if (nearZero(step)) writefln("WARNING: step size may be too small.") ; T searchRoot(T a, T b) { // search root by simple bisection T root ; int limit = 49 ; T gap = b - a ; while (!nearZero(gap) && limit--) { if (nearZero(f(a))) return a ; if (nearZero(f(b))) return b ; root = (b + a)/2.0L ; if (nearZero(f(root))) return root ; if (f(a) * f(root) < 0) b = root ; else a = root ; gap = b - a ; } return root ; } T dir = cast(T) (end > start ? 1.0 : -1.0) ; step = (end > start) ? abs(step) : - abs(step) ; for(T x = start ; x*dir <= end*dir ; x = x + step) if (f(x)*f(x + step) <= 0) { T r = searchRoot(x, x+ step) ; result[r] = f(r) ; } return result.keys.sort ; // reduce duplacated root, if any } real f(real x){ return x*x*x - 3*x*x + 2*x ; } void main(){ findroot(&f, -1.0L, 3.0L, 0.001L).report(&f) ; }
Output ( NB:smallest increment for real type in D is real.epsilon = 1.0842e-19 ):
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
Perl
sub f { my $x = shift; return ($x * $x * $x - 3*$x*$x + 2*$x); } my $step = 0.001; # Smaller step values produce more accurate and precise results my $start = -1; my $stop = 3; my $value = &f($start); my $sign = $value > 0; # Check for root at start print "Root found at $start\n" if ( 0 == $value ); for( my $x = $start + $step; $x <= $stop; $x += $step ) { $value = &f($x); if ( 0 == $value ) { # We hit a root print "Root found at $x\n"; } elsif ( ( $value > 0 ) != $sign ) { # We passed a root print "Root found near $x\n"; } # Update our sign $sign = ( $value > 0 ); }