# numerical.py
# Walker M. White (wmw2)
# November 5, 2015
"""Module to show while-loops and numerical computations.

This is one of the most powerful uses of a while loop: using it to run
a computation until it converges.  There are a lot of algorithms from
Calculus and Numerical Analysis that work this way.  In addition, this
is how we approach the clustering example from A5"""

def sqrt(c,err=1e-6):
    """Returns: the square root of c to with the given margin of error.
    
    Use Newton's Method to find the root of the polynomial f(x) = x^2-c
    
    Newton's Method produces a sequence where
    
        x_(n+1) = x_n - f(x_n)/f'(x_n) = x_n - (x_n*x_n-c)/(2x_n)
    
    which we simplify to
    
        x_(n+1) = x_n/2 + c/2 x_n
    
    We can stop this process and use x_n as an approximation of the square root
    of c.  The error when we do this is less than |x_n*x_n - c|.  So we loop
    until this error is less than err.
    
    Parameter c: The number to compute square root
    Precondition: c >= 0 is a number
    
    Parameter err: The margin of error (OPTIONAL: default is e-6)
    Precondition: err > 0 is a number"""
    x = c/2.0
    
    while abs(x*x-c) > err:
        # Compute next in Newton's Method sequence
        x = x/2.0+c/(2.0*x)
    
    return x


def pi(err):
    """Returns: the constant pi to with the given margin of error.
    
    We use the Power Series pi = 4(1 - 1/3 + 1/5 - 1/7 + ... + (-1 ** n)/(2n+1) + ...)
    
    We can approximate pi ~ 4(1 - 1/3 + 1/5 - 1/7 + ... + (-1 ** n)/(2n+1))
    In other words, we stop the infinite sum at some value of n.  The error when we 
    do this is < 4/(2n+1). So we loop until we reach a value of n with 4/(2n+1) <= error.
    
    Parameter err: The margin of error
    Precondition: err > 0 is a number."""
    value = 4       # The value of pi start with n = 0
    
    power = -1      # Alternate between 1 and -1
    n = 1           # The next value of n
    while 4.0/(2*n+1) > err:
        factor = (4.0*power)/(2*n+1)    # 4(-1 ** n)/(2n+1)
        value = value + factor          # Add factor
        power = -power
        n = n + 1
    
    return value