the case, on the contrary. Recall that the first-order algorithm $U_1(\tau)$ is all we need to program: As explained in Sec. \[sec3\], more accurate schemes can be implemented without extra programming. Let us consider the algorithm for the case of 1D. For 2D and 3D we simply repeat the steps described below two, respectively, three times. We have $$\begin{aligned}
U_1(\tau)&=&e^{\tau H_{1}} e^{\tau H_{2}}
\nonumber \\
&=&
\left\{\mathop{{\prod}'}_{i=1}^{n}
\exp\left[\beta_{i+1,i}\left(
\be^{{\phantom{T}}}_{i}\be^T_{i+1}-
\be^{{\phantom{T}}}_{i+1}\be^T_{i}
\right)\right]\right\}
\left\{\mathop{{\prod}'}_{i=1}^{n}
\exp\left[\beta_{i+1,i+2}\left(
\be^{{\phantom{T}}}_{i+1}\be^T_{i+2}-
\be^{{\phantom{T}}}_{i+2}\be^T_{i+1}
\right)\right]\right\}.
\label{U1D}\end{aligned}$$ where we used the block-diagonal structure of $H_1$ and $H_2$ (see Eqs.(\[h1m\]) and (\[h2m\])) to obtain an exact expression for $U_1(\tau)$ in terms of an ordered product of matrix exponentials. Each of these matrix exponentials only operates on a pair of elements of $\Psi(t)$ and leaves other elements intact. The indices of each of these pairs are given by the subscripts of $\be$ and $\be^T$. From Eq.(\[U1D\]) it is then clear what the program should do: Make loops over $i$ with stride 2. For each $i$ pick a pair of elements from $\Psi(t)$ according to the subscripts of $\be$ and $\be^T$, compute (or fetch from memory) the elements of the plane rotation (see Eq. (\[twobytwo\])), perform the plane rotation, i.e. multiply the $2\times2$ 