This is my submission for the June 2016 Community Challenge. It is R code that takes as input the number of squares on the board, the desired number of chutes and ladders, and the desired total sum of the chute and ladder deltas. It then samples uniformly at random from all feasible boards with these properties to return a board; each feasible board fitting the input requirements is given exactly the same probability of being selected.

To do this, the code defines a mathematical program that has a binary decision variable for each valid chute or ladder. The constraints are used to ensure a feasible board is returned (no more than one chute/ladder starting at a square, a chute cannot be chained to a ladder, and a ladder cannot be chained to a chute) and that the input requirements are met (desired number of chutes and ladders and desired sum of deltas). Each chute and ladder is assigned a random weight and the optimization model seeks the feasible board with minimum weight; this ensures we are selecting uniformly at random from all feasible boards.

library(lpSolve)

random.game.board <- function(squares=100, num.chute=9, num.ladder=9, delta=-50) {
  if (squares < 3) stop("Game board must have at least 3 squares")
  all.paired.squares <- expand.grid(from=seq_len(squares), to=seq_len(squares))
  chutes <- subset(all.paired.squares, from < squares & to < from)
  ladders <- subset(all.paired.squares, from > 1 & to > from)
  mod <- lp(objective.in = runif(nrow(chutes) + nrow(ladders)),
            const.mat =
              rbind(t(sapply(seq_len(squares), function(x) {
                      # Square x begins no more than one chute or ladder
                      as.numeric(c(chutes$from == x, ladders$from == x))
                    })),
                    t(sapply(seq_len(squares), function(x) {
                      # Square x does not begin a chute and end a ladder
                      as.numeric(c(chutes$from == x, ladders$to == x))
                    })),
                    t(sapply(seq_len(squares), function(x) {
                      # Square x does not end a chute and begin a ladder
                      as.numeric(c(chutes$to == x, ladders$from == x))
                    })),
                    # Number of chutes
                    rep(1:0, c(nrow(chutes), nrow(ladders))),
                    # Number of ladders
                    rep(0:1, c(nrow(chutes), nrow(ladders))),
                    # Delta of chutes and ladders
                    c(chutes$to - chutes$from, ladders$to - ladders$from)),
            const.dir=rep(c("<=", "="), c(3*squares, 3)),
            const.rhs=c(rep(1, 3*squares), num.chute, num.ladder, delta),
            all.bin=TRUE)
  if (mod$status != 0) stop("No feasible game boards")}
  board <- rbind(chutes, ladders)[mod$solution > 0.999,]
  board$delta <- board$to - board$from
  board <- board[order(board$from),]
  row.names(board) <- NULL
  board
}

This is an optimization problem with a huge number of binary variables (9702), but it is so loosely constrained that the open-source lpSolve package can solve it to optimality in about 1 second on my computer with 100 squares and 9 chutes and ladders, yielding a randomly selected board:

set.seed(144)
random.game.board()
#    from  to delta
# 1     3  75    72
# 2    10  43    33
# 3    12  38    26
# 4    18  41    23
# 5    26   4   -22
# 6    34  13   -21
# 7    50  69    19
# 8    53  25   -28
# 9    56  71    15
# 10   60  59    -1
# 11   64  22   -42
# 12   65  29   -36
# 13   66  85    19
# 14   72  84    12
# 15   74  23   -51
# 16   82  44   -38
# 17   89  57   -32
# 18   98 100     2

I would appreciate comments on any aspects of my code, though I would especially be interested in comments on:

1. Whether more efficient approaches exist to sample uniformly across all feasible boards
2. About approaches to vectorize the construction of the constraint matrix
3. About the error handling
4. Whether I could condense the three lines at the end of the function that order the rows, restore default row names, and return.