function [r] = biroot (x)
% Compute the square root of x >= 1 by the bisection algorithm.
a=1;  % left endpoint
b=x;  % right endpoint
while b-a > eps(a)
   c=(b+a)/2;  % midpoint
   if (c*c > x)
      b = c;
   else
      a = c;
   end
end
r=b;  % could use a, b, or c here - all within eps(a) of correct
end