/* 

  Landing on 25 in Picat.

  From BL "You Are Standing On The First Step Of An Infinitely Long Numbered Path -
  How likely will you step on step 25?"
  https://medium.com/puzzle-sphere/you-are-standing-on-the-first-step-of-an-infinitely-long-numbered-path-9e3173df242d
  """
  You have a fair coin which has the number 1 written on one side and the number 2 
  on the other. You throw the coin, and if it comes up N, you would then take N 
  steps to the right. For example, if you throw the coin and it comes up 2, you take 
  2 steps to the right to land on step number 3. You now repeat the exercise, throwing 
  the coin again and walking the number of steps that comes up on the coin. 
  If you throw the coin 24 times, you are certain to have landed on, or past, 
  step number 25.

  What is the probability that at some point you will land on step number 25?

  ...

  So the probability ... would be:
   u24 = 2/3 = 1/3 * (-1/2)^24 ~ 0.6666666865
  """

  The close form formula:
  Picat> X = 2/3 + (1/3 * (-1/2)**24)
  X = 0.666666686534882

  Cf my Gamble model gamble_landing_on_25.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.

/*

  var : len
  Probabilities:
  16: 0.2932000000000000
  17: 0.2445000000000000
  15: 0.2085000000000000
  18: 0.1158000000000000
  14: 0.0791500000000000
  19: 0.0375000000000000
  13: 0.0108500000000000
  20: 0.0090500000000000
  21: 0.0009500000000000
  12: 0.0003500000000000
  22: 0.0001500000000000
  mean = 16.2297

  var : tot
  Probabilities:
  25: 0.6677999999999999
  26: 0.3322000000000000
  mean = 25.3322

  var : p
  Probabilities:
  true: 0.6677999999999999
  false: 0.3322000000000000
  mean = 0.6678

*/
go ?=>
  reset_store,
  run_model(20_000,$model,[show_probs,mean]),
  nl,
  % show_store_lengths,nl,
  % fail,
  nl.
go => true.
  

f(A,Tot,Target) = Res =>
  if Tot >= Target then
    Res = [A,Tot]
  else
    C = uniform_draw([1,2]),
    Res = f(A++[C],Tot+C,Target)
  end.


model() =>
  Target = 25,

  [L,Tot] = f([],1,Target), % start at number 1
  Len = L.len,

  P = check(Tot == Target),

  add("len",Len),
  add("tot",Tot),
  add("p",P).