##
##  b r e n t D e k k e r . R
##


brentDekker <- function(f, a, b, ...,
                        maxiter = 100, tol = .Machine$double.eps^0.5)
# Brent and Dekker's root finding method,
# based on bisection, secant method and quadratic interpolation
{
    stopifnot(is.numeric(a), is.numeric(b),
              length(a) == 1, length(b) == 1)

    err <- try(fun <- match.fun(f), silent = TRUE)
    if (class(err) == "try-error") {
        stop("Argument function 'f' not known in parent environment.")
    } else {
        f <- function(x) fun(x, ...)
    }

	x1 <- a; f1 <- f(x1)
	if (f1 == 0) return(a)
	x2 <- b; f2 <- f(x2)
	if (f2 == 0) return(b)
	if (f1*f2 > 0.0)
	    stop("Root is not bracketed in [a, b].")

	x3 <- 0.5*(a+b)
	# Beginning of iterative loop
	niter <- 1
	while (niter <= maxiter) {
		f3 <- f(x3)
		if (abs(f3) < tol) {
		    x0 <- x3
		    break
		}

		# Tighten brackets [a, b] on the root
		if (f1*f3 < 0.0) b <- x3 else a <- x3
		if ( (b-a) < tol*max(abs(b), 1.0) ) {
		    x0 <- 0.5*(a + b)
		    break
	    }

		# Try quadratic interpolation
		denom <- (f2 - f1)*(f3 - f1)*(f2 - f3)
		numer <- x3*(f1 - f2)*(f2 - f3 + f1) + f2*x1*(f2 - f3) + f1*x2*(f3 - f1)
		# if denom==0, push x out of bracket to force bisection
		if (denom == 0) {
			dx <- b - a
		} else {
			dx <- f3*numer/denom
		}

		x <- x3 + dx
		# If interpolation goes out of bracket, use bisection.
		if ((b - x)*(x - a) < 0.0) {
			dx <- 0.5*(b - a)
			x  <- a + dx;
		}

		# Let x3 <-- x & choose new x1, x2 so that x1 < x3 < x2.
		if (x1 < x3) {
			x2 <- x3; f2 <- f3
		} else {
			x1 <- x3; f1 <- f3
		}

		niter <- niter + 1
		if (abs(x - x3) < tol) {
		    x0 <- x
		    break
	    }
		x3 <- x;
	}

    if (niter > maxiter) {
        warning("Maximum numer of iterations, 'maxiter', has been reached.")
        return(NA)
    } else {
        return(x0)
    }
}