# Random numbers with fixed sum

Your task is to write a program or a function that outputs n random numbers from interval [0,1] with fixed sum s.

# Input

n, n≥1, number of random numbers to generate

s, s>=0, s<=n, sum of numbers to be generated

# Output

A random n-tuple of floating point numbers with all elements from the interval [0,1] and sum of all elements equal to s, output in any convenient unambiguous way. All valid n-tuples have to be equally likely within the limitations of floating point numbers.

This is equal to uniformly sampling from the intersection of the points inside the n-dimensional unit cube and the n-1-dimensional hyperplane that goes through (s/n, s/n, …, s/n) and is perpendicular to the vector (1, 1, …, 1) (see red area in Figure 1 for three examples). Figure 1: The plane of valid outputs with n=3 and sums 0.75, 1.75 and 2.75

# Examples

n=1, s=0.8 → [0.8]
n=3, s=3.0 → [1.0, 1.0, 1.0]
n=2, s=0.0 → [0.0, 0.0]
n=4, s=2.0 → [0.2509075946818119, 0.14887693388076845, 0.9449661625992032, 0.6552493088382167]
n=10, s=9.999999999999 → [0.9999999999999,0.9999999999999,0.9999999999999,0.9999999999999,0.9999999999999,0.9999999999999,0.9999999999999,0.9999999999999,0.9999999999999,0.9999999999999]


# Rules

• Your program should finish under a second on your machine at least with n≤10 and any valid s.
• If you so wish, your program can be exclusive on the upper end, i.e. s<n and the output numbers from the half-open interval [0,1) (breaking the second example)
• If your language doesn't support floating point numbers, you can fake the output with at least ten decimal digits after the decimal point.
• Standard loopholes are disallowed and standard input/output methods are allowed.
• This is , so the shortest entry, measured in bytes, wins.
• Related. May 18, 2018 at 12:15
• When you say This is equal to uniformly sampling from the intersection - i can see a program choosing randomly from just the corners of that intersection. Would that be valid ? May 18, 2018 at 12:43
• @KevinCruijssen No, that's only true for s==0 or s==3. For all other values of s, the plane has nonzero area and you have to uniform-randomly choose a point on that plane. May 18, 2018 at 13:05
• Requiring the interval to be closed or half-closed (as opposed to open) is a theoretically unobservable requirement. Many random number generators give the output in (0,1). How to test that the output interval is [ 0,1) and not (0,1)? The value 0 "never" occurs anyway May 18, 2018 at 22:29
• Is it OK if our code uses rejection sampling, and so takes very long for test cases likes=2.99999999999, n=3? May we generate random reals in multiples of, say, 1e-9?
– xnor
May 19, 2018 at 0:11

# JavaScript (Node.js), 122 115 bytes

N=>W=S=>2*S>N?W(N-S).map(t=>1-t):(t=(Q=s=>n?[r=s-s*Math.random()**(1/--n),...r>1?[++Q]:Q(s-r)]:[])(S,n=N),Q?t:W(S))


Try it online!

# Python 2, 144128 119 bytes

from random import*
def f(n,s):
r=min(s,1);x=uniform(max(0,r-(r-s/n)*2),r);return n<2and[s]or sample([x]+f(n-1,s-x),n)


Try it online!

• -20 bytes, thanks to Kevin Cruijssen
• @LuisMendo Should be fixed now May 18, 2018 at 14:27
• they still seem not uniform
– l4m2
May 18, 2018 at 15:15
• @l4m2 I ran g(4, 2.0) 1,000 times to get 4,000 points and the results look like this which appears fairly uniform. May 18, 2018 at 15:45
• 122 bytes? May 18, 2018 at 15:47
• Still wrong
– l4m2
May 23, 2018 at 1:51

# Java 8, 194188196237 236 bytes

n->s->{double r[]=new double[n+1],d[]=new double[n],t;int f=0,i=n,x=2*s>n?1:0;for(r[n]=s=x>0?n-s:s;f<1;){for(f=1;i-->1;)r[i]=Math.random()*s;for(java.util.Arrays.sort(r);i<n;d[i++]=x>0?1-t:t)f=(t=Math.abs(r[i]-r[i+1]))>1?0:f;}return d;}


