function(k, n0, pstar, conf, rep)
{
	# last update 2/23/00
	# function to estimate critical values for comparisons with a standard
	# when mu0 is unknown, variances are unknown and unequal and
	# systems are independent
	# returns estimate of (g,h) with conf*100% upper confidence bound 
	# on each (confidence on h is actually 1 - (1-conf)/2)
	#
	# k = number of treatments excluding the standard
	# n0 = first-stage sample size
	# pstar = 1-alpha value (PCS); single value only in this version
	# rep = number of replications to use for estimate
	# conf = confidence level on upper bounds
	#
	# first get h 
	set.seed(13.)
	conf <- 1 - (1 - conf)/2
	Z <- matrix(rnorm(k * rep), ncol = k)
	Y <- matrix(rchisq(k * rep, n0 - 1.), ncol = k)
	C <- matrix(rchisq(rep, n0 - 1.), ncol = 1.)
	Cmat <- matrix(C, ncol = k, nrow = rep, byrow = F)
	denom <- sqrt((n0 - 1.) * (1./Y + 1./Cmat))
	H <- sort(apply(Z * denom, 1., max))
	Hstar <- quantile(H, pstar)
	upper <- ceiling(pstar * rep + qnorm(conf) * sqrt(pstar * (1. - pstar) *
		rep) + 0.5)
	#
	# then use h to get g
	Hupper <- H[upper]
	Z <- matrix(rnorm((k - 1.) * rep), ncol = k - 1.)
	Zk <- matrix(rnorm(rep), ncol = 1.)
	Y <- matrix(rchisq((k - 1.) * rep, n0 - 1.), ncol = k - 1.)
	C0 <- matrix(rchisq(rep, n0 - 1.), ncol = 1.)
	Ck <- matrix(rchisq(rep, n0 - 1.), ncol = 1.)
	Cmat <- matrix(Ck, ncol = k - 1., nrow = rep, byrow = F)
	denom <- sqrt((n0 - 1.) * (1./Y + 1./Cmat))
	G <- sort(apply(cbind(Zk * sqrt((n0 - 1.) * (1./C0 + 1./Ck)) + Hupper,
		Z * denom), 1., max))
	Gstar <- quantile(G, pstar)
	upper <- ceiling(pstar * rep + qnorm(conf) * sqrt(pstar * (1. - pstar) *
		rep) + 0.5)
	Gupper <- G[upper]
	result <- rbind(c(Hstar, Gstar), c(Hupper, Gupper))
	result
}