/* 

  Three dice product in Picat.

  From https://x.com/StudyMaths/status/1938559419719033125
  """
  Why is one of the boxes so much harder to fill than the rest?


  Roll a set of three dice 20 times and calculate their product.
  Record the results in the Carrol diagram

                 Odd product           Even product
                 -----------------------------------
  Less than 40   |0.09259259259259259  0.5185185185185185
                 |
  40 or more     |0.032407407407407460 0.35648148148148145
  """

  1. < 40: 0.6111111111111112  >40: 0.3888888888888889
  2. odd : 1/8                 even: 7/8                  

  Cf my Gamble model gamble_three_dice_product.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 : roll
  Probabilities (truncated):
  24: 0.0610000000000000
  36: 0.0590000000000000
  30: 0.0580000000000000
  12: 0.0560000000000000
  .........
  64: 0.0050000000000000
  216: 0.0040000000000000
  1: 0.0040000000000000
  125: 0.0030000000000000
  mean = 42.477

  var : < 40 odd
  Probabilities:
  false: 0.9010000000000000
  true: 0.0990000000000000
  mean = 0.099

  var : < 40 even
  Probabilities:
  true: 0.5130000000000000
  false: 0.4870000000000000
  mean = 0.513

  var : >= 40 odd
  Probabilities:
  false: 0.9660000000000000
  true: 0.0340000000000000
  mean = 0.034

  var : >= 40 even
  Probabilities:
  false: 0.6460000000000000
  true: 0.3540000000000000
  mean = 0.354

  var : < 40
  Probabilities:
  true: 0.6120000000000000
  false: 0.3880000000000000
  mean = 0.612

  var : >= 40
  Probabilities:
  false: 0.6120000000000000
  true: 0.3880000000000000
  mean = 0.388

  var : odd
  Probabilities:
  false: 0.8670000000000000
  true: 0.1330000000000000
  mean = 0.133

  var : even
  Probabilities:
  true: 0.8670000000000000
  false: 0.1330000000000000
  mean = 0.867


        Odd product  Even product
   < 40  0.0990        0.513000
  >= 40  0.0340        0.354000


*/
go ?=>
  reset_store,
  run_model(1_000,$model,[show_probs_trunc,mean]),
  nl,
  Store = get_store(),
  Lt40Odd = Store.get("< 40 odd").mean,
  Lt40Even = Store.get("< 40 even").mean,  
  Ge40Odd = Store.get(">= 40 odd").mean,
  Ge40Even = Store.get(">= 40 even").mean,  

  println("        Odd product  Even product"),
  printf(" < 40  %.4f        %.6f\n",Lt40Odd,Lt40Even),
  printf(">= 40  %.4f        %.6f\n",Ge40Odd,Ge40Even),  

  % show_store_lengths,nl,
  % fail,
  nl.
go => true.
  
  
model() =>
  N = 3,

  Roll = dice_n(6,N).prod,

  Lt40Odd  = check( (Roll <  40, Roll mod 2 == 1)),
  Lt40Even = check( (Roll <  40, Roll mod 2 == 0)),
  Ge40Odd  = check( (Roll >= 40, Roll mod 2 == 1)),
  Ge40Even = check( (Roll >= 40, Roll mod 2 == 0)),    

  Lt40 = check(Roll < 40),
  Ge40 = check(Roll >= 40),

  Odd  = check(Roll mod 2 == 1),
  Even = check(Roll mod 2 == 0),  

  add("roll",Roll),
  add("< 40 odd",Lt40Odd),
  add("< 40 even",Lt40Even),  
  add(">= 40 odd",Ge40Odd),
  add(">= 40 even",Ge40Even),  
  add("< 40",Lt40),
  add(">= 40",Ge40),  
  add("odd",Odd),
  add("even",Even).