% Random Iteration Algorithm
%  Set up for the Sierpinski triangle.


% Create the table for the IFS (number of functions x 6 numbers)
A=[0.5 0 0 0.5 0 0
   0.5 0 0 0.5 0.5 0
   0.5 0 0 0.5 0 0.5];
   
[NumFunc,N]=size(A);
if N~=6
  error('Number of entries in each row of matrix should be 6')   
end

%Randomly initialize a point:
x=randn(2,1);

%Number of iterations (number of points to plot):
NumIts=5000;

% The matrix Y will hold the points to plot at the end:
Y=zeros(2,NumIts);
Y(:,1)=x(:);

for j=1:500
  % Discard the first 100 points
  i=randi(NumFunc);
  x(1)=A(i,1)*Y(1,1)+A(i,2)*Y(2,1)+A(i,5);
  x(2)=A(i,3)*Y(1,1)+A(i,4)*Y(2,1)+A(i,6);
  Y(:,1)=x;
end


for j=1:NumIts-1
  % Select the function at random:
  i=randi(NumFunc);
  Y(1,j+1)=A(i,1)*Y(1,j)+A(i,2)*Y(2,j)+A(i,5);
  Y(2,j+1)=A(i,3)*Y(1,j)+A(i,4)*Y(2,j)+A(i,6);
end