#! /usr/bin/env python3
#
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

def vdc ( i ):

#*****************************************************************************80
#
## vdc() computes an element of the van der Corput sequence.
#
#  Discussion:
#
#    The van der Corput sequence is often used to generate a "subrandom"
#    sequence of points which have a better covering property
#    than pseudorandom points.
#
#    The van der Corput sequence generates a sequence of points in [0,1]
#    which never repeats.  The elements of the van der Corput sequence 
#    are strictly less than 1.
#
#    The van der Corput sequence writes an integer in a given base 2,
#    and then its digits are "reflected" about the decimal point.
#    This maps the numbers from 1 to N into a set of numbers in [0,1],
#    which are especially nicely distributed if N is one less
#    than a power of the base.
#
#    The generation is quite simple.  Given an integer I, the expansion
#    of I in base 2 is generated.  Then, essentially, the result R
#    is generated by writing a decimal point followed by the digits of
#    the expansion of I, in reverse order.  This decimal value is actually
#    still in base 2, so it must be properly interpreted to generate
#    a usable value.
#
#  Example:
#
#    I        I         van der Corput
#    decimal  binary    binary   decimal
#    -------  ------    ------   -------
#        0  =     0  =>  .0     = 0.0
#        1  =     1  =>  .1     = 0.5
#        2  =    10  =>  .01    = 0.25
#        3  =    11  =>  .11    = 0.75
#        4  =   100  =>  .001   = 0.125
#        5  =   101  =>  .101   = 0.625
#        6  =   110  =>  .011   = 0.375
#        7  =   111  =>  .111   = 0.875
#        8  =  1000  =>  .0001  = 0.0625
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    07 June 2022
#
#  Author:
#
#    John Burkardt
#
#  Reference:
#
#    John Halton,
#    On the efficiency of certain quasi-random sequences of points
#    in evaluating multi-dimensional integrals,
#    Numerische Mathematik,
#    Volume 2, pages 84-90, 1960.
#
#    John Hammersley,
#    Monte Carlo methods for solving multivariable problems,
#    Proceedings of the New York Academy of Science,
#    Volume 86, pages 844-874, 1960.
#
#    Johannes van der Corput,
#    Verteilungsfunktionen I & II,
#    Nederl. Akad. Wetensch. Proc.,
#    Volume 38, 1935, pages 813-820, pages 1058-1066.
#
#  Input:
#
#    integer I, the index of the element of the sequence.
#    I = 0 is allowed, and returns R = 0.
#
#  Output:
#
#    real R, the I-th element of the van der Corput sequence.
#

#
#  Isolate the sign, and only work with the integer part of I.
#
  if ( i < 0 ):
    s = -1
  else:
    s = +1

  t = abs ( int ( i ) )
