Algorithm for integer solutions of a circle?

Question

I am trying to search for integer solutions to the equation:

y^2 + x^2 = 2n^2

If I search this in wolfram alpha, they are all found almost immediately even for very large n. When I implemented a brute force approach it was very slow:

def psearch(n, count):
    for x in range(0, n):
        for y in range(0, n):
            if x*x + y*y == 2*n**2:
                print(x,y)
                count += 1
    return count

So I assume there is a much faster way to get all of the integer solutions to the equation above. How can I do this in python so that it will have much lower runtime?

Note: I have seen this question however it is about finding lattice points within a circle not the integer solutions to the equation of the circle. Also I am interested in finding the specific solutions not just the number of solutions.

Edit: I am still looking for something an order of magnitude faster. Here is an example: n=5 should have 12 integer solutions to find what those should be search this equation on Wolfram alpha.

Edit 2: @victor zen gave a phenomenal answer to the problem. Can anyone think of a way to optimize his solution further?

Answer 1

In your algorithm, you're searching for all possible y values. This is unnecessary. The trick here is to realize that

y^2 + x^2 = 2n^2

directly implies that

y^2 = 2n^2-x^2

so that means you only have to check that 2n^2-x^2 is a perfect square. You can do that by

y2 = 2*n*n - x*x 
#check for perfect square
y = math.sqrt(y2)
if int(y + 0.5) ** 2 == y2:
    #We have a perfect square.

Also, in your algorithm, you are only checking x values up to n. This is incorrect. Since y^2 will always be positive or zero, we can determine the highest x value we need to check by setting y^2 to its lowest value (i.e 0). Consequentially, we need to check all integer x values satisfying

x^2 <= 2n^2

which reduces to

abs(x) <= sqrt(2)*n.

Combine this with the optimization of only checking the top quadrant, and you have an optimized psearch of

def psearch(n):
    count = 0
    top = math.ceil(math.sqrt(2*n*n))
    for x in range(1, top):
        y2 = 2*n*n - x*x 
        #check for perfect square
        y = math.sqrt(y2)
        if int(y + 0.5) ** 2 == y2:
            count+=4
    return count

Answer 2

It is enough to search inside the first octant y>0, x<y (the four solutions (±n, ±n) are obvious and by symmetry a solution (x, y) yields 8 copies (±x, ±y), (±y, ±x)).

By monotonicity, for a given y there is at most one x. You can find it by following the circular arc incrementally, decreasing y then adjusting x. If you maintain the condition x²+y²≤2n² as tightly as possible, you get the code below which is optimized to use only elementary integer arithmetic (for efficiency, 2x is used instead of x).

    x, y, d= 2 * n, 2 * n, 0
    while y > 0:
        y, d= y - 2, d - y + 1
        if d < 0:
            x, d= x + 2, d + x + 1
        if d == 0:
            print(x >> 1, '² + ', y >> 1, '² = 2.', n, '²', sep='')

Here are all solutions for n between 1 and 100:

7² + 1² = 2.5²
14² + 2² = 2.10²
17² + 7² = 2.13²
21² + 3² = 2.15²
23² + 7² = 2.17²
28² + 4² = 2.20²
31² + 17² = 2.25²
35² + 5² = 2.25²
34² + 14² = 2.26²
41² + 1² = 2.29²
42² + 6² = 2.30²
46² + 14² = 2.34²
49² + 7² = 2.35²
47² + 23² = 2.37²
51² + 21² = 2.39²
56² + 8² = 2.40²
49² + 31² = 2.41²
63² + 9² = 2.45²
62² + 34² = 2.50²
70² + 10² = 2.50²
69² + 21² = 2.51²
68² + 28² = 2.52²
73² + 17² = 2.53²
77² + 11² = 2.55²
82² + 2² = 2.58²
84² + 12² = 2.60²
71² + 49² = 2.61²
79² + 47² = 2.65²
85² + 35² = 2.65²
89² + 23² = 2.65²
91² + 13² = 2.65²
92² + 28² = 2.68²
98² + 14² = 2.70²
103² + 7² = 2.73²
94² + 46² = 2.74²
93² + 51² = 2.75²
105² + 15² = 2.75²
102² + 42² = 2.78²
112² + 16² = 2.80²
98² + 62² = 2.82²
97² + 71² = 2.85²
113² + 41² = 2.85²
115² + 35² = 2.85²
119² + 17² = 2.85²
123² + 3² = 2.87²
119² + 41² = 2.89²
126² + 18² = 2.90²
119² + 49² = 2.91²
133² + 19² = 2.95²
137² + 7² = 2.97²
124² + 68² = 2.100²
140² + 20² = 2.100²