function [ r, seed ] = r8_uniform_01 ( seed ) %*****************************************************************************80 % %% 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 % % Licensing: % % This code is distributed under the GNU LGPL license. % % Modified: % % 21 September 2006 % % 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, integer SEED, the integer "seed" used to generate % the output random number. SEED should not be 0. % % Output, real R, a random value between 0 and 1. % % Output, integer SEED, the updated seed. This would % normally be used as the input seed on the next call. % i4_huge = 2147483647; if ( seed == 0 ) fprintf ( 1, '\n' ); fprintf ( 1, 'R8_UNIFORM_01 - Fatal error!\n' ); fprintf ( 1, ' Input SEED = 0!\n' ); error ( 'R8_UNIFORM_01 - Fatal error!' ); end seed = floor ( seed ); seed = mod ( seed, i4_huge ); if ( seed < 0 ) seed = seed + i4_huge; end k = floor ( seed / 127773 ); seed = 16807 * ( seed - k * 127773 ) - k * 2836; if ( seed < 0 ) seed = seed + i4_huge; end r = seed * 4.656612875E-10; return end