The other day, I found myself in a desperate need to sample uniformly from the Sierpiński triangle, but I couldn’t find any standard way to do this. It turns out that one can do this in a straightforward manner using barycentric coordinates. Barycentric coordinates are a very curious topic by themselves; given a triangle , we can specify any point on the plane with coordinates , additionally constrained by . This can be interpreted as “If we placed weight into vertex , weight into vertex and weight into vertex , the centre of mass would be in ”. We will only make use of a few facts holding in this coordinate system.

The Sierpiński triangle can be represented in barycentric coordinates as the set of points where each of the coordinates is written in binary (The sum of the three coordinates is then obviously , as expected.).

We may easily sample from the Sierpiński triangle in this representation as follows (Clearly a uniform sampling, as every point is chosen equiprobably):

For each , flip a 3-sided fair “coin” to set to one of , or .

Once we have obtained , we can readily convert to Cartesian coordinates using »standard formulae«. Let be the Cartesian coordinates of the vertices of the triangle we are computing barycentric coordinates with respect to. For an equilateral triangle (that is how the Sierpiński triangle is usually depicted), these will be . Let be our working precision (We will only approximate the points with many binary digits). Then the sampling algorithm is:

  • Generate as described above.
  • Convert binary to decimal :
  • By conversion formulate, the Cartesian coordinates of our sample will be:

Quite remarkably, the whole sampling process boils down to a simple loop.

from random import randint

def sample(precision=25):
    x = y = 0
    for i in range(1, precision):
        rand = randint(1, 3) # picking which coordinate gets 1
        v_i, w_i = int(rand==1), int(rand==2)

        x += (v_i + w_i/2) / (1 << i)
        y += w_i / (1 << i)
    return x, 3**0.5 * y / 2

You can try the sampling algorithm below:

Number of samples:
Point width (px):
Precision: