/* Create an Ulam sprial in SAS. Rick Wicklin: http://blogs.sas.com/content/iml/2014/01/27/ulam-spirals/ Creating the heat map requires using SAS/IML 13.1 */ proc iml; start SpiralFrame(n); if n=1 then return(1); if n=2 then return({4 3, 1 2}); if n=3 then return({9 8 7, 2 1 6, 3 4 5}); X = j(n,n,.); /* top of frame. 's' means 'start'. 'e' means 'end' */ s = n##2; e = s - n + 2; X[1,1:n-1] = s:e; /* right side of frame */ s = e - 1; e = s - n + 2; X[1:n-1,n] = (s:e)`; /* bottom of frame */ s = e - 1; e = s - n + 2; X[n,n:2] = s:e; /* left side of frame */ s = e - 1; e = s - n + 2; X[n:2,1] = (s:e)`; return( X ); finish; /* test the frame construction */ M2 = SpiralFrame(2); M4 = SpiralFrame(4); M6 = SpiralFrame(6); print M2, M4, M6; start Rot180(m); return( m[nrow(m):1,ncol(m):1] ); finish; /* Create N x N integer matrix with elements arranged in an Ulam spiral */ start UlamSpiral(n); X = SpiralFrame(n); /* create outermost frame */ k = 2; do i = n-2 to 2 by -2; /* iteratively create smaller frames */ r = k:n-k+1; /* rows (and cols) to insert smaller frame */ X[r, r] = SpiralFrame(i); /* insert smaller frame */ k = k+1; end; if mod(n,2)=1 then return(Rot180(X)); return(X); finish; U10 = UlamSpiral(10); /* test program by creating 10 x 10 matrix */ print U10[f=3.]; /* SAS/IML 13.1 Create heat map of prime numbers */ start primes_upto(n); v = 2:int(n/2); /* from Ian Wakeling */ return(setdif( 2:n, v`*v )); finish; N = 150; X = UlamSpiral(N); primes = primes_upto(N##2); E = element(X, primes); /* this image takes about 30 seconds, including loading ODS graphics */ ods graphics / width=800 height=800; call heatmapdisc(E, {white black}) displayoutlines=0 showlegend=0 title="Ulam Spiral";