function value = ignnbn ( n, p )

%*****************************************************************************80
%
%% IGNNBN generates a negative binomial random deviate.
%
%  Discussion:
%
%    This procedure generates a single random deviate from a negative binomial
%    distribution.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 March 2013
%
%  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.
%
%  Parameters:
%
%    Input, integer N, the required number of events.
%    0 <= N.
%
%    Input, real P, the probability of an event during a
%    Bernoulli trial.  0.0 < P < 1.0.
%
%    Output, integer VALUE, a random deviate from
%    the distribution.
%
  if ( n < 0 )
    fprintf ( 1, ' \n' );
    fprintf ( 1, 'IGNNBN - Fatal error!\n' );
    fprintf ( 1, '  N < 0.\n' );
    error ( 'IGNNBN - Fatal error!\n' );
  end

  if ( p <= 0.0 )
    fprintf ( 1, ' \n' );
    fprintf ( 1, 'IGNNBN - Fatal error!\n' );
    fprintf ( 1, '  P <= 0.0\n' );
    error ( 'IGNNBN - Fatal error!\n' );
  end

  if ( 1.0 <= p )
    fprintf ( 1, ' \n' );
    fprintf ( 1, 'IGNNBN - Fatal error!\n' );
    fprintf ( 1, '  1.0 <= P\n' );
    error ( 'IGNNBN - Fatal error!\n' );
  end
%
%  Generate Y, a random gamma (n,(1-p)/p) variable.
%
  r = n;
  a = p / ( 1.0 - p );
  y = gengam ( a, r );
%
%  Generate a random Poisson ( y ) variable.
%
  value = ignpoi ( y );

  return
end
