In[]:=
t1=SessionTime[];
In[]:=
genPTunder[lim_Integer?Positive]:=Module[{prim},prim=Join@@Table[If[CoprimeQ[m,n],{2mn,m^2-n^2,m^2+n^2},##&[]],{m,2,Floor@Sqrt@lim},{n,1+m~Mod~2,m,2}];​​Union@@(Range[lim~Quotient~Max@#]~KroneckerProduct~{Sort@#}&/@prim)]
In[]:=
data=Cases[genPTunder[Floor[
2*10^8
]],ele_/;(CoprimeQ@@(ele[[;;2]])&&Max[ele[[;;2]]]≤10^4)];
In[]:=
Length[data]
Out[]=
1788
In[]:=
%*2
Out[]=
3576
In[]:=
ListPlot[data[[All,;;2]]~Join~(data[[All,{2,1}]]),ImagePaddingAutomatic,AspectRatio1,GridLines{{10000},{10000}}]
Out[]=
In[]:=
SessionTime[]-t1
Out[]=
0.5351492