Conjectural coefficient pattern for Gaussian integer polynomials

Let $p$ be a prime and set $r=(p-1)^2$. Define

$$g_p(x)=\prod_{\substack{1\leq n,m\leq p-1\\(p,n^2+m^2)=1}}(x-(n+mi))=x^{(p-1)^2}+a_1x^{(p-1)^2-1}+\cdots+a_{r-1}x+a_r.$$

Conjectural coefficient pattern. If $p=4k+3>3$, then

$$g_p(x)\equiv 1+x^{2(p-1)}+x^{4(p-1)}+\cdots+x^{(p-1)^2}\equiv\frac{1-x^{p^2-1}}{1-x^{2(p-1)}}\pmod p.$$

If $p=4k+1>5$, then

$$g_p(x)\equiv 1+b_1x^{p-1}+b_2x^{2(p-1)}+\cdots+x^{(p-1)(p-3)}\pmod p,$$

where $b_j\equiv a_{j(p-1)}\pmod p$. The displayed pattern is based on computations; the source asks whether a closed formula for the coefficients $b_j$ exists.