Integral quaternary forms and theta functions The following question arises when I attempt to understand the modular parameterization of the elliptic curve $$E:y^2-y=x^3-x$$
In Mazur-Swinnerton-Dyer and Zagier's construction, a theta function associated with a positive definite quadratic form is induced:
$$\theta(q)=\sum_{x\in\mathbb{Z}^4}q^{\frac{1}{2}x^{T}Ax}$$
where $$A=\left(\begin{matrix}2 & 0 & 1 & 1\\
0 & 4 & 1 &2\\
1 & 1 & 10 & 1\\
1 & 2 & 1 & 20
\end{matrix}\right)$$
$A$ is a positive definite matrix of determinant $37^2$, and we have $37A^{-1}=K^TAK$ where $K$ is an integral matrix of determinant $\pm 1$.
Question: Suppose $A$ is a positive definite $4\times 4$ matrix with integral entries. All diagonal entries are even numbers. The determinant of $A$ is a square number $N^2$. Is it true that for every $N=p$ ($p>2$ is a prime number), there is at least one $A$ that $NA^{-1}=K^TAK$, where $K$ is an integral matrix of determinant $\pm 1$?
 A: The answer is true, using the following construction.
Let $B$ be the quaternion algebra of discriminant $p$ and let $O$ be a maximal order with an element $x$ satisfying $x^2 = -p$. The reduced norm is a quadratic form on $O$, with positive definite Gram matrix $A$ of determinant $p^2$. The matrix $A^{-1}$ then represents the reduced norm on the dual lattice $O^{\sharp}$. Recall that $\text{nrd}(\text{diff}(O)) = \text{discrd}(O) = p$ and since $O$ is maximal, $\text{diff}(O)$ is invertible and $O = \text{diff}(O) O^{\sharp}$ (see e.g.
Voight, John, Quaternion algebras,  ZBL07261776. Section 16.8), hence
$$
p O^{\sharp} \subseteq \text{diff}(O) O^{\sharp} = O
$$
which implies that $p \cdot (O^{\sharp} / O) = 0$ (therefore $O^{\sharp} / O$ is an abelian group of exponent $p$ and size $p^2$, so isomorphic to $(\mathbb{Z} / p \mathbb{Z})^2$). Next, we note that the matrix $A^{-1}$ is also the change of basis matrix between the chosen basis of $O$ and its dual, therefore $pA^{-1}$ is integral.
This is not enough, but so far we have not used the element $x$. The matrix $pA^{-1}$ is the matrix representing the norm form on the ideal $xO^{\sharp}$, since $x^2 = -p$. But $nrd(xO^{\sharp})=p$ is a maximal order with an element $x$ such that $x^2 = -p$. By a theorem of Ibukiyama (see reference below), it is isomorphic (hence isometric as lattices) to $O$.
This could also be seen directly (referring to some of the answers above - it is possible to do it using a finite number of cases) by explicitly writing down the maximal orders. I will list below explicit constructions, which are based on
Ibukiyama, Tomoyoshi, On maximal orders of division quaternion algebras over the rational number field with certain optimal embeddings, Nagoya Math. J. 88, 181-195 (1982). ZBL0473.12012. -
If $p \equiv 3 \bmod 4$, then $B = (-p,-1)$, and $O = \mathbb{Z}<(1+i)/2,j>$. In this case, we have
$$
A = \left( \begin{array}{cccc} 
2 & 1 & 0 & 0 \\ 
1 & \frac{p+1}{2} & 0 & 0 \\
0 & 0 & 2 & 1 \\
0 & 0 & 1 & \frac{p+1}{2}
\end{array} \right)
,
pA^{-1} = \left( \begin{array}{cccc} 
\frac{p+1}{2} & -1 & 0 & 0 \\ 
-1 & 2 & 0 & 0 \\
0 & 0 & \frac{p+1}{2} & -1 \\
0 & 0 & -1 & 2
\end{array} \right),
K = \left( \begin{array}{cccc} 
-1 & 1 & 0 & 0 \\ 
1 & 0 & 0 & 0 \\
0 & 0 & -1 & 1 \\
0 & 0 & 1 & 0
\end{array} \right)
$$
If $p \equiv 1 \bmod 4$, then choosing $q \equiv 3 \bmod 4$ such that
$\left( \frac{q}{p} \right) = -1$, we can find $c$ such that
$c^2 \equiv -p \bmod q$, and then $B = (-p,-q)$ and
$$
O = \mathbb{Z} \oplus \mathbb{Z} \frac{1+j}{2} \oplus
\mathbb{Z} \frac{i(1+j)}{2} \oplus \mathbb{Z} \frac{(c+i)j}{q}
$$
In this case, we compute that
$$
A = \left( \begin{array}{cccc} 
2 & 1 & 0 & 0 \\ 
1 & \frac{q+1}{2} & 0 & c \\
0 & 0 & \frac{p(q+1)}{2} & p \\
0 & c & p & \frac{2(p+c^2)}{q}
\end{array} \right)
,
pA^{-1} = \left( \begin{array}{cccc} 
\frac{(q+1)(c^2+p)}{2q} & -c-\frac{c^2+p}{q} & -c & \frac{c(q+1)}{2} \\ 
-c-\frac{c^2+p}{q} & 2(c^2 + \frac{c^2+p}{q}) & 2c & -c(q+1) \\
-c & 2c & 2 & -q \\
\frac{c(q+1)}{2} & -c(q+1) & -q & \frac{q(q+1)}{2}
\end{array} \right)
$$
and if $\{e_1, e_2, e_3, e_4 \}$ is the above basis for $O$, then we see that $ \{e_2 e_4, ce_1 - e_4, e_1, j e_2 \} $ is a basis for the resulting module, which written in terms of the original basis yields the matrix
$$
K = \left( \begin{array}{cccc}
0 & c & 1 & -\frac{q+1}{2} \\
-c & 0 & 0 & 1 \\
-1 & 0 & 0 & 0 \\
\frac{q+1}{2} & -1 & 0 & 0
\end{array}
 \right).
