function(x, y, m)
{
	n1 <- length(x)
	n2 <- length(y)
	nc <- n1 + n2
	a1 <- c(1:n1)
	a2 <- c((n1 + 1):nc)
	diff <- mean(x) - mean(y)
	ax <- c(1:nc)
	indx <- c(1:nc)
	for(i in 1:n1) {
		ax[i] <- x[i]
	}
	for(i in (n1 + 1):nc) {
		ax[i] <- y[i - n1]
	}
	ds <- c(1:m)
	for(i in 1:m) {
		ac <- sample(indx, nc, 
			replace = F)
		ac1 <- ac[1:n1]
		ac2 <- ac[(n1 + 1):nc]
		xs <- ax[ac1]
		ys <- ax[ac2]
		ds[i] <- mean(xs) - mean(ys)
	}
	prop <- 0 * jay(m, 1)
	for(i in 1:m) {
		if(abs(ds[i]) >= abs(diff))
			prop[i] <- 1
	}
	pvalue <- mean(prop)
	res <- list(diff, pvalue)
	names(res) <- c("DIFF", "P-Value")
	res
}
