This is not particularly fast but works: f = (-1 + Sqrt[3 + 6 #^2])/2 &; q = IntegerQ /@ (f /@ Range[1000000]); answer = Pick[Range[1000000], q] yields the triples : {{1, 11, 109, 1079, 10681, 105731} fanswer = f /@ answer {1, 13, 133, 1321, 13081, 129493} Dropping first case which is 0^2+1^2+2^2=1^2+2^2: s = StringTemplate["\!\(\*SuperscriptBox[\(``\), \(2\)]\)"]; func = Function[x, StringJoin @@ Riffle[TemplateApply[s, #] & /@ x, "+"]]; Grid[MapThread[{#1, "=", #2, "=", ToExpression[#2]} &, {func /@ trip, func /@ ftrip}]] yields: ![enter image description here][1] [1]: https://i.sstatic.net/CxI6I.png