%{ 
   Brent's Method
   John Bryan 
   2020
   Octave 5.1.0
%}


clc;
echo off;
clear all;
close all;
pkg load symbolic;
format long;
syms x;


function brentroot=BrentMethod(f,a,b,s,tolerance)
   if abs(f(a))<abs(f(b))
       temp=a;
       a=b;
       b=temp;
   endif
   c = a;
   mflag=1;
   while (((f(b)!=0) && (f(s)!=0)) && (abs(b-a)>tolerance))
	 if ((f(a) != f(c)) && (f(b) != f(c)))
	      % inverse quadrature interpolation
	      term1 = double((a*f(b)*f(c)) / ((f(a)-f(b))*(f(a)-f(c))));
	      term2 = double((b*f(a)*f(c)) / ((f(b)-f(a))*(f(b)-f(c))));
	      term3 = double((c*f(a)*f(b)) / ((f(c)-f(a))*(f(c)-f(b))));
	      s = double(term1 + term2 + term3);
	 else
	     %secant
	     s = double(b-(f(b)*(b-a)/(f(b)-f(a))));
	 endif
	 if (!(((3*a+b)/4<s)&&(s<b)) ||
	    (mflag&&(abs(s-b)>=abs(b-c)/2)) ||
	    ((!mflag)&&(abs(s-b)>=abs(c-d)/2)) ||
	    (mflag&&(abs(b-c)<abs(tolerance))) ||
	    ((!mflag)&&((abs(c-d))<abs(tolerance))))
	    %bisection 
	    s = double((a+b)/2);
	    mflag=1;
	 else
	    mflag=0;
	 endif
     d = c;
     c = b;
     if (f(a)*f(s) < 0)
	b = s;
     else
	a = s;
     endif
     if (abs(f(a)) < abs(f(b)))
	temp=a;
	a=b;
	b=temp;
     endif
   endwhile
   brentroot=s;
endfunction


f={@(x)(x+4)*(x+1),@(x)sin(x)};
a=double([-5,3*pi/4]);
b=double([-3,5*pi/4]);
tolerance=2e-8;
s=b;
for i=1:length(f)
   brentroot=BrentMethod(f{i},a(i),b(i),s(i),tolerance)
end