How many rounds to visit every place on a Monopoly style game board I see some kids playing on a circular gameboard trying to visit every field and I wonder how long it will take before they finish.
I can simulate it like below, but is there possibly an elegant solution that can compute or approximate this?
Rolling with a single 6 sided dice and moving on a Monopoly style board (36 fields) it takes on average about 12.55 rounds around the board with a variance of 13.8
The distribution seems to be approximately a Gumbel distribution (as this relates to the Coupon collectors problem which is also approximately Gumbel distributed). But there are some irregularities in the histogram which might have to deal with the moving around the board. We could use a model scheme like here but with an additional irregularity.

simulate = function(n = 4*9, d = 6) {
  fields = rep(0,n)
  position = 0
  counts = 0
  while(counts < n) {
    position = position + sample(1:d,1)
    if (fields[1 + position %% n] == 0) {
       fields[1 + position %% n] = 1
       counts = counts + 1
    }
  }
  return(position)
}


n = 36

set.seed(1)
x = replicate(5*10^4, simulate(n))
hist(x/n, breaks = seq((min(x)-4.5)/n,(max(x)+4.5)/n,4/n), xlab = "number of rounds", freq = FALSE)


mu = mean(x/n)
var = var(x/n)
scale = sqrt(var*6/pi^2)*1
location = mu - 0.57721 * scale
xs = seq(0,40,0.001)
lines(xs,evd::dgumbel(xs,location,scale), col = 2, lwd = 2)

 A: You can simulate or approximate with the coupon collector distribution
As a preliminary observation, I'll just note that you appear to be dealing with a simplified version of Monopoly where the player progresses through the board via the throw of a single die without any intervening "effects" that would cause special moves.  In the actual game you roll 2D6 for moves and it is possible to be "sent to jail" from certain board states.  My answer concerns your simplified version of the problem.
Consider the general case where you have a cycle of $N$ cells.  In theory, it is possible to analyse this process as a Markov chain, using a state vector $(x, i_1,...,i_N)$ where $x$ represents the position of the object and the binary values $i_1,...,i_N$ are indicator variables representing whether or not each of the $N$ cells have been visited.  If we let $\mathscr{S}$ denote the state space, it can be shown that there are $|\mathscr{S}| = 1 + 2^{N-1}$ possible states in the chain.$^\dagger$  The total number of distinct states that have been visited after $n$ moves is the sum of the indicator variables in the $n$th state of the chain, and the number of moves to cover all states is the "hitting time" for when this total number of states is $N$ (i.e., when all the visitation indicators are one).
While formal analysis using a Markov chain covers this type of process, with $N=36$ cells in the cycle there are $|\mathscr{S}| = 34,359,738,369$ (over thirty-four billion) possible states in the Markov chain.  This high number of states imposes computational difficulties on the analysis to the point where it probably is not a feasible method.  Your approach of using a simulation of the process seems like a reasonable fallback method to me.  So long as you use a sufficiently large number of simulations, you should be able to get a good estimate of the distribution of the number of rounds needed to visit all distinct cells in the cycle.

