/* 

  Sultan's dowry probabilities in Picat.

  From https://mathworld.wolfram.com/SultansDowryProblem.html
  """
  More specifically, the probability of obtaining the maximum dowry
  after waiting for x out of n daughters is
     P(n,x)=(1 + x*(H_(n-1)-H_x))  / n

  where H_n is a harmonic number (Havil 2003, p. 137) [...]

  the solution is therefore the smallest x such that
   H_x>=H_n-1. 	

  Solving,
    H_x=H_n-1 	
  (4)

  numerically and taking the ceiling function [x] then gives the solutions 
  0, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 5, 5, 5, ... (OEIS A054404) for n=1, 2, ... daughters. 
  """


  Cf my Gamble model gamble_sultans_dowry_probabilities.rkt

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

*/

import ppl_distributions, ppl_utils.
import util.
% import ordset.

main => go.

table
h_n(X) = harmonic_number(X).

prob(N,X) = Res =>
  Res = (1 + (X * (h_n(N-1) - h_n(X) ))) / N.

find_best(N) = Res =>
  OK = false,
  foreach(X in 1..N,break(OK != false))
    if h_n(X) >= h_n(N) - 1 then
      OK := X
    end
  end,
  Res = cond(OK == false,0,OK).

/*


*/
go ?=>
  println("prob(1..100):"),
  L = [prob(100,X) : X in 1..100],
  foreach(I in 1..L.len)
    println(I=L[I])
  end,
  nl,
  println(argmax=argmax(L)),
  println(find_best_1000=find_best(1000)),
  nl,
  println("For 0..20 daughters:"),
  foreach(N in 0..20)
    println(N=find_best(N)),
  end,
  nl.
go => true.