function [ an, bn ] = lcrg_anbn ( a, b, c, n )
%*****************************************************************************80
%
%% LCRG_ANBN computes the "N-th power" 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 values of AN and BN
% that describe the LCRG that is equivalent to N applications of the
% original LCRG.
%
% One use for such a facility would be to do random number computations
% in parallel. If each of N processors is to compute many random values,
% you can guarantee that they work with distinct random values
% by starting with a single value of SEED, using the original LCRG to generate
% the first N-1 "iterates" of SEED, so that you now have N "seed" values,
% and from now on, applying the N-th power of the LCRG to the seeds.
%
% If the K-th processor starts from the K-th seed, it will essentially
% be computing every N-th entry of the original random number sequence,
% offset by K. Thus the individual processors will be using a random
% number stream as good as the original one, and without repeating, and
% without having to communicate.
%
% 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-1) = a * SEED(N-2) + 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
% = AN * SEED + BN
%
% Thus, from any SEED, we can determine the result of N applications of the
% original LCRG directly if we can solve
%
% ( a - 1 ) * BN = ( a^N - 1 ) * b in modular arithmetic,
%
% and evaluate:
%
% AN = a^N
%
% Licensing:
%
% This code is distributed under the GNU LGPL license.
%
% Modified:
%
% 21 May 2008
%
% Author:
%
% John Burkardt
%
% Reference:
%
% Barry Wilkinson, Michael Allen,
% Parallel Programming:
% Techniques and Applications Using Networked Workstations and Parallel Computers,
% Prentice Hall,
% ISBN: 0-13-140563-2,
% LC: QA76.642.W54.
%
% 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.
%
% Output:
%
% integer AN, BN, the multiplier and added value for
% the LCRG that represent N applications of the original LCRG.
%
if ( n < 0 )
fprintf ( 1, '\n' );
fprintf ( 1, 'LCRG_ANBN - Fatal error!\n' );
fprintf ( 1, ' Illegal input value of N = %d\n', n );
error ( 'LCRG_ANBN - Fatal error!' );
end
if ( c <= 0 )
fprintf ( 1, '\n' );
fprintf ( 1, 'LCRG_ANBN - Fatal error!\n' );
fprintf ( 1, ' Illegal input value of C = %d\n', c );
error ( 'LCRG_ANBN - Fatal error!' );
end
if ( n == 0 )
an = 1;
bn = 0;
elseif ( n == 1 )
an = a;
bn = b;
else
%
% Compute A^N.
%
an = power_mod ( a, n, c );
%
% Solve
% ( a - 1 ) * BN = ( a^N - 1 ) mod B
% for BN.
%
am1 = a - 1;
anm1tb = ( an - 1 ) * b;
[ bn, ierror ] = congruence ( am1, c, anm1tb );
if ( ierror ~= 0 )
fprintf ( 1, '\n' );
fprintf ( 1, 'LCRG_ANBN - 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_ANGN - Fatal error!' );
end
end
return
end