function kelvin ( trial_num, p )

%*****************************************************************************80
%
%% kelvin() looks at Kelvin's fair results from a biased coin.
%
%  Discussion:
%
%    A biased coin comes up head with probability p.
%    To get a unbiased random value, toss the coin twice.
%    If you get TH, call it heads; if you get HT, call it tails.
%    If you get TT or HH, do another double toss.
%    On average, how many double tosses are necessary?
%
%  Example:
%
%     10,000 trials, average number of double tosses is 2.085.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    28 January 2018
%
%  Author:
%
%    Paul Nahin
%
%  Reference:
%
%    Paul Nahin,
%    Digital Dice: Computational Solutions to Practical Probability Problems,
%    Princeton, 2008,
%    ISBN: 978-0-691-15821-1.
%
%  Input:
%
%    integer trial_num: the number of trials.
%
%    real p: probability of heads.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'kelvin():\n' );
  fprintf ( 1, '  A biased coin comes up heads with probability %g.\n', p );
  fprintf ( 1, '  To get a unbiased random value, toss the coin twice.\n' );
  fprintf ( 1, '  If you get TH, call it heads; if you get HT, call it tails.\n' );
  fprintf ( 1, '  If you get TT or HH, do another double toss.\n' );
  fprintf ( 1, '  On average, how many double tosses are necessary?\n' );

  total_number_of_double_tosses = 0;

  for loop = 1 : trial_num

    number_of_double_tosses = 0;
    decision = 0;

    while ( decision == 0 )

      if ( rand ( ) < p )
        toss1 = 1;
      else
        toss1 = 0;
      end
      if ( rand ( ) < p )
        toss2 = 1;
      else
        toss2 = 0;
      end

      number_of_double_tosses = number_of_double_tosses + 1;

      if ( toss1 + toss2 == 1 )
        decision = 1;
      end    

    end

    total_number_of_double_tosses = total_number_of_double_tosses ...
      + number_of_double_tosses;

  end

  total_number_of_double_tosses = total_number_of_double_tosses / trial_num;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Average number of double tosses = %g\n', ...
    total_number_of_double_tosses );
  fprintf ( 1, '  Theoretical value               = %g\n', 1 / 0.48 );

  return
end
