Posons $f(x)=\det(A+xJ)$. On commence par démontrer que $f$ est un polynôme de degré inférieur ou égal à 1. En effet, dans la matrice $A+xJ$, retranchons la première colonne à toutes les autres colonnes. Alors le déterminant de $A+xJ$ est égal au déterminant d'une matrice dont la première colonne est constituée par des termes du type $a_{i,1}+x$ et tous les autres coefficients sont des constantes (ne dépendent pas de x). Si on développe ce déterminant par rapport à la première colonne, on trouve que $$f(x)=\sum_{i=1}^{2n} (-1)^{i+1}(a_{i,1}+x)\det(A_i)$$ où $A_i$ est une matrice à coefficients réels et donc $\det(A_i)$ est une constante (ne dépend pas de $x$).
On démontre ensuite que, pour tout $x\in\mathbb R$, $f(x)=f(-x)$. En effet, on a \begin{align*} f(x)&=\det(A+xJ)\\ &=\det ( (A+xJ)^T)\\ &=\det ( -A+xJ)\\ &=\det ( (-1)\times (A-xJ))\\ &=(-1)^{2n}\det(A-xJ)\\ &=\det(A-xJ)\\ &=f(-x). \end{align*} On conclut car les seuls polynômes de degré inférieur ou égal à $1$ qui sont pairs sont les polynômes constants.

Notons $\Delta_p$ ce déterminant (de taille $p+1$), et prouvons que $\Delta_p=\Delta_{p-1}$ pour tout $p\geq 1$. En effet, on effectue successivement les opérations suivantes : $$L_{p+1}-L_{p}\to L_{p+1},\ L_{p}-L_{p-1}\to L_{p},\dots, L_2-L_1\to L_2.$$ Clairement, on fait apparaître des zéros dans la première colonne, excepté sur la première ligne. Pour les autres colonnes, le coefficient du déterminant à la $i$-ème ligne et à la $j$-ème colonne vaut initialement $\binom{n+i-1}{j-1}$. Après les opérations, il vaut $$\binom{n+i-1}{j-1}-\binom{n+i-2}{j-1}=\binom{n+i-2}{j-2}$$ où la dernière égalité vient de la formule du triangle de Pascal. Autrement dit, on a prouvé que $$\Delta_p=\left| \begin{array}{ccccc} 1&*&\dots&\dots&\dots\\ 0&\binom{n}0&\binom n1&\dots&\binom n{p-1}\\ 0&\binom{n+1}0&\binom{n+1}1&\dots&\binom{n+1}{p-1}\\ 0&\vdots&\vdots&&\vdots\\ 0&\binom{n+p-1}0&\binom{n+p-1}1&\dots&\binom{n+p-1}{p-1} \end{array} \right|. $$ En développant par rapport à la première colonne, on trouve bien comme annoncé que $\Delta_p=\Delta_{p-1}$. Puisque $\Delta_0=1$, le déterminant recherché est égal à 1.

On a, en développant par rapport à la dernière colonne : \begin{eqnarray*} \det(A-xI_n)&=&\left|\begin{array}{ccccc} -x&\dots&\dots&0&a_0\\ 1&\ddots&&\vdots&a_1\\ 0&\ddots&\ddots&\vdots&\vdots\\ \vdots&\ddots&\ddots&-x&\vdots\\ 0&\dots&0&1&a_{n-1}-x \end{array}\right|\\ &=&(-x)^{n-1}(a_{n-1}-x)+\sum_{k=0}^{n-2}(-1)^{n+k-1}a_k\Delta_k \end{eqnarray*} où $\Delta_k$ est le déterminant suivant ($-x$ est répété $k$ fois sur la diagonale) : $$\Delta_k= \left| \begin{array}{cccccccc} -x&0&\dots&0&*&\dots&\dots&*\\ *&\ddots&\ddots&\vdots&\vdots&&&\vdots\\ \vdots&\ddots&\ddots&\vdots&\vdots&&&\vdots\\ *&\vdots&\vdots&-x&*&\dots&\dots&*\\ 0&\dots&\dots&0&1&*&\dots&*\\ \vdots& & & \vdots&0&\ddots&\ddots&\vdots\\ \vdots& & & \vdots&\vdots&\ddots&\ddots&*\\ 0&\dots&\dots&0 & 0&\dots&0&1 \end{array}\right|.$$ Ce déterminant $\Delta_k$ se calcule par blocs (on a une matrice triangulaire supérieure par blocs). De plus, chacun des blocs diagonaux est lui-même triangulaire (inférieur pour le premier, supérieur pour le second). On en déduit que $$\Delta_k=(-x)^k 1^{n-1-k}$$ et donc que \begin{eqnarray*} \det(A-xI_n)&=&(-x)^{n-1}(a_{n-1}-x)+\sum_{k=0}^{n-2}(-1)^{n+k-1}(-1)^k a_kx^k\\ &=&(-1)^n\left(x^n-\sum_{k=0}^{n-1}a_kx^k\right). \end{eqnarray*}

