On the product $\prod_{1\leq a<b\leq \frac{p-1}{2}}\,\left(a^2+b^2\right)$, where $p$ is prime
- Miklos Schweitzer 1962. Problem 4: https://artofproblemsolving.com/community/c7h226675p1257636
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.