S := {}: for x from 1 to 100 do for y from 1 to 100 do S := S union {x^2+y^2+2} od:od:S := sort(convert(S, list)): for i from 1 to 500 do if isprime(S[i]) then printf(`%d, `, S[i]) fi:od: # James A. Sellers, Feb 25 2003