Didn't have a chance to finish this one yet!
P(n) = div(3n^2 - n, 2)
Ps = map(P, 1:1000)
vec(collect(Iterators.product(ntuple(x -> Ps, 2)...)))
filter(x -> x[1] in Ps && x[2] in Ps, [
(P(b)-P(a), P(b)+P(a))
for a in 1:10000
for b in 1:10000
if a != b])