Algorithm or closed form solution that lets me directly compute the $n^{\text{th}}$ state of an LFSR? I have a large linear feedback shift register with a known state and taps.

Now I want to compute the state of the LFSR after $n$ iterations.

The problem: $n$ is in the order of $2^{256}$ so just brute forcing is not an option. Is there an algorithm or closed form solution that lets me directly compute the $n^{\text{th}}$ state?
If the initial state is $b_0,b_1,\dots,b_{k-1}$ and the recurrence relation is $b_k = \sum_{i = 0}^{k-1} a_ib_i$, then in linear-algebraic terms we have
$$ \begin{pmatrix}
b_1 \\ b_2 \\ \vdots \\ b_k
\end{pmatrix} = U \begin{pmatrix}
b_0 \\ b_1 \\ \vdots \\ b_{k-1}
\end{pmatrix}, $$
and more generally
$$ \begin{pmatrix}
b_{n} \\ b_{n+1} \\ \vdots \\ b_{n+k-1}
\end{pmatrix} = U^n \begin{pmatrix}
b_0 \\ b_1 \\ \vdots \\ b_{k-1}
\end{pmatrix}, $$
where of course
$$ U = \begin{pmatrix}
0 &amp; 1 &amp; 0 &amp; \cdots &amp; 0 \\
0 &amp; 0 &amp; 1 &amp; \cdots &amp; 0 \\
\vdots &amp; \vdots &amp; \vdots &amp; \ddots &amp; \vdots \\
0 &amp; 0 &amp; 0 &amp; \cdots &amp; 1 \\
a_0 &amp; a_1 &amp; a_2 &amp; \cdots &amp; a_{k-1}
\end{pmatrix}. $$
And $U^n$ can be computed efficiently by standard exponentiation algorithms.