# Sample from the discrete triangular distribution

Given an integer n >= 1 as input, output a sample from the discrete triangular distribution over the integers k, for 1 <= k <= n (1 <= k < n is also acceptable), defined by p(k) ∝ k.

E.g. if n = 3, then p(1) = 1/6, p(2) = 2/6, and p(3) = 3/6.

Your code should take constant expected time, but you are allowed to ignore overflow, to treat floating point operations as exact, and to use a PRNG (pseudorandom number generator).

You can treat your random number generator and all standard operations as constant time.

This is code golf, so the shortest answer wins.

• An explanation of the discrete triangular distribution with some test cases would be a lot clearer. Also reasonable quality isn't an objective criteria. Oct 4, 2022 at 11:42
• I suggest you remove the "reasonable quality" requirement for the PRNG. We usually assume that PRNG's are good enough; see for example here Oct 4, 2022 at 12:22
• Perhaps you should also explain what you intend by "your code should take constant expected time". Currently, answers seem to make different assumptions about whether a PRNG can be expected to run in constant time for all values of n in random_function(n), or whether this matters. Oct 4, 2022 at 13:30

# JavaScript (ES7), 37 bytes

-1 thanks to @Neil

n=>(Math.random()*n*-~n+1|0)**.5+.5|0


Try it online!

### How?

We essentially pick a random integer in $$\\left[1\dots\sum_{k=1}^{n}k\right]\$$ and then return the corresponding term in the sequence $$\1,2,2,3,3,3,4,4,4,4,\dots\$$ (this is A002024).

A more readable form of the formula is:

$$\left\lfloor\sqrt{2 \times\left\lfloor \operatorname{rand}()\times {n+1\choose 2}+1 \right\rfloor }+1/2\right\rfloor$$

$$\left\lfloor\sqrt{2 \times\lfloor \operatorname{rand}()\times n\times(n+1)/2+1 \rfloor }+1/2\right\rfloor$$

where $$\\operatorname{rand}()\$$ is assumed to return a random value drawn from the uniform distribution in the interval $$\[0,1[\$$.

• n=>(Math.random()*n*-~n+1)**.5+.5|0 avoids the second floor.
– Neil
Oct 4, 2022 at 15:55
• @Neil This may occasionally return n+1. Oct 4, 2022 at 18:34
• My bad, but at least +1 allows you to use |0 instead of &-1.
– Neil
Oct 4, 2022 at 22:06

# Python, 57 bytes

lambda n:max(r(1,n),r(n))
from random import*
r=randrange


Attempt This Online!

### Python, 60 bytes

lambda n:max(divmod(randrange(n,n*n),n))
from random import*


Attempt This Online!

## Picture:

   |  0  1  2  3  4  5
---+------------------
1  |  1  1  2  3  4  5
2  |  2  2  2  3  4  5
3  |  3  3  3  3  4  5
4  |  4  4  4  4  4  5
5  |  5  5  5  5  5  5


# R, 46 bytes

function(n)((1+4*runif(1)*n*(n+1))^.5-1)%/%2+1


Try it online!

runif(1) calculates a random number from the uniform distribution from zero to one. We then need to work-out which triangularly-increasing interval this lies in. So we multiply by the size of the full triangle (n*(n+1)/2) to get c, and solve the quadratic equation x(x+1)/2=c = x^2+x-2*c=0 (using the quadratic formula) to get the answer x.

• The question is a little unclear, but shouldn't sample(n,1,,1:n) work? Oct 4, 2022 at 12:04
• @Giuseppe - That was my original (now-deleted) answer, but it doesn't run in constant expected time. It takes much longer to run for large values of n, presumably due to the extra work of weighted sampling from bigger-and-bigger vectors... Oct 4, 2022 at 12:06
• @Giuseppe - try it Oct 4, 2022 at 12:08
• ah, I didn't test that. Good catch! Oct 4, 2022 at 12:35

# Python, 71 bytes

lambda x:x-abs((a:=r(x)-r(x+1))+(a>=0))
from random import*
r=randrange


Attempt This Online!

Surprisingly tricky. Uses the fact that the sum of 2 numbers form a rectangle, then using a>=0 to add 1 to negative numbers then abs merges both halves of the triangle.

• It seems like your explanation is out of date? Oct 4, 2022 at 13:16

# Vyxal, 7 6 bytes

-1 thanks to @emanresu A

›℅‹$℅x  Try it Online! • The : is unnecessary Jan 28 at 1:48 # Charcoal, 14 bytes Ｉ⌊⁺·⁵₂⊕‽Ｘ⁺·⁵Ｎ²  Try it online! Link is to verbose version of code. Explanation:  Ｎ Input n as a number ⁺·⁵ Plus literal number 0.5 Ｘ ² Squared ‽ Random element of implicit range ⊕ Incremented ₂ Square root ⁺·⁵ Plus literal number 0.5 ⌊ Floor Ｉ Cast to string Implicitly print  Example: For n=4, the implicit range is [0..20), which the remaining code maps 0..1 to 1, 2..5 to 2, 6..11 to 3 and 12..19 to 4, thus p(k)∝k as desired. • Neil - does sampling from a range execute in constant time in Charcoal for all range sizes? In R, for instance, it's slower to sample from a big range than from a small one... try it in R Oct 4, 2022 at 12:19 • @DominicvanEssen Certainly Charcoal calls the underlying API in constant time... – Neil Oct 4, 2022 at 13:00 • I also tried to port your approach, and I'm having difficulty (but maybe I just don't understand the explanation...). You explain: "12..19 maps to 4", but I calculate that floor(sqrt(12 + 1.5)) is equal to 3. Oct 4, 2022 at 13:19 • @DominicvanEssen Whoops, somehow I conflated the floor(sqrt(12+1)+0.5) into a floor(sqrt(12+1.5))... fixed now. – Neil Oct 4, 2022 at 15:49 # C (gcc), 43 bytes f(n){n=sqrt(2*(rand()/21e8*n*++n/2)+1)-.5;}  Try it online! # Pyt, 9 bytes Đř⇹0⇹Ř×ʀ↑  Try it online! Port of @loopy walt's answer # Raku, 24 bytes {Bag(1..*Z=>1..$_).pick}


Try it online!

A Bag object is a set with multiplicity, and it has a pick method that returns a random key, where each key is weighted by its multiplicity in the bag. Here it's initialized with a list of the numbers from 1 up to the input argument $_, each one weighted by itself. • 1 .. * is the infinite list of natural numbers. • 1 ..$_ is the list of natural numbers from 1 up to the input argument.
• Z=> zips those two lists together with the pair creation operator =>. It stops when it reaches the end of the shorter/noninfinite list, of course.