module montecarlo
	interface
		subroutine plainmc(f, a, b, N, Integral, err)
			implicit none
			real(8), external :: f
			real(8), intent(in) :: a(:), b(:)
			integer, intent(in) :: N
			real(8), intent(out) :: Integral, err
		end subroutine

		subroutine quasirandom(f, a, b, N, Integral, err)
			implicit none
			real(8), external :: f
			real(8), intent(in) :: a(:), b(:)
			integer, intent(in) :: N
			real(8), intent(out) :: Integral, err
		end subroutine
		
		recursive subroutine strata(f, a, b, N, acc, Integral, err)
			implicit none
			real(8), external :: f
			real(8), intent(in) :: a(:), b(:), acc
			integer, intent(in) :: N
			real(8), intent(out) :: err, Integral
		end subroutine
	end interface
end module

subroutine plainmc(f, a, b, N, Integral, err)
	implicit none
	real(8), external :: f
	real(8), intent(in) :: a(:), b(:)
	integer, intent(in) :: N
	real(8), intent(out) :: Integral, err
	real(8) :: x(size(a)), vol, s, s2, random(size(a)), sigma, avg, fx
	integer :: i, j

	s = 0.0
	s2 = 0.0
	do i = 1,N
		call random_number(random)
		x = a + random*(b-a)
		fx = f(x)
		s = s + fx
		s2 = s2 + fx**2
	enddo
	avg = s/N
	sigma = sqrt(s2/N - avg**2)

	vol = product(b-a)
	Integral = vol*avg
	err = vol*sigma/sqrt(real(N))
end subroutine plainmc

subroutine quasirandom(f, a, b, N, Integral, err)
	implicit none
	integer, parameter :: long = 10
	real(long), parameter :: offset = 2.71828182845904523536028747135266249775724709369995
	real(8), external :: f
	real(8), intent(in) :: a(:), b(:)
	integer, intent(in) :: N
	real(8), intent(out) :: Integral, err
	real(8) :: x(size(a)), vol, s, s2, random(size(a)), sigma, avg, fx
	integer :: i, j, d
	real(long) :: alpha(size(a))
	d = size(a)

	! Construct alpha of irational numbers: 
	do i = 1,d
		alpha(i) = sqrt(2*i+offset)
	enddo

	s = 0.0
	s2 = 0.0
	do i = 1,N
		random = (i*alpha)-floor(i*alpha)
		x = a + random*(b-a)
		fx = f(x)
		s = s + fx
		s2 = s2 + fx**2
	enddo
	avg = s/N
	sigma = sqrt(s2/N - avg**2)

	vol = product(b-a)
	Integral = vol*avg
	err = vol*sigma/sqrt(real(N))
end subroutine quasirandom


recursive subroutine strata(f, a, b, N, acc, Integral, err)
	!use montecarlo, only : plainmc
	implicit none
	real(8), external :: f
	real(8), intent(in) :: a(:), b(:), acc
	integer, intent(in) :: N
	real(8), intent(out) :: err, Integral
	real(8), dimension(size(a)) :: aR, bL, varL, varR, SL, SR, S2L, S2R, random, x
	integer, dimension(size(a)) :: NL, NR
	real(8) :: int1, int2, err1, err2, vol, avg, fx, sigma, s, s2
	integer :: i, d, k, d_big
	d = size(a)
	NL = 0;NR = 0;SL=0;SR=0;S2L=0;S2R=0;s=0;s2=0;

	! Estimate integral and error with Plain Monte Carlo:
	do k = 1,N
		call random_number(random)
		x = a + random*(b-a)
		fx = f(x)
		s = s+fx
		s2 = s2+fx**2
		do i = 1,d
			! Check wich half x lies in:
			if (x(i) < (a(i)+b(i))/2) then
				NL(i) = NL(i)+1
				SL(i) = SL(i)+fx
				S2L(i) = S2L(i)+fx**2
			else
				NR(i) = NR(i)+1
				SR(i) = SR(i)+fx
				S2R(i) = S2R(i)+fx**2
			endif
		enddo
	enddo

	! Calculate the integral and error from Plain Monte Carlo:
	vol = product(b-a)
	avg = s/N
	sigma = sqrt(s2/N - avg**2)
	Integral = vol*avg
	err = vol*sigma/sqrt(real(N))
	
	if (err >= acc) then
		! Calculate variance for the "two sides":
		varL = S2L/NL - (SL/NL)**2
		varR = S2R/NR - (SR/NR)**2
	
		! Find the dimension with bigest variance:
		d_big = maxloc(abs(varL-varR), 1)

		print *, d_big

		! Divide the volume in two along d_big:
		aR = a; bL = b
		aR(d_big) = (a(d_big)+b(d_big))/real(2)
		bL(d_big) = aR(d_big)

		! 
		call strata(f, a, bL, N, acc, int1, err1)
		call strata(f, aR, b, N, acc, int2, err2)
	
		Integral = int1+int2
		err = sqrt(err1**2 + err2**2)
	endif
end subroutine strata