Integer solutions to $x^2+y^2=N$? This is not homework; but I am completely lost on how to go about handling this especially for large N
if you want integer solutions, then the primes of the form $4k+3$ dividing $N$ must occur to an even power.  To see this you must show that

i) the product of a sum of two squares is a sum of two squares

ii) an odd prime $p$ is the sum of two squares iff $p\equiv1\mod 4$

i) $(a^2+b^2)(c^2+d^2)=(ac-bd)^2+(ad+bc)^2$ (think of the norms of complex numbers)

ii)wikipedia entry with some proofs
Ok, here's my attempt at this. The wikipedia page got better results though :(


define f(x,y) = x^2+y^2
define g(n) = N
rewrite the equation to look like f(x,y) = g(n)
Now use a pullback from category theory. You go backwards in the functions.
It starts like a pair (k,k), where k=k.
then left side of the pair tries to find x and y, while right side tries to find n.
the left side is more complicated. First consider a+b expression. You can split number N to two parts so that a+b = N (this N comes from left side of pair). Any splitting is ok. Think of range [0..N] and split it to [0..a] and [a..N]. Choose freely both a and b, but so that a+b=N.
Now you have chosen the a and b. These are numbers representing the area of a square. Now to find x and y, you just move from the area of square to the length of the side of the square. $x=\sqrt a$, $y=\sqrt b$ This will result in x : R, y : R. These are not in $N$, unless you choose a and b in certain way. So the process of choosing a and b fails if length of square's side is not integers.
Problem with algorithms solving this is that there are many different ways to choose the splitting, especially when the numbers are large. And then later in the process, many of those possible solutions (a,b) will fail to materialize.
There is one additional trick available which lets us narrow the solutions. The right side g(n)=N, can be rewritten as if it were a constant. Split $g : N \rightarrow N$ to $! : N \rightarrow 1$ and $g_1 : 1 \rightarrow N, g_1=N$. Now this might seem like technicality, but it gives us one big advantage - the g can only give us single element of $N$. So all values of n are in the same area inside N in right side. If we consider tuples (x,y,n) which are the solutions to the equation, large number of them are the same, and we only need to consider (x,y) as a shorthand notation for (x,y,N). This is only necessary because we chose to use a function g originally and did not regognize it to be a constant.
Now there is also intersection operation on powersets involved. To find out which (x,y,n) tuples are solutions of the equation, the projections (x,y,n) |-> (x,y) and (x,y,n) |-> (n) need to be "inverted" using the pullback. Left and right side of the equation gave us "areas" of values of x, y and n which are involved, now combining those areas gives us the solution. But then we end up with two (x,y,n) (powerset) tuples, and their intersection is needed to find the solutions of the equation. The pair (k,k) needs to be considered at this stage when both (x,y,n) tuples are combined - the pairs create equivalence class on the tuples which allows us to find several/all solutions to the equation. Jumping from left side to right side and back to left side of the equation might end up in different elements, and if we started from a solution to the equation, the new elements are again a solution to the equation.
Now going further than this likely requires breaking $\sqrt a$ and $\sqrt b$ to something smaller and discovering more of it's properties. (I will still need to read the wikipedia page properly to see how they did it)