#! /usr/bin/env python
#
def r8_random ( n, t, ix0, ix1 ):

#*****************************************************************************80
#
## R8_RANDOM is a portable pseudorandom number generator.
#
#  Discussion:
#
#    This random number generator is portable amoung a wide variety of
#    computers.  It generates a random number between 0.0 and 1.0
#    according to the algorithm presented by Bays and Durham.
#
#    The motivation for using this scheme, which resembles the
#    Maclaren-Marsaglia method, is to greatly increase the period of the
#    random sequence.  If the period of the basic generator is P,
#    then the expected mean period of the sequence generated by this
#    generator is given by
#
#      new mean P = sqrt ( pi * factorial ( N ) / ( 8 * P ) ),
#
#    where factorial ( N ) must be much greater than P in this
#    asymptotic formula.  Generally, N should be 16 to maybe 32.
#
#  Modified:
#
#    09 May 2016
#
#  Author:
#
#    Python version by John Burkardt.
#
#  Reference:
#
#    Carter Bays, Stephen Durham,
#    Improving a Poor Random Number Generator,
#    ACM Transactions on Mathematical Software,
#    Volume 2, Number 1, March 1976, pages 59-64.
#
#  Parameters:
#
#    Input, integer N, the number of random numbers in an auxiliary table.  
#
#    Input/output, real T(N+1), an array of random numbers, initialized
#    before first call by R8_RANDOM_INIT.
#
#    Input/output, integer IX0, IX1, two parameters used
#    as seeds for the random number generator.
#
#    Output, real VALUE, a random number between 0.0 and 1.0.
#
  from r8_rand import r8_rand
#
#  Pick J, a random index between 1 and N, and return T(J).
#
  j = int ( t[n] * abs ( n ) )
  t[n] = t[j]
  value = t[j]
#
#  Put a new value into T(J).
#
  r = 0.0
  t[j], ix0, ix1 = r8_rand ( r, ix0, ix1 )

  return value, t, ix0, ix1

def r8_random_init ( n, t, ix0, ix1 ):

#*****************************************************************************80
#
## R8_RANDOM_INIT initializes data for R8_RANDOM.
#
#  Discussion:
#
#    Before calling R8_RANDOM the first time, call R8_RANDOM_INIT
#    in order to initialize the T array.
#
#  Modified:
#
#    09 May 2016
#
#  Author:
#
#    Python version by John Burkardt.
#
#  Reference:
#
#    Carter Bays, Stephen Durham,
#    Improving a Poor Random Number Generator,
#    ACM Transactions on Mathematical Software,
#    Volume 2, Number 1, March 1976, pages 59-64.
#
#  Parameters:
#
#    Input, integer N, the number of random numbers in an auxiliary table.
#
#    Output, real T(N+1), an array of random numbers.
#
#    Input/output, integer IX0, IX1, two parameters used
#    as seeds for the random number generator.  On first call, these
#    might both be set to 0.
#
  import numpy as np
  from r8_rand import r8_rand

  r = 0.0

  t = np.zeros ( n + 1 )

  for i in range ( 0, n + 1 ):
    t[i], ix0, ix1 = r8_rand ( r, ix0, ix1 )

  return t, ix0, ix1

def r8_random_test ( ):

#*****************************************************************************80
#
## R8_RANDOM_TEST tests R8_RANDOM.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    09 May 2016
#
#  Author:
#
#    John Burkardt
#
  import numpy as np
  import platform

  n = 32
  t = np.zeros ( n + 1 )
  i_value = np.array ( [ 1, 2, 3, 4, 10, 100, 1000 ] )

  print ( '' )
  print ( 'R8_RANDOM_TEST:' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  R8_RANDOM is a random number generator.' )
  print ( '' )
  print ( '             I         R8_RANDOM' )
  print ( '' )
#
#  Start the sequence with IX0 = IX1 = 0.
#
  ix0 = 0
  ix1 = 0
  t, ix0, ix1 = r8_random_init ( n, t, ix0, ix1 )

  k = 0

  for i in range ( 1, 1001 ):

    r, t, ix0, ix1 = r8_random ( n, t, ix0, ix1 )

    if ( i == i_value[k] ):
      print ( '  %14d  %14.6g' % ( i, r ) )
      k = k + 1
#
#  Restart the sequence with IX0 = IX1 = 0.
#
  ix0 = 0
  ix1 = 0
  t, ix0, ix1 = r8_random_init ( n, t, ix0, ix1 )

  average = 0.0
  for i in range ( 0, 1000001 ):
    r, t, ix0, ix1 = r8_random ( n, t, ix0, ix1 )
    average = average + r
  average = average / 1000000.0
  print ( '' )
  print ( '       Average =  %14g  %14g' % ( average, 0.5 ) )
#
#  Restart the sequence with IX0 = IX1 = 0.
#
  ix0 = 0
  ix1 = 0
  [ t, ix0, ix1 ] = r8_random_init ( n, t, ix0, ix1 )

  variance = 0.0
  for i in range ( 0, 1000001 ):
    r, t, ix0, ix1 = r8_random ( n, t, ix0, ix1 )
    variance = variance + ( r - average ) ** 2
  variance = variance / 1000000.0
  print ( '       Variance = %14g  %14g' % ( variance, 1.0 / 12.0 ) )
#
#  Terminate.
#
  print ( '' )
  print ( 'R8_RANDOM_TEST:' )
  print ( '  Normal end of execution.' )
  return

if ( __name__ == '__main__' ):
  from timestamp import timestamp
  timestamp ( )
  r8_random_test ( )
  timestamp ( )