function [ value, ix0, ix1 ] = r8_rand ( r, ix0, ix1 )

%*****************************************************************************80
%
%% R8_RAND is a portable pseudorandom number generator.
%
%  Discussion:
%
%    This pseudo-random number generator is portable amoung a wide
%    variety of computers.  It is undoubtedly not as good as many
%    readily available installation dependent versions, and so this
%    routine is not recommended for widespread usage.  Its redeeming
%    feature is that the exact same random numbers (to within final round-
%    off error) can be generated from machine to machine.  Thus, programs
%    that make use of random numbers can be easily transported to and
%    checked in a new environment.
%
%    The random numbers are generated by the linear congruential
%    method described by Knuth in seminumerical methods (p.9),
%    addison-wesley, 1969.  Given the i-th number of a pseudo-random
%    sequence, the i+1 -st number is generated from
%      x(i+1) = (a*x(i) + c) mod m,
%    where here m = 2^22 = 4194304, c = 1731 and several suitable values
%    of the multiplier a are discussed below.  Both the multiplier a and
%    random number x are represented in real as two 11-bit
%    words.  The constants are chosen so that the period is the maximum
%    possible, 4194304.
%
%    In order that the same numbers be generated from machine to
%    machine, it is necessary that 23-bit integers be reducible modulo
%    2^11 exactly, that 23-bit integers be added exactly, and that 11-bit
%    integers be multiplied exactly.  Furthermore, if the restart option
%    is used (where r is between 0 and 1), then the product r*2^22 =
%    r*4194304 must be correct to the nearest integer.
%
%    The first four random numbers should be
%
%      0.0004127026,
%      0.6750836372,
%      0.1614754200,
%      0.9086198807.
%
%    The tenth random number is
%
%      0.5527787209.
%
%    The hundredth random number is
%
%      0.3600893021.
%
%    The thousandth number should be
%
%      0.2176990509.
%
%    In order to generate several effectively independent sequences
%    with the same generator, it is necessary to know the random number
%    for several widely spaced calls.  The I-th random number times 2^22,
%    where I=K*P/8 and P is the period of the sequence (P = 2^22), is
%    still of the form L*P/8.  In particular we find the I-th random
%    number multiplied by 2^22 is given by
%      I   =  0  1*p/8  2*p/8  3*p/8  4*p/8  5*p/8  6*p/8  7*p/8  8*p/8
%      RAND=  0  5*p/8  2*p/8  7*p/8  4*p/8  1*p/8  6*p/8  3*p/8  0
%    thus the 4*P/8 = 2097152 random number is 2097152/2**22.
%
%    Several multipliers have been subjected to the spectral test
%    (see Knuth, p. 82).  Four suitable multipliers roughly in order of
%    goodness according to the spectral test are
%      3146757 = 1536*2048 + 1029 = 2^21 + 2^20 + 2^10 + 5
%      2098181 = 1024*2048 + 1029 = 2^21 + 2^10 + 5
%      3146245 = 1536*2048 +  517 = 2^21 + 2^20 + 2^9 + 5
%      2776669 = 1355*2048 + 1629 = 5^9 + 7^7 + 1
%
%    In the table below log10(NU(I)) gives roughly the number of
%    random decimal digits in the random numbers considered I at a time.
%    C is the primary measure of goodness.  In both cases bigger is better.
%
%                     log10 nu(i)              c(i)
%         a       i=2  i=3  i=4  i=5    i=2  i=3  i=4  i=5
%
%      3146757    3.3  2.0  1.6  1.3    3.1  1.3  4.6  2.6
%      2098181    3.3  2.0  1.6  1.2    3.2  1.3  4.6  1.7
%      3146245    3.3  2.2  1.5  1.1    3.2  4.2  1.1  0.4
%      2776669    3.3  2.1  1.6  1.3    2.5  2.0  1.9  2.6
%     best
%      possible   3.3  2.3  1.7  1.4    3.6  5.9  9.7  14.9
%
%    Note that the original version of this routine used local static
%    variables IX0 and IX1.  In this revised version, IX0 and IX1 are
%    routine arguments.  To duplicate the behavior of the original code,
%    IX0 and IX1 should be set to zero before the first call.  
%    JVB, 04 May 2016.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    04 May 2016
%
%  Author:
%
%    MATLAB version by John Burkardt.
%
%  Input:
%
%    real R, determines the action.
%    * R = 0.0, the next random number of the sequence is generated.
%    * R < 0.0, the last generated number will be returned for
%    possible use in a restart procedure.
%    * R > 0.0, the sequence of random numbers will start with the
%    seed ( R mod 1 ).  This seed is also returned as the value of
%    R8_RAND provided the arithmetic is done exactly.
%
%    integer IX0, IX1, two parameters used
%    as seeds for the random number generator.  On first call, these
%    might both be set to 0.  Thereafter, they should simply have the
%    values output from the previous call.
%
%  Output:
%
%    real VALUE, a pseudo-random number between 0.0 and 1.0.
%
%    integer IX0, IX1, updated seeds for the random number generator.
%
  ia0 = 1029;
  ia1 = 1536;
  ia1ma0 = 507;
  ic = 1731;

  if ( r == 0.0 )
    iy0 = ia0 * ix0;
    iy1 = ia1 * ix1 + ia1ma0 * ( ix0 - ix1 ) + iy0;
    iy0 = iy0 + ic;
    ix0 = mod ( iy0, 2048 );
    iy1 = iy1 + ( iy0 - ix0 ) / 2048;
    ix1 = mod ( iy1, 2048 );
  end

  if ( 0.0 < r )
    ix1 = floor ( mod ( r, 1.0 ) * 4194304.0 + 0.5 );
    ix0 = mod ( ix1, 2048 );
    ix1 = ( ix1 - ix0 ) / 2048;
  end

  value = ix1 * 2048 + ix0;
  value = value / 4194304.0;

  return
end
