function [ a, b, c ] = newton ( )

%*****************************************************************************80
%
%% newton() simulates Newton's dice problem.
%
%  Discussion:
%
%    Estimate the probability that
%    A: a fair die rolled 6 times will show 6 at least once;
%    B: a fair die rolled 12 times will show 6 at least twice;
%    C: a fair die rolled 18 times will show 6 at least 3 times.
%
%  Licensing:
%
%    This code is distributed under the MIT license. 
%
%  Modified:
%
%    31 January 2018
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Paul Nahin,
%    Will You Be Alive 10 Years From Now?,
%    Princeton, 2014,
%    ISBN: 978-0691156804,
%    LC: QA273.25.N344
%
%  Output:
%
%    real a, b, c: the probabilities of cases A, B, and C.
%
  a = 0;
  b = 0;
  c = 0;
  trial_num = 10000;

  for sim = 1 : trial_num

    roll = randi ( [ 1, 6 ], 6, 1 );
    count = sum ( roll == 6 );
    if ( 1 <= count )
      a = a + 1;
    end

    roll = randi ( [ 1, 6 ], 12, 1 );
    count = sum ( roll == 6 );
    if ( 2 <= count )
      b = b + 1;
    end

    roll = randi ( [ 1, 6 ], 18, 1 );
    count = sum ( roll == 6 );
    if ( 3 <= count )
      c = c + 1;
    end

  end

  a = a / trial_num;
  b = b / trial_num;
  c = c / trial_num;

  return
end
