There is a method, algorithm really, that deserves to be better known. Given a symmetric matrix $H$ of integers, it provides a matrix $P$ with rational (or integer) entries and $\det P = 1,$ along with a diagonal matrix $D,$ such that
$$ P^T H P = D. $$
Since $\det P = 1 \;$ (and $P$ is usually upper triangular), it is not so hard to find $Q = P^{-1},$ after which
$$ Q^T D Q = H. $$
=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
$$ P^T H P = D $$
$$\left(
\begin{array}{rrr}
1 & 0 & 0 \\
- 1 & 1 & 0 \\
0 & - 1 & 1 \\
\end{array}
\right)
\left(
\begin{array}{rrr}
1 & 1 & 1 \\
1 & 2 & 2 \\
1 & 2 & 3 \\
\end{array}
\right)
\left(
\begin{array}{rrr}
1 & - 1 & 0 \\
0 & 1 & - 1 \\
0 & 0 & 1 \\
\end{array}
\right)
= \left(
\begin{array}{rrr}
1 & 0 & 0 \\
0 & 1 & 0 \\
0 & 0 & 1 \\
\end{array}
\right)
$$
$$ Q^T D Q = H $$
$$\left(
\begin{array}{rrr}
1 & 0 & 0 \\
1 & 1 & 0 \\
1 & 1 & 1 \\
\end{array}
\right)
\left(
\begin{array}{rrr}
1 & 0 & 0 \\
0 & 1 & 0 \\
0 & 0 & 1 \\
\end{array}
\right)
\left(
\begin{array}{rrr}
1 & 1 & 1 \\
0 & 1 & 1 \\
0 & 0 & 1 \\
\end{array}
\right)
= \left(
\begin{array}{rrr}
1 & 1 & 1 \\
1 & 2 & 2 \\
1 & 2 & 3 \\
\end{array}
\right)
$$
See, for example, reference for linear algebra books that teach reverse Hermite method for symmetric matrices
Illustrated here, with notation change $D$ = h2.
GP/PARI CALCULATOR Version 2.5.5 (released)
i686 running linux (ix86/GMP-5.1.2 kernel) 32-bit version
compiled: Sep 30 2013, gcc-4.8.1 (Ubuntu/Linaro 4.8.1-10ubuntu4)
(readline v6.3 enabled [was v6.2 in Configure], extended help enabled)
Copyright (C) 2000-2013 The PARI Group
PARI/GP is free software, covered by the GNU General Public License, and comes WITHOUT ANY WARRANTY WHATSOEVER.
Type ? for help, \q to quit.
Type ?12 for how to get moral (and possibly technical) support.
parisize = 4000000, primelimit = 500509
? h = [ 1,1,1; 1,2,2; 1,2,3]
%1 =
[1 1 1]
[1 2 2]
[1 2 3]
? ht = mattranspose(h)
%2 =
[1 1 1]
[1 2 2]
[1 2 3]
? ht - h
%3 =
[0 0 0]
[0 0 0]
[0 0 0]
? p1 = [1,-1,-1; 0,1,0; 0,0,1]
%4 =
[1 -1 -1]
[0 1 0]
[0 0 1]
? p1t = mattranspose(p1)
%5 =
[1 0 0]
[-1 1 0]
[-1 0 1]
? h1 = p1t * h * p1
%6 =
[1 0 0]
[0 1 1]
[0 1 2]
? p2 = [1,0,0; 0,1,-1; 0,0,1]
%7 =
[1 0 0]
[0 1 -1]
[0 0 1]
? p2t = mattranspose(p2)
%8 =
[1 0 0]
[0 1 0]
[0 -1 1]
? h2 = p2t * h1 * p2
%9 =
[1 0 0]
[0 1 0]
[0 0 1]
? p = p1 * p2
%10 =
[1 -1 0]
[0 1 -1]
[0 0 1]
? q = matadjoint(p)
%11 =
[1 1 1]
[0 1 1]
[0 0 1]
? qt = mattranspose(q)
%12 =
[1 0 0]
[1 1 0]
[1 1 1]
? qt * q
%13 =
[1 1 1]
[1 2 2]
[1 2 3]
? h
%14 =
[1 1 1]
[1 2 2]
[1 2 3]
?
? h2
%15 =
[1 0 0]
[0 1 0]
[0 0 1]
? qt * h2 * q
%16 =
[1 1 1]
[1 2 2]
[1 2 3]
?