Condition under which an odd prime is a sum of two squares
In additive number theory, Fermat's theorem on sums of two squares states that an odd prime "p" can be expressed as:
formula_1
with "x" and "y" integers, if and only if
formula_2
The prime numbers for which this is true are called Pythagorean primes.
For example, the primes 5, 13, 17, 29, 37 and 41 are all congruent to 1 modulo 4, and they can be expressed as sums of two squares in the following ways:
formula_3
On the other hand, the primes 3, 7, 11, 19, 23 and 31 are all congruent to 3 modulo 4, and none of them can be expressed as the sum of two squares. This is the easier part of the theorem, and follows immediately from the observation that all squares are congruent to 0 (if number squared is even) or 1 (if number squared is odd) modulo 4.
Since the Diophantus identity implies that the product of two integers each of which can be written as the sum of two squares is itself expressible as the sum of two squares, by applying Fermat's theorem to the prime factorization of any positive integer "n", we see that if all the prime factors of "n" congruent to 3 modulo 4 occur to an even exponent, then "n" is expressible as a sum of two squares. The converse also holds. This generalization of Fermat's theorem is known as the sum of two squares theorem.
History.
Albert Girard was the first to make the observation, characterizing the positive integers (not necessarily primes) that are expressible as the sum of two squares of positive integers; this was published in 1625. The statement that every prime "p" of the form formula_4 is the sum of two squares is sometimes called "Girard's theorem". For his part, Fermat wrote an elaborate version of the statement (in which he also gave the number of possible expressions of the powers of "p" as a sum of two squares) in a letter to Marin Mersenne dated December 25, 1640: for this reason this version of the theorem is sometimes called "Fermat's Christmas theorem."
Gaussian primes.
Fermat's theorem on sums of two squares is strongly related with the theory of Gaussian primes.
A Gaussian integer is a complex number formula_5 such that a and b are integers. The "norm" formula_6 of a Gaussian integer is an integer equal to the square of the absolute value of the Gaussian integer. The norm of a product of Gaussian integers is the product of their norms. This is the Diophantus identity, which results immediately from the similar property of the absolute value.
Gaussian integers form a principal ideal domain. This implies that Gaussian primes can be defined similarly as primes numbers, that is as those Gaussian integers that are not the product of two non-units (here the units are 1, −1, "i" and −"i"). 
The multiplicative property of the norm implies that a prime number p is either a Gaussian prime or the norm of a Gaussian prime. Fermat's theorem asserts that the first case occurs when formula_7 and that the second case occurs when formula_8 and formula_9 The last case is not considered in Fermat's statement, but is trivial, as formula_10
Related results.
The above point of view on Fermat's theorem is a special case of the theory of factorization of ideals in rings of quadratic integers. In summary, if formula_11 is the ring of algebraic integers in the quadratic field, then an odd prime number p, not dividing d, is either a prime element in formula_12 or the ideal norm of an ideal of formula_12 which is necessarily prime. Moreover, the law of quadratic reciprocity allows distinguishing the two cases in terms of congruences. If formula_11 is a principal ideal domain, then p is an ideal norm if and only 
formula_15
with a and b both integers.
In a letter to Blaise Pascal dated September 25, 1654 Fermat announced the following two results that are essentially the special cases formula_16 and formula_17 If p is an odd prime, then
formula_18
formula_19
Fermat wrote also:
 "If two primes which end in 3 or 7 and surpass by 3 a multiple of 4 are multiplied, then their product will be composed of a square and the quintuple of another square."
