Let's say we have two functions $h(s)$ and $g(s)$. We can easily simulate a stochastic integral, e.g.
$$t \mapsto \int\_0^t h(s) dB(s) \sim \mathcal{N}\bigg(0, \int\_0^t h(s)^2 ds \bigg). $$
What is the conditional distribution of stochastic integral of $g(s)$ with respect to $B(s)$ then?
$$ t \mapsto \int\_0^t g(s) dB(s) \bigg | \int\_0^t h(s) dB(s) \sim ?$$
I have an intuition that these two integrals might have a joint normal distribution with a covariance equal to $\int\_0^t h(s)g(s)ds$, but I need a rigorous derivation. Generally, I want to find a way to simulate a bunch of $n$ integrals with respect to the same Brownian motion:
$$ t \mapsto \bigg( \int\_0^t h\_1(s) dB(s), \int\_0^t h\_2(s) dB(s) , \dots , \int\_0^t h\_n(s) dB(s) \bigg)^T $$



Conditional stochastic integration
(From the context of your post, it is apparent that $B$ is assumed to be a (say standard) Brownian motion.)
The joint normal distribution follows from the way the stochastic integral is defined. However, if you are already convinced that
\begin{equation}
 I\_t(h):=\int\_0^t h(s) dB(s) \sim N\bigg(0, \int\_0^t h(s)^2 ds \bigg)
\end{equation}
for all $h\in L^2([0,t])$, then you have
\begin{equation}
 \sum\_1^n r\_i I\_t(g\_i)=I\_t\Big(\sum\_1^n r\_i g\_i\Big) \sim N\bigg(0, \int\_0^t \Big(\sum\_1^n r\_i g\_i(s)\Big)^2 ds \bigg) \tag{\*}
\end{equation}
for all real $r\_i$'s and all functions $g\_i\in L^2([0,t])$. So, all linear combinations of the $I\_t(g\_i)$'s are Gaussian and hence the $I\_t(g\_i)$'s are jointly Gaussian (think, e.g., of the joint characteristic function). Formula (\*) for $n=1,2$ also yields the covariances:
\begin{align\*}
 Var(I\_t(g\_1)+I\_t(g\_2))&=\int\_0^t (g\_1(s)+g\_2(s))^2 ds \\ 
 &=\int\_0^t g\_1(s)^2 ds+\int\_0^t g\_2(s)^2 ds
 +2\int\_0^t g\_1(s)g\_2(s)\,ds \\ 
 &=Var(I\_t(g\_1))+Var(I\_t(g\_2))+2\int\_0^t g\_1(s)g\_2(s)\,ds, 
\end{align\*}
so that
$Cov(I\_t(g\_1),I\_t(g\_2))=\int\_0^t g\_1(s)g\_2(s)\,ds$.


**Added in response to a comment by Dr\_Zaszuś:**
Thus, for any $g,h$ in $L^2([0,t])$, the pair $(I\_t(g),I\_t(h))$ has the bivariate normal distribution $N(\mu\_g,\mu\_h,\sigma^2\_g,\sigma^2\_h,\rho)$, where $\mu\_g=\mu\_h=0$, $\sigma\_g=\|g\|$, $\sigma\_h=\|h\|$, and $\rho=\frac{g\cdot h}{\|g\|\|h\|}$, where $g\cdot h:=\int\_0^t g(s)h(s)\,ds$ and $\|f\|:=\sqrt{f\cdot f}$. So, the conditional distribution of $I\_t(g)$ given $I\_t(h)=y$ is the normal distribution $N(\mu\_{g;y},(1-\rho^2)\sigma\_g^2)$ for any real $y$, where $\mu\_{g;y}:=\rho\sigma\_g y/\sigma\_h$.


