#!/net/python/bin/python
from numpy import *
import pylab as p

def HitOrMiss(f,N,xmin,xmax):
    frange = [f(x) for x in arange(xmin,xmax,0.01)]
    fmax = max( frange )
    fmin = min( frange )

    x = random.random(N)
    y = random.random(N)

    x = x * (xmax - xmin) -xmin
    y = y * (fmax - fmin) - fmin
    
    newx = x[ where(y < f(x)) ]

    return (xmax-xmin)*(fmax-fmin)*(1.0*len(newx)) / ( 1.0*N )

def error(measured,exact):
    return abs( (measured-exact)/exact )

def tableprint(f,Nlist,xmin,xmax,exact):
    print "N\tI_num\t\terror"
    for N in Nlist:
        num = HitOrMiss(f,N,xmin,xmax)
        print '%i\t' % N,
        print '%s\t\t' % round(num),
        print round(error(num,exact))

def f(x):
    return ( sin(8*x) )**4

def round(x, n=2):
    """
    rounds x to n significant digits, returns a string representation
    in scientific notation
    """
    if n < 1:
        raise ValueError("number of significant digits must be >= 1")
    return '%.*g' % (n,x)

if __name__ == "__main__":
    
    Nlist = [1.0]
    while True:
        if Nlist[-1]*2 > 200000:
            break
        else:
            Nlist.append(Nlist[-1]*2)
    
    tableprint(f,Nlist,0,2*pi,3*pi/4)

    Elist = []
    for N in Nlist:
        Elist.append( error(HitOrMiss(f,N,0,2*pi),3*pi/4 ) )
    Elist = array(Elist)
    
    p.loglog(Nlist,Elist)
    p.xlabel('N')
    p.ylabel('Error')
    p.title('Error in Hit-Or-Miss integration of ( sin(8x) )^8 for N points')
    p.show()