/* program from Ian Wakeling in response to article "The distribution of Pythagorean triples by angle" http://blogs.sas.com/content/iml/2015/04/15/pythagorean-triples-by-angle.html */ proc iml; /* Rick Wicklin code from blog - hops from one triple to another */ start GenerateTriples(depth=3); L1 = { 1 2 2, -2 -1 -2, 2 2 3}; L2 = { 1 2 2, 2 1 2, 2 2 3}; L3 = {-1 -2 -2, 2 1 2, 2 2 3}; NumTriples = sum(3##(0:Depth)); /* total number of triples */ triples = j(NumTriples,3,.); /* allocate array for results */ triples[1,] = {3 4 5}; /* parent triple */ k = 1; n = 2; do i = 1 to Depth; /* for each generation */ do j = 1 to 3##(i-1); /* generate 3##i children */ x = triples[k,]; triples[n,] = x*L1; triples[n+1,] = x*L2; triples[n+2,] = x*L3; k = k + 1; n = n + 3; end; end; return( triples ); finish; /* All at once method using large matrices */ a = 1000; /* sqrt of maximum hypotenuse */ m = repeat(t(1:a), 1, a); n = repeat(1:a, a); good = loc( (gcd(m, n)=1) # mod(m - n, 2) # (m > n) # ((m#m + n#n)<=(a#a)) ); m = m[good]; n = n[good]; t = (2#m#n) || (m#m - n#n) || (m#m + n#n); t1=t[,1]; /* sometimes the smallest side is in the 2nd column so flip these around */ t2=t[,2]; flipped = t1>t2; t[,1:2] = choose(flipped, t2, t1) || choose(flipped, t1, t2); t = t || atan(t[,1]/t[,2])#(180/constant('pi')); /* calculate the smallest angle in triangle */ call sort(t,{1 2 3}); create mytriples from t [colname={'X' 'Y' 'Z1' 'angle'}]; append from t; free /; /* Rick's method for comparison */ p = GenerateTriples(15); t = p[loc(p[,3]<=1E6),]; t1=t[,1]; t2=t[,2]; flipped = t1>t2; t[,1:2] = choose(flipped, t2, t1) || choose(flipped, t1, t2); call sort(t,{1 2 3}); create ricktriples from t [colname={'X' 'Y' 'Z2'}]; append from t; quit; data alltriples; merge mytriples ricktriples; by x y; run; proc sgplot data=alltriples; title "All 159,139 Pythagorean Triples"; histogram angle; xaxis values=(0 to 45 by 0.5); run; proc sgplot data=alltriples; where z2^=.; title "Triples found by the 'GenerateTriples' Module"; histogram angle; xaxis values=(0 to 45 by 0.5); run; proc sgplot data=alltriples; where z2=.; title "Triples not found by the 'GenerateTriples' Module"; histogram angle; xaxis values=(0 to 45 by 0.5); run;