function [xvect,nit,eta]=bisection_discont(a,b,toll,nmax,fun,cho,varargin)

% [xvect,nit]=bisection_discont(a,b,toll,nmax,fun,cho,varargin)
%
% Discontinuous bisection method (see for example annexe "Dichotomie discontinue" of http://utbmjb.chez-alice.fr/Polytech/MNBif/coursMNBif.pdf)
%
%  >> INPUT <<
%
%  a, b       two initial values
%  toll       tolerance required
%  nmax       maximum number of iterations
%  fun        function f
%  cho        optionnal integer in {1,2,3,4}
%               * 1      : research of change of values of f (default choice)
%               * 2 or 3 : research of change of sign of f
%                          if cho=2, we have f>=0 or f<0
%                          if cho=3, we have f>0  or f<=0
%               * 4      : research of roots of f, which is not necessarily
%                          continuous (at its roots)
%  varargin   additional variables of function fun
%             Warning, roots of function fun with respect to first
%             variable of function fun
%
%
%  >> OUTPUT <<
%
%  xvect      vector containing all the iterates, the expected value  is
%             xvect(end)
%  nit        number of iterations required to reach the
%             prescribed tolerance
%  eta        optionnal values of eta (see for example annexe "Dichotomie
%             discontinue" of http://utbmjb.chez-alice.fr/Polytech/MNBif/coursMNBif.pdf)
%
%  >> EXAMPLES <<
%
%  1)
%    a) 
%       [xvect,nit]=bisection_discont(-1/2,pi,eps,1000,'sign');nit,xvect(end)
%       [xvect,nit]=bisection_discont(-1/2,pi,0,1000,'sign');nit,xvect(end)
%    b)
%       a=0;
%       b=1;
%       eta=[1 0 0 1 0 1 1 0 0 1];
%       l=a+(b-a)*sum(eta./2.^(1:length(eta)));
%       [xvect,nit,eta]=bisection_discont(a,b,0,1000,@(x)   sign(x-l));xvect(end)-l,eta
%       [xvect,nit,eta]=bisection_discont(a,b,0,1000,@(x,l) sign(x-l),[],l);xvect(end)-l,eta
%       [xvect,nit,eta]=bisection_discont(a,b,0,1000,@(x)  (x-l~=0).*sign(x-l)+(x-l==0).*0./(x-l));xvect(end)-l,eta
%       [xvect,nit,eta]=bisection_discont(a,b,0,1000,@(x)  (x-l~=0).*sign(x-l)+(x-l==0).*3);xvect(end)-l,eta
%       [xvect,nit,eta]=bisection_discont(a,b,eps,1000,@(x)  (x-l~=0).*sign(x-l)+(x-l==0).*1);xvect(end)-l,eta
%  2)
%    a) 
%       [xvect,nit]=bisection_discont(-1/2,pi,eps,1000,@(x) (x>0)*x+(x<=0)*min([x^2*sin(1/x),0]),3);nit,xvect(end)
%    b)
%       a=0;
%       b=1;
%       eta=[1 0 0 1 0 1 1 0 0 1];
%       l=a+(b-a)*sum(eta./2.^(1:length(eta)));
%       [xvect,nit,eta]=bisection_discont(a,b,eps,1000,@(x) (x-l>0)*(x-l)+(x-l<=0)*min([(x-l)^2*sin(1/(x-l)),0]),3);xvect(end)-l,eta
%       [xvect,nit,eta]=bisection_discont(a,b,eps,1000,@(x) (x-l>0)*(x-l)+(x-l<0)*min([(x-l)^2*sin(1/(x-l)),0])+(x-l==0)*3,3);xvect(end)-l,eta
%       [xvect,nit,eta]=bisection_discont(a,b,eps,1000,@(x) (x-l>0)*(x-l)+(x-l>0)*min([(x-l)^2*sin(1/(x-l)),0])+(x-l==0)*(-3),3);xvect(end)-l,eta
%       [xvect,nit,eta]=bisection_discont(a,b,eps,1000,@(x) (x-l>0)*(x-l)+(x-l>0)*min([(x-l)^2*sin(1/(x-l)),0])+(x-l==0)*(0/(x-l)),3);xvect(end)-l,eta
%  3)
%    a) 
%       [xvect,nit]=bisection_discont(-1,pi/2,eps,1000,'atan',4);nit,xvect(end)
%    b)
%       a=0;
%       b=1;
%       eta=[1 0 0 1 0 1 1 0 0 1];
%       l=a+(b-a)*sum(eta./2.^(1:length(eta)));
%       [xvect,nit,eta]=bisection_discont(a,b,0,1000,@(x) atan(x-l),4);xvect(end)-l,eta
%       [xvect,nit,eta]=bisection_discont(a,b,0,1000,@(x) (x-l~=0).*atan(x-l)+(x-l==0).*0./(x-l),4);xvect(end)-l,eta
%       [xvect,nit,eta]=bisection_discont(a,b,eps,1000,@(x) (x-l~=0).*atan(x-l)+(x-l==0).*3,4);xvect(end)-l,eta
%
%
% (c) 2021 by      Jérôme BASTIEN,
%                  LIBM, Polytech, Université Lyon 1
%                  E-Mail : jerome.bastien@univ-lyon1.fr
%

if nargin<=5||isempty(cho);
    cho=1;
end
if cho==2||cho==3
    if feval(fun,a,varargin{:})==0&&feval(fun,b,varargin{:})==0
        error('f(a) et f(b) sont nulles');
    end
end
alpha=valf(a,cho,fun,varargin{:});
beta=valf(b,cho,fun,varargin{:});
if alpha==beta
    error('La fonction auxiliaire a la même valeur en a et b');
end
if isnan(alpha)||isnan(beta)
    error('La fonction auxiliaire n''est pas définie en a ou en b');
end
xvect=zeros(1,nmax);
if nargout>=3
    eta=xvect;
end
err=toll+1;
nit=0;
test=true;
while (nit<nmax&&err>toll&&test)
    nit=nit+1;
    c=(a+b)/2;
    xvect(nit)=c;
    try 
        fc=valf(c,cho,fun,varargin{:});
    catch
        fc=NaN;
    end
    if isnan(fc)||(fc~=alpha&&fc~=beta)
        test=false;
        if nargout>=3
            eta(nit)=1;
        end
    else
        if fc==alpha
            a=c;
            if nargout>=3
                eta(nit)=1;
            end
        else
            b=c;
        end
        
    end
    err=abs(b-a);
end
xvect=xvect(1:nit);
if nargout>=3
    eta=eta(1:nit);
end
end

% nested function
function y=valf(x,cho,fun,varargin)
try
    auxi=feval(fun,x,varargin{:});
catch
    auxi=NaN;
end
if cho==1
    y=auxi;
else
    auxi=sign(auxi);
    if cho==4
        y=auxi;
    else
        if isnan(auxi)
            y=NaN;
        else
            if cho==2
                if auxi>=0
                    y=1;
                else
                    y=0;
                end
            else
                if auxi>0
                    y=1;
                else
                    y=0;
                end
            end
        end
    end
end
end