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 $ABC$, we can specify any point $p$ on the plane with coordinates $p \equiv (\lambda_1, \lambda_2, \lambda_3)$, additionally constrained by $\lambda_1 + \lambda_2 + \lambda_3 = 1$. This can be interpreted as “If we placed weight $\lambda_1$ into vertex $A$, weight $\lambda_2$ into vertex $B$ and weight $\lambda_3$ into vertex $C$, the centre of mass would be in $p$”. 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 $S := \{(u, v, w) = (0.u_1u_2u_3\dots,\ 0.v_1v_2v_3\dots,\ 0.w_1w_2w_3\dots)\ \mid\ u_i + v_i + w_i = 1 \text{ for all } i\}$ where each of the coordinates is written in binary (The sum of the three coordinates is then obviously $0.111\dots = 1$, 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 $i$, flip a 3-sided fair “coin” to set $(u_i, v_i, w_i)$ to one of $(1, 0, 0)$, $(0, 1, 0)$ or $(0, 0, 1)$.

Once we have obtained $(u, v, w)$, we can readily convert to Cartesian coordinates using »standard formulae«. Let $(x_A, y_A), (x_B, y_B), (x_C, y_C)$ 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 $(0, 0), (1, 0), \left(\frac 1 2, \frac {\sqrt{3}} 2\right)$. Let $prec$ be our working precision (We will only approximate the points with $prec$ many binary digits). Then the sampling algorithm is:

• Generate $(u, v, w) = (0.u_1u_2u_3\dots u_{prec},\ 0.v_1v_2v_3\dots v_{prec},\ 0.w_1w_2w_3\dots w_{prec})$ as described above.
• Convert binary $(u, v, w)$ to decimal $(u', v', w')$:
• By conversion formulate, the Cartesian coordinates of our sample will be:

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

You can try the sampling algorithm below:

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