#! /usr/bin/env python
#
def r8_uni ( s1, s2 ):

#*****************************************************************************80
#
## R8_UNI returns a pseudorandom number between 0 and 1.
#
#  Discussion:
#
#    This function generates uniformly distributed pseudorandom numbers
#    between 0 and 1, using the 32-bit generator from figure 3 of
#    the article by L'Ecuyer.
#
#    The cycle length is claimed to be 2.30584E+18.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    19 August 2015
#
#  Author:
#
#    Original Pascal original version by Pierre L'Ecuyer
#    Python version by John Burkardt
#
#  Reference:
#
#    Pierre LEcuyer,
#    Efficient and Portable Combined Random Number Generators,
#    Communications of the ACM,
#    Volume 31, Number 6, June 1988, pages 742-751.
#
#  Parameters:
#
#    Input/output, integer S1, S2, two values used as the
#    seed for the sequence.  On first call, the user should initialize
#    S1 to a value between 1 and 2147483562;  S2 should be initialized
#    to a value between 1 and 2147483398.
#
#    Output, real VALUE, the next value in the sequence.
#
  k = s1 / 53668
  s1 = 40014 * ( s1 - k * 53668 ) - k * 12211
  if ( s1 < 0 ):
    s1 = s1 + 2147483563

  k = s2 / 52774
  s2 = 40692 * ( s2 - k * 52774 ) - k * 3791
  if ( s2 < 0 ):
    s2 = s2 + 2147483399

  z = s1 - s2
  if ( z < 1 ):
    z = z + 2147483562

  value = float ( z ) / 2147483563.0

  return value, s1, s2

def r8_uni_test01 ( ):

#*****************************************************************************80
#
## R8_UNI_TEST01 tests R8_UNI.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    08 July 2008
#
#  Author:
#
#    John Burkardt
#
  import platform

  print ( '' )
  print ( 'R8_UNI_TEST01' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  R8_UNI computes pseudorandom values.' )
  print ( '  Two seeds, S1 and S2, are used.' )

  s1 = 12345
  s2 = 34567

  print ( '' )
  print ( '      R                     S1            S2' )
  print ( '' )
  print ( '                  %12d  %12d' % ( s1, s2 ) )

  for i in range ( 0, 10 ):
    r, s1, s2 = r8_uni ( s1, s2 )
    print ( '  %14.6g  %12d  %12d' % ( r, s1, s2 ) )
#
#  Terminate.
#
  print ( '' )
  print ( 'R8_UNI_TEST01' )
  print ( '  Normal end of execution.' )
  return

def r8_uni_test02 ( ):

#*****************************************************************************80
#
## R8_UNI_TEST02 tests R8_UNI.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    08 July 2008
#
#  Author:
#
#    John Burkardt
#
  import numpy as np
  import platform

  n = 100000

  print ( '' )
  print ( 'R8_UNI_TEST02' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  Examine the average and variance of a' )
  print ( '  sequence generated by R8_UNI.' )
#
#  Start with known seeds.
#
  s1 = 12345
  s2 = 34567

  print ( '' )
  print ( '  Now compute %d elements.' % ( n ) )

  u = np.zeros ( n )

  u_avg = 0.0
  for i in range ( 0, n ):
    u[i], s1, s2 = r8_uni ( s1, s2 )
    u_avg = u_avg + u[i]

  u_avg = u_avg / float ( n )

  u_var = 0.0
  for i in range ( 0, n ):
    u_var = u_var + ( u[i] - u_avg ) * ( u[i] - u_avg )

  u_var = u_var / float ( n - 1 )

  print ( '' )
  print ( '  Average value = %g' % ( u_avg ) )
  print ( '  Expecting       %g' % ( 0.5 ) )

  print ( '' )
  print ( '  Variance =      %g' % ( u_var ) )
  print ( '  Expecting       %g' % ( 1.0 / 12.0 ) )
#
#  Terminate.
#
  print ( '' )
  print ( 'R8_UNI_TEST02' )
  print ( '  Normal end of execution.' )
  return

def r8_uni_test03 ( ):

#*****************************************************************************80
#
## R8_UNI_TEST03 tests R8_UNI.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    08 July 2008
#
#  Author:
#
#    John Burkardt
#
  import platform

  print ( '' )
  print ( 'R8_UNI_TEST03' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  Show how the seeds used by R8_UNI,' )
  print ( '  which change on each step, can be reset to' )
  print ( '  restore any part of the sequence.' )

  s1_save = 12345
  s2_save = 34567

  s1 = s1_save
  s2 = s2_save

  print ( '' )
  print ( '  Begin sequence with following seeds:' )
  print ( '' )
  print ( '  S1 = %d' % ( s1 ) )
  print ( '  S2 = %d' % ( s2 ) )
  print ( '' )
  print ( '        I      R                     S1            S2' )
  print ( '' )

  for i in range ( 0, 10 ):

    r, s1, s2 = r8_uni ( s1, s2 )
    print ( '  %8d  %14.6g  %12d  %12d' % ( i, r, s1, s2 ) )

    if ( i == 5 ):
      s1_save = s1
      s2_save = s2

  s1 = s1_save
  s2 = s2_save

  print ( '' )
  print ( '  Restart the sequence, using the seeds' )
  print ( '  produced after step 5:' )
  print ( '' )
  print ( '  S1 = %d' % ( s1 ) )
  print ( '  S2 = %d' % ( s2 ) )
  print ( '' )
  print ( '         I      R                     S1            S2' )
  print ( '' )

  for i in range ( 0, 10 ):

    r, s1, s2 = r8_uni ( s1, s2 )
    print ( '  %8d  %14.6g  %12d  %12d' % ( i, r, s1, s2 ) )
#
#  Terminate.
#
  print ( '' )
  print ( 'R8_UNI_TEST03' )
  print ( '  Normal end of execution.' )
  return

if ( __name__ == '__main__' ):
  from timestamp import timestamp
  timestamp ( )
  r8_uni_test01 ( )
  r8_uni_test02 ( )
  r8_uni_test03 ( )
  timestamp ( )