function(x, y, d, m)
{
	n1 <- length(x)
	n2 <- length(y)
	a1 <- mean(x)
	a2 <- mean(y)
	sq1 <- var(x)
	sq2 <- var(y)
	tst <- (a1 - a2 - d)/sqrt(sq1/n1 + 
		sq2/n2)
	da <- a1 - a2
	tsts <- c(1:m)
	for(i in 1:m) {
		xs <- sample(x, n1, replace
			 = T)
		ys <- sample(y, n2, replace
			 = T)
		as1 <- mean(xs)
		as2 <- mean(ys)
		vs1 <- var(xs)
		vs2 <- var(ys)
		tsts[i] <- (as1 - as2 - da)/
			sqrt(vs1/n1 + vs2/n2
			)
	}
	prop <- 0 * jay(m, 1)
	for(i in 1:m) {
		if(abs(tsts[i]) >= abs(tst)
			)
			prop[i] <- 1
	}
	pvalue <- sum(prop)/m
	res <- list(da, tst, pvalue)
	names(res) <- c("Diff. Means", 
		"t-Diff.", "P-Value")
	res
}
