Solving matrix equation in $GL_n(\mathbb Z)$ I was doing some computation for finding Bravis groups. In the process I came to solve some kind of matrix equation. But I am not very familiar with algorithms related to it. I am hoping if someone can help me by providing a any results or references for the following question.

Given a matrix $M\in GL_n(\mathbb Z)$ how to solve the following matrix equation that is how to find matrices $N \in GL_n(\mathbb Z)$ such that
  $$N^{tr}MN=M$$

Any help will be greatly appreciated. Thanks in advance.
 A: I assume that the considered relation is $N^TMN=M$. We can do the job when $n\leq 3$; I am very skeptical about the case $n\geq 4$.
The following is a generic (random) case. 
Let $M=\begin{pmatrix}-1& -16& -3\\-23& 19& 5\\-4& 4& 1\end{pmatrix},N=[n_{i,j}]$. The system in the $n$ unknowns $(n_{i,j})$.
$\det(N)=1,N^TM=MAdjoint(N)$ 
has degree $2$; over $\mathbb{C}$, the solutions depend on $1$ parameter. More precisely, a Grobner basis has $8$ relations of which $7$ are independent linear and the last one has degree $2$. Thus, using the linear equations, we obtain 
$N=\begin{pmatrix}1+5a_1& -7+348a_1+7a_2& -1+49a_1+a_2\\6-329a_1-6a_2& 14-623a_1-13a_2& 1-42a_1-a_2\\-42+2273a_1+42a_2& -49+2273a_1+49a_2&a_2\end{pmatrix}$
and the quadratic Diophantine equation $15911a_1^2+618a_1a_2+6a_2^2-721a_1-14a_2+8=0$.
The  following website works about these equations
https://www.alpertron.com.ar/JQUAD.HTM 
It gives a sequence of solutions $(x_n,y_n)=(a_1,a_2)$ in the form 
$x_{n+1}=Px_n+Qy_n+K,y_{n+1}=Rx_n+Sy_n+L$; unfortunately, often, it does not give any initial couple! Fortunately, we know a particular solution $N=I_3$, that is associated to $a_1=0,a_2=1$. The software gives us
$P=-305, Q=-6, K= 7, R= 15911, S= 313, L= -364$.
We obtain
$N=\begin{pmatrix}6& -16& -3\\-17& 54& 10\\89& -275& -51\end{pmatrix},\begin{pmatrix}41& -135& -25\\-130& 438& 81\\670& -2249& -416\end{pmatrix},\cdots$.
