Sample from the discrete triangular distribution

Question

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.

Answer 1

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$$

leading to:

$$\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[\$.

Answer 2

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

Answer 3

Vyxal, 7 6 bytes

-1 thanks to @emanresu A

›℅‹$℅x

Try it Online!

Answer 4

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.

Answer 5

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.

Answer 6

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.

Answer 7

C (gcc), 43 bytes

f(n){n=sqrt(2*(rand()/21e8*n*++n/2)+1)-.5;}

Try it online!

Answer 8

Pyt, 9 bytes

Đř⇹0⇹Ř×ʀ↑

Try it online!

Port of @loopy walt's answer

Answer 9

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.