IntBrownFunc <- function(K,m){	
	#The goal is to calculate the k-1 fold integration of a two sided-Brownian motion starting from 0, for k in 1...K 
	# m is the index of the partial sum approximating a Brownian Bridge. It also measures the precision of the approximation.
	# IntschauderFunc calculates the k-1 fold integration of any shauder function h_pi
	# grid over which we would like to evaluate the Brownian motion
	
	grid <- seq(0, 1, 2^{- m})
	L.g <- length(grid)
	# Approximation of the k-1 fold integration of a two sided Brownian motion with   #  a partial sum
	
	Zv <- rnorm(2^{m + 1} - 1, 0, 1)
	Z <- rnorm(1, 0, 1)
	
	IntUm <- matrix(0, nrow=K,ncol=L.g)
	IntBrowm <- matrix(0, nrow=K,ncol=L.g)

	
	for(y in 2:L.g) {
		for(k in 1:K){
		   for(n in 0:m) {
			for(j in 0:(2^n - 1)) {
				 IntUm[k,y] <- IntUm[k,y] + Zv[2^n + j] * IntSchauderFunc(k, n, j, grid[y])
			}
		  }
	 }
  }
	
	for(k in 1:K){	
	IntBrowm[k,] <- IntUm[k,] +  Z*( (grid)^{k})/factorial(k)
   	}	
	
	IntBrowm	
}