# Christopher Ross

ToB:=proc(a,k)
	local i, j, t:
	if((a >= 2^(k)) or (a < 0)) then:
		return(failure):
	end if:
	if(k = 0) then: return []: end if:
	j:=a:
	if(j >= 2^(k-1)) then: return [1 ,op(ToB(a-2^(k-1),k-1))];
	else: return [0,op(ToB(a,k-1))]: end if:
end:
Rule1:=proc(a)
	local i, T, r, Li:
	r:=ToB(a,8):
	Li:=[seq(ToB(7-i,3),i=0..7)];
	for i from 1 to nops(Li) do:
		T[Li[i]]:=r[i]:
	end do:
	op(T);
end:
ApRule:=proc(Li,a)
	local Li1, i, T:
	Li1:=[0,op(Li),0]:
	T:=Rule1(a):
	return([seq( T[[Li1[i],Li1[i+1],Li1[i+2]]] ,i=1..nops(Li1)-2)]);
end:
Iter:=proc(L,a,t)
	local i, j, M, pts, md:
	with(plots):
	if(2*t >= nops(L)) then: return (failure): end if: 
	M:=L:
	md:=nops(L)/2:
	for i from 1 to t do:
		for j from 1 to nops(M) do:
			if(M[j] = 0) then: 
pts[i][j]:=polygonplot([[md,0],[md,0],[md,0] ,[md,0]],color=WHITE,axes=none,scaling=CONSTRAINED):
			else: 
pts[i][j]:=polygonplot([[j,-i],[j,-(i+1)],[j+1,-(i+1)] ,[j+1,-i]],color=BLACK,axes=none,scaling=CONSTRAINED):
			end if:
		end do:
		printf("\n"):
		M:=ApRule(M,a):
	end do:
	display(seq(seq(pts[i][j],j=1..nops(M)),i=1..t));
end:	

Iter([0$50,1,0$50],30,50);