function [ value, t, ix0, ix1 ] = r8_random ( n, t, ix0, ix1 )

%*****************************************************************************80
%
%% R8_RANDOM is a portable pseudorandom number generator.
%
%  Discussion:
%
%    This random number generator is portable amoung a wide variety of
%    computers.  It generates a random number between 0.0 and 1.0
%    according to the algorithm presented by Bays and Durham.
%
%    The motivation for using this scheme, which resembles the
%    Maclaren-Marsaglia method, is to greatly increase the period of the
%    random sequence.  If the period of the basic generator is P,
%    then the expected mean period of the sequence generated by this
%    generator is given by
%
%      new mean P = sqrt ( pi * factorial ( N ) / ( 8 * P ) ),
%
%    where factorial ( N ) must be much greater than P in this
%    asymptotic formula.  Generally, N should be 16 to maybe 32.
%
%  Modified:
%
%    09 May 2016
%
%  Author:
%
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Carter Bays, Stephen Durham,
%    Improving a Poor Random Number Generator,
%    ACM Transactions on Mathematical Software,
%    Volume 2, Number 1, March 1976, pages 59-64.
%
%  Input:
%
%    integer N, the number of random numbers in an auxiliary table.  
%
%    real T(N+1), an array of random numbers, initialized
%    before first call by R8_RANDOM_INIT.
%
%    integer IX0, IX1, two parameters used
%    as seeds for the random number generator.
%
%  Output:
%
%    real VALUE, a random number between 0.0 and 1.0.
%
%    real T(N+1), an updated array of random numbers.
%
%    integer IX0, IX1, updated seeds for the random number generator.
%

%
%  Pick J, a random index between 1 and N, and return T(J).
%
  j = floor ( t(n+1) * abs ( n ) + 1.0 );
  t(n+1) = t(j);
  value = t(j);
%
%  Put a new value into T(J).
%
  r = 0.0;
  [ t(j), ix0, ix1 ] = r8_rand ( r, ix0, ix1 );

  return
end
 
