# Set show to True to see guesses.
show = False

def iter_solve(guess, done, update, iteration_limit=32):
    """Return the result of repeatedly applying UPDATE, 
    starting at GUESS, until DONE yields a true value 
    when applied to the result.  Causes error if more than
    ITERATION_LIMIT applications of UPDATE are necessary."""

    while not done(guess):
        if iteration_limit <= 0:
            raise ValueError("failed to converge")
        guess, iteration_limit = update(guess), iteration_limit-1
        if show:
            print("Guess {}".format(guess))
    return guess

def newton_solve(func, deriv, start, tolerance):
    """Return x such that |FUNC(x)| < TOLERANCE, given initial
    estimate START and assuming DERIV is the derivatative of FUNC."""
    def close_enough(x):
        return abs(func(x)) < tolerance
    def newton_update(x):
        return x - func(x) / deriv(x)

    return iter_solve(start, close_enough, newton_update)

def square_root(a):
     return newton_solve(lambda x: x*x - a, lambda x: 2 * x, 
                         a/2, 1e-5)

def logarithm(a, base = 2):
     return newton_solve(lambda x: base**x - a, 
                         lambda x: x * base**(x-1),
                         1, 1e-5)

def iter_solve2(guess, done, update, state=None):
    """Return the result of repeatedly applying UPDATE, 
    starting at GUESS and STATE, until DONE yields a true value 
    when applied to the result.  Besides a guess, UPDATE 
    also takes and returns a state value, which is also passed to
    DONE."""
    while not done(guess, state):
        guess, state = update(guess, state)
        if show:
            print("Guess: {}".format(guess))
    return guess

def secant_solve(func, start0, start1, tolerance):
    def close_enough(x, state):
        return abs(func(x)) < tolerance
    def secant_update(xk, xk1):
        return (xk - func(xk) * (xk - xk1) 
                                / (func(xk) - func(xk1)), 
                xk)
    return iter_solve2(start1, close_enough, secant_update, start0)

# Integration

def integrate_trapezoidal(f, low, high, step):
    """An approximation to the definite intregral of F from
    LOW to HIGH, computed by adding the areas of trapezoids of
    height STEP."""
    
    area = 0
    while low + step < high:
        area += (f(low) + f(low + step)) * step * 0.5
        low += step
    # Before returning, take care of the case where the final value
    # of low is less than high.
    return area + (f(low) + f(high)) * (high - low) * 0.5


# Every floating-point addition or multiplication produces a result that may
# be off by as much as 0.5 units in the last place (ulps), where a unit
# (on our machines) is a power of two.  So, in place of low += step, it
# may be better to use multiplication:

def integrate_trapezoidal2(f, low, high, step):
    """An approximation to the definite intregral of F from
    LOW to HIGH, computed by adding the areas of trapezoids of
    height STEP."""
    
    area = 0
    k = 0
    while True:
        x = low + k * step
        if x + step >= high:
            return area + (f(x) + f(high)) * (high - x) * 0.5
        area += (f(x) + f(x + step)) * step * 0.5
        k += 1


# The loop above redundantly calculates f(x).  If you consider adjacent
# iterations of the loop, you'll see that we add f(x + step) in one step
# and then f(x) in the next: the same value.  So we can group terms and
# save time.

def integrate_trapezoidal3(f, low, high, step):
    """An approximation to the definite intregral of F from
    LOW to HIGH, computed by adding the areas of trapezoids of
    height STEP."""
    
    area0 = 0
    k = 1
    while True:
        x = low + k * step
        if x + step >= high:
            return (f(low) + 2 * area0 + f(x)) * step * 0.5 \
                   + (f(x) + f(high)) * (high - x) * 0.5
        area0 += f(x)
        k += 1