Roots of a function

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 example, use f(x)=x^3-3x^2+2x.

Ada

with Ada.Text_Io; use Ada.Text_Io;

procedure Roots_Of_Function is
   package Real_Io is new Ada.Text_Io.Float_Io(Long_Float);
   use Real_Io;
   
   function F(X : Long_Float) return Long_Float is
   begin
      return (X**3 - 3.0*X*X + 2.0*X);
   end F;
   
   Step  : constant Long_Float := 1.0E-6;
   Start : constant Long_Float := -1.0;
   Stop  : constant Long_Float := 3.0;
   Value : Long_Float := F(Start);
   Sign  : Boolean := Value > 0.0;
   X     : Long_Float := Start + Step;
 
begin
   if Value = 0.0 then
      Put("Root found at ");
      Put(Item => Start, Fore => 1, Aft => 6, Exp => 0);
      New_Line;
   end if;
   while X <= Stop loop
      Value := F(X);
      if (Value > 0.0) /= Sign then
         Put("Root found near ");
         Put(Item => X, Fore => 1, Aft => 6, Exp => 0);
         New_Line;
      elsif Value = 0.0 then
         Put("Root found at ");
         Put(Item => X, Fore => 1, Aft => 6, Exp => 0);
         New_Line;
      end if;
      Sign := Value > 0.0;
      X := X + Step;
   end loop;
end Roots_Of_Function;

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.") ;
}

bool 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

Maple

 f := x^3-3*x^2+2*x;
 roots(f,x);

outputs:

 [[0, 1], [1, 1], [2, 1]]

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).

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.

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 );
}