Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Having seen a twitter post speed by about Ulam Spirals I started to read up.  As the story goes in 1963 Stanislaw Ulam was bored at conference and started scribbling numbers in a spiral. What he discovered was a strange diaginal pattern of Prime Numbers. I was a little bored as well, so i decided to write up some R to plot these patterns.

Here is the main function. I am sure there are more elegant ways of calculating a Ulam spiral but here is my hack. First we calculate the diagonals and then fill in the gaps.


Ulam.Spiral<-function(N){
if (even(N)){
cat(sprintf("Error: function only accepts odd integers because it a poorly written and fragile piece of code.\n"))
}else{
m <- matrix(NA, nrow=(N), ncol=N)
top.left<-c(1,1)
bottom.right<-c(N,N)
top.right<-c(1,N)
bottom.left<-c(N,1)
n<-N
a<-N
m[median(1:N),median(1:N)]<-1

while(a>=3){
# This is an adaptation of a Euler Problem 28 solution. It calculates the diaginals.
m[bottom.right,bottom.right]<- a^2 #bottom right
m[top.left,top.left]<- a^2 - 2*a + 2 # top left
m[top.right,top.right]<- a^2 - 3*a + 3 # top right
m[bottom.left,bottom.left]<- a^2 - a + 1 #bottom left

#Both Horizontal Rows
m[top.left,(top.left+1):(top.right-1)]<-seq( m[top.left,top.left]-1,m[top.left,top.right]+1) #Fill in top row
m[bottom.left,(bottom.left+1):(bottom.right-1)]<-seq(m[bottom.left,bottom.left]+1, m[bottom.left,bottom.right]-1) #Fill in bottom row
#Left Vertical Rows
m[(top.left+1):(bottom.left-1), top.left] <- seq(m[top.left,top.left]+1, m[bottom.left,top.left]-1)#Left Hand Row
#Right Verical Row
m[(top.right+1):(bottom.right-1),top.right] <-seq(m[top.right,top.right]-1 , m[top.right,top.right] - (a-2) )
#drop down on square and repeat
a<-a-2
top.left<-c(top.left+1,top.left+1)
bottom.right<-c(bottom.right-1,bottom.right-1)
top.right<-c(top.right+1,top.right-1)
bottom.left<-c(bottom.left-1,bottom.left+1)
}
return(m)
}
}

Here is the output for a 5X5 spiral

     [,1] [,2] [,3] [,4] [,5]
[1,]   17   16   15   14   13
[2,]   18    5    4    3   12
[3,]   19    6    1    2   11
[4,]   20    7    8    9   10
[5,]   21   22   23   24   25

The code is on git hub. https://github.com/jofusa/ulam-spirals-R

Many Thanks to http://librestats.wordpress.com/2011/08/20/prime-testing-function-in-r/ for a useful IsPrime() function and saving me some time.

Here is an example Plot- Prime Number are in red: 