_.py

from collections import defaultdict
from itertools import combinations

l = defaultdict(list)
for i, j in combinations(xrange(1, 10**3), 2):
	k = i*i + j*j
	l[k].append((i, j))
	if len(l[k]) > 2:
		# if k is a perfect square -- approximation
		if int(k**.5)**2 == k:
			print int(k**.5)
			print l[k]
			break