# generate binary trees with the specification T = E + ZT^2
# the number of such trees of size n is binom(2n,n)/(n+1)

genTree := proc(n)
	# handle union
	if n=0 then
	   return E;
	end if;

	#set up
	Bn := binomial(2*n,n)/(n+1);
	roll := rand(Bn+1);
	x := roll()/Bn;
	
	# now work with n-1 in the while loop for the product T^2
	k := 0;
	s := (binomial(2*k,k)/(k+1) * binomial(2*(n-k-1),n-k-1)/(n-k))/Bn;
	while x > s do
	      k := k+1;
	      s := s + (binomial(2*k, k)/(k+1) * binomial(2*(n-k-1), n-k-1)/(n-k))/Bn;
	end do;
	return [Z, genTree(k), genTree(n-k-1)];
end proc;