function roots ()
% test program for root finding routine.
  fprintf('  %s\n', 'Bisection method');
  fprintf('  %8s  %16s  %16s  %6s\n', 'function', 'xend', 'error', 'nmbf');
  a0 = 1; b0 = 2;
  [xend, nmbf] = findroot_bisect ('func0', a0, b0);
  fprintf('  %8s  %16.7e  %16.7e  %6d\n', 'func0', xend, xend-sqrt(2), nmbf);
  a0 = 1; b0 = 2;
  [xend, nmbf] = findroot_bisect ('func1', a0, b0);
  fprintf('  %8s  %16.7e  %16.7e  %6d\n', 'func1', xend, xend-sqrt(2), nmbf);
  a0 = 1; b0 = 2;
  [xend, nmbf] = findroot_bisect ('func2', a0, b0);
  fprintf('  %8s  %16.7e  %16.7e  %6d\n', 'func2', xend, xend-sqrt(2), nmbf);
  a0 = 1; b0 = 2;
  [xend, nmbf] = findroot_bisect ('func3', a0, b0);
  fprintf('  %8s  %16.7e  %16.7e  %6d\n', 'func3', xend, xend-sqrt(2), nmbf);
  a0 = 1; b0 = 2;
  [xend, nmbf] = findroot_bisect ('func4', a0, b0);
  fprintf('  %8s  %16.7e  %16.7e  %6d\n', 'func4', xend, xend-sqrt(2), nmbf);
  a0 = 1; b0 = 2;
  [xend, nmbf] = findroot_bisect ('func5', a0, b0);
  fprintf('  %8s  %16.7e  %16.7e  %6d\n', 'func5', xend, xend-sqrt(2), nmbf);
  a0 = -1; b0 = 2;
  [xend, nmbf] = findroot_bisect ('func6', a0, b0);
  fprintf('  %8s  %16.7e  %16.7e  %6d\n', 'func6', xend, xend-0, nmbf);
  fprintf('  %s\n', 'Fractional linear method');
  fprintf('  %8s  %16s  %16s  %6s\n', 'function', 'xend', 'error', 'nmbf');
  a0 = 1; b0 = 2;
  [xend, nmbf] = findroot ('func0', a0, b0);
  fprintf('  %8s  %16.7e  %16.7e  %6d\n', 'func0', xend, xend-sqrt(2), nmbf);
  a0 = 1; b0 = 2;
  [xend, nmbf] = findroot ('func1', a0, b0);
  fprintf('  %8s  %16.7e  %16.7e  %6d\n', 'func1', xend, xend-sqrt(2), nmbf);
  a0 = 1; b0 = 2;
  [xend, nmbf] = findroot ('func2', a0, b0);
  fprintf('  %8s  %16.7e  %16.7e  %6d\n', 'func2', xend, xend-sqrt(2), nmbf);
  a0 = 1; b0 = 2;
  [xend, nmbf] = findroot ('func3', a0, b0);
  fprintf('  %8s  %16.7e  %16.7e  %6d\n', 'func3', xend, xend-sqrt(2), nmbf);
  a0 = 1; b0 = 2;
  [xend, nmbf] = findroot ('func4', a0, b0);
  fprintf('  %8s  %16.7e  %16.7e  %6d\n', 'func4', xend, xend-sqrt(2), nmbf);
  a0 = 1; b0 = 2;
  [xend, nmbf] = findroot ('func5', a0, b0);
  fprintf('  %8s  %16.7e  %16.7e  %6d\n', 'func5', xend, xend-sqrt(2), nmbf);
  a0 = -1; b0 = 2;
  [xend, nmbf] = findroot ('func6', a0, b0);
  fprintf('  %8s  %16.7e  %16.7e  %6d\n', 'func6', xend, xend-0, nmbf);

  function [xend, nmbf] = findroot_bisect (fname, a0, b0)
% rootfinding with use of bisection
  a = a0;
  b = b0;
  fa = feval(fname,a);
  fb = feval(fname,b);
  nmbf = 2;
  if fa==0
   c = a; fc = fa;
  elseif fb==0
   c = b; fc = fb;
  else
   if signum(fa)==signum(fb)
    error ('error in findroot -- no bracket for root')
   end
   c = bisect(a,b);
   fc = feval(fname,c);
   nmbf = nmbf+1;
  end
  while (c~=a&c~=b)
   if fc==0
    a = c; fa = fc;
    b = c; fb = fc;
   elseif signum(fc)==signum(fa)
    a = c; fa = fc;
    c = bisect(a,b);
    fc = feval(fname,c);
    nmbf = nmbf+1;
   elseif signum(fc)==signum(fb)
    b = c; fb = fc;
    c = bisect(a,b);
    fc = feval(fname,c);
    nmbf = nmbf+1;
   else
    error ('error in findroot -- f generates NaN?')
   end
  end
  xend = c;

  function [xend, nmbf] = findroot (fname, a0, b0)
