docs: public static double root(double eps, int maxIterations, double a, double b, Function1D f)
      {
          /*
           * Code has few comments, taken fro algorithum descriptoin http://en.wikipedia.org/wiki/Brent%27s_method#Algorithm ,
           * which is from Brent's book (according to comments, I would like to get the book either way)
           * 
           */
     
          double fa = f.f(a);
          double fb = f.f(b);
          
          if(abs(fa-0) < 2*eps)
              return a;
          if(abs(fb-0) < 2*eps)
              return b;
          
          if(fa * fb >= 0)
              throw new ArithmeticException("The given search interval does not appear to contain the root ");
          
          if(abs(fa) < abs(fb)) //swap
          {
              double tmp = a;
              a = b;
              b = tmp;
              
              tmp = fa;
              fa = fb;
              fb = tmp;
          }
          
          
          double c = a;
          double fc = fa;
          boolean mflag = true;
          double s;
          double d = 0;//inital value dosnt matter, and will not be used
          
          double fs;
          
          do
          {
              if(fa != fc && fb != fc)//inverse quadratic interpolation
              {
                  s = a*fb*fc/( (fa-fb)*(fa-fc) ) + b*fa*fc/( (fb-fa)*(fb-fc) ) + c*fa*fb/( (fc-fa)*(fc-fb) );
              }
              else//secant rule
              {
                  s = b - fb*(b-a)/(fb-fa);
              }
              
              
              //Determin wethor or not we must use bisection
              
              boolean cond1 = (s - ( 3 * a + b) / 4 ) * ( s - b) >= 0;
              boolean cond2 = mflag && (abs(s - b) >= (abs(b - c) / 2));
              boolean cond3 = !mflag && (abs(s - b) >= (abs(c - d) / 2));
              boolean cond4 = mflag && (abs(b-c) < 2*eps);
              boolean cond5 = !mflag && abs(c-d) < 2*eps;
              
              if(cond1 || cond2 || cond3 || cond4 || cond5)//Bisection must be used
              {
                  s = (a+b)/2;
                  mflag = true;
              }
              else
                  mflag = false;
              
              fs = f.f(s);
              d = c;
              c = b;
              
              //adjust the interval accordingly
              if(fa*fs < 0)
              {
                  b = s;
                  fb = fs;
              }
              else
              {
                  a = s;
                  fa = fs;
              }
              
              if(abs(fa) < abs(fb))//swap
              {
                  double tmp = a;
                  a = b;
                  b = tmp;
                  
                  tmp = fa;
                  fa = fb;
                  fb = tmp;
                  
              }
  
          }
          while( fb != 0.0 && fs != 0.0 && abs(b-a) > 2*eps && maxIterations-- > 0);
          
          
          return b;
      }