function [ x, seed ] = gamma_sample ( a, b, c, seed )

%*****************************************************************************80
%
%% GAMMA_SAMPLE samples the Gamma PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 October 2004
%
%  Author:
%
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    J H Ahrens and U Dieter,
%    Generating Gamma Variates by a Modified Rejection Technique,
%    Communications of the ACM,
%    Volume 25, Number 1, January 1982, pages 47 - 54.
%
%    J H Ahrens and U Dieter,
%    Computer Methods for Sampling from Gamma, Beta, Poisson and
%      Binomial Distributions.
%    Computing, Volume 12, 1974, pages 223 - 246.
%
%    J H Ahrens, K D Kohrt, and U Dieter,
%    Algorithm 599,
%    ACM Transactions on Mathematical Software,
%    Volume 9, Number 2, June 1983, pages 255-257.
%
%  Parameters:
%
%    Input, real A, B, C, the parameters of the PDF.
%    0.0 < B,
%    0.0 < C.
%
%    Input, integer SEED, a seed for the random number generator.
%
%    Output, real X, a sample of the PDF.
%
%    Output, integer SEED, an updated seed for the random number generator.
%
  a1 =   0.3333333;
  a2 = - 0.2500030;
  a3 =   0.2000062;
  a4 = - 0.1662921;
  a5 =   0.1423657;
  a6 = - 0.1367177;
  a7 =   0.1233795;
  e1 = 1.0;
  e2 = 0.4999897;
  e3 = 0.1668290;
  e4 = 0.0407753;
  e5 = 0.0102930;
  euler = 2.71828182845904;
  q1 =   0.04166669;
  q2 =   0.02083148;
  q3 =   0.00801191;
  q4 =   0.00144121;
  q5 = - 0.00007388;
  q6 =   0.00024511;
  q7 =   0.00024240;
%
%  Allow C = 0.
%
  if ( c == 0.0 )
    x = a;
    return
  end
%
%  C < 1.
%
  if ( c < 1.0 )

    while ( 1 )

      [ u, seed ] = r8_uniform_01 ( seed );
      t = 1.0 + c / euler;
      p = u * t;

      [ s, seed ] = exponential_01_sample ( seed );

      if ( p < 1.0 )
        x = exp ( log ( p ) / c );
        if ( x <= s )
          break
        end
      else
        x = -log ( ( t - p ) / c );
        if ( ( 1.0 - c ) * log ( x ) <= s )
          break
        end
      end

    end

    x = a + b * x;
    return
%
%  1 <= C.
%
  else

    s2 = c - 0.5;
    s = sqrt ( c - 0.5 );
    d = sqrt ( 32.0 ) - 12.0 * sqrt ( c - 0.5 );

    [ t, seed ] = normal_01_sample ( seed );
    x = ( sqrt ( c - 0.5 ) + 0.5 * t )^2;

    if ( 0.0 <= t )
      x = a + b * x;
      return
    end

    [ u, seed ] = r8_uniform_01 ( seed );

    if ( d * u <= t^3 )
      x = a + b * x;
      return
    end

    r = 1.0 / c;

    q0 = ( ( ( ( ( ( ...
           q7   * r ...
         + q6 ) * r ...
         + q5 ) * r ...
         + q4 ) * r ...
         + q3 ) * r ...
         + q2 ) * r ...
         + q1 ) * r;

    if ( c <= 3.686 )
      bcoef = 0.463 + s - 0.178 * s2;
      si = 1.235;
      co = 0.195 / s - 0.079 + 0.016 * s;
    elseif ( c <= 13.022 ) 
      bcoef = 1.654 + 0.0076 * s2;
      si = 1.68 / s + 0.275;
      co = 0.062 / s + 0.024;
    else
      bcoef = 1.77;
      si = 0.75;
      co = 0.1515 / s;
    end

    if ( 0.0 < sqrt ( c - 0.5 ) + 0.5 * t )

      v = 0.5 * t / s;

      if ( 0.25 < abs ( v ) )
        q = q0 - s * t + 0.25 * t * t + 2.0 * s2 * log ( 1.0 + v );
      else
        q = q0 + 0.5 * t * t * ( ( ( ( ( ( ...
               a7   * v ...
             + a6 ) * v ...
             + a5 ) * v ...
             + a4 ) * v ...
             + a3 ) * v ...
             + a2 ) * v ...
             + a1 ) * v;
      end

      if ( log ( 1.0 - u ) <= q )
        x = a + b * x;
        return
      end

    end

    while ( 1 )

      [ e, seed ] = exponential_01_sample ( seed );

      [ u, seed ] = r8_uniform_01 ( seed );

      u = 2.0 * u - 1.0;
      if ( u < 0.0 )
        t = bcoef - si * e;
      else
        t = bcoef + si * e;
      end

      if ( -0.7187449 <= t )

        v = 0.5 * t / s;

        if ( 0.25 < abs ( v ) )
          q = q0 - s * t + 0.25 * t * t + 2.0 * s2 * log ( 1.0 + v );
        else
          q = q0 + 0.5 * t * t * ( ( ( ( ( ( ...
               a7   * v ...
             + a6 ) * v ...
             + a5 ) * v ...
             + a4 ) * v ...
             + a3 ) * v ...
             + a2 ) * v ...
             + a1 ) * v;
        end

        if ( 0.0 < q )

          if ( 0.5 < q )
            w = exp ( q ) - 1.0;
          else
            w = ( ( ( ( ...
                    e5   * q ...
                  + e4 ) * q ...
                  + e3 ) * q ...
                  + e2 ) * q ...
                  + e1 ) * q;
          end

          if ( co * abs ( u ) <= w * exp ( e - 0.5 * t * t ) )
            x = a + b * ( s + 0.5 * t )^2;
            return
          end

        end

      end

    end

  end

  return
end