%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% zerof.m %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function [xnext,work]=zerof(x,f,work);
% reverse communication bisection method for finding a zero 
% *** if no sign change, the routine tries to find one,
% *** but not necessarily in the required range
% *** if no good starting points are avaialble
% *** it may be better to use a minimizer for |f| (not f^2 !)
% *** such as http://solon.cma.univie.ac.at/~neum/software/ls/
% *** to locate first good starting points, and possibly sign changes 
% 
% usage: 
% execute the following to find a zero in [range(1),range(2)]:
%
% x=first_point;work=second_point;	% first two points
%					% must be in range
% out=0;
% nf=0;
% for times=1:nfmax,			% max number of function values
%   if x<range(1) | x>range(2),
%     out=out+1;
%     f=100^out*work(4);		% dummy function value
%   else
%     nf=nf+1;				% number of function values used
%     f=function(x);			% calculate function value
%   end;
%   [x,work]=zerof(x,f,work);		% new guess
%   done=(work(1)<=accuracy); 		% accuracy test
%   if done | out>3, break; end;
% end;
% if ~work(2), 
%   disp('no sign change found; check function value'); 
%   x=work(3);f=work(4);
% end;
% 
% then x is the zero (if work(2)) or the best point found (otherwise)
% IMPORTANT: do not change work outside of zerof!!!
%
% accuracy=0; works without problems (but is slightly slower)
% 
% method: find sign change with safeguarded Muller steps;
% then bisect with safeguarded secant steps
% is slow on triple roots 
% (worst case factor, ~0.25 in 4 iterations, is achieved)
%
% source:     http://solon.cma.univie.ac.at/~neum/software/zeros/zerof.m
% send bug reports to: Arnold Neumaier (neum@cma.univie.ac.at)
% more software at:    http://solon.cma.univie.ac.at/~neum/software.html
%
function [xnext,work]=zerof(x,f,work);

qmin=0.707;

prt=0; 	% for debug only

% check for exact zero
if f==0, xnext=x;work=[0,1]; return; end;

if length(work)==1,
  % initialization
  xnext=work(1);
  work(1)=inf;	% accuracy 
  work(2)=-1;	% sign change indicator
  work(3)=x;	% first point
  work(4)=f;	% first function value
  work(5)=NaN;	% second point
  work(6)=NaN;	% second function value
  work(7)=NaN;	% third point
  work(8)=NaN;	% third function value
  work(9)=1;	% number of calls to zerof.m
  if prt>1, 
    disp('************************************************') 
    disp(' f1         x1               x2               f2') 
  elseif prt==1, 
    disp('***********************************') 
    disp('log10(fbest)  improv    sign change')
  end; 

  return;
end;

sign_change=work(2);
x1=work(3);f1=work(4);
x2=work(5);f2=work(6);
x3=work(7);f3=work(8);	% currently unused
nf=work(9)+1;

% place current point
if sign_change<0, 
  % assign second point
  x2=x;f2=f;
  sign_change=(f1*f2<=0);
  % enforce abs(f1)<=abs(f2)
  if abs(f1)>abs(f2),
    x=x1;f=f1;x1=x2;f1=f2;x2=x;f2=f;
  end;
  q=0;
  % use mirroring step
  xnext = x1 + (x1 - x2);    
  acc=inf;
elseif sign_change==0, 
  % no bracket exists yet 
  sign_change=(f1*f<=0);

  if sign_change,
    disp(['first sign change at nf = ',num2str(nf)]);
  else
    if abs(f)<=abs(f1), 
      % best point into x1
      x3=x;f3=f; x=x1;f=f1; x1=x3;f1=f3;
    end;    
    if abs(x-x1)>abs(x-x2); 
      % farthest away point into x2
      x3=x;f3=f; x=x2;f=f2; x2=x3;f2=f3;
    end;   
    % compute new point xnext and accuracy 
    f10=(f1-f)/(x1-x); 
    f20=(f2-f)/(x2-x); 
    f120=(f10-f20)/(x1-x2); 
    f12=(f1-f2)/(x1-x2); 
    df1=f12+f10-f20;
    d1=df1^2;
    d2=4*f1*f120;
    den=df1+sign(df1)*sqrt(max(0,d1-d2));
    if d2>d1,
      % truncated parabolic optimization step
      % sign change safeguards not yet added
      xnext=x1+0.5*(x2-x1-f12/f120);
    elseif den~=0,
      % Muller step
      xnext=x1-2*f1/den; 
    else
      xnext=0.5*(x1+x2);
    end;
  end;
  acc=abs(x2-x1);
  fac=(xnext-x1)/(x2-x1);
  if fac<-0.01,
    % 10-fold extrapolation
    xnext=x1+10*(x1-x2);
  elseif fac<0.01,
    % tiny step
    xnext=x1-0.01*sign(fac)*(x2-x1);
    % tiny extrapolation
    % xnext=x1+0.01*(x1-x2);
  end;
%  if abs(f)<abs(f2), 
    % exchange with second point to ensure that best info is present
  if ~(x==x1|x==x2), 
    % exchange with second point to ensure that new info is present
    x3=x;f3=f; x=x2;f=f2; x2=x3;f2=f3;
  end;   

else 
  % update bracket
  if f1*f>0, q=(x2-x)/(x2-x1);x3=x1;f3=f1;x1=x;f1=f;
  else       q=(x1-x)/(x1-x2);x3=x2;f3=f2;x2=x;f2=f;
  end;
  if q<qmin & sign_change~=3,	% forces two consecutive mean steps 
    sign_change=1; 
    if sign_change==3, disp(3); end;
  else
    sign_change=sign_change+1;
  end;
end;

if prt>1, 
  disp(sprintf('   q=%7.1e;    sign_change = %2.0f',q,sign_change))
  disp(sprintf('%9.2e  %15.8e  %15.8e  %9.2e',f1,x1,x2,f2)) 
elseif prt==1, 
  disp([log10(abs(f1)),q,sign_change])
end; 

if sign_change,
  % compute new point xnext and accuracy
  if sign_change==1, 
    % secant step
    if abs(f1)<abs(f2), 
      xnext = x1 - (x1 - x2)*f1/(f1 - f2);
      acc=abs(xnext-x1); 
    else 
      xnext = x2 - (x2 - x1)*f2/(f2 - f1);
      acc=abs(xnext-x2);
    end;
  elseif sign_change==2, 
    % safeguarded secant extrapolation step 
    xmid=0.5*(x1+x2);
    xsec = x + (x - x3)*abs(f/(f - f3));
    if x<xmid, xnext=min(xsec,xmid);
    else       xnext=max(xsec,xmid); 
    end;
    acc=abs(xnext-x);
  else
    % safeguarded geometric mean step
    if x1*x2>0, xnext=x1*sqrt(x2/x1);
    elseif x1==0, xnext=0.1*x2;
    elseif x2==0, xnext=0.1*x1;
    else          xnext=0;
    end;
    acc=max(abs(xnext-[x1,x2]));
  end;
end;

work(1)=acc;
work(2)=sign_change;
work(3)=x1;work(4)=f1;	% best point if no sign change
work(5)=x2;work(6)=f2;
work(7)=x3;work(8)=f3;
work(9)=nf;