Simulate a tournament with N teams Based on a programming problem I saw recently, I am wondering if it is possible to simulate the following situation more efficiently:

*

*There are $N$ teams playing in a tournament.

*Each team plays every other team $K$ times.

*Each team has $1/2$ chance of winning every game it plays.

*There are no draws: in each game, one game wins and one game loses.

*All game outcomes are independent.

*After the tournament is complete, vector $S$ contains $N$ integers; $S[i]$ gives the number of wins by team $i$.

I am looking for a probability distribution with parameters $(N, K)$ with a typical realization $S$. I can simulate $S$ by drawing $K * N * (N - 1)$ Bernoulli random variables (see code below).
Question: Is there a known distribution (as a function of $(N, K)$) that I can sample from to get realizations $S$?.
Here is a sample Python code for $N = 8$ and $K = 3$.
from itertools import combinations
import random
random.seed(0)

n_teams = 8
k_playoffs = 3

def simulate_naive(n_teams=n_teams, k_playoffs=k_playoffs):
  # store score for each team
  # eventually, the score[i] contains the number of times
  # team i has won
  scores = [0] * n_teams
  # iterate over all pairs of teams
  for i, j in combinations(range(n_teams), 2):
    # each pair of teams plays k times
    for k in range(k_playoffs):
      # within each game, each team has 50/50 chance of winning
      if random.random() < 0.5: # condition is met 50% of the time
        scores[i] += 1
      else:
        scores[j] += 1
  return scores

example_scores = simulate_naive()
print(example_scores)
# [8, 9, 13, 4, 14, 11, 12, 13]


First thing to note is that there are $KN(N-1)/2$ Bernoulli trials to simulate, not $KN(N-1)$ since we are dealing with unique pairs of $N$ teams.  For instance, if there were only 2 teams then clearly only $K$ games would be played, not $2K$.
Now, among each pair of teams, we can simulate all $K$ games by using a $Binomial(K,.5)$ distribution instead of performing $K$ separate $Bernoulli(.5)$ simulations.  Unfortunately, I do not think we could do better due to the following reason.  Let $S_{ij}$ denote the number of games that Team $i$ won against Team $j$.  As mentioned before, we can simulate this by setting $S_{ij} \sim Binomial(K,.5)$.  Now we do not need to simulate $S_{ji}$ after simulating $S_{ij}$ since the number of games Team $j$ won against Team $i$ is $K-S_{ij}$.  Therefore, we could simulate $N(N-1)/2$ random variables ($\{S_{ij}: i<j; i,j \in \{1,\cdots,N\}\}$).
Now it may be asked if we can simulate less.  For instance, when $N=3$ we could simulate $S_{12}, S_{13},$ and $S_{23}$.  Let $S_i$ denote the number of games Team $i$ won in total.  Now the number of games that Team 1 won is $S_1 = S_{12} + S_{13}$.  Hence, one may be led to opine that it is possible to simulate $S_1$ as a $Binomial(2K,.5)$ distribution.  Then we would need only simulate $S_{23}$, saving a single simulation call.  Unfortunately, it would be impossible to solve for $S_2=(K-S_{12})+S_{23}$ and $S_3=(K-S_{13}) + (K-S_{23})$, since given $S_1 \ne 0$ we would not be able to recover $S_{12}$ and $S_{13}$.
I recommend initiating an $N\times N$ matrix with zeroes along the main diagonal.  Call this $\boldsymbol{X} = (x_{ij})$, with $x_{ij}$ denoting the element in the $i$th row and $j$th column.  Next simulate the $N(N-1)/2$ random variables $S_{ij} \sim Binomial(K,.5)$ for $i<j$ and set $x_{ij}=S_{ij}$.  Thus you will have a strictly upper-triangular matrix.  Next define a strictly lower triangular $N\times N$ matrix, $\boldsymbol{Y}$, where each element equals $K$.  Then define the matrix $\boldsymbol{Z} = \boldsymbol{Y} - \boldsymbol{X}^{\prime} + \boldsymbol{X}$.  The row sums of $\boldsymbol{Z}$ will give you $S_1, \cdots S_N$.
Here is some R code which performs quite quickly for $N$ as large as $10,000$.  For larger $N$, there may be storage issues with the matrices and the $N(N-1)/2$ length vector of Binomial random variables.
library(Matrix)
set.seed(555)
K = 10
N = 3
X = diag(rep(0,N))
Y = X
S = rbinom(choose(N,2),K,.5)
up.index=which(upper.tri(X))
X[up.index] = S
low.index=which(lower.tri(X))
Y[low.index] = K
Z = Y-t(X)+X
team.points = rowSums(Z)

