/* 

  Infinite dice in Picat.

  https://medium.com/mathadam/this-math-puzzle-looks-hard-fa8ec6e8ed5
  """
  This Math Puzzle Looks Hard
  Infinitely-sided dice produce a surprising result.

  Here’s the challenge. You have three infinitely-sided dice. When you roll one 
  of these dice, you get a Real Number between 0 and 1. So when you roll three of 
  these dice, you get a Real Number between 0 and 3. Capish?

  Here’s the question. You roll the three dice. You square the outcome of each die.

  What are the odds that the sum of the three squares will be less than or equal to 1?
  
  ...

  We need the volume of that chunk (an eight) of the sphere, and the volume of the cube. 
  Divide one by the other, and we’re done.

   1/8*(Volume of Sphere)   1/8*(Pi*4/3 * 1^3)
   ---------------------- = ------------------ = Pi / 6 ~ 0.52
     Volume of Cube              1^3

  The odds are slightly in favour of the result being less than 1.

  Isn’t it awesome that those infinitely-sided dice — which would each be spherical — 
  produce π in this way?
  """

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

main => go.

/*

  pi/6 = 0.523599
  var : p
  Probabilities:
  true: 0.5240500000000000
  false: 0.4759500000000000
  mean = [true = 0.52405,false = 0.47595]

  var : s
  Probabilities (truncated):
  2.872128301291969: 0.0000500000000000
  2.837978069161114: 0.0000500000000000
  2.789398842788187: 0.0000500000000000
  2.788923774895088: 0.0000500000000000
  .........
  0.002598930807224: 0.0000500000000000
  0.002228332843047: 0.0000500000000000
  0.002223421418779: 0.0000500000000000
  0.002046151330282: 0.0000500000000000
  mean = 1.00071


*/
go ?=>
  println("pi/6"=(pi/6)),
  reset_store,
  run_model(20_000,$model,[show_probs_trunc,mean]),
  nl,
  % fail,
  nl.
go => true.
  
  
model() =>
  [A,B,C] = uniform_n(0,1,3),

  S = A**2 + B**2 + C**2,
  P = check(S < 1),

  add("p",P),
  add("s",S).