symmetric integer matrices Suppose I have a symmetric positive definite matrix $M$ with integer entries. I want to decide whether $M = A A^t,$ with $A$ likewise integral. I assume that decision problem is NP-complete, as is the question of finding the $A$ even if an oracle tells you such an $A$ exists. Can someone provide a reference (I would very much like to be wrong about the hardness of the problem...)
EDIT A remark: this question is equivalent to finding a collection of integral vectors (the columns of $A$) with prescribed distances (by the parallelogram law, the inner products give us the distances). If we require $A$ to be a $0-1$ matrix, I am pretty sure that this can encode knapsack, so is NP-complete. It seems that as per Will Jagy and Gerhard Paseman, this question (via Hasse-Minkowski) might only be as hard as factoring (which is generally conjectured to be less than NP-complete), but I haven't yet completely understood what is entailed in the Hasse-Minkowski approach...
Further EDIT
In fact, the Hasse local-to-global principle works fine for small dimensions, since the class number of identity equals one in that case, and one can enumerate solutions by the 
Smith-Minkowski-Siegel mass formula. This apparently works only in dimension at most eight. This gives the oracle (the wiki article cited seems to imply that the right hand side can be computed in polynomial time, though I am none-too-certain of this), so this gives the required oracle in small dimensions, though not obviously an algorithm for finding solutions. In dimensions greater than eight we seem to be sunk.
 A: Let me just describe what I actually do in solving $$  M = A N A^t,$$ where $M,N$ are symmetric positive definite and integral, and I want $A$ integral.
If the determinants are small enough, I check whether $ \det M / \det N$ is an integral square. If I am sure that is false, I quit. ( I often solve  $B G B^t = k F,$  where the positive integer $k$ makes the determinant ratio an integral square).
If it is worth continuing, I take the largest entry $\bar{m} = \max M_{ii}$ in the diagonal entries of $M,$ and use that to get explicit bounds on the entries of $A.$ This is easy enough, the set  of row vectors $X$ with $ X N X^t = \bar{m}$ is a smooth ellipsoid, Lagrange multipliers give bounds on each entry of $X.$ Then I run a multiple loop. Each time I get some $ X N X^t = M_{jj}$ I save that $X$ as a vector in list number $j.$ Note that there is little benefit to running this search for $X$ vectors for separate diagonal entries, it is better to just run it once for $\bar{m} = \max M_{ii}.$
Finally, I have a bunch of lists of vectors, where $L_j$ is a list of all vectors that can serve as row $j.$ I run through $L_1.$ For each one, I run through $L_2.$ If the inner product $ X_1 N X_2^t = M_{1,2},$ I then go through all possible $L_3.$ And so on. If all inner product terms work out I print out the matrix $A.$ For what I do, I usually want all possible $A,$ so I do not stop after one solution.
I suppose in general what I do would be called a backtracking algorithm. The thing that is worth emphasizing is that it is just not possible to vary an entire square matrix at once, and there is no reason to do so.
EDIT: I had forgotten this. The program I wrote for Kap, very much on his instructions, was for the "oracle" part of your question. The surprise is that I wrote it for arbitrary dimension! I then put in features to take an integer sextuple and reproduce the related 3 by 3 positive symmetric matrix. What it does is compute the principal minor determinants in order, screw around with those in a recipe that goes back to Minkowski, Hasse, and Witt. It is in Mathematica. Anyway, you are welcome to it. 
A: A minor observation: $\left( \begin{smallmatrix} N & 0 \\ 0 & N \end{smallmatrix} \right)$ is of the form $A A^T$ if and only $N$ is of the form $a^2+b^2$. Specifically, if $N=a^2+b^2$ then  $\left( \begin{smallmatrix} N & 0 \\ 0 & N \end{smallmatrix} \right) = \left( \begin{smallmatrix} a & b \\ -b & a \end{smallmatrix} \right) \left( \begin{smallmatrix} a & -b \\ b & a \end{smallmatrix} \right)$. The converse is left as an exercise.
So this problem is at least as hard as determining whether or not an integer is a sum of two squares.
We discussed the complexity of determining whether an integer is a sum of two squares here. Nobody really knew the answer, but it seemed that it might be as hard as factoring.
A: This is a comment not an answer (don't have enough points to comment). If you de-symmetrize and
 consider $M=PQ^T$. Then you get a problem obviously harder than factoring if you choose $\det P, \det Q$ to be large
 primes. There should be efficient probabilistic algorithm
 for generating large pseudo-primes in determinant form so this may be
 a practical one way function. Also by choosing the size of the
 matrix to be moderately large say 100 , one can keep the matrix
 entries small and this may be advantageous. There seems to be no
 analogue of little Fermat (and high dimensional RSA), but there may
 be some other ways to use this.
