RootKTSMM(k,x0,x1,epsm)

function [Root Ze] = RootKTSMM(k,x0,x1,epsm)
% algorithm MM strictly following the J. of Complexity paper
% On a class of Omnibus Algorithms for Zero Finding
% A. Kozek & A. Trzmielak-Stanislawska, 
% Journal of Complexity 5 (1989) pp. 80-95
% we follow notation of that paper adding fun at the name of the function
% used in the paper.

CTR = 2;  

% CTR is a treshold variable controlling when to switch from the Brent algorithm to bisection
%  CTR = 2 was found optimal on a broad set of test functions, see
%  Kozek, Andrzej S. and Trzmielak-Stanislawska, Anna (1988)
%  On dependence of {MR}- and {MM}-algorithms upon the value of
%        switching control variable {CTR}},
%  J. Comput. Appl. Math. 23(1) pp. 109-115 

if (nargin == 0)
    x0 = -1;
    x1 = 0.8;
    k = @try4kts3;
end;
if (nargin < 4)
    epsm = m_eps();
end;

f0 = k(x0);
f1 = k(x1);
if (f0*f1>0)
    Ze = false;
    Root = 0;
    display('f(x0)*f(x1)>0');
    return;
end;
if abs(f1)<= abs(f0)
    b = x1;
    fb = f1;
    a = x0;
    fa = f0;
    c = x0;
else
    a = x1;
    fa = f1;
    b = x0;
    b = f0;
    c = x1;
end;

% moreover
d = a;
ext = 0;
ctr = 0;

fc = fa;
whiledo = true;

while whiledo
    if ((ext<=1)&&(ctr<=CTR))
        x = funw(funl(a,b,fa,fb),b,c);
    elseif ((ext == 2)&&(ctr<=CTR))
        x = funw(funr(a,b,d,fa,fb,fd),b,c);
    else
        x = funm(b,c);
    end;

    fx = k(x);

    [a1,b1,c1,d1,fa1,fb1,fc1,fd1,ext1] = functionMiu(a,b,c,fa,fb,fc,x,fx,ext,CTR);
    if fc1*fb1>0
        display('f(x0)*f(x1)>0');
    else
        a=a1;
        b=b1;
        c=c1;
        d=d1;
        fa=fa1;
        fb=fb1;
        fc=fc1;
        fd=fd1;
        ext=ext1;
    end;
    if ext>CTR
        ctr = ctr+1;
    end;
    db = fundelta(b);
    whiledo = ((abs(b-c)-2.*db)>0)||(abs(b-c)>=epsm);
end;
Ze = (fb.*fc <= 0);
Root = x;
o=1;


function l = funl(a,b,fa,fb)
dfba = fb - fa;
if (dfba == 0)
    if (fa == 0)
        l = NaN;
    else
        l = b;
    end;
else
    l = b-fb.*(b-a)./dfba;
end;

function r = funr(a,b,d,fa,fb,fd)
alpha = (fb-fd).*fa./(b-d);
beta = (fa-fd).*fb./(a-d);
dba = beta-alpha;
if (dba == 0)
    if (alpha == 0)
        r = b;
    else
        r = NaN;
    end;
else
    r = b-beta.*(b-a)./dba;
end;

function d = fundelta(b)
delta0 = 1e-8;
delta1 = 1e-5;
d = delta0+delta1.*abs(b);

function m = funm(b,c)
m=(b+c)./2;

function h = funh(b,c)
h = b+sign(c-b).*fundelta(b);

function w = funw(lambda,b,c)
h = funh(b,c);
m = funm(b,c);
if ((lambda - h).*(lambda-m)<0)
    w = lambda;
else
    insidebm = (lambda-b).*(lambda-m);
    db = fundelta(b);
    if ((insidebm <=0) && (abs(lambda - b)<= db))
        w =  funh(b,c);
    else
        w = m;
    end;
end;

function [a1,b1,c1,d1,fa1,fb1,fc1,fd1,ext1] = functionMiu(a,b,c,fa,fb,fc,x,fx,ext,CTR)
absfx = abs(fx);
absfc = abs(fc);
A = ((absfx<=absfc)&&(sign(fx).*sign(fc)<0));
B = ((absfx<=absfc)&&(sign(fx).*sign(fc)>0));
C = ((absfx>absfc)&&(sign(fx).*sign(fc)<0));
D = ((absfx>absfc)&&(sign(fx).*sign(fc)>0));

if A
    b1 = x;
    fb1 = fx;
    c1 = c;
    fc1 = fc;
    a1 = b;
    fa1 = fb;
    d1 = a;
    fd1 = fa;
elseif B
    b1 = x;
    fb1 = fx;
    c1 = b;
    fc1 = fb;
    a1 = b;
    fa1 = fb;
    d1 = a;
    fd1 = fa;
elseif C
    b1 = c;
    fb1 = fc;
    c1 = x;
    fc1 = fx;
    a1 = x;
    fa1 = fx;
    d1 = b;
    fd1 = fb;
elseif D
    b1 = b;
    fb1 = fb;
    c1 = x;
    fc1 = fx;
    a1 = x;
    fa1 = fx;
    d1 = a;
    fd1 = fa;
else % it can only be that fx==0 or fc == 0
    if fx == 0
        b1 = x;
        fb1 = fx;
        c1 = x;
        fc1 = fx;
        a1 = x;
        fa1 = fx;
        d1 = x;
        fd1 = fx;
    else
        b1 = c;
        fb1 = fc;
        c1 = c;
        fc1 = fc;
        a1 = c;
        fa1 = fc;
        d1 = c;
        fd1 = fc;
    end
end
if (A||C||((b1 == funm(b,c))&&(ext<=CTR)))
    ext1 = ext+1;
else
    ext1 = 0;
end;



function meps = m_eps()
eps = 1.0;
tol1 = 1+eps;
while tol1 > 1.0
    eps = eps/2;
    tol1 = 1+eps;
end
meps = sqrt(eps);