+49 bytes (188 → 196 and 196 → 237) to fix the speed of test cases nearing 1, as well as fix the algorithm in general.

Try it online

Explanation:

Uses the approach in this StackoverFlow answer, inside a loop as long as one of the items is still larger than 1.
Also, if 2*s>n, s will be changed into n-s, and a flag is set to indicate we should use 1-diff instead of diff in the result-array (thanks for the tip @soktinpk and @l4m2).

n->s->{              // Method with integer & double parameters and Object return-type
double r[]=new double[n+1]
//  Double-array of random values of size n+1
d[]=new double[n],
//  Resulting double-array of size n
t;          //  Temp double
int f=0,           //  Integer-flag (every item below 1), starting at 0
i=n,           //  Index-integer, starting at n
x=             //  Integer-flag (average below 0.5), starting at:
2*s>n?       //   If two times s is larger than n:
1           //    Set this flag to 1
:            //   Else:
0;          //    Set this flag to 0
for(r[n]=s=        //  Set both the last item of r and s to:
x>0?          //   If the flag x is 1:
n-s          //    Set both to n-s
:             //   Else:
s;           //    Set both to s
f<1;){         //  Loop as long as flag f is still 0
for(f=1;         //   Reset the flag f to 1
i-->1;)      //   Inner loop i in range (n,1] (skipping the first item)
r[i]=Math.random()*s;
//    Set the i'th item in r to a random value in the range [0,s)
for(java.util.Arrays.sort(r);
//   Sort the array r from lowest to highest
i<n;         //   Inner loop i in the range [1,n)
;d[i++]=     //     After every iteration: Set the i'th item in d to:
x>0?       //      If the flag x is 1:
1-t       //       Set it to 1-t
:          //      Else:
t)        //       Set it to t
f=(t=Math.abs( //    Set t to the absolute difference of:
r[i]-r[i+1]))
//     The i'th & (i+1)'th items in r
>1?          //    And if t is larger than 1 (out of the [0,1] boundary)
0           //     Set the flag f to 0
:            //    Else:
f;}         //     Leave the flag f unchanged
return d;}         //  Return the array d as result

• Time out for test(10, 9.99);
– l4m2
May 18, 2018 at 15:33
• @l4m2 Yeah, noticed the same with 10, 9.0 right after I edited to fix the n=10, s=9.999999999999 test case.. Not sure if there is a fix in Java while still holding its uniformly randomness.. Will have to think about it for a while. For now I'll edit it to state it times out. May 18, 2018 at 16:52
• If n-s<1 you can call f(n,n-s) and flip every number over 1/2 (i.e. replace x with 1-x) like l4m2 did. This might solve the issue for numbers where s is close to n. May 19, 2018 at 2:58
• @soktinpk Thanks for the tip. It's actually s+s>n instead of n-s<1, but when I looked at the other JavaScript answers it indeed made sense. Everything is fixed now, including another bug that was still present. Bytes went up quite a bit, but every works now. Will work the byte-count down from here. :) May 19, 2018 at 11:31
• I don't know of a general proof but I believe this algorithm works because an N-dimensional hypercube can be cut into N N-dimensional hyperpyramids.
– Neil
Jul 22, 2018 at 10:48

# JavaScript (Node.js), 153 bytes

(n,s)=>s+s>n?g(n,n-s).map(r=>1-r):g(n,s)
g=(n,s)=>{do(a=[...Array(n-1)].map(_=>Math.random()*(s<1?s:1))).map(r=>t-=r,t=s);while(t<0|t>=1);return[...a,t]}


Try it online!

## C++11, 284 267 bytes

-17 bytes thanks to Zacharý
Uses C++ random library, output on the standard output

#include<iostream>
#include<random>
typedef float z;template<int N>void g(z s){z a[N],d=s/N;int i=N;for(;i;)a[--i]=d;std::uniform_real_distribution<z>u(.0,d<.5?d:1-d);std::default_random_engine e;for(;i<N;){z c=u(e);a[i]+=c;a[++i]-=c;}for(;i;)std::cout<<a[--i]<<' ';}


