<s>[INST] Fermat's infinite descent for finding the squares that sum to a prime Fermat's theorem on sum of two squares states that an odd prime $p = x^2 + y^2 \iff p \equiv 1 \pmod 4$
Applying the descent procedure I can get to $a^2 + b^2 = pc$ where $c \in \mathbb{Z} \gt 1$
I want $c = 1$, so how do I proceed from here? How do I apply the procedure iteratively?
Example:
$$
p = 97
$$
$$97 \equiv 1 \pmod 4 \implies \left(\frac{-1}{97}\right) = 1 \implies x^2 \equiv -1 \pmod {97}$$ has a solution
$$x^2 + 1 \equiv 0 \pmod {97}$$
$$x^2 + 1 = 97m$$
We find an $x,m$ that solves the equation.
$$x = 75, m = 58$$
Now, we pick an $a,b$ such that $\frac{-m}{2} \leq a,b \leq \frac{m}{2}$
$$a \equiv x \pmod m = 17$$
$$b \equiv y \pmod m = 1$$
Observations:


*

*$ a^2 + b^2 \equiv x^2 + 1 \equiv 0 (\mod m)$

*$ (a^2 + b^2) = mc$

*$ (x^2 + 1) = mp$
Plugging in $a,b,m$ for 2, we get $c = 5$

By this identity, we know that
$(a^2 + b^2)(c^2 + d^2) = (ac + bd)^2 + (ad - bc)^2$
**$(a^2 + b^2)(x^2 + 1^2) = (ax + b)(a - bx) = m^2pc$
Dividing ** by $m^2$, $pc = (\frac{ax+b}{m})^2 + (\frac{a-bx}{m})^2$
Plugging in $a,b,m,p,c$ we get that $22^2 + (-1)^2 = 97*5$
So we have two squares that add up to 5 times our $p$. How do we turn the 5 into a 1? What is the next step in the descent?
 [/INST] $$22^2+1=97\times 5$$
Use the idendity $$(2a+b)^2+(a-2b)^2=5a^2+5b^2$$
The numbers $a,b$ obvioulsy solve the equation $a^2+b^2=97$ , if 
$$2a+b=22\ ,\ a-2b=1$$
Solving this system gives $a=9\ ,\ b=4$
</s>