Square roots of integral matrices Define $n\times n$ matrices $L , M$  as follows: $\,$ L$_{ij}$= $1$ if $\,$i + j $\geq$n + 1
and is zero otherwise; $\,$ $\qquad$ $\qquad$ $\quad$M$_{ij}$= $\min(i,j)$ for $1 \leq i , j \leq n$ .  It is straightforward to check that $M = L^2$ .  For example, when $n=3$ we have
$\qquad$ $\qquad$ $\qquad$$\begin{bmatrix}1 & 1 & 1\\1 & 2 & 2\\1 & 2 & 3\\ \end{bmatrix}$ = 
$\begin{bmatrix}0 & 0 & 1\\0 & 1 & 1\\1 & 1 & 1\\ \end{bmatrix}^2$ .
[Remarks: $\det(M) = 1$ ; $M$ is positive definite but $L$ is not  ; the eigenvalues of $M$ are distinct ( to see this, it's easier to work with $M^{-1}$ which is tridiagonal ). ]
General theory then predicts that $L$ can be written as a polynomial in $M$.  A priori, this polynomial might be expected to have algebraic numbers as coefficients, but the examples below suggest that the situation may be nicer than that. First, when $n=2$ we find by inspection that  $L = M - I.$ Next, for $n=3$, a calculation shows that $L = M^2 - 5M + 2I$. And yet again for $n=4$, we 
find $L = 2M^3 - 19M^2 + 21M - 5I.$
In each case (so far) the polynomial turned out to have integral coefficients.
Questions: (1) Is it always true, for arbitrary $n$ , that $L$ can be expressed as a polynomial in $M$ with integral coefficients?
(2) Presumably, matrices like $M$ above which have an integral square root are
the exception rather than the rule.  Is it easy to show that almost all ( in any suitable sense) matrices in $\text{SL}(n,Z)$ do not possess a square root in $\text{GL}(n,Z)$ whenever $n \geq 2$ ?
 A: If $L$ and $M$ are $n \times n$ matrices, $g(M) = L$ for polynomial $g(z) = c_0 + \ldots + c_{n-1} z^{n-1}$ translates to the set of linear equations $\sum_{k=0}^{n-1} c_k (M^k)_{ij} = L_{ij}$ for $i,j = 1\ldots n$. 
If $L$ and $M$ have rational entries and solutions exist, then rational solutions exist.  The only surprise is that the solutions are in integers.
If $L$ has minimal polynomial $f(\lambda)$ (in your case these appear to always be the characteristic polynomial), then $g(L^2) = L$ if and only if 
$f(\lambda)$ divides $g(\lambda^2) - \lambda$.  If $f$ has degree $n$ and $g(\lambda) = \sum_{j=0}^{n-1} c_j \lambda^j$, we take $a_{kj}$ so that
$\lambda^{2k} \equiv \sum_{j=0}^{n-1} a_{kj} \lambda^j \mod f(\lambda)$.
Then $c_j$ are obtained by solving $A c = (0,1,0,\ldots,0)^T$ where $A$ is the matrix of $a_{kj}$, $k,j = 0 \ldots n-1$.  In your case it appears that $A$ always has determinant $\pm 1$, which will cause $c$ to have integer entries.  I have verified that the determinant is $\pm 1$ for $n$ up to $200$.  That seems not to be the case for "most" matrices in $SL(n,\mathbb Z)$, and the polynomials $g$ in those cases will tend to have non-integer entries.
The characteristic polynomial $f_n(\lambda)$ of $M$ appears to satisfy a recurrence: $f_{n+4}(\lambda) + (1 - 2 \lambda^2) f_{n+2}(\lambda) + \lambda^4 f_n(\lambda) = 0$.  I'm not sure if that will help.
A: PART 1. From About positive semidefiniteness of one matrix , we can deduce
$Spectrum(L)=\dfrac{(-1)^n}{2\cos(2k\pi/(2n+1))},k=1\cdots,n$ and 
$Spectrum(M)=\dfrac{1}{4(\cos(2k\pi/(2n+1)))^2},k=1\cdots,n$.
Thus, $M$ admits $n$ distinct eigenvalues, $\det(M)=1$ and there is a unique polynomial $P=p_0+\cdots +p_{n-1}x^{n-1}\in \mathbb{Q}[x]$ of degree $< n$ s.t. $L=P(M)$. We show that the $p_i$ are integers. Let $(e_i)$ be the canonical basis of $\mathbb{R}^n$.


*

*Since $\det(M)=1$, $M^{-1}=Q(M)$ where $Q\in\mathbb{Z}[x]$ has degree $<n$.

*The key is that $N=M^{-1}=[n_{i,j}]$, has $n$ distinct eigenvalues, is defined by: if $i<n$, then $n_{i,i}=2$, $n_{n,n}=1$, $n_{i,i+1}=n_{i+1,i}=1$ and the other entries are $0$. Moreover, there is a unique polynomial $R=r_0+\cdots +r_{n-1}x^{n-1}\in \mathbb{Q}[x]$ of degree $< n$ s.t. $L=R(N)$.
Since $e_n$ is the first column of $L$, $r_0,\cdots,r_{n-1}$ is a solution of the linear system (S) $r_0e_1+\cdots+r_{n-1}N^{n-1}e_1=e_n$. The matrix of (S) is $C=[C_0,\cdots,C_{n-1}]$ where $C_i$ is the first column of $N^i$. Since $C$ is a triangular matrix with diagonal $(1,-1,1,\cdots)$, $\det(C)=\pm 1$ and (S) has a unique solution that necessarily is $(r_i)$ and the $(r_i)$ are integers.


*Finally, $L=R(Q(M))=P(M)$ where $R,Q\in\mathbb{Z}[x]$ and, consequently,  $P\in\mathbb{Z}[x]$.


EDIT. PART 2. About $G=SL_n(\mathbb{Z})$. $G$ is spanned by the products of the transvections $\{I_n+\lambda E_{i,j};\lambda\in\mathbb{Z},i\not= j\}$. A random matrix $A\in G$ is obtained by the product of (say $30$ such matrices) -randomly choose the $(i,j)'s$ and the $\lambda's$, using a probability with $\mathbb{Z}$ as support-. 
Firstly, note that when $n\rightarrow\infty$, the probability that $A$ admits at least one $<0$ eigenvalue tends to $1$ and the prob. that it has only simple eigenvalues tends also to $1$. Then, when $n\rightarrow\infty$, the prob. that $A$ has no square root in $GL_n(\mathbb{R})$ tends to $1$.
Now $n$ is fixed; we are waiting for the release of $A\in G$ without $<0$ eigenvalues; with a strong probability ($<1$), such an  $A$ has simple eigenvalues $(\lambda_i)$; their square roots are $(\pm\sqrt{\lambda_i})$, calculated with the principal $\log$. $A$ admits at most $2^n$ square roots in $GL_n(\mathbb{R})$. If $A$ admits a square root in $G$, then, a necessary condition  is: 
(*) there are $\epsilon_i\in \{\pm 1\}$ s.t. $\sum_i\epsilon_i\sqrt{\lambda_i}\in\mathbb{Z}$.
Let $P$ be the unitary polynomial that admits as roots the previous $2^n$ sums; note that $P\in\mathbb{Z}[x]$. Then (*) $\Leftrightarrow $ $P$ has a linear factor. This condition can be written as algebraic equations. That shows that the probability $p$ that $A$ has a square root in $G$ is very small. Remark that $p>0$ because $M$ is such a matrix and a probability with support $\mathbb{Z}$ is nuclear, that is the prob. to obtain any element is $>0$.
