'''Bisection method for cubic equations of the form 
f(x) = x**3 + c[0]*x**2 + c[1]*x + c[2] = 0.

Input:
    c = [c[0], c[1], c[2]] ... list of coefficients.
    tol ... absolute accuracy for f(xi), computation stops if |f(xi)| <= tol.
    maxk ... maximal number of iterations, program terminates if maxk iterations
    are reached.
    
Output: 
    xi... root of the function f.
    interval... array of the interval boundaries of the iteration.
    
Typical call of program:
    from python06_ex7a import bisection
    c = [-3, 1, 1]
    bisection(c)

'''
import numpy as np

def bisection(c, tol = 1e-6, maxk = 100):
    '''We chose some default arguments for tol and maxk. That means we can call
    bisection with only c as an input argument.
    Example: bisection(c) is the same as bisection(c, 1e-6, 100)
    '''
    
    
    '''search for a, b with f(a)*f(b) < 0:'''  
    n = 0
    p = c[2]
    q = c[2]
    
    while p*q >= 0:
       n += 1   # Shortcut for n = n+1
       p = n**3 + c[0]*n**2 + c[1]*n + c[2]
       q = -n**3 + c[0]*n**2 - c[1]*n + c[2]

    a = -float(n)
    b = -a
    
     
    '''iteration of the bisection'''
    for k in range(1,maxk+1):
        y = (a + b)/2
        fy = y**3 + c[0]*y**2 + c[1]*y + c[2]
        
        '''We stop if the error is small enough'''
        if abs(fy) <= tol:
            break
            
        elif fy > tol:
            b = y
            
        else:
            a = y
        print 'interval = [',a,b,']'
    
    print 'solution = ', y
    return