Iterate sum of two squares

Question

I'm looking for an efficient way to iterate over the infinite non-decreasing sequence defined by a^2+b^2 where a, b are both positive integers. What I mean by iterate here, is that given an existing list of n entries, the algorithm should efficiently (I'm hoping for O(log(m)), where m=a^2+b^2) find the next entry.

The start of this list is: 1^2+1^2=2, 1^2+2^2=5, 2^2+2^2=8, 1^2+3^2=10, 2^2+3^2=13, 1^2+4^2=17, ...

Here's the python code used to generate these entries (the first 100 are correct):

xs=[]
for i in range(1, 100):
    for j in range(i, 100):
        xs.append((i**2+j**2, i, j))
xs.sort()

I've looked at the start of the list and can't see any pattern at all. Anyone know an algorithm to do this?

[edit] Upon some searching, I found Cornacchia's algorithm which requires computing a quadratic residue. However I'm still hoping for something better, as we already know the previous numbers during iteration.

Answer 1

The isSumOfSquares function returns True if n can be written as a sum of squares greater than zero and False otherwise, based on an algorithm from Edsgar Dijkstra's 1976 book The Discipline of Programming: x sweeps downward from the square root of n, decreasing when the sum of the squares is too large, and y sweeps upward from 1, increasing when the sum of the squares is too small.

from math import sqrt

def isSumOfSquares(n):
    x, y = int(sqrt(n)), 1
    while x >= y:
        if x*x + y*y < n: y = y + 1
        elif x*x + y*y > n: x = x - 1
        else: return True
    return False

filter(isSumOfSquares, range(1,100))

This returns the list [2, 5, 8, 10, 13, 17, 18, 20, 25, 26, 29, 32, 34, 37, 40, 41, 45, 50, 52, 53, 58, 61, 65, 68, 72, 73, 74, 80, 82, 85, 89, 90, 97, 98]. I discussed a similar question at my blog.