# -*- coding: utf-8 -*-
"""
Created on Sun Mar 16 22:14:48 2014

@author:Paul Hofstra
"""
from __future__ import division
from math import sin
import sys

EPS = sys.float_info.epsilon

def dekker(fun, x1, x2, args=(), xtol=1e-8, disp=False):
    """see wiki Brent
       c is last iterate of b
       a, b bracket solution abs(fb) < abs(fa)
       if f(s) has same sign as fa, then old b beomes new contrapoint
       the solution is within 6*EPS + 2 * xtol
    """
    a, b   = float(x1), float(x2)
    fa, fb = fun(a, *args), fun(b, *args)
    if fa * fb > 0.0: raise ValueError("f(%f)=%f - f(%f)=%f" % (a, b, fa, fb))
    c, fc = a, fa
    feval = 2
    while 1:
        if abs(fa) < abs(fb):
            a, fa, b, fb, c, fc = b, fb, a, fa, b, fb
        xm = (a - b) / 2  # bisection step
        tol1 = 2 * EPS * abs(b) + xtol
        if abs(xm) <= tol1 or fb == 0:
            return b
        # now calculate secant step - use points b and c
        s = fb / fc
        p = (c - b) * s
        q = 1 - s
        if p > 0:
            q = -q
        p = abs(p)
        d = p / q if p < (q * xm) else xm
        c, fc = b, fb
        # do not take steps smaller than tol1
        b += (d if abs(d) > tol1 else sign(tol1, xm))
        fb = fun(b, *args)
        feval += 1
        if disp:
            msg = "b" if abs(d) == abs(xm) else "s"
            fmt = "%3d a=%7.4g b=%12.9g c=%7.4g fb=%10.3g xm=%8.2g %s"
            print fmt % (feval, a, b, c, fb, xm, msg)
        if fa * fb > 0:     # if new iterate has same sign as a - old b take s place of a
            a, fa = c, fc

        
class RootNotBracketedException(Exception): pass
class TooManyIterationsException(Exception): pass

def __sign(x, y):
    if y > 0.0  :
        return x
    elif y < 0.0:
        return -x
    else:
        return 0.0
        
def __sss(f):
    return "+" if f > 0 else "-"
    
def brent(fun, x1, x2, args=(), xtol=1e-12, disp=False, maxiter=1000):
    """
    Bisection: k = log2((b-a)/tol) function eval.
    Dekker   : upper bound: 2**k function eval.
    Brent    : theoretical upper bound: (k+1)**2 ? 2
                practical upper bound: 3*(k+1)
    Example:
        tol = 1e-6
        b-a = 1
        Then k = 20
        Upperbounds
            Bisection: 20 (upperbound is also number iterations)
            Dekker   : 2**20
            Brent    : 63 (practical bound)

    """
    if disp: print "BRENT"
    a,  b  = float(x1), float(x2)
    fa, fb = fun(a, *args), fun(b, *args)
    ncalls = 2
    if fa * fb > 0.0:
        raise RootNotBracketedException(a, b, fa, fb)
    c, fc = b, fb
    for iters in range(maxiter):
        if disp: msg = ""
        if fc * fb > 0.0:
            if disp: msg += "c=a "
            c, fc = a, fa
            e = d = b - a
        if abs(fc) < abs(fb):
            if disp: msg += "swap "
            a, b, c, fa, fb, fc  = b, c, b, fb, fc, fb
        tol1 = 2 * EPS * abs(b) + xtol # like in Brent's original paper
        xm = (c - b) / 2
        if abs(xm) <= tol1 or fb == 0:
            return b, iters, ncalls, 0
        if abs(e) >= tol1 and abs(fa) > abs(fb):
            # try quadratic interpolation
            s = fb / fa
            if a == c:
                p = 2 * xm * s
                q = 1 - s
            else:
                q = fa / fc
                r = fb / fc
                p = s * (2 * xm * q * (q - r) - (b - a) * (r - 1))
                q = (q - 1) * (r - 1) * (s - 1)
            if p > 0:
                q = -q
            p = abs(p)
            if 2 * p < abs(e * q) and 2 * p < 3 * xm * q - abs(tol1 * q):
                if disp: msg += ("LINEAR" if a == c else "QUAD  ")
                e = d
                d = p / q
            else:
                e = d = xm
                if disp: msg += "Bisect"
        else:
            e = d = xm
            if disp:
                msg += "Bisect"
        a, fa = b, fb
        if abs(d) > tol1:
            b += d
        else:
            b += __sign(tol1, xm)
        fb = fun(b, *args)
        ncalls += 1
        if disp:
            print "%15s  a %8.5g %s    b %8.5g %s    c %8.5g %s   " \
            "fb %12.5g   xm %8.5g" % \
            (msg, a, __sss(fa), b, __sss(fb), c, __sss(fc), fb, xm)
    raise TooManyIterationsException(maxiter)
    
if __name__ == "__main__":
    import math
        
    def f(x):
        return math.exp(x) - math.exp(0.8)
        
    def g(x):
        return sin(x) * x ** 3
        
    def __dirty(x):
        """zero close to 1.0"""
        D = 1e-2
        epsi = 1e-50
        if x > 1 - D:
            res = ((x - (1 - D)) / D) ** 2
        else:
            res = 0
        return res - epsi
        
    def brent_example(x, delta=1e-3,a=0.0, b=1.0):
        if a + delta <= x <= b:
            return 2 ** (x / delta)
        if x == a:
            return -(b - a - delta) / delta * 2 ** (b / delta)
        return 1
        
    disp = True
    fun = brent_example
    xmin = 0
    xmax = 1
    precision = 8
    _xtol = 10 ** -(precision + 1)
    fmt = "{:10." + str(precision) + "f}"
    
        
    print "Brent"
    res =  brent(fun, xmin, xmax, xtol=_xtol, disp=disp)
    print res[0], res[1]
 
    print "Dekker"
    res =  dekker(fun, xmin, xmax, xtol=_xtol, disp=disp)
    print res


# your code goes here