Let's say that you have a 20x4 binary matrix mapping people to tables. Rows are people, columns are iterations of the process. Each row contains some permutation of the numbers [0, 3] to indicate the order in which each person traverses the tables. Each column contains exactly five elements corresponding to each table, meaning that at each iteration, all tables are filled evenly. Here is a sample starting configuration:
np.array([
[0, 1, 2, 3], # Person 0 goes to tables 0, 1, 2, 3, in that order
[0, 1, 2, 3],
[0, 1, 2, 3],
[0, 1, 2, 3],
[0, 1, 2, 3],
[1, 2, 3, 0], # Person 5 goes to tables 1, 2, 3, 0, in that order
[1, 2, 3, 0],
[1, 2, 3, 0],
[1, 2, 3, 0],
[1, 2, 3, 0],
[2, 3, 0, 1],
[2, 3, 0, 1],
[2, 3, 0, 1],
[2, 3, 0, 1],
[2, 3, 0, 1],
[3, 0, 1, 2],
[3, 0, 1, 2],
[3, 0, 1, 2],
[3, 0, 1, 2],
[3, 0, 1, 2],
])
This is a valid configuration. Any configuration with 5 of each element in a column and unique elements in each row is a valid configuration. One way to move to another valid configuration is to swap any of the rows: two people will swap schedules, but the constraint on the rows and columns will not change. Similarly, you can swap any two columns: the order of visits will change, but everything will remain valid. Shuffling entire rows and columns like this is not very interesting: groups of five people will end up sticking together for all the steps, even if their paths will be somewhat randomized.
More generally, you can swap a single pair of elements in a given row or column. That will break the validity of the configuration, so you will have to swap other elements a bunch of times to restore it. Here is an example of what happens if you swap the first and last element of the first row:
3, 1, 2, 0
Now you need to find another row that starts with 3 and swap the 3 and the 0 in that row. That will restore the first column, but will possibly invalidate the column that now contains the 3. So you find another 3 elsewhere in the new column, and swap with 0 in the row that contains the 3. You keep doing that until the 3 ends up in the last column for some row, and order is restored.
You can implement something like a Fisher-Yates shuffle for each row, with the added step of resolving the conditions as you go. Here is a sample implementation. I'm using numpy arrays for storage and indexing convenience, but you can do this with lists just as well:
P = 20
M = 4
N = P // M
idx = np.arange(M)
# Start with a valid configuration, same as the array written out above
paths = ((idx + idx[:, None]) % M).repeat(N, axis=0)
for i in range(P):
# Shuffle each row
for j in range(M - 1):
swap = j + np.random.choice(M - j)
if swap == j:
continue
n0 = paths[i, j] # First number to swap
n1 = paths[i, swap] # Second number to swap
# fix up the consequences
ix = swap # Index of the other column
while j != swap:
# Swap the elements in the previous row
paths[i, [j, ix]] = paths[i, [ix, j]]
# Find a new row with n1 at column j
i = np.random.choice(np.flatnonzero(paths[:, j] == n1))
# Find the location of n0 in the new row
ix = j
j = np.flatnonzero(paths[i] == n0)[0]
paths[i, [j, ix]] = paths[i, [ix, j]]
# Check the results:
assert (np.sort(paths, axis=1) == idx).all(None)
assert (np.sum(paths, axis=0) == N * M * (M - 1) // 2).all()
Comments in the code should help you follow along. I can't prove that the solution is unbiased, but it seems pretty good (and fast) from the cursory tests I ran. The assertions at the end can be moved out of the loop to check the final result, or removed entirely if you are confident in the implementation.
With the conditions as you've specified, there are 6 possible paths for each person starting at a given table (3 other tables to visit). Since there are 5 people starting at each table, it works out that in almost all cases, one or more pairs of folks starting at a given table will still end up on the same path. I'm not sure if this just happens because it's likely (~91% for M=4, N=5 without priors [see here]), or because there is something in the conditions that effectively requires it.