Is it possible to compute factorials by converting to matrix multiplications? An $n$-th term of the Fibonacci sequence can be computed by a nice trick by converting the recurrence relation in a matrix form. Then we compute $M^n$ in $O(\log n)$ steps using exponentiation by squaring.
Would it be possible to use such a trick to compute factorials? If not, can it be proved? I figured out how to compute any polynomial in $n$ using this approach, but for factorials I wasn't able to express the factorial's recurrence relation as a linear transformation.
 A: The Fibonacci recurrence relation:
$$F_{n+1} = F_n + F_{n-1}$$
is linear. By that I mean it gives the next term as a linear combination of previous terms. Other recurrences having this feature are:
$$A_{n+1} = 3A_n + A_{n-1}$$
$$B_{n+1} = -4B_n$$
$$C_{n+1} = C_n + C_{n-1} + C_{n-2}$$
so any of those could be represented by matrices. Here they are!
$$
\begin{aligned}
M_A &= \begin{bmatrix}3 & 1\\1& 0\end{bmatrix}\\
M_B &= \begin{bmatrix}-4\end{bmatrix}
\end{aligned}
$$
And the last one is a challenge to the reader.
The Fibonacci matrix
$$\begin{bmatrix}1&1\\0&1\end{bmatrix}$$
maps a vector
$$\begin{bmatrix}F_n\\F_{n-1}\end{bmatrix} \rightarrow \begin{bmatrix}F_{n+1}\\F_n\end{bmatrix}$$
and these work analogously.
But the factorial recurrence relation:
$$n! = n(n-1!)$$
is not a linear combination (with constant coefficients) of previous arguments, so we can't represent it as a matrix. That is, we can't represent it by the same matrix, no matter what $n$ we're trying to compute. When the coefficents are constant, we can, which lets us do exponentiation by squaring, etc.
Matrices are a convenient shorthand for representing linear functions on vectors. (They're not only a shorthand, but often I find it useful to think of them this way).
A: Well, this is probably not the answer you are looking for, but it might be interesting.
Let $D_k$ be the $k \times k$ matrix with zeros everywhere except on the superdiagonal, where it has the values $1, 2, \dots, k - 1$. Then the factorial $(k-1)!$ is equal to the first element in the vector
$$D_k^{k-1} e_k$$
where $e_k$ is the unit vector with a one at position $k$. Equivalently, the element at position $(1,k)$ in $D_k^{k-1}$ is $(k-1)!$.
Example
$$D_5 = \left(
\begin{array}{ccccc}
 0 & 1 & 0 & 0 & 0 \\
 0 & 0 & 2 & 0 & 0 \\
 0 & 0 & 0 & 3 & 0 \\
 0 & 0 & 0 & 0 & 4 \\
 0 & 0 & 0 & 0 & 0
\end{array}
\right)$$
$$D_5^4 = \left(
\begin{array}{ccccc}
 0 & 0 & 0 & 0 & 24 \\
 0 & 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0
\end{array}
\right)$$
Explanation
The matrix $D_k$ acts on the space of polynomials of degree $\leq k$ as differentiation. Then one can use the following characterization of the factorial:
$$n!= \frac{d^n}{dx^n} x^n$$
and get the result above.
Observations
The matrix $D_k$ can also be used to generate an upper triangular Pascal matrix $U_k$ through $U_k = e^{D_k}$, using the matrix exponential.
A closer look at $D_k$
When you put $D_k$ on Jordan normal form. For example with $D_5$ above you get $D_5 = S^{-1}JS$:
$$D_5 = \underbrace{\left(
\begin{array}{ccccc}
 1 & 0 & 0 & 0 & 0 \\
 0 & 1 & 0 & 0 & 0 \\
 0 & 0 & 2 & 0 & 0 \\
 0 & 0 & 0 & 6 & 0 \\
 0 & 0 & 0 & 0 & 24
\end{array}
\right)}_{=S^{-1}}
\underbrace{\left(
\begin{array}{ccccc}
 0 & 1 & 0 & 0 & 0 \\
 0 & 0 & 1 & 0 & 0 \\
 0 & 0 & 0 & 1 & 0 \\
 0 & 0 & 0 & 0 & 1 \\
 0 & 0 & 0 & 0 & 0
\end{array}
\right)}_{=J}
\underbrace{
\left(
\begin{array}{ccccc}
 1 & 0 & 0 & 0 & 0 \\
 0 & 1 & 0 & 0 & 0 \\
 0 & 0 & \frac{1}{2} & 0 & 0 \\
 0 & 0 & 0 & \frac{1}{6} & 0 \\
 0 & 0 & 0 & 0 & \frac{1}{24}
\end{array}
\right)}_{=S}
$$
where $J$ is one big Jordan block (which is not surprising since $D_k$ is nilpotent), and $S^{-1}$ is a diagonal matrix with factorials on the diagonal!
Investigating this a little but further reveals that this is not so strange. $D_k$ has only one eigenvalue, namely 0, with algebraic multiplity $k$ and geometric multiplicity 1. The eigenvector belonging to this eigenvalue is $e_1$.
When we put something on Jordan normal form, we use generalized eigenvectors, which are vectors $v$ that satisfy $(D_k - \lambda I)^n v = 0$ for some $k$. In our case $\lambda = 0$, so we just look at $D_k^n v = 0$. This gives us that $e_k$ is a generalized eigenvector. Remember our original formulation?
$$D_k^{k-1}e_k = (k-1)! e_1$$
which explains why we have factorials in $S^{-1}$.
