#! /usr/bin/env python3
#
def asa183_test ( ):

#*****************************************************************************80
#
## ASA183_TEST tests the ASA183 library.
#
#  Licensing:
#
#    This code is distributed under the MIT license. 
#
#  Modified:
#
#    19 August 2015
#
#  Author:
#
#    John Burkardt
#
  import platform

  print ( '' )
  print ( 'ASA183_TEST' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  Test asa183().' )

  r8_random_test01 ( )
  r8_random_test02 ( )
  r8_random_test03 ( )

  r8_uni_test01 ( )
  r8_uni_test02 ( )
  r8_uni_test03 ( )
#
#  Terminate.
#
  print ( '' )
  print ( 'ASA183_TEST:' )
  print ( '  Normal end of execution.' )

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 MIT 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.
#
#  Input:
#
#    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:
#
#    Output, real VALUE, the next value in the sequence.
#
#    integer S1, S2, S3, updated seed values.
#
  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 MIT 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 MIT 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 MIT 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

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 MIT 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.
#
#  Input:
#
#    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.
#
#    integer S1, S2, updated seed values.
#
  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 MIT 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 MIT 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 MIT 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

def timestamp ( ):

#*****************************************************************************80
#
## TIMESTAMP prints the date as a timestamp.
#
#  Licensing:
#
#    This code is distributed under the MIT license. 
#
#  Modified:
#
#    06 April 2013
#
#  Author:
#
#    John Burkardt
#
  import time

  t = time.time ( )
  print ( time.ctime ( t ) )

  return None

if ( __name__ == '__main__' ):
  timestamp ( )
  asa183_test ( )
  timestamp ( )