#
#  Carry out the computation.
#
  base_inv = 0.5

  r = 0.0

  while ( t != 0 ):
    d = ( t % 2 )
    r = r + d * base_inv
    base_inv = base_inv / 2.0
    t = ( t // 2 )
#
#  Recover the sign.
#
  r = r * s

  return r

def vdc_test ( ):

#*****************************************************************************80
#
## vdc_test() tests vdc().
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    07 June 2022
#
#  Author:
#
#    John Burkardt
#
  print ( '' )
  print ( 'vdc_test():' )
  print ( '  vdc() returns the I-th element of a van der Corput sequence.' )
  print ( '' )
  print ( '    I          vdc(I)' )
  print ( '' )

  for i in range ( -10, 11 ):
    r = vdc ( i )
    print ( '  %3d  %14.8f' % ( i, r ) )

  return

def vdc_base ( i, b ):

#*****************************************************************************80
#
## vdc_base() computes an element of the van der Corput sequence in any base.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    07 June 2022
#
#  Author:
#
#    John Burkardt
#
#  Input:
#
#    integer I, the index of the element of the sequence.
#    I = 0 is allowed, and returns R = 0.
#
#    integer B, the base of the sequence.  The standard sequence
#    uses B = 2, and this function expects 2 <= B.
#
#  Output:
#
#    real R, the I-th element of the van der Corput sequence.
#

#
#  2 <= B.
#
  if ( b < 2 ):
    print ( '' )
    print ( 'vdc_base(): Fatal error!' )
    print ( '  2 <= B is required.' )
    raise Exception ( 'vdc_base(): Fatal error!' )
#
#  B should be an integer.
#
  b = int ( b )
#
#  Isolate the sign, and only work with the integer part of I.
#
  if ( i < 0 ):
    s = -1
  else:
    s = +1

  t = abs ( int ( i ) )
#
#  Carry out the computation.
#
  base_inv = 1.0 / b

  r = 0.0

  while ( t != 0 ):
    d = ( t % b )
    r = r + d * base_inv
    base_inv = base_inv / b
    t = ( t // b )
#
#  Recover the sign.
#
  r = r * s

  return r

def vdc_base_test ( ):

#*****************************************************************************80
#
## vdc_base_test() tests vdc_base().
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    07 June 2022
#
#  Author:
#
#    John Burkardt
#
  print ( '' )
  print ( 'vdc_base_test():' )
  print ( '  vdc_base() returns the I-th element of a van der Corput' )
  print ( '  sequence in base B.' )
  print ( '' )
  print ( '    I          vdc_base(I,2)   vdc_base(I,3)   vdc_base(I,5)' )
  print ( '' )

  for i in range ( -10, 11 ):
    r2 = vdc_base ( i, 2 )
    r3 = vdc_base ( i, 3 )
    r5 = vdc_base ( i, 5 )
    print ( '  %3d         %14.8f  %14.8f  %14.8f' % ( i, r2, r3, r5 ) )

  return

def vdc_inverse ( r ):

#*****************************************************************************80
#
## vdc_inverse() inverts an element of the van der Corput sequence.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    07 June 2022
#
#  Author:
#
#    John Burkardt
#
#  Input:
#
#    real R, the I-th element of the van der Corput sequence.
#    |R| < 1.0
#
#  Output:
#
#    integer I, the index of the element of the sequence.
#
  if ( r < 0.0 ):
    s = -1
  else:
    s = +1

  t = abs ( r )

  if ( 1.0 <= t ):
    print ( '' )
    print ( 'vdc_inverse(): Fatal error!' )
    print ( '  |R| < 1.0 is required.' )
    raise Exception ( 'vdc_inverse(): Fatal error!' )

  i = 0
  p = 1

  while ( t != 0.0 ):
    t = t * 2.0
    d = t - ( t % 1.0 )
    i = i + d * p
    p = p * 2
    t = ( t % 1.0 )
#
#  Recover the sign.
#
  i = i * s

  return i

def vdc_inverse_test ( ):

#*****************************************************************************80
#
## vdc_inverse_test() tests vdc_inverse().
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    07 June 2022
#
#  Author:
#
#    John Burkardt
#
  print ( '' )
  print ( 'vdc_inverse_test():' )
  print ( '  vdc_inverse() inverts an element of a van der Corput sequence.' )
  print ( '' )
  print ( '    I        R=vdc(I)  vdc_inverse(R)' )
  print ( '' )

  for i in range ( -10, 11 ):
    r = vdc ( i )
    i2 = vdc_inverse ( r )
    print ( '  %3d  %14.8f  %3d' % ( i, r, i2 ) )

  return

def vdc_sequence ( i1, i2 ):

#*****************************************************************************80
#
## vdc_sequence() computes a sequence of elements of the van der Corput sequence.
#
#  Discussion:
#
#    This function could be rewritten to take advantage of MATLAB's
#    vectorization capabilities.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    07 June 2022
#
#  Author:
#
#    John Burkardt
#
#  Input:
#
#    integer I1, I2, the indices of the first and last elements.
#
#  Output:
#
#    real R(|I2+1-I1|), elements I1 through I2 of the van der
#    Corput sequence.
#
  import numpy as np

  n = abs ( i2 - i1 ) + 1
  r = np.zeros ( n )

  if ( i1 <= i2 ):
    i3 = +1
  else:
    i3 = -1

  j = 0
#
#  The syntax of the Python RANGE function strikes again!
#
  for i in range ( i1, i2 + i3, i3 ):
#
#  Isolate the sign, and only work with the integer part of I.
#
    if ( i < 0 ):
      s = -1
    else:
      s = +1

    t = abs ( int ( i ) )
#
#  Carry out the computation.
#
    base_inv = 0.5

    r[j] = 0.0

    while ( t != 0 ):
      d = ( t % 2 )
      r[j] = r[j] + d * base_inv
      base_inv = base_inv / 2.0
      t = ( t // 2 )
#
#  Recover the sign.
#
    r[j] = r[j] * s

    j = j + 1   

  return r

def vdc_sequence_test ( ):

#*****************************************************************************80
#
## vdc_sequence_test() tests vdc_sequence().
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    07 June 2022
#
#  Author:
#
#    John Burkardt
#
  print ( '' )
  print ( 'vdc_sequence_test():' )
  print ( '  vdc_sequence() returns elements I1 through I2 of ' )
  print ( '  a van der Corput sequence.' )
  print ( '' )

  i1 = 7
  i2 = 7
  n = abs ( i2 - i1 ) + 1
  r = vdc_sequence ( i1, i2 )
  print ( '  R=vdc_sequence(  7,  7):' )
  print ( r )

  i1 = 0
  i2 = 8
  n = abs ( i2 - i1 ) + 1
  r = vdc_sequence ( i1, i2 )
  print ( '  R=vdc_sequence(  0,  8):' )
  print ( r )

  i1 = 8
  i2 = 0
  n = abs ( i2 - i1 ) + 1
  r = vdc_sequence ( i1, i2 )
  print ( '  R=vdc_sequence(  8,  0):' )
  print ( r )

  i1 = -3
  i2 = +5
  n = abs ( i2 - i1 ) + 1
  r = vdc_sequence ( i1, i2 )
  print ( '  R=vdc_sequence( -3,  5):' )
  print ( r )

  i1 = 100
  i2 = 105
  n = abs ( i2 - i1 ) + 1
  r = vdc_sequence ( i1, i2 )
  print ( '  R=vdc_sequence(100,105):' )
  print ( r )

  return

def van_der_corput_test ( ):

#*****************************************************************************80
#
## van_der_corput_test() test van_der_corput().
#
#  Licensing:
#
#    This code is distributed under the MIT license. 
#
#  Modified:
#
#    07 June 2022
#
#  Author:
#
#    John Burkardt
#
  import platform

  print ( '' )
  print ( 'van_der_corput_test():' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  Test van_der_corput().' )

  vdc_test ( )
  vdc_inverse_test ( )
  vdc_sequence_test ( )
  vdc_base_test ( )
#
#  Terminate.
#
  print ( '' )
  print ( 'van_der_corput_test():' )
  print ( '  Normal end of execution.' )
  return

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