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.
On
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)!$.
$$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)$$
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.
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.
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}$.
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).