x^2 + y^2 = 19451945.  I'll factor the left side using complex numbers: (x+yi) * (x-yi) = 19451945.  So then this problem is equivalent to finding the ways to factor 19451945 into Gaussian integer conjugates.

The regular prime factorization of 19451945 is 5*73*137*389.  All four factors are of the form 4k+1, which in Gaussian integers are all semiprimes: 5=(2+i)*(2-i), 73=(8+3i)*(8-3i), 137=(11+4i)*(11-4i), and 389=(17+10i)*(17-10i).

To generate all the basic Gaussian factorizations of 19451945 it is sufficient to take all 16 ways of choosing one of the pairs of each of the four conjugate pairs, but this will count everything exactly twice since for each of the 16 ways, taking all the conjugates will produce an equivalent value.  So then I'll just fix one factor, and I'll let that be 17+10i.

Then we have eight Gaussian products to evaluate:

(2+i)*(8+3i)*(11+4i)*(17+10i) = -518+4372i

(2-i)*(8+3i)*(11+4i)*(17+10i) = 3149+3088i

(2+i)*(8-3i)*(11+4i)*(17+10i) = 2437+3676i

(2-i)*(8-3i)*(11+4i)*(17+10i) = 4403+256i

(2+i)*(8+3i)*(11-4i)*(17+10i) = 2363+3724i

(2-i)*(8+3i)*(11-4i)*(17+10i) = 4397+344i

(2+i)*(8-3i)*(11-4i)*(17+10i) = 4229+1252i

(2-i)*(8-3i)*(11-4i)*(17+10i) = 3539-2632i

Now to get the ordered pairs Ady wants, all that's left to do is to take the absolute values of the real and imaginary parts and put them in descending order creating the solution set:

{(4372,518), (3149,3088), (3676,2437), (4403,256), (3724,2363), (4397,344), (4229,1252), (3539,2632)}.

Note: If you put just the number 19451945 into Wolfram Alpha, one of the things it generates is the answer to this problem, all eight distinct ways to express 19451945 as a sum of two squares.