function ix = genmul ( n, p, ncat )

%*****************************************************************************80
%
%% genmul() generates a multinomial random deviate.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    03 September 2018
%
%  Author:
%
%    Original FORTRAN77 version by Barry Brown, James Lovato.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Luc Devroye,
%    Non-Uniform Random Variate Generation,
%    Springer, 1986,
%    ISBN: 0387963057,
%    LC: QA274.D48.
%
%  Input:
%
%    integer N, the number of events, which will be
%    classified into one of the NCAT categories.
%
%    real P(NCAT-1).  P(I) is the probability that an event
%    will be classified into category I.  Thus, each P(I) must be between
%    0.0 and 1.0.  Only the first NCAT-1 values of P must be defined since
%    P(NCAT) would be 1.0 minus the sum of the first NCAT-1 P's.
%
%    integer NCAT, the number of categories.
%
%  Output:
%
%    integer IX(NCAT), a random observation from
%    the multinomial distribution.  All IX(i) will be nonnegative and their
%    sum will be N.
%
  if ( n < 0 )
    fprintf ( 1, ' \n' );
    fprintf ( 1, 'GENMUL - Fatal error!\n' );
    fprintf ( 1, '  N < 0\n' );
    error ( 'GENMUL - Fatal error!\n' );
  end

  if ( ncat <= 1 )
    fprintf ( 1, ' \n' );
    fprintf ( 1, 'GENMUL - Fatal error!\n' );
    fprintf ( 1, '  NCAT <= 1\n' );
    error ( 'GENMUL - Fatal error!\n' );
  end

  do i = 1, ncat - 1

    if ( p(i) < 0.0 )
      fprintf ( 1, ' \n' );
      fprintf ( 1, 'GENMUL - Fatal error!\n' );
      fprintf ( 1, '  Some P(i) < 0.\n' );
      error ( 'GENMUL - Fatal error!\n' );
    end

    if ( 1.0 < p(i) )
      fprintf ( 1, ' \n' );
      fprintf ( 1, 'GENMUL - Fatal error!\n' );
      fprintf ( 1, '  Some 1 < P(i).\n' );
      error ( 'GENMUL - Fatal error!\n' );
    end

  end

  ptot = 0.0;
  for i = 1 : ncat - 1
    ptot = ptot + p(i);
  end

  if ( 0.99999 < ptot )
    fprintf ( 1, ' \n' );
    fprintf ( 1, 'GENMUL - Fatal error!\n' );
    fprintf ( 1, '  1 < Sum of P().\n' );
    error ( 'GENMUL - Fatal error!\n' );
  end
%
%  Initialize variables.
%
  ntot = n;
  ptot = 1.0;
  ix = zeros ( ncat, 1 );
%
%  Generate the observation.
%
  for icat = 1 : ncat - 1
    prob = p(icat) / ptot;
    ix(icat) = ignbin ( ntot, prob );
    ntot = ntot - ix(icat);
    if ( ntot <= 0 )
      return
    end
    ptot = ptot - p(icat);
  end

  ix(ncat) = ntot;

  return
end