% rootfinding with use of fractional linear interpolation
  a = a0;
  b = b0;
  fa = feval(fname,a);
  fb = feval(fname,b);
  nmbf = 2;
  if (fa==0)
   c = a; fc = fa;
  elseif (fb==0)
   c = b; fc = fb;
  else
   if signum(fa)==signum(fb)
    error ('error in findroot -- no bracket for root')
   end
   c = bisect(a,b);
   fc = feval(fname,c);
   nmbf = nmbf+1;
  end
  while (c~=a&c~=b)
   if (fc==0)
    a = c; fa = fc;
    b = c; fb = fc;
   elseif (signum(fc)==signum(fa))
    t0 = fraclin(a,fa,b,fb,c,fc);
    if (abs(b-c)<2*abs(b-a)/3&min(c,b)<t0&t0<max(c,b))
     xnew = t0;
    else
     xnew = bisect(c,b);
    end
    a = c; fa = fc;
    c = xnew;
    fc = feval(fname,c);
    nmbf = nmbf+1;
   elseif (signum(fc)==signum(fb))
    t0 = fraclin(a,fa,b,fb,c,fc);
    if (abs(c-a)<2*abs(b-a)/3&min(a,c)<t0&t0<max(a,c))
     xnew = t0;
    else
     xnew = bisect(a,c);
    end
    b = c; fb = fc;
    c = xnew;
    fc = feval(fname,c);
    nmbf = nmbf+1;
   else
    error ('error in findroot -- f generates NaN?')
   end
  end
  xend = c;

  function c = bisect(a,b)
% Returns a bisection point for the interval (a,b) or (b,a).
% If a and b are both non-zero, of the same sign, and not too far
% apart on a logarithmic scale, returns c=(a+b)/2.  Other cases ...
% see the code.  The intent is that bisect(a,b) divides the interval
% so that there are about equal number of floating point values in
% either half.  Some preference towards zero.
  eps = 2.0^(-52); % IEEE machine precision
  tiny = 2.0^(-1022); % IEEE smallest positive normal number
  if (signum(a)*signum(b)<0)
   c = 0;
  elseif (signum(a)*signum(b)>0)
   t0 = min(abs(a),abs(b));
   t1 = max(abs(a),abs(b));
   if (t1<4*t0|t1<4*tiny)
    c = (a+b)/2;
   else
    c = signum(a)*sqrt(t0)*sqrt(t1);
   end
  elseif (a==0|b==0)
   t0 = a+b;
   if (abs(t0)>=1/sqrt(eps))
    c = signum(t0);
   elseif (abs(t0)>=sqrt(eps))
    c = t0*sqrt(eps);
   elseif (abs(t0)>=sqrt(tiny))
    c = t0*abs(t0);
   elseif (abs(t0)>=4*tiny)
    c = signum(t0)*tiny;
   elseif (abs(t0)>eps*tiny)
    c = signum(t0)*eps*tiny;
   else
    c = 0;
   end
  else
%  must be NaN
   c = (a+b)/2;
  end

  function c = fraclin(x0,y0,x1,y1,x2,y2)
% fractional linear interpolation
  dx0 = x0-x2;
  dx1 = x1-x2;
  xnew0 = x2+dx0*dx1*y2*(y1-y0)/(dx1*y1*(y2-y0)-dx0*y0*(y2-y1));
  xnew1 = x2+dx0*dx1*(y1-y0)/(dx1*(y2-y0)-dx0*(y2-y1));
  if (min(min(x0,x1),x2)<=xnew0&xnew0<=max(max(x0,x1),x2))
   c = xnew0;
  elseif (min(min(x0,x1),x2)<=xnew1&xnew1<=max(max(x0,x1),x2))
   c = xnew1;
  else
   % this can only occur because of roundoff error.  Use the median.
   c = min(min(max(x0,x1),max(x1,x2)),max(x2,x0));
  end

  function s = signum(x)
% returns the sign of x; returns 0 if x==0.
  if x<0
    s = -1;
  elseif x==0
    s = 0;
  elseif x>0
    s = 1;
  else
    s = 0/0;
  end

  function y = func0(x)
% test function for rootfinding routine.
% func0 is smooth and has a simple root at x=sqrt(2).
  y = (x-sqrt(2))*sqrt(1+x^2);

  function y = func1(x)
% test function for rootfinding routine.
% func1 is smooth and has a triple root at x=sqrt(2).
  y = (x-sqrt(2))^3;

  function y = func2(x)
% test function for rootfinding routine.
% func2 is smooth and has a simple root at x=sqrt(2); the
% asymptotic behavior may pose problems for a Newton-like method.
  y = atan(x-sqrt(2));

  function y = func3(x)
% test function for rootfinding routine.
% func3 is smooth and has a simple root at x=sqrt(2); the
% asymptotic behavior may pose problems for a Newton-like method.
  y = atan(256*(x-sqrt(2)));

  function y = func4(x)
% test function for rootfinding routine.
% func4 is discontinuous with sign change at x=sqrt(2).
  if x<0|x*x<2
    y = -1;
  else
    y = 1;
  end

  function y = func5(x)
% test function for rootfinding routine.
% func5 has a vertical asymptote and sign change at x=sqrt(2).
  t = (x-sqrt(2))*sqrt(1+x^2);
  if t~=0
    y = 1/t;
  else
    y = 0;
  end

  function y = func6(x)
% test function for rootfinding routine.
% func6 is discontinous with sign change at 0.
  if x<0
    y = -1;
  else
    y = 1;
  end