Prove a matrix of binomial coefficients over $\mathbb{F}_p$ satisfies $A^3 = I$. (This problem is problem $1.16$ in Stanley's Enumerative Combinatorics Vol. 1).

Let $p$ be a prime, and let $A$ be the matrix $A = \left[\binom{j+k}{k} \right]_{j,k = 0}^{p-1}$, taken over the field $\mathbb{F}_p$. Show that $A^3 = I$, the identity matrix. (Note that $A$ vanishes below the main antidiagonal, i.e. $A_{jk} = 0$ if $j + k \geq p$).  Moreover, how many eigenvalues of $A$ are equal to $1$?

In the case of $p = 5$, for instance, we have $$A = \left(\begin{array}{ccccc}
1 & 1 & 1 & 1 & 1\\
1 & 2 & 3 & 4 & 0\\ 
1 & 3 & 1 & 0 & 0\\
1 & 4 & 0 & 0 & 0\\
1 & 0 & 0 & 0 & 0\\
\end{array} \right). $$
Somewhat surprisingly, we have $$A^2 = \left(\begin{array}{ccccc}
0 & 0 & 0 & 0 & 1\\
0 & 0 & 0 & 4 & 1\\ 
0 & 0 & 1 & 3 & 1\\
0 & 4 & 3 & 2 & 1\\
1 & 1 & 1 & 1 & 1\\
\end{array} \right), $$
that is the matrix flips horizontally and vertically; this can also be seen as reflecting across the main antidiagonal.  In all other cases I've tested, this held as well; the strategy I came up with for proving $A^3 = I$ is to show that this flipping occurs, and that the flipping happens again when squaring $A^2$.  This would show $A^4 = A$.  Since $A$ is clearly full rank, this would imply $A^3 = I$.  
Numerically, this flipping is captured in the statement $$ 
(A^2)_{i,j} = \sum\limits_{k = 0}^{p-1}\binom{i + k}{k} \binom{k + j}{j} \equiv \binom{2p - i - j - 2}{p - i - 1} \pmod p.
$$
I haven't been able to make any sense of this equivalence, combinatorially or otherwise; I've also made no progress on the eigenvalue part of the question, but admittedly I haven't spent as much time thinking about it.  Any help---even if it has nothing to do with the strategy that I've tried---is greatly appreciated.  

EDIT: I just figured out how to prove $A^3 = I$, and am now working on the eigenvalue problem.  Here's the proof for $A^3 = I$.  
Lemma: $A_{p - i - 1,j} \equiv (-1)^{j}\binom{i}{j} \pmod{p}$.
Proof of Lemma: We note that $A_{p - i - 1,j} = \binom{p - i + j - 1}{j}$ is a polynomial in $p$ with integer coefficients.  Thus, when viewed $\pmod{p}$, we have that it is equivalent to its constant term.  The constant term is $$ \frac{(-i + j - 1)(-i + j - 2)\cdots(-i)}{j!} = (-1)^j\frac{i(i-1)\cdots(i - (j - 1))}{j!} = (-1)^{j}\binom{i}{j} $$ as claimed.
We now have \begin{align*}
A^2_{p-i-1,p-j-1} &= \sum\limits_k A_{p-i-1,k}A_{k,p-j-1} \\
&= \sum\limits_k A_{p-i-1,k}A_{p-j-1,k} &\text{ due to symmetry of } A \\
&= \sum\limits_k (-1)^{2k}\binom{i}{k}\binom{j}{k} &\text{ by the Lemma} \\
&= \sum\limits_k \binom{i}{k}\binom{j}{j - k} \\
&= \binom{i + j}{j} &\text{ by Vandermonde's Identity} \\
&= A_{i,j}.
\end{align*}
Similarly, we have \begin{align*} 
A^4_{i,j} &= \sum\limits_k A^2_{i,k} A^2_{k,j} \\
&= \sum\limits_k A_{p - i - 1,p - k - 1} A_{p - k - 1, p - j -1}\\
&= \sum\limits_k A_{p - i - 1,p - k - 1} A_{p - j - 1, p - k -1} &\text{ by symmetry of }A \\
&= \sum\limits_k (-1)^{p - k - 1}\binom{i}{p - k - 1}(-1)^{p - k -1}\binom{j}{p - k -1} &\text{ by the Lemma} \\
&= \sum\limits_k \binom{i}{p - k - 1}\binom{j}{p - k - 1} \\
&= \sum\limits_k \binom{i}{k} \binom{j}{k} &\text{ by redefining }k \\
&= \binom{i + j}{j} &\text{ by Vandermonde}\\
&= A_{i,j}. 
\end{align*}
Thus $A^4 = A$.  Since $A$ is full rank, this implies that $A^3 = I$, as desired.
 A: Notice that ${k+j\choose j} \equiv (-1)^j{p-k-1 \choose j} \pmod p$ for $0\le k+j \le p-1$ and ${k+j\choose j} \equiv 0 \pmod p$ for  $k+j\ge p$, $j\le p-1$ (since $j!$ is coprime to $p$). Therefore,
\begin{align}
\sum_{j = 0}^{p-1} {k+j\choose j} z^j = \sum_{j = 0}^{p-k-1} (-1)^j{p-k-1 \choose j} z^j = (1-z)^{p-k-1} \tag{1} 
\end{align}
in $\mathbb F_p(z)$. 
Then for $x\in \mathbb{F}_p$, $x\neq 1$, using that $(1-x)^{p-1} = 1$ in $\mathbb F_p$,
$$
\sum_{j=0}^{p-1} {k+j\choose j}  x^j = \frac{1}{(1-x)^{k}},\ k = 0,\dots,p-1.\tag{2}$$
Consider now the vectors $e(x) = (1,x,x^2,\dots,x^{p-1})^{\top}$, $x\in \mathbb{F}_p$. For $x\notin \{0,1\}$, from (2) we have 
$$
e(x) \overset{A}{\longrightarrow} e((1-x)^{-1}) \overset{A}{\longrightarrow}  e(-(1-x) x^{-1})\overset{A}{\longrightarrow} e(x).\tag{3}
$$
Further, 
$$
e(0) \overset{A}{\longrightarrow} e(1) \overset{A}{\longrightarrow} (0,\dots,0,1)^\top \overset{A}{\longrightarrow} e(0),
$$
where the middle arrow follows from plugging $z=1$ to (1).
The vectors $e(x)$, $x\in \mathbb F_p$, are linearly independent. Therefore, $A^3 = I$. 
Concerning the eigenvalue question, there are two cases depending on whether the equation $x = (1-x)^{-1}$, equivalently, $(2x-1)^2 = -3$, is solvable in $\mathbb F_p$. (A pair of) solutions exist iff $p= 1\pmod 6$.


*

*There is no solution, so $p=6k-1$ or $p=2$. Then all elements except $0,1$ are divided into the groups of three: $x,(1-x)^{-1}, (1-x)x^{-1}$. Hence, noting that $(0,\dots,0,1)^\top = -e(0)-\dots - e(p-1)$, a vector has eigenvalue $1$ iff its decomposition in the base $\{e(0), \dots,e(p-1)\}$ has the same coefficients before $e(x), e((1-x)^{-1}), e((1-x)x^{-1})$ and zero coefficients before $e(0)$ and $e(1)$. Consequently, the space of vectors with eigenvalue $1$ has dimension $(p-2)/3 = 3k-1$ (or $0$ for $p=2$), and this is the answer.

*There are solutions $\omega$, $\omega^{-1}$, then $p=6k+1$. As above, elements except $0,1,\omega,\omega^{-1}$ are divided into threes; a vector has eigenvalue $1$ iff its decomposition in the base $\{e(0), \dots,e(p-1)\}$ has the same coefficients before $e(x), e((1-x)^{-1}), e((1-x)x^{-1})$ and zero coefficients before $e(0)$ and $e(1)$. Consequently, the space of vectors with eigenvalue $1$ has dimension $(p-4)/3 + 2= 3k+1$, and this is the answer.
