  Consider the fractional Brownian Motion (fBM) $B^H=\{B^H(t): t \in [0,1] \}$
with Hurst index $H\in (0,1)$. We construct a probability space supporting both
$B^H$ and a fully simulatable process $\hat B_{\epsilon}^H $ such that
$$\sup_{t\in [0,1]}|B^H(t)-\hat B_{\epsilon}^H(t)| \le \epsilon$$ with
probability one for any user specified error parameter $\epsilon>0$. When
$H>1/2$, we further enhance our error guarantee to the $\alpha$-H\"older norm
for any $\alpha \in (1/2,H)$. This enables us to extend our algorithm to the
simulation of fBM driven stochastic differential equations $Y=\{Y(t):t
\in[0,1]\}$. Under mild regularity conditions on the drift and diffusion
coefficients of $Y$, we construct a probability space supporting both $Y$ and a
fully simulatable process $\hat Y_{\epsilon}$ such that $$\sup_{t\in
[0,1]}|Y(t)-\hat Y_{\epsilon}^H(t)| \le \epsilon$$ with probability one. Our
algorithms enjoy the tolerance-enforcement feature, i.e., the error bounds can
be updated sequentially. Thus, the algorithms can be readily combined with
other simulation techniques like multilevel Monte Carlo to estimate expectation
of functionals of fBMs efficiently.
