考虑把未知的位置替换成变量，如样例一的 \begin{bmatrix}0&1\\0&2\end{bmatrix} 就可以变成 \begin{bmatrix}x&1\\y&2\end{bmatrix}，而它的三次方就是 \begin{bmatrix}x^3+2xy+2y&x^2+2x+y+4\\x^2y+y^2+2xy+4y&xy+4y+8\end{bmatrix}。关于矩阵 p 次幂的计算，可以由 A^2_{i,j}=\sum _{k}A_{i,k}\times A_{k,j} 类比出 A^{p}_{i,j}=\prod _{k_{1,2,\dots,p-1}} A_{i,k_1}\times A_{k_1,k_2}\times\dots\times A_{k_{p-1},j}。
注意到，\sum_{i=1}^{p-1}i^{k}\bmod p 只有当 k=p-1 时才不为零（并且刚好是 p-1）。
考虑证明。找到 p 的原根 r，则原式等价于 \sum_{i=0}^{p-2}r^{ik}=\sum_{i=0}^{p-2}(r^k)^i。当 k=p-1 时 r^k\equiv 1，则原式 \equiv p-1。反之，可以运用等比数列求和公式得出原式 \equiv \frac{r^{k(p-1)}-1}{r^k-1}，根据费马小定理，r^{k(p-1)}\equiv 1\pmod p，则分子是 0，整个分数就是 0，得证。
有了这个强有力的结论，正解就呼之欲出了。
先把没有未知数的贡献用矩阵快速幂算出来。可以发现矩阵的某个位置的某一项中所有未知数的指数都是 p-1 时这一项的贡献才不为零。p=2 的情况当作 corner case，很好处理。其余情况只有当 k_{1,2,\dots,p-1}=i 或者 k_{1,2,\dots,p-1}=j 时才可能会出现未知数的指数是 p-1 的情况，只有两种情况有贡献，这部分很好做。当 p=3 时，k_1=j,k_2=i，A_{i,j} 的指数是 2，也符合条件，要把它看成一个 corner case。贡献记得乘上选其它未知数的方案数。
时间复杂度是 O(n^3\log p)，瓶颈在于矩阵快速幂。
代码。