import numpy
from matplotlib.pyplot import *
                
def integrale(N):
    x = numpy.linspace(0,numpy.pi,N)
    f = numpy.sin(x)
    return abs(f.sum()*numpy.pi/N-2.0)
                

e = integrale(100)
                

e = integrale(1000)
                

def integrale2(N):
    dx = numpy.pi/N
    def f(i,j):
        return numpy.sin(i*dx)*numpy.cos(j*dx)
    f = numpy.fromfunction(f,(N,N))
    return abs(f.sum()*(numpy.pi/N)**2)
                

e=integrale2(100)
                

e=integrale2(1000)
                

import numpy.random
import random
import math
import numpy
                

def integration1d(fonction,a,b,N):
    x = a+(b-a)*numpy.random.random_sample(N)
    p = 1.0/(b-a)
    f = fonction(x)
    moyenne = f.sum()/(N*p)
    g = f*f
    variance = g.sum()*1.0/(N*p*p)-moyenne*moyenne
    return (moyenne,math.sqrt(variance/N)*1.96)
                 

def f(x):
    return numpy.sin(x)
(integrale,intervalle) = integration1d(f,0,math.pi,1000)
                 

def integration1d_preferentiel(fonction,a,b,N):
    delta = (b-a)*1.0/3
    c1 = a+delta
    c2 = a+2*delta
    somme = 0.0
    somme2 = 0.0
    for i in range(N):
        p = random.randint(1,4)
        if p==1:
            x = random.uniform(a,c1)
            f = fonction(x)*4
        elif p==2 or p==3:
            x = random.uniform(c1,c2)
            f = fonction(x)*2
        else:
            x = random.uniform(c2,b)
            f = fonction(x)*4
        somme += f
        somme2 += f*f
    moyenne = somme*delta/N
    variance = somme2*delta*delta/N-moyenne*moyenne
    return (moyenne,math.sqrt(variance/N)*1.96) 
                   

(integrale,intervalle) = integration1d_importance(f,0,math.pi,1000)
                   

def integration2d(fonction,a,b,c,d,N):
    x = a+(b-a)*numpy.random.random_sample(N)
    y = c+(d-c)*numpy.random.random_sample(N)
    p = 1.0/((b-a)*(d-c))
    somme = 0.0
    somme2 = 0.0
    for i in range(N):
        f = fonction(x[i],y[i])
        somme += f
        somme2 += f*f
    moyenne = somme/(p*N)
    variance = somme2/(p*p*N)-moyenne*moyenne
    return (moyenne,math.sqrt(variance/N)*1.96)
                  

def f(x,y):
    if x*x+y*y <= 1:
        return 1
    else:
        return 0

(integrale,intervalle) = integration2d(f,-1,1,-1,1,10000)
                   

def integrationDisque(fonction,R,N):
    theta = 2*numpy.pi*numpy.random.random_sample(N)
    r = R*numpy.sqrt(numpy.random.random_sample(N))
    p = 1.0/(numpy.pi*R*R)
    somme = 0.0
    somme2 = 0.0
    for i in range(N):
        f = fonction(r[i],theta[i])
        somme += f
        somme2 += f*f
    moyenne = somme/(p*N)
    variance = somme2/(p*p*N)-moyenne*moyenne
    return (moyenne,math.sqrt(variance/N)*1.96)
                       

def f(r,theta):
    return 1-r*r
    
(integrale,intervalle) = integrationDisque(f,1.0,10000)
                        

def volume(dim,R1,R2,N):
    somme = 0.0
    somme2 = 0.0
    R12 = R1*R1
    R22 = R2*R2
    p = 1.0/math.pow(2*R2,dim)
    for i in range(N):
        r2 = 0.0
        for d in range(dim):
            x = random.uniform(-R2,R2)
            r2 += x*x
        if r2>= R12 and r2<=R22:
            f = 1.0
        else:
            f = 0.0
        somme += f
        somme2 += f*f
    moyenne = somme/(p*N)
    variance = somme2/(p*p*N)-moyenne*moyenne
    return (moyenne,variance,math.sqrt(variance/N)*1.96)
                      

(m,v,e) = volume(3,0,1,10000)
                      

(m,v,e) = volume(5,0,1,10000)
                      

(m,v,e) = volume(5,0,1,10000*90)
                      

(m,v,e) = volume(5,0.9,1,10**5)
                       

(m,v,e) = volume(7,0,1,10**7)
                       

(m,v,e) = volume(7,0.9,1,10**7)