Prime numbers are strange beasts. They exhibit properties of both randomness and regularity. Recently I watched an excellent nine-minute video on the *Numberphile* video blog that shows that if you write the natural numbers in a spiral pattern (called the Ulam spiral), then there are certain lines in the pattern that are very likely to contain prime numbers.

The image to the left shows the first 22,500 natural numbers arranged in the Ulam spiral pattern within a 150 x 150 grid. Cells that contain prime numbers are colored black. Cells that contain composite numbers are colored white. You can see certain diagonal lines with slope ±1 that contain a large number of prime numbers. There are conjectures in number theory that explain why certain lines along the Ulam spiral have a greater density of prime numbers than other lines. (The diagonal lines in the spiral correspond to quadratic equations.)

I don't know enough about prime numbers to blog about the mathematical properties of the Ulam spiral, but as soon as I saw the Ulam spiral explained, I knew that I wanted to generate it in SAS. The Ulam spiral packs the natural numbers into a square matrix in a certain order. This post describes how to construct the Ulam spiral in the SAS/IML matrix language.

### The Ulam spiral construction

The Ulam spiral is simple to describe. On a gridded piece of paper, write down the number 1. Then write successive integers as you spiral out in a counter-clockwise fashion, as shown in the accompanying figure.

Although you can stop writing numbers at any time, for the purpose of this post let's assume you want to fill an *N* x *N* array with numbers. Notice that when *N* is even, the last number in the array (*N*^{2}=4, 16, 36,...), appears in the upper left corner of the array, whereas for odd *N*, the last number (9, 25, 49,...) appears in the lower right corner. I found this asymmetry hard to deal with, so I created an algorithm that fills the *N* x *N* matrix so that the *N*^{2} term is always in the upper left. For odd values of *N*, the algorithm rotates the matrix after inserting the elements. I've previously shown that it is easy to write SAS/IML statements that rotate elements in a matrix.

My algorithm constructs the spiral iteratively. Notice that if you have constructed the spiral correctly for an
(*N*–2) x (*N*–2) array, you can construct the *N* x *N* array by adding two vertical columns (first and last) and two horizontal rows (first and last). I call the outer rows and columns a "frame" because they remind me of a picture frame. Given a value for *N*, you can figure out the starting and ending values of each row and column of the frame. You can use the SAS/IML index creation operator to create each row and column, as shown in the following program:

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; |

You can see that the M4 matrix fits inside the M6 frame. Similarly, the M2 matrix fits inside the M4 frame.

After writing the code that generates the frames, the rest of the construction is easy. Start by creating the frame of size *N*. Decrease *N* by 2 and iteratively create a frame of size *N*, taking care to insert the (*N*–2) x (*N*–2) array into the interior of the existing *N* x *N* array. After the array is filled, rotate the array by 180 degrees if *N* is odd. The following SAS/IML statements implement this algorithm:

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.]; |

You might not have much need to generate Ulam spirals in your work, but this exercise demonstrates several important principles of SAS/IML programming:

- It is often convenient to calculate the first and last element of an array and to use the color index operator (:) to generate the vector.
- Notice that I generated the largest frame first and inserted the smaller frames inside of it. That is more efficient than generating the smaller frames and then using concatenation to grow the matrix.
- In the same way, sometimes an algorithm is simpler or more efficient to implement if your DO loops run in reverse order.

The heat map at the beginning of this article is created by using the new HEATMAPDISC subroutine in SAS/IML 13.1. I will blog about heat maps in a future post. If you have SAS/IML 13.1 and want to play with Ulam spirals yourself, you can download the SAS code used in this blog post.

## 1 Comment

So in writing a field of numbers in a spiral on a surface... one is also giving a two dimensional structure depth. The spiral becomes a one at one end and forty-two at the other, in other words, you are adding a dimension. The spiral can now be viewed much like the rifling in a gun barrel. The pattern seen after circling the primes will be more understandable after one untwists the spiral now that it has a third dimension to be observed in. What do you get... (starting with a square spiral) a well shaded perfect vanishing point 3d hallway I would guess, but perhaps a cube, or the always possible, petunia??? I can only see this in my head, if you can afford to be so kind, will you guys do the cool stuff!?

P.S. these guys put a number of spiraling primes thru a single point over time and got and image of shapes (a cross section of a 4 dimensional object)... http://gizmodo.com/5892424/five-years-of-every-nba-shot-attempt-visualized