function seed_lcrg = lcrg_seed ( a, b, c, n, seed )

%% LCRG_SEED computes the N-th seed of a linear congruential generator.
%
%  Discussion:
%
%    We are considering a linear congruential random number generator.
%    The LCRG takes as input an integer value called SEED, and returns
%    an updated value of SEED, 
%
%      SEED(out) = a * SEED(in) + b, mod c.
%
%    and an associated pseudorandom real value
%
%      U = SEED(out) / c.
%
%    In most cases, a user is content to call the LCRG repeatedly, with
%    the updating of SEED being taken care of automatically.
%
%    The purpose of this routine is to determine the value of SEED that
%    would be output after N successive applications of the LCRG.  This
%    allows the user to know, in advance, what the 1000-th value of
%    SEED would be, for instance.  Obviously, one way to do this is to
%    apply the LCRG formula 1,000 times.  However, it is possible to
%    do this in a more direct and efficient way.
%
%    One use for such a facility would be to do random number computations
%    in parallel.  If each processor is to compute 1,000 values, you can
%    guarantee that they work with distinct random values by starting the
%    first processor with SEED, the second with the value of SEED after 
%    1,000 applications of the LCRG, and so on.
%
%    To evaluate the N-th value of SEED directly, we start by ignoring 
%    the modular arithmetic, and working out the sequence of calculations
%    as follows:
%
%      SEED(0) =     SEED.
%      SEED(1) = a * SEED      + b
%      SEED(2) = a * SEED(1)   + b = a^2 * SEED + a * b + b
%      SEED(3) = a * SEED(2)   + b = a^3 * SEED + a^2 * b + a * b + b
%      ...
%      SEED(N) = a * SEED(N-1) + b = a^N * SEED 
%                                    + ( a^(n-1) + a^(n-2) + ... + a + 1 ) * b
%
%    or, using the geometric series,
%
%      SEED(N) = a^N * SEED + ( a^N - 1) / ( a - 1 ) * b
%
%    Therefore, we can determine SEED(N) directly if we can solve
%
%      ( a - 1 ) * BN = ( a^N - 1 ) * b in modular arithmetic,
%
%    and evaluated:
%
%      AN = a^N
%
%    Using the formula:
%
%      SEED(N) = AN * SEED + BN, mod c
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license. 
%
%  Modified:
%
%    15 November 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer A, the multiplier for the LCRG.
%
%    Input, integer B, the added value for the LCRG.
%
%    Input, integer C, the base for the modular arithmetic.  For 32 bit
%    arithmetic, this is often 2^31 - 1, or 2147483647.  It is required
%    that 0 < C.
%
%    Input, integer N, the "index", or number of times that the LCRG
%    is to be applied.  It is required that 0 <= N.
%
%    Input, integer SEED, the starting value of SEED.  It is customary
%    that 0 < SEED.
%
%    Output, integer SEED_LCRG, the value of SEED that would be output
%    if the LCRG were applied to the starting value N times.
%
  if ( n < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'LCRG_SEED - Fatal error!\n' );
    fprintf ( 1, '  Illegal input value of N = %12d\n', n );
    error ( 'LCRG_SEED - Fatal error!' )
  end

  if ( c <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'LCRG_SEED - Fatal error!\n' );
    fprintf ( 1, '  Illegal input value of C = %12d\n', c );
    error ( 'LCRG_SEED - Fatal error!' )
  end

  if ( n == 0 )
    seed_lcrg = mod ( seed, c );
    if ( seed_lcrg < 0 )
      seed_lcrg = seed_lcrg + c;
    end
    return
  end
%
%  Get A^N.
%
  an = power_mod ( a, n, c );
%
%  Solve ( a - 1 ) * BN = ( a^N - 1 ) for BN.
%
%  The LCRG I have been investigating uses B = 0, so this code
%  has not been properly tested yet.
%
  [ bn, ierror ] = congruence ( a-1, c, ( an - 1 ) * b );

  if ( ierror ~= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'LCRG_SEED - Fatal error!\n' );
    fprintf ( 1, '  An error occurred in the CONGRUENCE routine.\n' );
    fprintf ( 1, '  The error code was IERROR = %d\n', ierror );
    error ( 'LCRG_SEED - Fatal error!' );
  end
%
%  Set the new SEED.
%
  value2 = an * seed + bn;

  value2 = mod ( value2, c );
%
%  Guarantee that the value is positive.
%
  if ( value2 < 0 )
    value2 = value2 + c;
  end

  seed_lcrg = value2;

  return
end