1 #=  
 2     Brent algorithm 
 3     Julia 1.1.1
 4     John Bryan 2020
 5 =#
 6 
 7 
 8 using Printf
 9 using Pkg
10 
11 
12 function brent(f,a,b,s,tolerance)::Float64
13    if abs(f(a))< abs(f(b))
14        temp = a;
15        a = b;
16        b = temp;
17    end
18    c = a;
19    d = c;
20    mflag=true;
21    while (((f(b)!=0) && (f(s)!=0)) && (abs(b-a)>tolerance))
22          # inverse quadratic interpolation
23          if ((f(a) != f(c)) && (f(b) != f(c)))
24               term1 = (a*f(b)*f(c)) / ((f(a)-f(b))*(f(a)-f(c)));
25               term2 = (b*f(a)*f(c)) / ((f(b)-f(a))*(f(b)-f(c)));
26               term3 = (c*f(a)*f(b)) / ((f(c)-f(a))*(f(c)-f(b)));
27               s = term1 + term2 + term3;
28          else
29              # secant linear interpolation
30              s = b-(f(b)*(b-a)/(f(b)-f(a)));
31          end
32          if (!(((3*a+b)/4<s)&&(s<b)) ||
33             (mflag&&(abs(s-b)>=abs(b-c)/2)) ||
34             ((!mflag)&&(abs(s-b)>=abs(c-d)/2)) ||
35             (mflag&&(abs(b-c)<abs(tolerance))) ||
36             ((!mflag)&&((abs(c-d))<abs(tolerance))))
37             # bisection 
38             s = (a+b)/2;
39              mflag=true;
40          else
41              mflag=false;
42          end
43      d = c;
44      c = b;
45      if (f(a)*f(s) < 0)
46         b = s;
47      else
48         a = s;
49      end
50      if (abs(f(a)) < abs(f(b)))
51         temp = a;
52         a = b;
53         b = temp;
54      end
55    end
56    convert(Float64,s);
57    return s;
58 end
59 
60 
61 f(x) = (x+4)*(x-1);
62 g(x) = sin(x);
63 a = [-5,3*pi/4];
64 b = [-3,5*pi/4];
65 s = b;
66 fg = [f,g]
67 tolerance = 2e-8;
68 for i=1:length(fg)
69    @printf("\n%lf\n", brent(fg[i],a[i],b[i],s[i],tolerance));
70 end