function [P,PM]=birthdays(k,trials)
% What is the probability of having two students in a group of k share
% the same birthday?
% format P=birthdays(k,trials)
% input N = number of individuals in group
%       trials (optional) = number of Monte Carlo trials
%       P = exact probability
%       PM = probability based on Monte Carlo simulation
% the birthday problem
% Exact solution & Monte Carlo simulation
% uses Gallagher's ties program for checking birthdays
% Written by E Gallagher, Eugene.Gallagher@umb.edu
% http://alpha.es.umb.edu/faculty/edg/files/edgwebp.html
% Written: 8/5/01, Revised 1/16/11;
P=(365^k - exp(gammaln(366)-gammaln(366-k)))/365^k;

if nargin < 2 & nargout ==2;
    trials=100000; % about 5 seconds
end
if nargout==2
    tally=0;
    for i=1:trials           
        draw=randi(365,k,1);
        if length(unique(draw))<k  % unique will drop shared birthdays
             tally=tally+1;
        end
    end
    PM=tally/trials;disp(P);
end