xing <- function() {
  # Perfect simulation of the time taken for a BM to hit +/- 1 starting at 0
  # Zaeem Burq 23.3.06, Owen Jones 24.3.06
  accepted <- FALSE
  while (!accepted) {
    X <- rgamma(1, 1.088870, 0.810570)
    Y <- runif(1)*1.243707*dgamma(X, 1.088870, 0.810570)
    sqrt2piX3 <- sqrt(2*pi*X^3)
    N <- max(ceiling(0.275*X), 3)
    K <- (1+2*(-N:N))
    fN0 <- sum((-1)^(-N:N)*K*exp(-K^2/(2*X)))/sqrt2piX3
    N <- N + 1
    fN1 <- fN0 + (-1)^N*((1-2*N)*exp(-(1-2*N)^2/(2*X))
        + (1+2*N)*exp(-(1+2*N)^2/(2*X)))/sqrt2piX3;
    while (sign((Y - fN0)*(Y - fN1)) == -1) {
        fN0 <- fN1
        N <- N + 1
        fN1 <- fN0 + (-1)^N*((1-2*N)*exp(-(1-2*N)^2/(2*X))
            + (1+2*N)*exp(-(1+2*N)^2/(2*X)))/sqrt2piX3;
    }
    if (Y <= fN1) accepted <- TRUE
  }
  return(X)
}

BM_xing <- function(n, scale = 1) {
  # n co-ordinates of BM observed when it crosses levels scale*Z
  # Zaeem Burq 23.3.06, Owen Jones 24.3.06
  t <- rep(0, n+1)
  for (i in 1:n) t[i+1] <- t[i] + xing()
  Wt <- c(0, cumsum(sign(runif(n, -1, 1))))
  if (scale != 1) {
      t <- scale^2*t
      Wt <- scale*Wt
  }
  return(matrix(c(t, Wt), ncol = 2))
}