Sum of two squares - Number of steps in Fermat descent If a prime $p$ can be written as the sum of two squares, then one can construct this representation via Fermat descent if we know an $x$ such that $x^2 \equiv −1 \mod p$. Is there a possibility to say how much steps the Fermat descent will have without just going through the descent?
I have also posted this question on MathStackExchange a week ago (here), but received no answer. I hope it is ok to post it here, too.
EDIT: For example, let $p=1553$. Then $x=339$ and $339^2 + 1^2 = 74 \cdot 1553$. With the descent, we next compute numbers $x_2$, $y_2$ with $x_2^2+y_2^2 = k \cdot 1553$ with $k < 74$. In this case, $x_2=-142$, $y_2=5$ and $x_2^2+y_2^2=13 \cdot 1553$. Going on, in the next step we get $(-9)^2 + (-55)^2 = 2 \cdot 1553$ and finally $(-32)^2 + 23^2 = 1 \cdot 1553$. So here we have three steps: $74 \mapsto 13, 13 \mapsto 2, 2 \mapsto 1$.
 A: In each iteration of Fermat descent, the multiple of $p$ written as a sum of two squares is at least divided by 2. So if you start with $x>0$ satisfying $x\equiv -1\mod p$, the procedure ends in at most $\log_2(\frac{x^2+1}{p})$ iterations.   
A: Let $l$ be the square root of $-1$: $l^2+1\equiv 0 \mod p$. Then from factorization $x^2+y^2\equiv(x+ly)(x-ly)$ we see that a set
$$\{(x,y)\in\mathbb{Z}^2:x^2+y^2\equiv 0\mod p\}$$ is the union of two latticies
$$\Lambda_\pm=\{(x,y)\in\mathbb{Z}^2:x\pm ly\equiv 0\mod p\}.$$
Each step of Fermat descent is almost the same (but not exactly the same) as a step from Gauss reduction of given basis $(0,p)$, $(\mp l,1)$. It means that it works like "nearest integer continued fraction" algorithm. But vectors $(x_k,y_k)$ may jump between $\Lambda_+$ and $\Lambda_-$, numbers $x_k$, $y_k$ may be reminders or swapped reminders and so on. If you'll need more details on Fermat descent then it will be necessary to take care about all these technical difficulties.
