function sierpinski

n = 100000;
X = zeros(n,1);
Y = zeros(n,1);

for i = 2:n
    r = rand;
    
    if(r < 0.33)
        X(i) = 0.5*X(i-1);
        Y(i) = 0.5*Y(i-1);
    elseif(r > 0.33 && r < 0.66)
        X(i) = 0.5*X(i-1) + 0.25;
        Y(i) = 0.5*Y(i-1) + 0.5*sqrt(3)/2;
    else
        X(i) = 0.5*X(i-1) + 0.5;
        Y(i) = 0.5*Y(i-1);
    end
end
figure
plot(X,Y,'.k');
axis equal;