To call, you just need to do that :

g<2>(0.0);


Where the template parameter ( here, 2 ) is N, and the actual parameter ( here, 0.0 ) is S

• I think you can remove the space between <z> and u May 23, 2018 at 12:34
• I got it down further: typedef float z;template<int N>void g(z s){z a[N],d=s/N;int i=N;for(;i;)a[--i]=d;std::uniform_real_distribution<z>u(.0,d<.5?d:1-d);std::default_random_engine e;for(;i<N;){z c=u(e);a[i]+=c;a[++i]-=c;}for(;i;)std::cout<<a[--i]<<' ';}. A newline doesn't have to be the separator between items May 23, 2018 at 13:36
• Suggest get rid of d completely by changing d=s/N to s/=N Suggest reworking the second loop for(z c;i<N;a[++i%N]-=c)a[i]+=c=u(e); instead of for(;i<N;){z c=u(e);a[i]+=c;a[++i]-=c;} (note the added %N in order to make the program calculate the first number correctly) Sep 7, 2018 at 19:21

# C, 132127125118110 107 bytes

-2 bytes thanks to @ceilingcat

i;f(s,n,o,d)float*o,s,d;{for(i=n;i;o[--i]=d=s/n);for(;i<n;o[++i%n]-=s)o[i]+=s=(d<.5?d:1-d)*rand()/(1<<31);}


Try it online!

• Unfortunately, this answer does not meet the challenge specification. The output random numbers are not restricted to [0,1], and their joint distribution is not uniform. Jul 16, 2019 at 21:24
• @Nitrodon hey, could you please provide an input for which the output is not restricted to [0,1]? I tried a couple different examples and they all seemed to be correct, unless I misunderstood the objective. Jul 17, 2019 at 12:16
• With the RNG state on TIO, and keeping your n=4, the values s=3.23 and s=0.89 gives outputs outside of the range. More to the point, the distribution of X-s/n should depend on s, but it doesn't. Jul 17, 2019 at 18:12
• @Nitrodon Oops, my bad. I was converting some parts from the C++ answer above and forgot to add a thing. It should be fixed now? Also golfed a few bytes in the process. Jul 18, 2019 at 0:01

# Clean, 221 201 bytes

Clean, code-golf, or random numbers. Pick two.

import StdEnv,Math.Random,System._Unsafe,System.Time
g l n s#k=toReal n
|s/k>0.5=[s/k-e\\e<-g l n(k-s)]
#r=take n l
#r=[e*s/sum r\\e<-r]
|all((>)1.0)r=r=g(tl l)n s

g(genRandReal(toInt(accUnsafe time)))


Try it online!

Partial function literal :: (Int Real -> [Real]). Will only produce new results once per second.
Accurate up to at least 10 decimal places.

# R, 99 bytes (with gtools package)

f=function(n,s){if(s>n/2)return(1-f(n,n-s))
while(any(T>=1)){T=gtools::rdirichlet(1,rep(1,n))*s}
T}


Try it online!

