function 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.
%
%  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
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'NEWTON:\n' );
  fprintf ( 1, '  Estimate the probability that\n' );
  fprintf ( 1, '  A: a fair die rolled 6 times will show 6 at least once;\n' );
  fprintf ( 1, '  B: a fair die rolled 12 times will show 6 at least twice;\n' );
  fprintf ( 1, '  C: a fair die rolled 18 times will show 6 at least 3 times.\n' );

  a = 0;
  b = 0;
  c = 0;
  m = 10000;

  for sim = 1 : m

    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 / m;
  b = b / m;
  c = c / m;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  A: %6.4g  %6.4g\n', a, 0.6651 );
  fprintf ( 1, '  B: %6.4g  %6.4g\n', b, 0.6187 );
  fprintf ( 1, '  C: %6.4g  %6.4g\n', c, 0.5973 );

  return
end
