nrep <- 100000
ris <- matrix(0,nrow=nrep,ncol=2)
for (i in 1:nrep)
{
   s1 <- 0; n1 <- 0
   s2 <- 0; n2 <- 0
   while (s1<=1 | s2<=1) {
      u1 <- runif(1)
      u2 <- 1 - u1
      if (s1<=1) {
         s1 <- s1 + u1
         n1 <- n1 + 1
      }
      if (s2<=1) {
         s2 <- s2 + u2
         n2 <- n2 + 1
      }
   }
   ris[i,1] <- n1
   ris[i,2] <- n2
}
out <- apply(ris,1,mean)

xm <- mean(out)
se <- sqrt(var(out)/nrep)
cil <- xm - 1.96*se
cih <- xm + 1.96*se
cat(xm,cil,cih,(cih-cil),"\n")

sum(apply(ris,1,max))