n = 10;
N = 5000000;
M = 500;
halfwidth = 4;
distribution='Normal';
denominator='root_j_factorial';
realcomplex='complex';
%denominator='one';

figure(1);
colormap gray;
hold off;
pdf = zeros(2*M+1, 2*M+1);
%maxp = 0;
for j = 1:N
    if mod(j, 10000) == 0 
        j 
    end

    if strcmp(distribution, 'Cauchy' )
        x = rand(n+1,1);
        a = tan(pi*(x-1/2));
	if strcmp(realcomplex, 'complex')
            y = rand(n+1,1);
            b = tan(pi*(x-1/2));
	    a = a+b*i;
	end
    elseif strcmp(distribution, 'Normal')
        a = randn(n+1,1);
	if strcmp(realcomplex, 'complex')
            b = randn(n+1,1);
	    a = a+b*i;
	end
    elseif strcmp(distribution, 'Uniform(-1,1)')
        a = 2*(rand(n+1,1)-0.5);
    elseif strcmp(distribution, 'Uniform(0,1)')
        a = rand(n+1,1);
	if strcmp(realcomplex, 'complex')
	    a = a+b*i;
	end
    elseif strcmp(distribution, 'Bernoulli(0,1)')
        a = randi(2,n+1,1)-1;
	a(1) = 1;
    elseif strcmp(distribution, 'Bernoulli(-1,1)')
        a = 2*randi(2,n+1,1)-3;
    end

    if strcmp(denominator, 'one')
    elseif strcmp(denominator, 'j_factorial')
      jfact = cumprod(1:n);
      jfact = jfact(end:-1:1);
      a(1:n) = a(1:n)./jfact';
%      for i = 1:n
%          a(i) = a(i)/factorial(n+1-i);
%      end
    elseif strcmp(denominator, 'root_j_factorial')
      jfact = cumprod(1:n);
      rootjfact = sqrt(jfact(end:-1:1));
      a(1:n) = a(1:n)./rootjfact';
    end

    r = roots(a);
    % sanity check: see if the roots are roots
%    for k = 1:n
%        x = r(k);
%	p = a(n+1);
%	for i = 1:n
%	    p = p + a(i)*x^(n+1-i);
%	end
%    end
%    maxp = max(abs(p),maxp);
    x = real(r);
    y = imag(r);
    % figure out which bin the root falls in...
    ii = min(2*M+1,max(1,ceil((2*M+1)*(x+halfwidth)/(2*halfwidth))));
    jj = min(2*M+1,max(1,ceil((2*M+1)*(y+halfwidth)/(2*halfwidth))));
    % increment that bin
    for k = 1:numel(ii)
        pdf(ii(k),jj(k)) = pdf(ii(k),jj(k)) + 1;
    end
end
% display sanity check...
%maxp
% trim off the edges
pdf(1,:) = 0;
pdf(end,:) = 0;
pdf(:,end) = 0;
pdf(:,1) = 0;
if strcmp(realcomplex,'real')
    % extract real roots from complex roots
    pdfreal = zeros(2*M+1, 2*M+1);
    pdfreal(:,M+1) = pdf(:,M+1);
    pdf = pdf - pdfreal;
    % renormalize complex and real roots separately
    max(max(pdf))
    max(max(pdfreal))
    pdf = pdf/max(max(pdf));
    pdfreal = pdfreal/max(max(pdfreal));
    pdf = pdf + pdfreal;
elseif strcmp(realcomplex,'complex')
    max(max(pdf))
    pdf = pdf/max(max(pdf));
end
% recombine into a single empirical density

% plot empirical density
pdf = pdf';
%pdf(1:end,:) = pdf(end:-1:1,:);
imagesc([-halfwidth halfwidth], [-halfwidth halfwidth], pdf);
axis equal;
axis tight;
axis xy;

%figure(2);
%colormap gray;
%hold off;
%for j = 1:10000
%    x = rand(n+1,1);
%    a = tan(pi*(x-1/2));
%    r = roots(a);
%    x = real(r);
%    y = imag(r);
%    plot(x,y,'k.');
%    hold on;
%end
%xlim ( [-5 5] );
%ylim ( [-5 5] );