Approximation as an occupancy problem: To develop an approximation, we can first recognise that the coverage induced by the moves are "auto-regressive" in the sense that they are affected by the closeness of position between consecutive moves.  It is clear that coverage induced by moves that are far away from each other (i.e., separated by many other moves) are only going to be weakly correlated.  Since there are a large number of moves made to obtain full coverage, one average, the coverage induced by two separate moves is going to be only very weakly correlated.  Consequently, one reasonable approach to approximating the distribution would be to analyse an analogous "occupancy problem" where the positions covered on each move are independent uniform values over the cells.  This gives rise to the well-known coupon-collector distribution for the number of moves required to obtain full coverage (see O'Neill 2022).
To facilitate this approximation, suppose we let let $M_N$ denote the true number of moves under your cyclical movement process and let $\tilde{M}_N$ denote the number of moves required to cover all cells if the positions at the end of each move are taken to be IID uniform over all the cells in the cycle.  The probability mass function for the latter is:
$$\mathbb{P}(\tilde{M}_N = N+t) 
= \text{CoupColl}(t|N,1)
= \frac{1}{N^{N+t}} \cdot N! \cdot S(N+t-1, N-1),$$
where $S( \cdot, \cdot)$ denotes the Stirling numbers of the second kind.  This distribution is available using the dnegocc function in the occupancy package.  So long as the number of moves is large, this ought to give you a reasonably good approximation to the distribution of the number of moves required to get full coverage.
To get the corresponding distribution for the number of  rounds (i.e., the number of full or partial cycles) we would then treat each move as contributing a random distance using your die-roll mechanism.  If we let $X_1,X_2,X_3,... \sim \text{IID U} \{ 1,...,6 \}$ denote the outcomes of a sequence of dice rolls then the total distance travelled in order to obtain full coverage of the cycle would be approximated by $\tilde{D}_N = \sum_{k=1}^{\tilde{M}_N} X_k$.  We can form an approximating distribution for the total distance by applying the central limit theorem to form a mixture distribution that is a coupon-collector mixture of normal approximations.

Simulation in R: In the code below I generate simulations of your process using some alternative code and then plot the number of rounds required to visit all cells in the cycle.  (I have defined a "round" as including the last partial round.)  I will represent the state with a value POSITION and a vector INDICATORS, with the position taken over the cell values $0,...,N-1$.
#Create vector of simulations (of distance to full coverage)
N <- 36
S <- 10^5

#Generate the simulations
set.seed(1)
SIMULATIONS <- data.frame(Distance = rep(0, S), Moves = rep(0, S))
for (i in 1:S) {
  
  #Create starting state
  INDICATORS <- rep(0, N)
  POSITION   <- 0
  DISTANCE   <- 0
  MOVES      <- 0
  
  #Simulate process
  while (sum(INDICATORS) < N) {
    X <- sample(1:6, size = 1)
    POSITION <- (POSITION + X) %% N
    DISTANCE <- DISTANCE + X
    MOVES    <- MOVES + 1
    INDICATORS[POSITION+1] <- 1 }
  SIMULATIONS$Distance[i] <- DISTANCE
  SIMULATIONS$Moves[i]    <- MOVES }

#Add number of rounds
SIMULATIONS$Rounds <- ceiling(SIMULATIONS$Distance/N)



$^\dagger$ If we denote the cells as $1,...,N$ then all states after the starting state have a free value $x=1,...,N$ and $N-1$ free indicator variables (one of the indicators is fixed by the fact that the present position $x$ implies a positive indicator for that cell).  Taking the starting state as a special state and then applying the multiplication rule of counting gives the stated size.  (Note that if the movement mechanism has a maximum move length that is less than the number of cells then we may be able to eliminate some states as impossible, which would give a slight reduction in the state space.)
A: My best idea so far is still simulation. Here is a vectorized approach that can complete a large number of simulations pretty quickly. The idea is to use a state vector initialized at $S=2^c-1$ for each replication, where $c$ is the number of positions on the board. $\mod\!\!{(S,2^{p+1})}<2^p$ when position $p\in0...c-1$ has not been visited yet. When a position is visited for the first time, update $S$ to $S-2^p$. When all positions have been visited, $S=0$
rdist <- function(n, cycle = 36L, d = 6L) {
  state <- rep(2^cycle - 1, n)
  fields <- integer(n)
  totdist <- integer(n)
  k <- 0
  
  while(k < n) {
    fields <- fields + sample(d, length(fields), TRUE)
    pos <- fields %% cycle
    a <- 2^pos
    state <- state - ((state %% (2*a)) >= a)*a
    idx <- which(bln <- (state == 0))
    if (length(idx)) {
      totdist[(k + 1L):(k + length(idx))] <- fields[idx]
      k <- k + length(idx)
      fields <- fields[-idx]
      state <- state[-idx]
    }
  }
  totdist
}

system.time(totdist <- rdist(1e5))
#>    user  system elapsed 
#>    1.50    0.13    1.63

