r <- 1
v <- c(1,1)
inA <- function(x) sum((x-v)^2)<r^2
N <- 100000; K <- 1000; T <- 0
for (k in 1:K){
	X <- matrix(rnorm(2*N),N)
	C <- sum(as.integer(apply(X,1,inA)))
	T <- T+C
}
cat(sprintf("P(A)=%g\n",T/N/K))