$$
Referring to Will Jagy's wondering in the first answer, the reason that this does not work for most lattices is that they do not correspond to a maximal order (as in the answer by few_reps, the quotient of the lattices is actually cyclic) and there are only one or two (depending on $p \bmod 4$ isomorphism classes of maximal orders which contain a root of $-p$. For example, for $p = 37$, there are only two isomorphism classes of maximal orders in the quaternion algebra, and only one of them contains a root of $-p$.
A: It turns out that the OP will be satisfied with just one form for each of these square discriminants that maps to itself. I should already say that this thing reminds me of Watson transformations. However, few_reps has shown that some of the forms of the genus interchange. It seems this mapping permutes the genus, and one might need to search for a very long time to find a case when this permutation is a derangement. 
I  have two (infinite sets of) examples that suggest a derangement is going to be hard to find.  If prime $q \equiv 3 \pmod 4,$ make a quaternary form out of two copies of the binary $x^2 + xy + \left( \frac{q+1}{4} \right) y^2.$ This works, so half the primes are finished.
Next, if $p = 6k-1,$ take matrix
$$
\left(
\begin{array}{cccc}
2 & 1 & 1 & 1 \\
1 & 2 & 0 & 1 \\
1 & 0 & 4k & 2k \\
1 & 1 & 2k & 4k
\end{array}
\right)
$$
with determinant $p^2.$
The inverse times $p$ is
$$
\left(
\begin{array}{rrrr}
4k & -2k & -1 & 0 \\
-2k & 4k & 1 & -1 \\
-1 & 1 & 2 & -1 \\
0 & -1 & -1 & 2
\end{array}
\right)
$$
I will need to check for the explicit change of variables matrix, but it looks good. 
Got it, in 
$$
K =
\left(
\begin{array}{rrrr}
0 & 0 & 1 & 0 \\
0 & 0 & 0 & -1 \\
-1 & 0 & 0 & 0 \\
0 & 1 & 0 & 0
\end{array}
\right)
$$
