# Boltzmann generator for binary trees with the specification T = E + ZT^2
# the number of such trees of size n is binom(2n,n)/(n+1)
# the generating function for such trees is T(x) = (1-sqrt(1-4x))/(2x)

BoltzmannTree := proc(x)
	      # set up	
	      u := evalf(rand()/10^12);	
	      	   # maple rand() gives a 12 digit nonneg int	

	      # handle union		
	      if u < 1/((1-sqrt(1-4*x))/(2*x)) then
	      	 return E;
	      end if;		      

	      # handle product
	      return [Z, BoltzmannTree(x), BoltzmannTree(x)];
end proc;