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). 



Given a positive-definite integral unimodular Gram matrix, how to find a basis of the associated lattice (over $\mathbf Q$)?
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$.] 


