How to generate a list of random numbers between -1 and 1 with sum 0 in Julia?

I'm currently working on implementing a simple genetic algorithm in Julia, and I need a way to generate n random numbers between -1 and 1 with a sum of 0. Ideally, the distribution would be uniformly random, but it would still work as long as it's a reasonable approximation. A search on Google returns these two posts (1, 2) which are related but not exactly what I'm looking for.

I've reduced the problem to generating n numbers in the range [0, 1] with sum n/2, but I'm not sure where to go from there.

Simulate u_1, ..., u_n between -1 and 1. Set x1 = (u_1-u_n)/2, x2 = (u_2-u_1)/2, x3 = (u_3-u_2)/2, ... xn = (u_n-u_{n-1})/2. Then each number x is between -1 and 1, and the numbers xs sum to 0.

n = 10
u = 2*rand(n) .- 1
x = [(u[1]-u[n])/2; diff(u) ./ 2]
sum(x)

EDIT

Here is the DRS algorithm. The call DRS(n, s, a, b) generates n random numbers between a and b that sum to s, and they are uniform.

# adapted from Roger Stafford's Matlab implementation
# https://nl.mathworks.com/matlabcentral/fileexchange/9700-random-vectors-with-fixed-sum
# Roger Stafford (2023). Random Vectors with Fixed Sum (https://www.mathworks.com/matlabcentral/fileexchange/9700-random-vectors-with-fixed-sum), MATLAB Central File Exchange. Retrieved March 23, 2023.
function DRS(n, s, a, b)
  if s < n*a || s > n*b || a >= b
    error("invalid values")
  end
  # Rescale to a unit cube: 0 <= x(i) <= 1
  s = (s - n*a) / (b-a)
  # Construct the transition probability table, t.
  # t(i,j) will be utilized only in the region where j <= i + 1.
  k = max(min(Int(floor(s)), n-1), 0) # Must have 0 <= k <= n-1
  s = max(min(s, k+1), k) # Must have k <= s <= k+1
  s1 = 
展开收缩
 # s1 & s2 will never be negative
  s2 = [i - s for i in k+n:-1:k+1]
  w = zeros(n, n+1)
  w[1, 2] = prevfloat(typemax(Float64)) # Scale for full 'double' range
  t = zeros(n-1, n)
  tiny = eps(Float64) # The smallest positive 'double' 
  for i = 2:n
    tmp1 = w[i-1, 2:i+1] .* s1[1:i]/i
    tmp2 = w[i-1, 1:i] .* s2[n-i+1:n]/i
    w[i, 2:i+1] = tmp1 + tmp2
    tmp3 = w[i, 2:i+1] .+ tiny # In case tmp1 & tmp2 are both 0,
    tmp4 = s2[n-i+1:n] .> s1[1:i] # then t is 0 on left & 1 on right
    t[i-1, 1:i] = (tmp2 ./ tmp3) .* tmp4 + (1 .- tmp1 ./ tmp3) .* (.!tmp4)
   end
  # Derive the polytope volume v from the appropriate element in the bottom row of w.
  v = n^(3/2) * (w[n, k+2] / prevfloat(typemax(Float64))) * (b - a)^(n - 1)
  # Now compute the matrix x.
  x = zeros(n)
  rt = rand(n - 1) # For random selection of simplex type
  rs = rand(n - 1) # For random location within a simplex
  j = k + 1 # For indexing in the t table
  sm = 0
  pr = 1 # Start with sum zero & product 1
  for i = (n-1):(-1):1  # Work backwards in the t table
    e = (rt[n-i] <= t[i, j]) # Use rt to choose a transition
    sx = rs[n-i] ^ (1/i) # Use rs to compute next simplex coord.
    sm = sm + (1 - sx) * pr * s / (i + 1) # Update sum
    pr = sx * pr # Update product
    x[n-i] = sm + pr * e # Calculate x using simplex coords.
    s = s - e
    j = j - e # Transition adjustment
  end
  x[n] = sm + pr * s # Compute the last x
  # Randomly permute the order in x and rescale.
  rp = rand(n) # Use rp to carry out a random permutation
  p = sortperm(rp) # 
  return (b - a) * x
 .+ a
end