/* 

  Birthday paradox (birthday problem) in Picat.

  http://en.wikipedia.org/wiki/Birthday_paradox
  """
  In probability theory, the birthday problem, or birthday paradox, pertains to the 
  probability that in a set of randomly chosen people some pair of them will have 
  the same birthday. In a group of 23 (or more) randomly chosen people, there is 
  more than 50  probability that some pair of them will both have been born on the 
  same day of the year. For 57 or more people, the probability is more than 99%, 
  tending toward 100  as the pool of people grows. The mathematics behind this 
  problem leads to a well-known cryptographic attack called the birthday attack.
  """

  Also see my (Swedish) essay 
  "Simulering av sammanträffanden - I" (Simulating coincidences)
  
  http://www.hakank.org/sims/coincidence_simulating.html


  This Picat model was created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Picat page: http://www.hakank.org/picat/

*/


% import util.
import mip.


main => go.


go =>

  N = 365, % number of days of a year
  prob(N, ProbDup),
  % println(probDup=ProbDup),
  
  % when do we reach 50%?
  Found = false,
  foreach(I in 1..ProbDup.len, Found ==false)
     println(I=ProbDup[I]),
     if ProbDup[I] >= 0.5 then
        printf("We reach 50%% for I=%d\n", I),
        Found := true
     end
  end,

  nl.



% How many items are needed to get a 50% probability of duplications
% where the attribute ("birthday") is in 1..N?
go2 => 
  foreach(N in 1..10000) 
    prob(N, ProbDup),
    % when do we reach 50%?
    Found = false,
    foreach(I in 1..N, Found ==false)
       if ProbDup[I] >= 0.5 then
          println(N=I),
          Found := true
       end
    end
  end,
  nl.


go3 => 
  foreach(N in 1..10000) 
    P = prob2(N),
    printf("%5d= %2.6f ~ %5d\n",N,P,round(P))
  end,
  nl.

prob(N, ProbDup) =>
  Len = ceiling(2*sqrt(N)),
  ProbNoDup = new_list(Len), % probablity of no duplicate
  ProbNoDup[1] := 1.0,
  foreach(I in 2..Len) 
    ProbNoDup[I] := ProbNoDup[I-1] * (N - I +1.0)/N
  end,
  ProbDup = [1.0-ProbNoDup[I] : I in 1..Len].


% From https://en.wikipedia.org/wiki/Birthday_problem
prob2(N) = 1/2+sqrt(1/4 + 2*log(2)*N).