Values for integral of square of standard Brownian process I am trying to generate values in a table for the following function:
$$
W = \int_0^1 [B(t)]^2 dt 
$$
Where $B(t)$ is a standard Brownian motion.
Example: $W_{0.05} = 1.656$, $W_{0.025} = 2.135$.
This is taken from a table C.6 in 'Survival Analysis' by Klein & Moeschberger. (They don't give a reference as to how values were computed).
Am using R for this:
BM2 <- function (x) (rnorm(n=1, mean=0, sd=x))^2

But of course the standard integrate(BM2, lower=0, upper=1)
doesn't work for this: evaluation of function gave a result of wrong length. 
Looked at this similar question is similar but it's clearly not as simple as:
$$
W_t = \frac{t^3}{3}
$$
I'm guessing the function is integrated by Ito's method but can't seem to find any worked examples similar to the formula above.
Any help or pointers (particularly in R which doesn't seem to have an existing function for this) would be greatly appreciated, particularly if the numeric values given above could be reproduced.
UPDATE 1/29:
Thanks @ThePawn. Question corrected. I like the resampling method and it does produce reasonable approximations. At first I thought this was due to the sample size or loop iterations being too small. Function doesn't appear to 'vecorize' easily in R so I rewrote it in C++. Using samples >10,000 or dt <0.0001 doesn't tend to improve things much (which I should probably have guessed) and gets rather slow. The method doesn't seem to give precision beyond 1 significant figure (although this is certainly good enough for a practical application):
require(rcpp)
require(inline)

src <- "
arma::colvec b1 = Rcpp::as<arma::vec>(B);
arma::colvec w1 = Rcpp::as<arma::vec>(B);
double nb1 = b1.size();
double dt1 = as<double>(dt);
double dtrt1 = pow(dt1, 0.5); // square root of dt
int lt1 = as<int>(lt);

for (int i=0; i<lt1; i++){
RNGScope scope; // reset random number generator each time
arma::colvec rn1 = rnorm(nb1, 0, dtrt1);
b1 += rn1;
w1 += pow(b1, 2) * dt1;
}

return Rcpp::wrap(w1);
"

f1 <- inline::cxxfunction(signature(B="numeric",lt="integer",dt="double"),
 src, plugin = "RcppArmadillo")

n <- 1e4 # sample size
B <- rep(0, n) # hold result
dt <- 1e-5
T <- seq(0, 1, dt)
lt <- length(T)

B <- f1(B,lt,dt)

sum(B >2.135)/n # approx = 0.025
sum(B >1.656)/n # approx = 0.05

@probabilityislogic - not sure what limiting process I'm using, open to suggest there.  Will try implementing the formula you suggest but looks like it will be another approximate solution...

Am still wondering if there's an exact method for this - or do you think the values were generated from a resampling approach like the above?
 A: It makes no sense to index $W$ with respect to $t$ because $t$ is the integration variable. To clarify, you could express $W$ as a function of $\omega$ which is going to represent the state in $\Omega$ (have a look at wikipedia for a rigourous definition of a stochastic process). 
Because $W$ is a random variable, it would make sense to compute its cumulative distribution function $P(W\le u)$.
A basic Monte Carlo method would consist in discretizing the time space using a time step $\Delta t$, simulate a Brownian motion at those times and integrate. 
For example (using R):
nmc = 1000;
dt = 0.01;
T = seq(0, 1, dt);
W = rep(0, nmc);
B = rep(0, nmc);
for (j in 1:length(T))
{
    B = B + rnorm(nmc, 0, sqrt(dt)); # brownian motion at time delta_t * j
    W = W + B^2 * dt;
}
cdfW = function(u) { return(mean(W < u)) }
mean(W) # 0.4983

Where I checked at the end that the mean is close to the theoretical one that is
$$E[W] = \int_0^1 E[B_t^2]dt = \int_0^1 tdt = \frac{1}{2}$$
You will need to study the behaviour of the empirical cumulative distribution function as $dt \rightarrow 0$ if you want to compute an accurate approximation.
For more properties concerning the integral of a squared brownian motion, I found this link.
I think another way of doing would consist in calculating all the moment of $W$, deduce its characteristic function and invert it using a Fourier transform for example. It seems to me that the moments can be computed in closed form:
$$E[W^k] = \int_{t_1 = 0}^{1} \ldots \int_{t_k = 0}^{1}E[B_{t_1}^2 \ldots B_{t_k}^2]dt_1 \ldots dt_k$$
The terms $E[B_{t_1}^2 \ldots B_{t_k}^2]$ can be calculated by ranking the different $t_i$ and using the properties of the Brownian motion. 
For example, with $k=2$ and $s \le t$, $$E[B_t^2B_s^2] = E[B_s^4] + E[B_s^2(B_t-B_s)^2]+2E[B_s^3(B_t-B_s)]$$ That is $$E[B_t^2B_s^2] = 3s^2+(t-s)s$$.
If it is too cumbersome to calculates all moments, you can still use a few to improve the convergence of your empirical distribution function by using variance reduction techniques.
A: P. Erdös and M. Kac gave the exact expression for the distribution function of W on pp. 292-293 (see III) of the paper, reproduced here:

Erdös, P.; Kac, M. On certain limit theorems of the theory of probability. Bull. Amer. Math. Soc. 52, (1946). 292–302. 
The formula is far from simple, so I do not know if it is useful
for you.
A: One way to figure this out is to use the following limit formula
$$ W(s)=\int_{0}^{s}[B(t)]^2dt= \lim_{n\to\infty}\Delta t\sum_{i=0}^{n-1}[ B(i\Delta t)]^{2}$$
where $\Delta t=\frac{s}{n}$ (this holds because squared brownian motion is continuous wrt time).  You use the properties of the sum , which suggest a generalised chi-square distribution for $W(s)$ - has a representation as a weighted sum of non-central  chi squares, and would be difficult to work with.
