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()
}