FUNCTION BIRTHDAY_CLASH, n=n, m=m, seed=seed, dp=dp

; ----------------------------------------------------------
;+
; NAME:
;       BIRTHDAY_CLASH
;
; PURPOSE:
;       Calculate probability of no shared birthdays by Monte Carlo
;
; AUTHOR:
;       Simon Vaughan
;
; CALLING SEQUENCE:
;       p = BIRTHDAY_CLASH(n=30,m=1000,seed=seed)
;
; INPUTS:
;       n      - (integer) number of people in group
;       m      - (integer) number of groups to simulate
;
; OPTIONAL INPUTS:
;       seed   - (long integer) seed value for random numbers
;
; OUTPUTS:
;       p      - (float) probability estimate
;
; OPTIONAL OUTPUTS:
;       dp     - (float) error on p
;
; DETAILS:
;       Find the probability that, out of a group of N people, 
;       none of them share a birthday. This is a classic
;       problem in elementary probability theory.
;
;       This example function uses a Monte Carlo scheme whereby 
;       random birthdays are allocated and then checked for 
;       coincidences. The fraction of groups that show at least 
;       one coincidence is the estimate of the probability
;       of at least one match. Then:
;
;          Prob(0 matches) = 1 - Prob(at least 1 match)
;
;        An estimate of the error on p is returned using
;        the binomial formula:
;
;          var[p] = p*(1-p)/N
;
;        This comes from the formula for the variance on
;        number of successes in N independent trials if the
;        chance of success (see: variance of binomial dist.)
;
; EXAMPLE USAGE:
;        p = BIRTHDAY_CLASH(n=30,m=10000,seed=seed,dp=dp)
;        print,p,' +/-',dp
;
; HISTORY:
;       08/11/2007  -- v1.0 -- first working version
;-
; ----------------------------------------------------------

; watch out for errors
  on_error,2

; ----------------------------------------------------------
; Check the arguments and input

; if N is not set use default
  if (n_elements(n) eq 0) then n=30

; if M not set use default
  if (n_elements(m) eq 0) then m=10000

; ----------------------------------------------------------
; Main part of function

; number of classes with at least one matching birthday

  matches = 0L

; arrays of days in the year - range: [0,364]

  day = make_array(365,/integer,value=0)

; loop over M independent simulated groups

  for i=0L,m-1L do begin

; array to keep track of birthdays 

      day[*] = 0

; generate a N random birthdays in the interval [0,364]

      birthday = floor(randomu(seed,n) * 365)
    
; flag which days have birthdays

      day[birthday] = 1

; are there N (independent) days flagged?
; If not, some values of BIRTHDAY must be the same

      if (total(day) lt N) then matches = matches + 1

  endfor                       

; fraction of classes with matches

  p = 1.0 - float(matches)/float(m)

; error on this from binomial formula

  dp = sqrt(p*(1.0-p)/m)

; ----------------------------------------------------------
; return result to user

  return,p

END