samp <- sample(1:60000, 1)
    sampleimage <- mnistdata$X[samp,]
    allimages <- mnistdata$X[1:60000, ]
    distances <- apply(allimages, 1, eucdist, y=sampleimage)

...


eucdist <- function(x, y) {
    dist <- sqrt(sum((x-y)^2))
}