k            <– 10000
perms        <- replicate(k, sample(c(g1, g2)))
my.mean.difs <- apply(perms[1:length(g1), ], 2, mean) -
                apply(perms[(length(g1)+1):(length(g1)+length(g2)), ], 2, mean)

p            <- sum(my.mean.difs < (mean(g1) - mean(g2))) / k