import random,math
def radapt(a,b,F,acc=1e-3,eps=1e-3,old_N=0,old_sum_f=0,old_sum_f2=0):
	N=8; sum_f=0; sum_f2=0
	N_left=0; sum_f_left=0; sum_f2_left=0
	N_right=0; sum_f_right=0; sum_f2_right=0
	for i in range(N):
		rnd=1.0-random.random()
		x=a+rnd*(b-a)
		f=F(x)
		sum_f  += f; sum_f2 += f*f
		if x<(a+b)/2.0:
			N_left  += 1; sum_f_left  += f; sum_f2_left  += f*f
		else:
			N_right += 1; sum_f_right += f; sum_f2_right += f*f
	sum_f_tot  = sum_f  + old_sum_f
	sum_f2_tot = sum_f2 + old_sum_f2
	N_tot = N + old_N
	average_f = sum_f_tot/N_tot
	average_f2 = sum_f2_tot/N_tot
	variance = math.sqrt(average_f2-average_f*average_f)
	integral = average_f*(b-a)
	error = variance*(b-a)/math.sqrt(N_tot)
	tolerance = acc + abs(integral)*eps
	if error < tolerance :
		return integral
	else :
		left_integral = \
radapt(a,(a+b)/2.0,F,acc/1.414,eps,N_left,sum_f_left,sum_f2_left)
		right_integral = \
radapt((a+b)/2.0,b,F,acc/1.414,eps,N_right,sum_f_right,sum_f2_right)
		return left_integral+right_integral

def ccradapt(a,b,F,acc=1e-3,eps=1e-3):
	return radapt(0,math.pi,\
	lambda x: F((a+b)/2.0+(b-a)/2.0*math.cos(x))*math.sin(x)*(b-a)/2.0,\
	acc,eps,0,0,0)

import sys
sys.setrecursionlimit(10000)
a=0;b=1;acc=1e-5;eps=0
print("1) Variance as error estimate")
print("---- integral from 0 to 1 of x ------"," acc=",acc," eps=",eps)
npoints=0
def F(x) : 
	global npoints
	npoints+=1
	return x
integral=radapt(a,b,F,acc,eps)
print("exact =",1/2)
print("result=",integral," npoints=",npoints)
print("---- integral from 0 to 1 of sqrt(x) ------"," acc=",acc," eps=",eps)
npoints=0
def F(x) : 
	global npoints
	npoints+=1
	return math.sqrt(x) 
integral=radapt(a,b,F,acc,eps)
print("exact =",2/3)
print("result=",integral," npoints=",npoints)
print("---- clenshaw-curtis integral from 0 to 1 of 1/sqrt(x) ------"," acc=",acc," eps=",eps)
npoints=0
def F(x) : 
	global npoints
	npoints+=1
	return 1/math.sqrt(x) 
integral=ccradapt(a,b,F,acc,eps)
print("exact =",2)
print("result=",integral," npoints=",npoints)