Let $(X_i)_{i\in \mathbb{N}_0}$ be a sequence of iid standard Normal distributed numbers. In dimension one define
$$B_t^H = X_0 t + \sum_{k=0}^\infty X_{2^k+\lfloor 2^k t\rfloor}\frac{|2^k t-\lfloor 2^k t + \tfrac{1}{2}\rfloor|}{(2^k)^H}, \ t\in [0,1].$$
For $H=\tfrac{1}{2}$ this is the standard Brownian motion and its paths are continuous but nowhere differentiable (almost surely).
When $t$ is a dyadic rational ($t \in \{\tfrac{i}{2^m} : i,m \in \mathbb{N}_0\}$) the above series has only $m$ summands.
That means, sampling at these positions can be done in linear time.
As this is a dense subset and numbers in a computer are represented anyway by binaries this should do the trick.
In dimension two take $(X_{i,j})_{i,j\in \mathbb{N}_0}$ and define
$$B_{s,t}^H = X_{0,0} \tfrac{s+t}{2} + \sum_{k=0}^\infty \frac{(\tilde X +\dot X) a(2^k s,2^k t)+(\tilde X -\dot X)b(2^k s,2^k t)}{(2^k)^H}$$
for $ \ (s,t)\in [0,1]^2$, where
$$\tilde X = X_{2^k+\lfloor 2^k s\rfloor,2^k+\lfloor 2^k t\rfloor},\\
\dot X = X_{2^k+(\lfloor 2^k s\rfloor+ c(2^k s,2^k t)) \text{mod}\, 2^k,2^k+ (\lfloor 2^k t\rfloor+ c(2^k t,2^k s)) \text{mod}\, 2^k},\\
c(x,y)=\tfrac{\text{sign}(\{x\}+\{y\}-1)+\text{sign}(\{x\}-\{y\})}{2},\\
a(x,y)=1 -||\{x\}+\{y\}-1|-|\{x\}-\{y\}||,\\
b(x,y)=1- |\{x\}+\{y\}-1|-|\{x\}-\{y\}| $$
and $\{\cdot\}$ denotes the fractional part.