function key_lcrg = lcrg_key ( a, b, c, n, key )

%*****************************************************************************80
%
%% lcrg_key() computes the N-th key of a linear congruential generator.
%
%  Discussion:
%
%    We are considering a linear congruential random number generator.
%    The LCRG takes as input an integer value called key, and returns
%    an updated value of key, 
%
%      key(out) = a * key(in) + b, mod c.
%
%    and an associated pseudorandom real value
%
%      U = key(out) / c.
%
%    In most cases, a user is content to call the LCRG repeatedly, with
%    the updating of key being taken care of automatically.
%
%    The purpose of this routine is to determine the value of key 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
%    key 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 key, the second with the value of key after 
%    1,000 applications of the LCRG, and so on.
%
%    To evaluate the N-th value of key directly, we start by ignoring 
%    the modular arithmetic, and working out the sequence of calculations
%    as follows:
%
%      key(0) =     key.
%      key(1) = a * key      + b
%      key(2) = a * key(1)   + b = a^2 * key + a * b + b
%      key(3) = a * key(2)   + b = a^3 * key + a^2 * b + a * b + b
%      ...
%      key(N) = a * key(N-1) + b = a^N * key 
%                                    + ( a^(n-1) + a^(n-2) + ... + a + 1 ) * b
%
%    or, using the geometric series,
%
%      key(N) = a^N * key + ( a^N - 1) / ( a - 1 ) * b
%
%    Therefore, we can determine key(N) directly if we can solve
%
%      ( a - 1 ) * BN = ( a^N - 1 ) * b in modular arithmetic,
%
%    and evaluated:
%
%      AN = a^N
%
%    Using the formula:
%
%      key(N) = AN * key + BN, mod c
%
%  Licensing:
%
%    This code is distributed under the MIT license. 
%
%  Modified:
%
%    15 November 2004
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer A, the multiplier for the LCRG.
%
%    integer B, the added value for the LCRG.
%
%    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.
%
%    integer N, the "index", or number of times that the LCRG
%    is to be applied.  It is required that 0 <= N.
%
%    integer key, the starting value of key.  It is customary
%    that 0 < key.
%
%  Output:
%
%    integer key_LCRG, the value of key that would be output
%    if the LCRG were applied to the starting value N times.
%
  if ( n < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'LCRG_key - Fatal error!\n' );
    fprintf ( 1, '  Illegal input value of N = %12d\n', n );
    error ( 'LCRG_key - Fatal error!' )
  end

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

  if ( n == 0 )
    key_lcrg = mod ( key, c );
    if ( key_lcrg < 0 )
      key_lcrg = key_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_key - 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_key - Fatal error!' );
  end
%
%  Set the new key.
%
  value2 = an * key + bn;

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

  key_lcrg = value2;

  return
end