We wish to sample uniformly from the set $$\\tilde{\mathcal A}=\{w_1,\ldots,w_n: \forall i, 0. I'll divide all the $$\w_i\$$ by $$\s\$$ and instead sample from $$\\mathcal A=\{w_1,\ldots,w_n: \forall i, 0.

If $$\s=1\$$, this is easy: it corresponds to sampling from the $$\Dirichlet(1, 1, \ldots, 1)\$$ distribution (which is the uniform over the simplex). For the general case $$\s\neq 1\$$, we use rejection sampling: sample from the Dirichlet distribution until all entries are $$\<1/s\$$, then multiply by $$\s\$$.

The trick to mirror when $$\s>n/2\$$ (which I think l4m2 was the first to figure out) is essential. Before I saw that, the number of iterations in the rejection sampler exploded for the last test case, so I spent a lot of time trying to sample efficiently from well-chosen truncated Beta distributions, but it is not necessary in the end.

# Wolfram Language (Mathematica), 92 90 bytes

If[2#2>#,1-#0[#,#-#2],While[Max[v=Differences@Sort@Join[{0,#2},RandomReal[#2,#-1]]]>1];v]&


Try it online!

Un-golfed code:

R[n_, s_] := Module[{v},
If[s <= n/2,             (* rejection sampling for s <= n/2:                        *)
While[
v = Differences[Sort[
Join[{0},RandomReal[s,n-1],{s}]]];         (* trial randoms that sum to s *)
Max[v] > 1           (* loop until good solution found                          *)
];
v,                     (* return the good solution                                *)
1 - R[n, n - s]]]      (* for s > n/2, invert the cube and rejection-sample       *)


Here is a solution that works in 55 bytes but for now (Mathematica version 12) is restricted to n=1,2,3 because RandomPoint refuses to draw points from higher-dimensional hyperplanes (in TIO's version 11.3 it also fails for n=1). It may work for higher n in the future though:

RandomPoint[1&~Array~#~Hyperplane~#2,1,{0,1}&~Array~#]&


Try it online!

Un-golfed code:

R[n_, s_] :=
RandomPoint[                           (* draw a random point from *)
Hyperplane[ConstantArray[1, n], s],  (* the hyperplane where the sum of coordinates is s *)
1,                                   (* draw only one point *)
ConstantArray[{0, 1}, n]]            (* restrict each coordinate to [0,1] *)


# Haskell, 122217 208 bytes

import System.Random
r p=randomR p
(n#s)g|n<1=[]|(x,q)<-r(max 0$s-n+1,min 1 s)g=x:((n-1)#(s-x)$q)
g![]=[]
g!a|(i,q)<-r(0,length a-1)g=a!!i:q![x|(j,x)<-zip[0..]a,i/=j]
n%s=uncurry(!).(n#s<$>).split<$>newStdGen


Try it online!

Sometimes the answers are slightly off due, I assume, to floating point error. If it's an issue I can fix it at a cost of 1 byte. I'm also not sure how uniform this is (pretty sure it's fine but I'm not all that good at this kind of thing), so I'll describe my algorithm.

Basic idea is to generate a number x then subtract it from s and recur until we have n elements then shuffle them. I generate x with an upper bound of either 1 or s (whichever is smaller) and a lower bound of s-n+1 or 0 (whichever is greater). That lower bound is there so that on the next iteration s will still be less than or equal to n (derivation: s-x<=n-1 -> s<=n-1+x -> s-(n-1)<=x -> s-n+1<=x).

EDIT: Thanks to @michi7x7 for pointing out a flaw in my uniformity. I think I've fixed it with shuffling but let me know if there's another problem

EDIT2: Improved byte count plus fixed type restriction

• Chaining uniform samples will never lead to a uniform distribution (the last coordinate is almost always bigger than 0.99 in your example) May 18, 2018 at 18:05
• @michi7x7 I see your point. What if I shuffled the order of the list after generating it? I should've taken more stats classes May 18, 2018 at 18:35
• The numbers doesn't look very uniform. Here, 8 results are >0.99, 1 is 0.96, and the last is 0.8. This is what it looks like. May 18, 2018 at 20:33
• @user1472751 there are several good answers here: stackoverflow.com/q/8064629/6774250 May 19, 2018 at 11:02
• Uniformness still has some problem, see here - there are too many zeroes generated (plot of sorted values from 1000%500)
– Angs
May 29, 2018 at 6:45

# Haskell, 188 bytes

import System.Random
import Data.List
n!s|s>n/2=map(1-)<$>n!(n-s)|1>0=(zipWith(-)=<<tail).sort.map(*s).(++[0,1::Double])<$>mapM(\a->randomIO)[2..n]>>= \a->if all(<=1)a then pure a else n!s


Ungolfed:

n!s
|s>n/2       = map (1-) <$> n!(n-s) --If total more than half the # of numbers, mirror calculation |otherwise = (zipWith(-)=<<tail) --Calculate interval lengths between consecutive numbers . sort --Sort . map(*s) --Scale . (++[0,1::Double]) --Add endpoints <$> mapM(\a->randomIO)[2..n] --Calculate n-1 random numbers
>>= \a->if all(<=1)a then pure a else n!s   --Retry if a number was too large


Try it online!