function probs = BinProbs(n,p,eps)
% return vector for Bin(n,p) probs for terms greater than
%  eps times maximum term
% probs(k+1) contains b(n,p,k) = P{Bin(n,p) = k} for k =0,1,...,n


% Use fact that b(n,p,k+1)/b(n,p,k) = (n-k)*p/(k+1)/(1-p)

if (nargin <3)
	eps = .000001;
end


ratio = p/(1-p);
probs = zeros(n+1,1);

% the biggest term is somewhere near the expected value
modek = ceil(n*p);
probs(modek+1) = 1;  % rescale soon; start by getting correct ratios

k = modek;
while(k<n)
	probs(k+2) = probs(k+1)*ratio*(n-k)/(k+1);
	if(probs(k+2) < eps)
% 	fprintf(1,'Stop calculation at k > %d\n',k+1);
		break;
	end;
	k = k+1;
end;

k = modek;
while(k>0)
	probs(k) = probs(k+1)*k/ratio/(n-k+1);
	if(probs(k) < eps)
% 	fprintf(1,'Stop calculation at k < %d\n',k);
		break;
	end;
	k = k-1;
end;

probs = probs/sum(probs); % rescale to get probs summing to 1
