chisqhdr <- function(p,n){
  theta0 <- function(d) d/tanh(d/(n-2))
  pr <- function(d) pchisq(theta0(d)+d,n)-pchisq(theta0(d)-d,n)
  d0 <- (dchisq(qchisq(p/2,n),n)+dchisq(qchisq(1-p/2,n),n))/2
  d1 <- d0
  while (pr(d1) > p) d1 <- d0/2
  while (pr(d0) < p) d0 <- 2*d0
  maxiter <- 100
  iter <- 0
  pp <- 0
  while ((abs(pp-p)>0.0001)&&(iter<maxiter)){
    iter <- iter+1
    dm <- (d0+d1)/2
    if ((pr(d1)-p)*(p-pr(dm)) > 0) d0 <- dm
    if ((pr(dm)-p)*(p-pr(d0)) > 0) d1 <- dm
  }
low <- theta0(d0)-d0  
high <- theta0(d0)+d0  
return(c(low,high))
}