$\mnr$ est la somme directe du sous-espace vectoriel des matrices symétriques et des matrices antisymétriques. Soit $(A_1,\dots,A_p)$ et $(B_1,\dots,B_q)$ une base respective de l'espace vectoriel des matrices symétriques et antisymétriques. $(A_1,\dots,A_p,B_1,\dots,B_q)$ forme une base de $\mnr$, et il suffit de calculer le déterminant dans cette base. Mais $\phi(A_i)=A_i$ tandis que $\phi(B_j)=-B_j$. On a donc $\det(\phi)=(-1)^q$. Il suffit ensuite de se souvenir que $p=\frac{n(n+1)}{2}$, ou $q=\frac{n(n-1)}{2}$.

Trouvons d'abord une condition nécessaire. Puisque $M\in \mathcal M_n(\mathbb Z)$, $\det(M)\in \mathbb Z$. D'autre part, si $M^{-1}\in \mathcal M_n(\mathbb Z)$, son déterminant est un entier et donc $\det(M^{-1})=\frac{1}{\det M}\in\mathbb Z$. Ceci entraine que $\det(M)=\pm 1$.
Réciproquement, si $\det(M)=\pm 1$, alors $M$ est inversible. De plus, toutes les entrées de sa comatrice, qui sont obtenues comme déterminants de sous-matrices de $M$, sont des entiers. De la formule de Cramer, on déduit que $M^{-1}=\frac{1}{\det(M)}\ ^t\!\textrm{comat}(M)$ est une matrice à coefficients dans $\mathbb Z$.

Calculons le déterminant de cette famille (de $(n+1)$ vecteurs dans un espace de dimension $n+1$) par rapport à la base canonique $(1,X,\dots,X^n)$. On a $$(X-z_i)^n=\sum_{j=0}^n \binom{n}{j}(-1)^{n-j}z_i^{n-j}X^j.$$ Le déterminant recherché est donc \begin{eqnarray*} \Delta&=&\left| \begin{array}{cccc} \binom{n}{0}(-z_0)^n&\binom{n}{0}(-z_1)^n&\dots&\binom{n}{0}(-z_n)^n\\ \binom{n}{1}(-z_0)^{n-1}&\binom{n}{1}(-z_1)^{n-1}&\dots&\binom{n}{1}(-z_n)^{n-1}\\ \vdots&\vdots&\vdots&\vdots\\ \binom{n}{n}&\binom{n}{n}&\dots&\binom{n}{n}\\ \end{array}\right|\\ &=&\binom{n}{0}\binom{n}{1}\dots\binom{n}{n}\left| \begin{array}{cccc} (-z_0)^n&(-z_1)^n&\dots&(-z_n)^n\\ (-z_0)^{n-1}&(-z_1)^{n-1}&\dots&(-z_n)^{n-1}\\ \vdots&\vdots&\vdots&\vdots\\ 1&1&\dots&1\\ \end{array}\right|. \end{eqnarray*} On reconnaît un déterminant de Vandermonde, qui est non-nul puisque les $z_i$ sont supposés tous distincts. La famille considérée est effectivement une base de $\mathbb C_n[X]$.