%  R8_UNIFORM_01 returns a unit pseudorandom R8
%  Discussion:
%
%    This routine implements the recursion
%
%      seed = 16807 * seed mod ( 2**31 - 1 )
%      r8_uniform_01 = seed / ( 2**31 - 1 )
%
%    The integer arithmetic never requires more than 32 bits,
%    including a sign bit.
%
%    If the initial seed is 12345, then the first three computations are
%
%      Input     Output      R8_UNIFORM_01
%      SEED      SEED
%
%         12345   207482415  0.096616
%     207482415  1790989824  0.833995
%    1790989824  2035175616  0.947702
%
%  Author:   John Burkardt
%
%  Reference:
%
%    Paul Bratley, Bennett Fox, Linus Schrage,
%    A Guide to Simulation,
%    Springer Verlag, pages 201-202, 1983.
%
%    Pierre L'Ecuyer,
%    Random Number Generation,
%    in Handbook of Simulation,
%    edited by Jerry Banks,
%    Wiley Interscience, page 95, 1998.
%
%    Bennett Fox,
%    Algorithm 647:
%    Implementation and Relative Efficiency of Quasirandom
%    Sequence Generators,
%    ACM Transactions on Mathematical Software,
%    Volume 12, Number 4, pages 362-376, 1986.
%
%    Peter Lewis, Allen Goodman, James Miller,
%    A Pseudo-Random Number Generator for the System/360,
%    IBM Systems Journal,
%    Volume 8, pages 136-143, 1969.
%
%  Parameters:
%
%    Input/output, integer SEED, the "seed" value, which should NOT be 0.
%    On output, SEED has been updated.
%
%    Output, double precision R8_UNIFORM_01, a new pseudorandom variate,
%    strictly between 0 and 1.
 function R8 = r8_uniform_01(seed)
 
      if ( seed == 0 )
        disp('R8_UNIFORM_01 - Fatal error!')
        disp('  Input value of SEED = 0.')
        stop
      end

      k = seed / 127773;

      seed = 16807 * ( seed - k * 127773 ) - k * 2836;

      if ( seed < 0 ) 
        seed = seed + i4_huge;
      end 
%
%  Although SEED can be represented exactly as a 32 bit integer,
%  it generally cannot be represented exactly as a 32 bit real number!
%
      R8 = double(seed) * 4.656612875e-10;
return