 function [Pest, Psim] = birthdayprob( N, K, Nt);
%function [Pest, Psim] = birthdaypob( N, K, Nt);
%compute all birthday probabilities.
%end the computation at 10*floor(N/30);
%K = 10*floor(N/30);
kall =[0:K-1];
%Pnohit =(N)_k/N^k;

Pnohit = cumsum( log(1 -kall/N) );
Pnohit = exp(Pnohit);
Pnoapp = exp( -(kall+1).*(kall)/(2*N) );
% a very good approximation!
X = [1:K]'*[1 1 1 1];
P = [ Pnohit;1-Pnohit; Pnoapp; 1-Pnoapp]';

% note that we are plotting both tails
semilogy( X, P,'-*');grid

break
%Birthday Simulation

edges = [0:1/N:1];
% use built-in histogram
knt = histc(rand(K,Nt), edges);
%whos
% OK -- 23 = max(sum(knt)),min(sum(knt));
maxknt = max(knt);length(maxknt);
maxnhits = 6;tot = zeros(maxnhits,1);%col
for nhit=1:maxnhits
    tot(nhit) = length(find(maxknt >= nhit));
end
semilogy(tot/Nt, '*-')
Psim = tot(2)/Nt;
Pest = exp( - K*(K-1)/(2*N));
% process each column of knt(each col is 366-tuple
% to determine if 