/* Monte Carlo Simulation of Pi */
% 18 Nov 2002. modified: 19 Nov 2002.

:-write('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%'),nl,
   write(' A Monte Carlo Simulation of Pi'),nl,
   write('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%'),nl,
   write('type sim_pi(mc_trial,1000,P). or so.'),nl.

mc_trial:-
  % By `P is random(11)', with uniform probability,
  % we get a value of P as an element of [0, ..., 10].
  X is random(11) - 5,
  Y is random(11) - 5,
  Z is X^2 + Y^2,
  update_counter,
  (Z =< 25 -> update_pi;true).
update_pi:-
  retract(pi(P)),
  P1 is P + 1,%write(pi(P1)),
  assert(pi(P1)).
update_counter:-
  retract(n(N)),
  N1 is N + 1,%write(n(N1)),
  assert(n(N1)).
pi(0).
n(0).
n_init:-
  abolish(sim_result,2),
  abolish(pi,1),
  abolish(n,1),
  assert(pi(0)),
  assert(n(0)).

sim_pi(_,N,abend):-
  N > 1000,
  write('It will go beyond the self-preservation. Please try a smaller.').
sim_pi(G,N,Pi):-
  G=mc_trial,
  N =< 1000,
  n_init,
  length(X,N),
  forall(
    nth1(K,X,_),
    (
     G,
     sim_result(Pi),
     assert(sim_result(K,Pi)),
     trim_stacks
    )
  ),
  sim_result(N,Pi).
sim_result(Pi):-
  pi(A),
  n(N),
  Pi is A * 400 / N.