library(ggplot2)

num_rows <- 200
divisors <- 2:100

# Create the grid on which to make the raster plot of the triangle.
# Since even and odd rows of the triangle don't have their squares
# align, break them up into four pieces each (only need to break them
# into halves to get the alignment right; break into quarters so that
# the horizontal and vertical scales are equal).
grid_x <- seq(-num_rows/2 + 0.25, num_rows/2 - 0.25, 0.5)
grid_y <- seq(0.25, num_rows - 0.25, 0.5)
grid <- as.data.frame(expand.grid(x=grid_x, y=grid_y))

# Calculate the grid indices for each row of the triangle;
# broken into even and odd grid cells for each of the two
# rows of grid cells that comprise a triangle row.
list_odds1 <- list()
list_odds2 <- list()
list_evens1 <- list()
list_evens2 <- list()

for (n in 0:(num_rows-1)) {
  grid_indices_y <- which((grid$y > n) & (grid$y < n+1))
  grid_indices_x <- which((grid$x > -(n+1)/2) & (grid$x < (n+1)/2))
  grid_indices <- intersect(grid_indices_x, grid_indices_y)
  
  # Since each entry of Pascal's triangle is broken into 4 squares, the length of the grid_indices
  # vector is 4 times the length of pascal_row.
  len <- length(grid_indices)/2
  odds1 <- seq(1, len-1, 2)
  odds2 <- seq(len+1, 2*len-1, 2)
  evens1 <- seq(2, len, 2)
  evens2 <- seq(len+2, 2*len, 2)
  
  list_odds1[[n+1]] <- grid_indices[odds1]
  list_odds2[[n+1]] <- grid_indices[odds2]
  list_evens1[[n+1]] <- grid_indices[evens1]
  list_evens2[[n+1]] <- grid_indices[evens2]
}


for (divisor in divisors) {
  print(sprintf("Starting modulo %d", divisor))
  
  # (Re-)initialise the column of the grid data frame holding the Pascal's triangle values
  # modulo the divisor:
  grid$pascal <- NA
  
  for (n in 0:(num_rows-1)) {
    if (n == 0) {
      new_row <- 1
      prev_row <- c(1, 0)
    } else {
      new_row <- rep(0, n+1)
      new_row[1] <- 1
      for (k in 2:(n+1)) {
        new_row[k] <- (prev_row[k-1] + prev_row[k]) %% divisor
      }
      prev_row <- c(new_row, 0)
    }
    
    pascal_row <- (new_row == 0) * 1
    
    grid$pascal[list_odds1[[n+1]]] <- pascal_row
    grid$pascal[list_odds2[[n+1]]] <- pascal_row
    grid$pascal[list_evens1[[n+1]]] <- pascal_row
    grid$pascal[list_evens2[[n+1]]] <- pascal_row
  }
  
  pascal_raster <- ggplot() + coord_equal() + scale_x_continuous(label=function(x) format(x + num_rows/2)) + scale_y_reverse() + guides(fill=F)
  pascal_raster <- pascal_raster + ggtitle(sprintf("Entries of Pascal's triangle divisible by %d", divisor))
  pascal_raster <- pascal_raster + ylab("n") + xlab("k")
  pascal_raster <- pascal_raster + geom_raster(data=grid, aes(x=x, y=y, fill=pascal))
  
  out_file <- sprintf("images/%03d.png", divisor)
  png(filename=out_file, width=500, height=500, units="px")
  print(pascal_raster)
  dev.off()
}