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])
```

© Daniel Marvin. Last modified: July 22, 2021.