2 added 502 characters in body

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.

1

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.