In other words, if "p, q" are of the form 20"k" + 3 or 20"k" + 7, then "pq" = "x"2 + 5"y"2. Euler later extended this to the conjecture that
formula_20
formula_21
Both Fermat's assertion and Euler's conjecture were established by Joseph-Louis Lagrange. This more complicated formulation relies on the fact that formula_22 is not a principal ideal domain, unlike formula_23 and formula_24
Algorithm.
There is a trivial algorithm for decomposing a prime of the form formula_8 into a sum of two squares: For all n such formula_26, test whether the square root of formula_27 is an integer. If this is the case, one has got the decomposition.
However the input size of the algorithm is formula_28 the number of digits of p (up to a constant factor that depends on the numeral base). The number of needed tests is of the order of formula_29 and thus exponential in the input size. So the computational complexity of this algorithm is exponential.
A Las Vegas algorithm with a probabilistically polynomial complexity has been described by
Stan Wagon in 1990, based on work by Serret and Hermite (1848), and Cornacchia (1908).
The probabilistic part consists in finding a quadratic non-residue, which can be done with success probability formula_30 and then iterated if not successful. Conditionally this can also be done in deterministic polynomial time if the generalized Riemann hypothesis holds as explained for the Tonelli–Shanks algorithm.
Description.
Given an odd prime formula_31 in the form formula_32, first find formula_33 such that formula_34.
This can be done by finding a quadratic non-residue modulo formula_31, say formula_36, and letting
formula_37. 
Such an formula_38 will satisfy the condition since quadratic non-residues satisfy formula_39.
Once formula_38 is determined, one can apply the Euclidean algorithm with formula_31 and formula_38. Denote the first two remainders that are less than the square root of formula_31 as formula_44 and formula_45. Then it will be the case that formula_46.
Example.
Take formula_47. A possible quadratic non-residue for 97 is 13, since formula_48. so we let formula_49.
The Euclidean algorithm applied to 97 and 22 yields:
formula_50
formula_51
formula_52
formula_53
The first two remainders smaller than the square root of 97 are 9 and 4; and indeed we have formula_54, as expected.
Proofs.
Fermat usually did not write down proofs of his claims, and he did not provide a proof of this statement. The first proof was found by Euler after much effort and is based on infinite descent. He announced it in two letters to Goldbach, on May 6, 1747 and on April 12, 1749; he published the detailed proof in two articles (between 1752 and 1755). Lagrange gave a proof in 1775 that was based on his study of quadratic forms. This proof was simplified by Gauss in his "Disquisitiones Arithmeticae" (art. 182). Dedekind gave at least two proofs based on the arithmetic of the Gaussian integers. There is an elegant proof using Minkowski's theorem about convex sets. Simplifying an earlier short proof due to Heath-Brown (who was inspired by Liouville's idea), Zagier presented a non-constructive one-sentence proof in 1990.
And more recently Christopher gave a partition-theoretic proof.
Euler's proof by infinite descent.
Euler succeeded in proving Fermat's theorem on sums of two squares in 1749, when he was forty-two years old. He communicated this in a letter to Goldbach dated 12 April 1749. The proof relies on infinite descent, and is only briefly sketched in the letter. The full proof consists in five steps and is published in two papers. The first four steps are Propositions 1 to 4 of the first paper and do not correspond exactly to the four steps below. The fifth step below is from the second paper.
For the avoidance of ambiguity, zero will always be a valid possible constituent of "sums of two squares", so for example every square of an integer is trivially expressible as the sum of two squares by setting one of them to be zero.
1. "The product of two numbers, each of which is a sum of two squares, is itself a sum of two squares."
This is a well-known property, based on the identity
formula_55
due to Diophantus.
2. "If a number which is a sum of two squares is divisible by a prime which is a sum of two squares, then the quotient is a sum of two squares."
(This is Euler's first Proposition).
Indeed, suppose for example that formula_56 is divisible by formula_57 and that this latter is a prime. Then formula_58 divides
formula_59
Since formula_57 is a prime, it divides one of the two factors. Suppose that it divides formula_61. Since
formula_55
(Diophantus's identity) it follows that formula_57 must divide formula_64. So the equation can be divided by the square of formula_57. Dividing the expression by formula_66 yields:
formula_67
and thus expresses the quotient as a sum of two squares, as claimed.
On the other hand if formula_57 divides formula_69, a similar argument holds by using the following variant of Diophantus's identity:
formula_70
3. "If a number which can be written as a sum of two squares is divisible by a number which is not a sum of two squares, then the quotient has a factor which is not a sum of two squares." (This is Euler's second Proposition).
Suppose formula_36 is a number not expressible as a sum of two squares, which divides formula_72. Write the quotient, factored into its (possibly repeated) prime factors, as formula_73 so that formula_74. If all factors formula_75 can be written as sums of two squares, then we can divide formula_72 successively by formula_77, formula_78, etc., and applying step (2.) above we deduce that each successive, smaller, quotient is a sum of two squares. If we get all the way down to formula_36 then formula_36 itself would have to be equal to the sum of two squares, which is a contradiction. So at least one of the primes formula_75 is not the sum of two squares.
4. "If formula_44 and formula_45 are relatively prime positive integers then every factor of formula_56 is a sum of two squares."
(This is the step that uses step (3.) to produce an 'infinite descent' and was Euler's Proposition 4. The proof sketched below also includes the proof of his Proposition 3).
Let formula_85 be relatively prime positive integers: without loss of generality formula_72 is not itself prime, otherwise there is nothing to prove. Let formula_36 therefore be a "proper" factor of formula_72, not necessarily prime: we wish to show that formula_36 is a sum of two squares. Again, we lose nothing by assuming formula_90 since the case formula_91 is obvious.
Let formula_92 be non-negative integers such that formula_93 are the closest multiples of formula_36 (in absolute value) to formula_85 respectively. Notice that the differences formula_96 and formula_97 are integers of absolute value strictly less than formula_98: indeed, when formula_90 is even, gcdformula_100; otherwise since gcdformula_101, we would also have gcdformula_102.
Multiplying out we obtain
formula_103
uniquely defining a non-negative integer formula_104. Since formula_36 divides both ends of this equation sequence it follows that formula_106 must also be divisible by formula_36: say formula_108. Let formula_109 be the gcd of formula_110 and formula_111 which by the co-primeness of formula_85 is relatively prime to formula_36. Thus formula_114 divides formula_115, so writing formula_116, formula_117 and formula_118, we obtain the expression formula_119 for relatively prime formula_120 and formula_121, and with formula_122, since
formula_123
Now finally, the "descent" step: if formula_36 is not the sum of two squares, then by step (3.) there must be a factor formula_125 say of formula_126 which is not the sum of two squares. But formula_127 and so repeating these steps (initially with formula_128 in place of formula_129, and so on "ad infinitum") we shall be able to find a strictly decreasing infinite sequence formula_130 of positive integers which are not themselves the sums of two squares but which divide into a sum of two relatively prime squares. Since such an infinite descent is impossible, we conclude that formula_36 must be expressible as a sum of two squares, as claimed.
5. "Every prime of the form formula_4 is a sum of two squares."
(This is the main result of Euler's second paper).
If formula_133, then by Fermat's Little Theorem each of the numbers formula_134 is congruent to one modulo formula_31. The differences formula_136 are therefore all divisible by formula_31. Each of these differences can be factored as
formula_138
Since formula_31 is prime, it must divide one of the