/* 

  Multiplying dice in Picat.

  From 
  https://enigmaticcode.wordpress.com/2025/10/23/braintwister-96-multiplying-dice/
  """
  BrainTwister #96: Multiplying dice

  From New Scientist #3566, 25th October 2025 [link] [link]
 
    I roll two standard six-sided dice, but instead of adding the values as usual, 
    I multiply them. What is the average value of this product?

    What is the probability of a result that is greater than the average value?

    What is the probability of a greater-than-average result if I roll three dice 
    instead of two dice?

  """


  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.

main => go.

/*

  n = 2
  var : dice prod
  Probabilities (truncated):
  6: 0.1089000000000000
  12: 0.1078000000000000
  4: 0.0819000000000000
  30: 0.0591000000000000
  .........
  1: 0.0292000000000000
  25: 0.0278000000000000
  36: 0.0276000000000000
  16: 0.0275000000000000
  mean = 12.277

  var : ds mean
  Probabilities:
  12.25: 1.0000000000000000
  mean = 12.25

  var : p
  Probabilities:
  false: 0.6371000000000000
  true: 0.3629000000000000
  mean = [false = 0.6371,true = 0.3629]


  n = 3
  var : dice prod
  Probabilities (truncated):
  24: 0.0737000000000000
  12: 0.0698000000000000
  36: 0.0537000000000000
  30: 0.0535000000000000
  .........
  1: 0.0046000000000000
  27: 0.0045000000000000
  125: 0.0043000000000000
  216: 0.0042000000000000
  mean = 42.6181

  var : ds mean
  Probabilities:
  42.875: 1.0000000000000000
  mean = 42.875

  var : p
  Probabilities:
  false: 0.6430000000000000
  true: 0.3570000000000000
  mean = [false = 0.643,true = 0.357]


*/
go ?=>
  member(N,2..3),
  println(n=N),
  reset_store,
  run_model(10_000,$model(N),[show_probs_trunc,mean]),
  nl,
  % show_store_lengths,nl,
  fail,
  nl.
go => true.
  
  
model(N) =>
  Dice = dice_n(6,N),
  DiceProd = Dice.prod,
  if N == 3 then
    Ds = [D1*D2*D3 : D1 in 1..6, D2 in 1..6, D3 in 1..6]
  else
    Ds = [D1*D2 : D1 in 1..6, D2 in 1..6]
  end,
  DsMean = Ds.mean,
  P = check(DiceProd > DsMean),
   
  add("dice prod",DiceProd),
  add("ds mean",DsMean),  
  add("p",P).