Relations between $m,x,y$ for prime $p=x^2+y^2$ and $m^2 \equiv -1 \pmod p$ Things we know from elementary number theory:

*

*If a prime $p=x^2+y^2$, then $p \equiv 1 \pmod 4$

*Given the above, $(-1 \mid p) = 1$, and $\exists m : [m^2= -1 \pmod p \land 0<m< \frac12 (p-1)]$
I was looking at the relationship between $m$ and $x,y$ in these situations. Assuming WLOG that $x>y>0$, the following hold:

*

*If $y=1$, then $m=x$ (trivially).

*If $x=y+1$, then $m=x+y$. (Substituting $y=x-1$ shows this.)


After moving numbers around a bit, it seems that in general, $m=ay+bx$ holds, where $a$ and $b$ are the absolute values of the minimal Bezout coefficients for $x$ and $y$.
I haven't worked toward actually proving this, but I suspect this is a well-known result, given its simplicity. Alas, Approach $0$ and Google have let me down in terms of finding info.
Does anyone have a reference (or link to another question) with info on this sort of problem? Is this a lemma/theorem I should know the name of?
Edit: variables reversed, oops.
Edit 2: For absolute clarity, given the definitions above, the following seems to hold true empirically:
$$ux+vy=1 \implies |uy| +|vx| =m$$
But why?
Edit 3: Further tinkering yields further "coincidences." So far as I can tell, for the $4k+1$ primes, $m$ and $-m$ are the unique residues where $-m = m'$. That is, the two square roots of $-1$ mod these primes are additive inverses--which is expected--but also multiplicative inverses. Though now that I consider it, this makes sense: $m^2 \equiv -1 \implies -m^2 \equiv 1 \implies m(-m) \equiv 1 \pmod p$.
But importantly this gives me something else to search for.
 A: Proposition. Let $p\equiv 1~({\rm mod}~4)$ be a prime. Assume that $x>y>0$ are integers such that $$p=x^2+y^2\qquad (1)$$ and $$ax+by=1,\qquad (2)$$ where $a,b$ are integers. Then $m:=|bx-ay|$ satisfies $$m^2+1\equiv 0~({\rm mod~}p).$$ Furthermore, if $a,b$ are the Bézout coefficients for $x,y$, then one has $0<m\leq \frac {p-1}2$, where equality holds only when $p=5$.
Proof. $$m^2+1=(bx-ay)^2+1$$
$$=b^2(p-y^2)-2abxy+a^2(p-x^2)+1$$
$$=-(ax+by)^2+1+p(a^2+b^2)$$
$$=p(a^2+b^2)\equiv 0~({\rm mod~}p),$$ where one has used (1) and (2).
Now by construction $m>0$, so it remains to prove that $m\leq \frac{p-1}2,$ where equality holds only when $p=5$. Note that one has $$a^2+b^2\leq \frac p4,$$ which is clear when $y=1$, since $a=0,b=1,$ and $p\geq 5$. If $y>1$ (or even $y=1$), one uses the result on the Bézout coefficients to see that $|a|\leq \frac y2$ and $|b|\leq \frac x2,$ hence $a^2+b^2\leq \frac {x^2+y^2}4=\frac p4.$ Substituting this into the beginning part of the proof, one has $$m^2+1=p(a^2+b^2)\leq \frac {p^2}4$$
$$\Rightarrow m^2<\frac{p^2}4$$
$$\Rightarrow m<\frac p2$$
$$\Rightarrow m\leq \frac{p-1}2.$$ If $m=\frac{p-1}2,$ then, $p$ being odd, $$m^2+1\equiv 0~({\rm mod ~}p)$$
$$\Leftrightarrow (p-1)^2+4\equiv 0~({\rm mod~}p)$$
$$\Leftrightarrow p=5,$$ which proves the case for equality. QED
Remarks.

*

*Given $x$ with $x^2\equiv -1~({\rm mod~}p),$ there is an algorithm with a polynomial complexity by Stan Wagon (1990), based on work by Serret and Hermite (1848) and Cornacchia (1908) (see wikipedia page here).


*Gauss has shown that if $p=4k+1$ then $p=\alpha^2+\beta^2$ with $\alpha=\langle \frac 12{{2k}\choose{k}}\rangle$ and $\beta=\langle(2k)!\alpha\rangle,$ where $\langle n\rangle$ denotes the residue of $n$ mod $p$ with $|\langle n\rangle|<\frac p2.$ (See the following reference.)
Reference
Wagon, Stan (1990), “Editor’s Corner: The Euclidean Algorithm Strikes Again”, American Mathematical Monthly, 97 (2): 125.
