function [value,it]=Brent(a,b,machep,t,f)
% Se calcula un c en [a,b] tal que f(c)=0.
%   En [a,b] debe darse un cambio de signo en f(x).
%   La precisión de la raíz es 6*machep*abs(c)+2*t.
%
 sa = a; sb = b;  fa = f(sa); fb = f(sb);
 c = sa; fc = fa; e  = sb-sa; d  = e; it=0;
 while 1, it=it+1;
    if abs(fc)<abs(fb), sa=sb; sb=c; c=sa; fa=fb; fb=fc; fc=fa; end
    tol = 2.0*machep*abs(sb) + t;
    m   = 0.5*(c-sb);
    if abs(m)<=tol || fb==0.0, break, end
    if abs(e)<tol  || abs(fa)<=abs(fb), e = m; d = e;
    else
       s = fb/fa;
       if sa==c, p=2.0*m*s; q=1.0-s;             % Interpolación lineal
       else
          q = fa/fc; r = fb/fc;                  % Interpolación cuadrática
          p = s*(2.0*m*q*(q-r)-(sb-sa)*(r-1.0)); %    inversa
          q = (q-1.0)*(r-1.0)*(s-1.0);
       end
       if 0.0<p, q=-q;
       else
          p=-p;
       end
       s = e; e = d;
       if 2.0*p<3.0*m*q-abs(tol*q) && p<abs(0.5*s*q)
          d = p/q;                               % Interpolación aceptable
       else
          e = m; d = e;                          % Bisección
       end
    end
    sa = sb; fa = fb;                            % Completa paso
    if tol<abs(d), sb=sb+d;
    elseif 0.0<m
       sb = sb + tol;
    else
       sb = sb - tol;
    end
    fb = f(sb);
    if (0.0<fb&&0.0<fc) || (fb<=0.0&&fc<=0.0)
       c = sa; fc = fa; e = sb-sa; d = e;
    end
 end
 value = sb;
end
