function [ value, t ] = r4_random ( t, n )

%*****************************************************************************80
%
%% R4_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:
%
%    25 September 2011
%
%  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.
%
%  Parameters:
%
%    Input, integer N.  The absolute value of N is the number
%    of random numbers in an auxiliary table.  Note though that abs(N)+1 is
%    the number of items in array T.  If N is positive and differs from its
%    value in the previous invocation, then the table is initialized for
%    the new value of N.  If N is negative, abs(N) is the number of items
%    in an auxiliary table, but the tables are now assumed already to
%    be initialized.  This option enables the user to save the table T at
%    the end of a long computer run and to restart with the same sequence.
%    Normally, this function would be called at most once with negative N.
%    Subsequent invocations would have N positive and of the correct magnitude.
%
%    Input/output, real T(abs(N)+1), an array of random numbers
%    from a previous invocation of this function.  Whenever N is positive
%    and differs from the old N, the table is initialized.  The first
%    abs(N) numbers are the table discussed in the reference, and the
%    last value is Y.  This array may be saved in order to restart a sequence.
%
%    Output, real VALUE, a random number between 0.0 and 1.0.
%
  persistent floatn
  persistent nold

  if ( isempty ( floatn ) )
    floatn = -1.0;
    nold = -1;
  end

  if ( n ~= nold )

    nold = abs ( n );
    floatn = nold;

    if ( n < 0 )

      dummy = r4_rand ( t(nold+1) );

    else

      for i = 1 : nold
        t(i) = r4_rand ( 0.0 );
      end
      t(nold+1) = r4_rand ( 0.0 );

    end

  end

  j = r4_aint ( t(nold+1) * floatn + 1.0 );
  t(nold+1) = t(j);
  r4_random = t(j);
  t(j) = r4_rand ( 0.0 );

  return
end