function gs ( trial_num, n )

%*****************************************************************************80
%
%% gs(): the Gamow-Stern elevator problem.
%
%  Discussion:
%
%    A building has 7 floors, and there are n elevators, each of which
%    is at a randomly chosen floor.
%    A person on floor 2 requests an elevator, wishing to go up.
%    What is the probability that the first elevator to arrive is going down?
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    23 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.
%
%    integer N, the number of elevators.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'gs():\n' );
  fprintf ( 1, '  A building has 7 floors, and there are n elevators, each of which\n' );
  fprintf ( 1, '  is at a randomly chosen floor.\n' );
  fprintf ( 1, '  A person on floor 2 requests an elevator, wishing to go up.\n' );
  fprintf ( 1, '  What is the probability that the first elevator to arrive is going down?\n' );

  G = 1 / 6;

  total_going_down = 0;
  elevator = zeros ( n, 4 );

  for loop = 1 : trial_num

    for j = 1 : n

      decision(j) = rand ( );

      if ( decision(j) < 0.5 )
        elevator(j,1) = 0;
      else
        elevator(j,1) = 1;
      end

      elevator(j,2) = rand ( );
      if ( elevator(j,2) < G )
        if ( elevator(j,1) == 0 )
          elevator(j,3) = G + elevator(j,2);
        else
          elevator(j,3) = G - elevator(j,2);
        end
        elevator(j,4) = 1;
      else
        if ( elevator(j,1) == 0 )
          elevator(j,3) = elevator(j,2) - G;
        else
          elevator(j,3) = ( 11 / 6 ) - elevator(j,2);
        end
        elevator(j,4) = 0;
      end

    end

    min = elevator(1,3);
    index = 1;
    for k = 2 : n
      if ( elevator(k,3) < min )
        min = elevator(k,3);
        index = k;
      end
    end

    if ( elevator(index,4) == 0 )
      total_going_down = total_going_down + 1;
    end

  end

  total_going_down = total_going_down / trial_num;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Estimated probability of down elevator = %g\n', total_going_down );
  fprintf ( 1, '  Theoretical probability is %g\n', 0.5 + 0.5 * ( 2/3 ) ^ n );

  return
end
