baluev_3.1-3

import sys, math, random, ROOT

input_x_a = 0.01
input_x_b = 1
input_y_a = 0
input_y_b = 1

def f(x,y):
    #return 1 / ( math.pow(x,2) * ( y + math.sin(x)/2 ) ) - 1/(math.pow(x,2) + y)
    return 1 / ( math.pow(x,2) * ( y + math.sin(x)/2 ) ) - 1/(math.pow(x,2) + y) + 1/( math.pow(x,2) * ( y + 0.01/2 ) )


def single(f, x_limits, y_limits):
    x_a, x_b = x_limits
    y_a, y_b = y_limits
    x = ROOT.gRandom.Uniform(x_a, x_b)
    y = ROOT.gRandom.Uniform(y_a, y_b)
    return (x_b - x_a) * (y_b - y_a) * f(x, y)

outerN = 10
innerN = 1000

sumtot = 0
sum2tot = 0
ntot = 0

for i in range(outerN):
    sumloc, sum2loc, nloc = [0,0,0]
    for j in range(innerN):
        r = (single(f, (input_x_a, input_x_b), (input_y_a, input_y_b)))
        sumloc += r
        sum2loc += r**2

    sumtot += sumloc
    sum2tot += sum2loc
    ntot += 1000

    run_average = 525.027 + sumloc / innerN
    run_variance = (sum2loc / innerN) - (run_average)**2
    run_sigma = math.sqrt( run_variance / innerN )

    average = 525.027 + sumtot / ntot
    variance = (sum2tot / ntot) - (sumtot / ntot)**2
    sigma = math.sqrt( variance / ntot )

    print run_average, " +/- ", run_sigma,
    print average, " +/- ", sigma

def f_for_drawing(x, p):
    return f(x[0], x[1])

tf2 = ROOT.TF2("f", f_for_drawing, input_x_a, input_x_b, input_y_a, input_y_b, 1)
tf2.Draw("surf1z")
sys.stdin.readline()