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

def f_x(x):
    return 1./(x**2 + (cos(x))**2)

def w_x(x,a):
    A = a/(1.-exp(-a*pi))
    return A*exp(-a*x)

def y_x(x,a):
    return (1.-exp(-a*x))/(1-exp(-a*pi))

def x_y(y,a):
    return -(1./a) * log(1-y*(1-exp(-a*pi)))

a_vals = [0.1, 0.3, 0.5, 0.6, 0.7, 0.75, 0.8, 0.81, 0.82, 0.83, 0.84, 0.85, 0.86, 0.87, 0.88, 0.89, 0.9, 1.0]
N = 1E6

print 'a\tmean\tvar_f\tvar_I'
for a in a_vals:
    #generate non-uniform random numbers
    r = x_y( random.random(N),a )
    I = f_x(r)/w_x(r,a)
    
    m1 = I.mean()
    m2 = (I**2).mean()
    var = sqrt(m2 - m1**2)
    print '%.5f\t%.5f\t%.5f%.5f' % (a,m1,var,var/sqrt(N))