Return the number of bisections from the bisect() function.

[octave.git] / bisect.m

diff --git a/bisect.m b/bisect.m
old mode 100755 (executable)
new mode 100644 (file)
index 56e099e..3570e57
--- a/bisect.m
+++ b/bisect.m
@@ -1,6 +1,4 @@
-#!/usr/bin/octave --silent
-
-function root = bisect(f, a, b, epsilon)
+function [root,iterations] = bisect(f, a, b, epsilon)
   ## Find a root of the function `f` on the closed interval [a,b].
   ##
   ## It is assumed that there are an odd number of roots within [a,b].
@@ -11,6 +9,26 @@ function root = bisect(f, a, b, epsilon)
   ##
   ## If `f` has more than one root on [a,b], the smallest root will be
   ## returned. This is an implementation detail.
+  ##
+  ##
+  ## INPUTS:
+  ##
+  ##   * ``f`` - The function whose root we seek.
+  ##
+  ##   * ``a`` - The "left" endpoint of the interval in which we search.
+  ##
+  ##   * ``b`` - The "right" endpoint of the interval in which we search.
+  ##
+  ##   * ``epsilon`` - How close we should be to the actual root before we
+  ##     halt the search and return the current approximation.
+  ##
+  ## OUTPUTS:
+  ##
+  ##   * ``root`` - The root that we found.
+  ##
+  ##   * ``iterations`` - The number of bisections that we performed
+  ##     during the search.
+  ##
 
   ## Store these so we don't have to recompute them.
   fa = f(a);
@@ -20,16 +38,19 @@ function root = bisect(f, a, b, epsilon)
   ## interval.  
   if (fa == 0)
     root = a;
+    iterations = 0;
     return;
   end
 
   if (fb == 0)
     root = b;
+    iterations = 0;
     return;
   end
 
   ## Bisect the interval.
   c = (a+b)/2;
+  iterations = 1;
   
   ## If we're within the prescribed tolerance, we're done.
   if (b-c < epsilon)
@@ -38,8 +59,10 @@ function root = bisect(f, a, b, epsilon)
   end
 
   if (has_root(fa,f(c)))
-    root = bisect(f,a,c,epsilon);
+    [root,subiterations] = bisect(f,a,c,epsilon);
+    iterations = iterations + subiterations;
   else
-    root = bisect(f,c,b,epsilon);
+    [root,subiterations] = bisect(f,c,b,epsilon);
+    iterations = iterations + subiterations;
   end
 end