Find the remainder $R$ of $(1^2+1)(2^2+1)...(p^2+1)$ divided by $p$

Find the remainder $R$ of $(1^2+1)(2^2+1)...(p^2+1)$ divided by $p$, with $p$ being a prime number greater than $3$.

For $p \equiv 1 \mod 4$, there exists an integer $j$ such that $p\mid j^2+1$ (since $-1$ is a quadratic residue of $p$), therefore $R=0$.

For $p \equiv 3 \mod 4$, how can we find $R$? Does $R$ depend on the value of $p$?

If $p=2$ or $p\equiv 1\pmod{4}$, then $R$ is obviously $0$ because $-1$ is a quadratic residue modulo $p$. We assume from now on that $p\equiv 3\pmod{4}$.

Let $f(x)$ denote the polynomial $x^p-x$ over the Galois field $K=\operatorname{GF}(p)$. Note that $f$ factorizes into $$f(x)=\prod_{t\in K}(x-t).$$ Let $E$ be the extension $K[X]/(X^2+1)$ of $K$, which is just the Galois field $\operatorname{GF}(p^2)$. (This is where we use the assumption that $p\equiv 3\pmod{4}$, since $X^2+1$ is irreducible over $K$.) Write $i$ for the image of $X$ under the canonical projection $K[X]\to K[X]/(X^2+1)=E$.

Now, we have $$R=\prod_{k=1}^p(k^2+1)=\prod_{k=1}^p(i-k)(-i-k)=\left(\prod_{k=1}^p(i-k)\right)\left(\prod_{k=1}^p(-i-k)\right).$$ By the definition of $f$, we get $$R=f(i)f(-i)=\left(i^p-i\right)\left((-i)^p-(-i)\right)=\left(i^{p-1}-1\right)\left((-i)^{p-1}-1\right).$$ Since $\frac{p-1}{2}$ is odd, we have $$R=\left((-1)^{\frac{p-1}{2}}-1\right)\left((-1)^{\frac{p-1}{2}}-1\right)=\big((-1)-1\big)^2=4.$$ This is an equality in $K$. Translating this result back to $\Bbb{Z}$, we conclude that $$R=\begin{cases} 0&\text{if}\ p=2 \vee p\equiv 1\pmod{4},\\ 1&\text{if}\ p=3,\\ 4&\text{if}\ p>3 \wedge p\equiv 3\pmod{4}. \end{cases}$$

We can also use the symmetric polynomials.

Let $A=(1^2+1)(2^2+1) \cdots ((p-1)^2+1)(p^2+1)$. Then, $A \equiv (1^2+1)(2^2+1) \cdots ((p-1)^2+1) \mod p$, since $p^2+1\equiv 1 \mod p$.

Consider the symmetric polynomial $\displaystyle\prod_{i=1}^{\frac{p-1}2}(1+x_i)=\sum_{j=0}^{\frac{p-1}2}e_j(x)$, where $e_j$ is the $j$-th elementary symmetric polynomial. Then $\displaystyle A\equiv(\sum_{j=0}^{\frac{p-1}2}e_j(y))^2$, where $y_i$ belongs to the set $QR:=\{1^2,\cdots,(\frac{p-1}2)^2\}$.

Now notice that $QR$ consists of roots of $x^{(p-1)/2}-1$ (over $\mathbb F_p$). So $e_j(y)$ is $(-1)^j$ times the coefficient of $x^{\frac{p-1}2-j}$ in $x^{(p-1)/2}-1$. Thus $e_0(y)=1,\,e_{(p-1)/2}(y)=(-1)^{\frac{p+1}2}$, and $e_j(y)=0$ otherwise.

Therefore $A\equiv(1+(-1)^{(p+1)/2})^2\pmod p$.

If $ p $ is not $ 3 $ modulo $ 4 $ the given expression is trivially $ 0 $, so assume $ p \equiv 3 \pmod{4} $.

Consider the polynomial

$$ P(x) = \prod_k (x + k) $$

where the product is taken over all $ k $ which are quadratic residues mod $ p $ (excluding $ 0 $). First, observe that if $ r $ is a quadratic residue mod $ p $, then

$$ P(r) = \prod_k (r + k) = r^{(p-1)/2} \prod_k (1 + k r^{-1}) = \prod_k (1 + k) = P(1) = c $$

In addition, it follows by pairing each non-identity quadratic residue with its inverse that $ P(0) = 1 $. The only polynomial mod $ p $ of degree $ \leq (p-1)/2 $ which satisfies these conditions is $ (c-1) x^{(p-1)/2} + 1 $ (it's easy to check that it has the required values, and uniqueness is easy to see.) Since $ P $ has leading coefficient $ 1 $, it follows $ c = 2 $, and so $ P(1) = 2 $. The expression given in the question is $ P(1)^2 $, which is therefore always equal to $ 4 $ mod $ p $.

Supposing $$(x-1)(x-2)...(x-(p-1))\equiv x^{p-1}\pm1\pmod{p}$$ The left-hand side is $$(x^2-1^2)(x^2-2^2)...(x^2-({p-1\over2})^2)\pmod p$$ Replace $x^2$ by $y$; and then replace $y$ by $-1$. The left-hand side is now a square-root of your product $\pmod p$, give or take a sign. The right-hand side is $1\pm1$.
Can you fill in the details?