idoneal.py

import heapq
import math
from fractions import gcd

D = 33
n = 10000

mx = int( math.sqrt( n ) )
my = int( math.sqrt( n / D ) )

heap = [ (x ** 2, x, 0) for x in xrange( 0, mx + 1 ) ]

heapq.heapify( heap )

lastValue = 0
lastX = None
lastY = None
duplicates = set( [0, 1] )

while len( heap ) > 0:
    (s, x, y) = heap[0]
    if s > n:
        break
    
    if s % 2 == 1:
        if s == lastValue:
            duplicates.add( s )
        elif lastValue not in duplicates:
           # Check GCD *after* uniqueness
           if gcd( lastX, D * lastY ) == 1:
              print lastValue, lastX, lastY
        lastValue = s
        lastX = x
        lastY = y
                
    if y + 1 <= my:
        heapq.heapreplace( heap, ( x ** 2 + D * (y + 1) ** 2, x, y+1 ) )
    else:
        heapq.heappop( heap )

Thanks for this programme. I am studying it very hard.