Ulam Spirals in R and ggplot

November 29, 2011
By

(This article was first published on Command-Line Worldview, and kindly contributed to R-bloggers)

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

Created by Pretty R at inside-R.org

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:

To leave a comment for the author, please follow the link and comment on his blog: Command-Line Worldview.

R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...



If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Tags: , , , ,

Comments are closed.