Given a positive-definite integral unimodular Gram matrix, how to find a basis of the associated lattice (over $\mathbf Q$)? Let $G$ be a $n\times n-$symmetric matrix with integral coefficients and determinant $1$ (i.e. unimodular) such that the associated quadratic form is positive-definite.
I am interested in having an algorithm to find a rational basis of a lattice $L$ such that $G$ is the Gram matrix of $L$. Concretely, this consists in finding a square matrix $M$ (with rational coefficients) such that $G$ factors as $^tM \cdot M$.
Note that:
1) The unimodular condition and Hasse-Minkowski theory predict that this is indeed possible (this is a remark in Serre's course in arithmetics, §1.3.6 p. 86 in the French edition).
2) I demand the matrix $M$ to be square. (For non-square $M$, this is much easier: first by the Gram-Schmidt process, one can assume that $M$ is diagonal (with positive rational entries); then write each diagonal entry as a sums of  (at most 4) squares). 
 A: Here is a method.
Let $\mathrm{V}$ be the bilinear space $(\mathbf Q^n, b)$ with Gram Matrix $G$ in the canonical basis. 
Let $M_0$ be the lattice $\mathbf Z^n$ on $\mathbf Q^n$. Find a primitive vector $v_0\in M_0$ such that $b(v_0,v_0)$ is a square, say $b(v_0,v_0)=m^2$ in $M_0$ (note that this can be done quite efficiently : use LLL to get a basis of vectors with small norms, take the two* first vectors, study the lattice they generate).
( * Edit 1: in fact you need to take the first 4 vectors ... indeed, there exists (many) $3\times 3$ positive definite Gram matrices that don't represent any square over $\mathbf Q$ ... The simplest example is $7.I_3$.)
( Edit 2: If you don't use this trick, you can arrange so as to find a $v_0$ such that $m=2^s$ ... this will allow you to find an orthonormal basis inside $\mathbf Z[\frac 12]^n$, whose existence is predicted by the Strong Approximation Theorem.)
Let $K_1$ be the lattice made of vectors $w\in M_0$ such that $(v,w)\equiv 0$ mod. $m$.
Let $L_1$ be the lattice $K_1+\mathbf Z\cdot\frac{v_0}{m}$. Then $L_1$ is integral, unimodular, and has a vector of length $1$, namely $w_0:=\frac{v_0}{m}$. 
Let $M_1$ be the orthogonal complement of $w_0$ in $L_1$. Then $M_1$ is an $n-1$-dimensional unimodular lattice (contained in the $n-1$-dimensional subspace $<w_0>^\perp$ of $\mathrm{V}$). 
Go on to find $v_1\in M_1$, construct $M_2$ and so on ... you're done (the matrix $M$ you were looking for is the inverse of the matrix whose columns are $(w_0,\dots,w_{n-1})$. 
[The classification of unimodular lattices shows that at least when $n\leq 7$, you can always arrange so as to have $(v_0,v_0)=1$ (so that $M_1\subset M_0$), and for $n\leq 24$, you can find some $v_0$ with $(v_0,v_0)=4$.] 
