function[y]=birthn(n,m)
y=0;
for i=1:100
  y=y+birthn1(n,m);
end
y=y/100;
if m==2
    p1=365:-1:365-n+1;
    p=1-prod(p1)/365^n;
    disp('the theorical value is ');
    disp(p);
end
end

 
function[y]=birthn1(n,m)
%n个人中有m个人生日相同的概率
k=0;
for i=1:1000
    A=floor(rand(1,n)*365+1);
    B=tabulate(A);
    x=max(B(:,2));
    if x>=m
        k=k+1;
    end
end
y=k/1000;
end