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 & 1 & 0 & \cdots & 0 \\ 0 & 0 & 1 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & 1 \\ a_0 & a_1 & a_2 & \cdots & a_{k-1} \end{pmatrix}.$$ And $U^n$ can be computed efficiently by standard exponentiation algorithms.
• Addition: it is common to use a primitive polynomial for the LFSR; in that case, the LFSR repeats after $2^k-1$ steps, and it follows that the $n^\text{th}$ state is also the $(n\bmod(2^k-1))^\text{th}$ state, and we can compute $U^{(n\bmod(2^k-1))}$. In the context of the question, it will help when the polynomial is of degree less than about 256. – fgrieu Aug 4 '16 at 1:36