On the product $\prod_{1\leq aLet $p$ be an odd prime natural number.  Define 
$$M(p):=\prod_{1\leq a<b\leq \frac{p-1}{2}}\,\left(a^2+b^2\right)\,.$$
Prove that, over the field $\mathbb{F}_p=\mathbb{Z}/p\mathbb{Z}$ (i.e., in modulo $p$),
$$M(p)= \left\{
\begin{array}{ll}
0&\text{if }p\equiv1\pmod{4}\,,\\
+1&\text{if }p\equiv3\pmod{16}\text{ or }p\equiv15\pmod{16}\,,\\
-1&\text{if }p\equiv7\pmod{16}\text{ or }p\equiv11\pmod{16}\,.
\end{array}
\right.$$
The first equality (i.e., that $M(p)=0$ for $p\equiv 1\pmod{4}$) is trivial.  For now (see this thread), we know that $M(p)=\pm1$ for $p\equiv 3\pmod{4}$.  Under my answer in that thread, jpvee found that the claim above holds for prime natural numbers $p$ such that $p\leq  1000$.  Does anybody know whether the claim is true, and if it so, how to prove the claim?
 A: *

*Miklos Schweitzer 1962. Problem 4: https://artofproblemsolving.com/community/c7h226675p1257636
A: Let's assume that $p\equiv3\pmod4$ is a prime. As you already did in the linked thread, we can use the field $K=\Bbb{F}_p[i]=\{a+bi\mid a,b\in\Bbb{F}_p\}$.
By the properties of the norm map $N:K^*\to\Bbb{F}_p^*, N(a+bi)=a^2+b^2$, we know that every element of $F_p^*$ is the image of exactly $p+1$ elements. This is because the norm map is a surjective homomorphism between cyclic groups of respective orders $p^2-1$ and $p-1$.
To get a handle on the constraints on $a,b\in\Bbb{F}_p^2\setminus\{(0,0)\}$ we need to do the following. We have the dihedral group $G=D_4$ of order eight acting as obvious symmetries of the norm form: $(a,b)\mapsto (\pm a,\pm b)$ or $(a,b)\mapsto (\pm b,\pm a)$. It is clear that each orbit of $G$ on $K^*$ contains a unique pair $(a,b)$ such that $0\le a\le b\le(p-1)/2$. 
Here the orbits of $(a,b)$ such that $0<a<b\le(p-1)/2$ have full size $8$, but the points $(a,0),0<a\le (p-1)/2$, and $(a,a),0<a\le(p-1)/2$ have orbits of size four only as their stabilizers in $G$ have order two.
So if we write $N_1(x)$ for the number of solutions of $x=a^2+b^2,0<a<b\le(p-1)/2$, $N_2(x)$ for the number of solutions of $x=a^2,0<a\le(p-1)/2$, and $N_3(x)$ for the number of solutions of $x=a^2+a^2$, then we have, for all $x\in\Bbb{F}_p^*$, the identity
$$
p+1=8N_1(x)+4N_2(x)+4N_3(x).
$$
Clearly $N_2(x)=1$ or $0$ according to whether $x$ is a quadratic residue or not. Similarly $N_3(x)=1$ or $0$ according to whether $x$ is two times a quadratic residue or not.
I am fairly sure that the formula follows from these observations and the facts that 


*

*the product of all the elements of $\Bbb{F}_p^*$ is $-1$,

*the product of all the quadratic reside of $\Bbb{F}_p^*$ is equal to $+1$,

*the product of all the elements that are twice a square is equal to $2^{(p-1)/2}=\left(\dfrac2p\right)$ where the Legendre symbol is $+1$ when $p\equiv 7\pmod8$ and $-1$ when $p\equiv3\pmod8$.


The way this shows up in your product is that the factor $x$ occurs in it
$$
N_1(x)=\frac{p+1}8-\frac{N_2(x)+N_3(x)}2
$$
times. If two is a quadratic non-residue, then $N_2(x)+N_3(x)=1$ for all $x$
and $N_1(x)=(p-3)/8$ for all $x$ meaning that your product is $(-1)^{(p-3)/8}$. If two is a quadratic residue, then $N_2(x)+N_3(x)=2$ when $x$ is a QR and zero otherwise. In this case your product is $(-1)^{(p+1)/8}$ because the extra factor from quadratic residues is equal to one.
