I tried to build in as many of the properties as possible as I was digging through numbers. My for loops create squares of x-y and y-z. I then realized that I could check x+y before entering the z loop. That sped things up quite a bit.
bound = 10**6
def issquare(n):
if n**.5%1 == 0:
return True
else:
return False
def squaretest(x,y,z):
if issquare(x+y) == False:
return False
elif issquare(x+z) == False:
return False
elif issquare(y+z) == False:
return False
elif issquare(x-z) == False:
return False
elif issquare(x-y) == False:
return False
elif issquare(y-z) == False:
return False
else:
return True
squares = []
for n in range(1, bound+1):
if issquare(n) == True:
squares.append(n)
for x in xrange(3, bound):
if x%1000 == 0: print x
for yi in squares:
y = x - yi
if y < 0:
break
elif issquare(x+y) == False:
continue
for zi in squares:
z = y - zi
if z < 0:
break
elif squaretest(x,y,z) == True:
print x, y, z, x+y+z