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

#*****************************************************************************80
#
## R8_RANDOM returns a pseudorandom number between 0 and 1.
#
#  Discussion:
#
#    This function returns a pseudo-random number rectangularly distributed
#    between 0 and 1.   The cycle length is 6.95E+12.  (See page 123
#    of Applied Statistics (1984) volume 33), not as claimed in the
#    original article.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    19 August 2015
#
#  Author:
#
#    FORTRAN77 original version by Brian Wichman, David Hill.
#    Python version by John Burkardt.
#
#  Reference:
#
#    Brian Wichman, David Hill,
#    Algorithm AS 183: An Efficient and Portable Pseudo-Random
#    Number Generator,
#    Applied Statistics,
#    Volume 31, Number 2, 1982, pages 188-190.
#
#  Parameters:
#
#    Input/output, integer S1, S2, S3, three values used as the
#    seed for the sequence.  These values should be positive
#    integers between 1 and 30,000.
#
#    Output, real VALUE, the next value in the sequence.
#
  s1 = ( ( 171 * s1 ) % 30269 )
  s2 = ( ( 172 * s2 ) % 30307 )
  s3 = ( ( 170 * s3 ) % 30323 )
 
  value = float ( s1 ) / 30269.0 \
        + float ( s2 ) / 30307.0 \
        + float ( s3 ) / 30323.0

  value = ( value % 1.0 )

  return value, s1, s2, s3

def r8_random_test01 ( ):

#*****************************************************************************80
#
## R8_RANDOM_TEST01 tests R8_RANDOM.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    19 August 2015
#
#  Author:
#
#    John Burkardt
#
  import platform

  print ( '' )
  print ( 'R8_RANDOM_TEST01' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  R8_RANDOM computes pseudorandom values.' )
  print ( '  Three seeds, S1, S2, and S3, are used.' )

  s1 = 12345
  s2 = 34567
  s3 = 56789

  print ( '' )
  print ( '      R                 S1        S2        S3' )
  print ( '' )
  print ( '                  %8d  %8d  %8d' % ( s1, s2, s3 ) )

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

def r8_random_test02 ( ):

#*****************************************************************************80
#
## R8_RANDOM_TEST02 tests R8_RANDOM.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    19 August 2015
#
#  Author:
#
#    John Burkardt
#
  import platform
  import numpy as np

  n = 100000

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

  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, s3 = r8_random ( s1, s2, s3 )
    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_RANDOM_TEST02' )
  print ( '  Normal end of execution.' )
  return

def r8_random_test03 ( ):

#*****************************************************************************80
#
## R8_RANDOM_TEST03 tests R8_RANDOM.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    19 August 2015
#
#  Author:
#
#    John Burkardt
#
  import platform

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

  s1_save = 12345
  s2_save = 34567
  s3_save = 56789

  s1 = s1_save
  s2 = s2_save
  s3 = s3_save

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

  for i in range ( 0, 10 ):

    r, s1, s2, s3 = r8_random ( s1, s2, s3 )
    print ( '  %8d  %14.6g  %8d  %8d  %8d' % ( i, r, s1, s2, s3 ) )

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

  s1 = s1_save
  s2 = s2_save
  s3 = s3_save

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

  for i in range ( 0, 10 ):

    r, s1, s2, s3 = r8_random ( s1, s2, s3 )
    print ( '  %8d  %14.6g  %8d  %8d  %8d' % ( i, r, s1, s2, s3 ) )
#
#  Terminate.
#
  print ( '' )
  print ( 'R8_RANDOM_TEST03' )
  print ( '  Normal end of execution.' )
  return

if ( __name__ == '__main__' ):
  from timestamp import timestamp
  timestamp ( )
  r8_random_test01 ( )
  r8_random_test02 ( )
  r8_random_test03 ( )
  timestamp ( )