# include <cstdlib>
# include <fstream>
# include <iomanip>
# include <iostream>
# include <iomanip>
# include <ctime>
# include <cmath>

using namespace std;

# include "van_der_corput.H"

//
//  These variables, which define the van der corput base and seed,
//  are accessible to the user via some routines.
//
static int van_der_corput_BASE = 2;
static int van_der_corput_SEED = 1;

//****************************************************************************80

char digit_to_ch ( int i )

//****************************************************************************80
//
//  Purpose:
//
//    DIGIT_TO_CH returns the base 10 digit character corresponding to a digit.
//
//  Example:
//
//     I     C
//   -----  ---
//     0    '0'
//     1    '1'
//   ...    ...
//     9    '9'  
//    10    '*'
//   -83    '*'
//
//  Modified:
//
//    16 June 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int I, the digit, which should be between 0 and 9.
//
//    Output, char DIGIT_TO_CH, the appropriate character '0' through '9' or '*'.
//
{
  char c;

  if ( 0 <= i && i <= 9 )
  {
    c = '0' + i;
  }
  else
  {
    c = '*';
  }

  return c;
}
//****************************************************************************80

int get_seed ( void )

//****************************************************************************80
//
//  Purpose:
//
//    GET_SEED returns a random seed for the random number generator.
//
//  Modified:
//
//    17 November 2004
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, int GET_SEED, a random seed value.
//
{
# define I_MAX 2147483647

  time_t clock;
  int i;
  int ihour;
  int imin;
  int isec;
  int seed;
  struct tm *lt;
  time_t tloc;
//
//  If the internal seed is 0, generate a value based on the time.
//
  clock = time ( &tloc );
  lt = localtime ( &clock );
//
//  Hours is 1, 2, ..., 12.
//
  ihour = lt->tm_hour;

  if ( 12 < ihour )
  {
    ihour = ihour - 12;
  }
//
//  Move Hours to 0, 1, ..., 11
//
  ihour = ihour - 1;

  imin = lt->tm_min;

  isec = lt->tm_sec;

  seed = isec + 60 * ( imin + 60 * ihour );
//
//  We want values in [1,43200], not [0,43199].
//
  seed = seed + 1;
//
//  Remap ISEED from [1,43200] to [1,IMAX].
//
  seed = ( int ) 
    ( ( ( double ) seed )
    * ( ( double ) I_MAX ) / ( 60.0 * 60.0 * 12.0 ) );
//
//  Never use a seed of 0.
//
  if ( seed == 0 )
  {
    seed = 1;
  }

  return seed;
# undef I_MAX
}
//****************************************************************************80

int i4_log_10 ( int i )

//****************************************************************************80
//
//  Purpose:
//
//    I4_LOG_10 returns the whole part of the logarithm base 10 of an I4.
//
//  Discussion:
//
//    It should be the case that 10^I4_LOG_10(I) <= |I| < 10^(I4_LOG_10(I)+1).
//    (except for I = 0).
//
//    The number of decimal digits in I is I4_LOG_10(I) + 1.
//
//  Examples:
//
//        I    I4_LOG_10(I)
//
//        0     0
//        1     0
//        2     0
//
//        9     0
//       10     1
//       11     1
//
//       99     1
//      100     2
//      101     2
//
//      999     2
//     1000     3
//     1001     3
//
//     9999     3
//    10000     4
//    10001     4
//
//  Modified:
//
//    17 June 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int I, the integer.
//
//    Output, int I4_LOG_10, the whole part of the logarithm of abs ( I ).
//
{
  int ten_pow;
  int value;

  i = abs ( i );

  ten_pow = 10;
  value = 0;

  while ( ten_pow <= i )
  {
    ten_pow = ten_pow * 10;
    value = value + 1;
  }

  return value;
}
//****************************************************************************80

char *i4_to_s ( int i )

//****************************************************************************80
//
//  Purpose:
//
//    I4_TO_S converts an I4 to a string.
//
//  Examples:
//
//    INTVAL  S
//
//         1  1
//        -1  -1
//         0  0
//      1952  1952
//    123456  123456
//   1234567  1234567
//
//  Modified:
//
//    13 March 2004
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int I, an integer to be converted.
//
//    Output, char *I4_TO_S, the representation of the integer.
//
{
  int digit;
  int j;
  int length;
  int ten_power;
  char *s;

  length = i4_log_10 ( i );

  ten_power = ( int ) pow ( ( double ) 10, ( double ) length );

  if ( i < 0 )
  {
    length = length + 1;
  }
//
//  Add one position for the trailing null.
//
  length = length + 1;

  s = new char[length];

  if ( i == 0 )
  {
    s[0] = '0';
    s[1] = '\0';
    return s;
  }
//
//  Now take care of the sign.
//
  j = 0;
  if ( i < 0 )
  {
    s[j] = '-';
    j = j + 1;
    i = abs ( i );
  }
//
//  Find the leading digit of I, strip it off, and stick it into the string.
//
  while ( 0 < ten_power )
  {
    digit = i / ten_power;
    s[j] = digit_to_ch ( digit );
    j = j + 1;
    i = i - digit * ten_power;
    ten_power = ten_power / 10;
  }
//
//  Tack on the trailing NULL.
//
  s[j] = '\0';
  j = j + 1;

  return s;
}
//****************************************************************************80

double i4_to_van_der_corput ( int seed, int base )

//****************************************************************************80
//
//  Purpose:
//
//    I4_TO_VAN_DER_CORPUT computes an element of a 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 (theoretically) never repeats.  Except for SEED = 0, the
//    elements of the van der Corput sequence are strictly between 0 and 1.
//
//    The van der Corput sequence writes an integer in a given base B,
//    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.
//
//    Hammersley suggested generating a set of N nicely distributed
//    points in two dimensions by setting the first component of the
//    Ith point to I/N, and the second to the van der Corput 
//    value of I in base 2.  
//
//    Halton suggested that in many cases, you might not know the number 
//    of points you were generating, so Hammersley's formulation was
//    not ideal.  Instead, he suggested that to generated a nicely
//    distributed sequence of points in M dimensions, you simply
//    choose the first M primes, P(1:M), and then for the J-th component of
//    the I-th point in the sequence, you compute the van der Corput
//    value of I in base P(J).
//
//    Thus, to generate a Halton sequence in a 2 dimensional space,
//    it is typical practice to generate a pair of van der Corput sequences,
//    the first with prime base 2, the second with prime base 3.
//    Similarly, by using the first K primes, a suitable sequence
//    in K-dimensional space can be generated.
//
//    The generation is quite simple.  Given an integer SEED, the expansion
//    of SEED in base BASE is generated.  Then, essentially, the result R
//    is generated by writing a decimal point followed by the digits of
//    the expansion of SEED, in reverse order.  This decimal value is actually
//    still in base BASE, so it must be properly interpreted to generate
//    a usable value.
//
//  Example:
//
//    BASE = 2
//
//    SEED     SEED      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
//
//  Modified:
//
//    25 February 2003
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    J H Halton,
//    On the efficiency of certain quasi-random sequences of points
//    in evaluating multi-dimensional integrals,
//    Numerische Mathematik,
//    Volume 2, pages 84-90, 1960.
// 
//    J M Hammersley,
//    Monte Carlo methods for solving multivariable problems,
//    Proceedings of the New York Academy of Science,
//    Volume 86, pages 844-874, 1960.
//
//    J G van der Corput,
//    Verteilungsfunktionen I & II,
//    Nederl. Akad. Wetensch. Proc.,
//    Volume 38, 1935, pages 813-820, pages 1058-1066.
//
//  Parameters:
//
//    Input, int SEED, the index of the desired element.
//    SEED should be nonnegative.
//    SEED = 0 is allowed, and returns R = 0.
//
//    Input, int BASE, the van der Corput base, which is usually 
//    a prime number.  BASE must be greater than 1.
//
//    Output, double VAN_DER_CORPUT, the SEED-th element of the van 
//    der Corput sequence for base BASE.
//
{
  double base_inv;
  int digit;
  double r;

  if ( base <= 1 )
  {
    cout << "\n";
    cout << "I4_TO_VAN_DER_CORPUT - Fatal error!\n";
    cout << "  The input base BASE is <= 1!\n";
    cout << "  BASE = " << base << "\n";
    exit ( 1 );
  }

  if ( seed < 0 ) 
  {
    cout << "\n";
    cout << "I4_TO_VAN_DER_CORPUT - Fatal error!\n";
    cout << "  SEED < 0.";
    cout << "  SEED = " << seed << "\n";
    exit ( 1 );
  }

  r = 0.0;

  base_inv = 1.0 / ( ( double ) base );

  while ( seed != 0 )
  {
    digit = seed % base;
    r = r + ( ( double ) digit ) * base_inv;
    base_inv = base_inv / ( ( double ) base );
    seed = seed / base;
  }

  return r;
}
//****************************************************************************80

void i4_to_van_der_corput_sequence ( int seed, int base, int n, 
  double r[] )

//****************************************************************************80
//
//  Purpose:
//
//    I4_TO_VAN_DER_CORPUT_SEQUENCE: next N elements of a van der Corput sequence.
//
//  Modified:
//
//    27 February 2003
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    J H Halton,
//    On the efficiency of certain quasi-random sequences of points
//    in evaluating multi-dimensional integrals,
//    Numerische Mathematik,
//    Volume 2, pages 84-90, 1960.
//
//    J G van der Corput,
//    Verteilungsfunktionen I & II,
//    Nederl. Akad. Wetensch. Proc.,
//    Volume 38, 1935, pages 813-820, pages 1058-1066.
//
//  Parameters:
//
//    Input, int SEED, the index of the first desired element.
//    SEED should be nonnegative.
//    SEED = 0 is allowed, and returns R = 0.
//
//    Input, int BASE, the van der Corput base, which is typically a
//    prime number.
//
//    Input, int N, the number of elements desired.
//
//    Output, double R[N], the SEED-th through 
//    (SEED+N-1)-th elements of the van der Corput sequence for 
//    the given base.
//
{
  double base_inv;
  int digit;
  int i;
  int seed2;

  if ( base <= 1 ) 
  {
    cout << "\n";
    cout << "I4_TO_VAN_DER_CORPUT_SEQUENCE - Fatal error!\n";
    cout << "  The input base BASE is <= 1!\n";
    cout << "  BASE = " << base << "\n";
    exit ( 1 );
  }

  if ( seed < 0 ) 
  {
    cout << "\n";
    cout << "I4_TO_VAN_DER_CORPUT_SEQUENCE - Fatal error!\n";
    cout << "  The input base SEED is < 0!\n";
    cout << "  SEED = " << seed << "\n";
    exit ( 1 );
  }

  for ( i = 0; i < n; i++ )
  {
    r[i] = 0.0;
    seed2 = seed + i;
    base_inv = 1.0E+00 / ( ( double ) base );

    while ( seed2 != 0 )
    {
      digit = seed2 % base;
      r[i] = r[i] + ( ( double ) digit ) * base_inv;
      base_inv = base_inv / ( ( double ) base );
      seed2 = seed2 / base;
    }
  }

  return;
}
//****************************************************************************80

double r8_epsilon ( void )

//****************************************************************************80
//
//  Purpose:
//
//    R8_EPSILON returns the R8 round off unit.
//
//  Discussion:
//
//    R8_EPSILON is a number R which is a power of 2 with the property that,
//    to the precision of the computer's arithmetic,
//      1 < 1 + R
//    but 
//      1 = ( 1 + R / 2 )
//
//  Modified:
//
//    01 July 2004
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double R8_EPSILON, the double precision round-off unit.
//
{
  double r;

  r = 1.0;

  while ( 1.0 < ( double ) ( 1.0 + r )  )
  {
    r = r / 2.0;
  }

  return r;
}
//****************************************************************************80

void timestamp ( void )

//****************************************************************************80
//
//  Purpose:
//
//    TIMESTAMP prints the current YMDHMS date as a time stamp.
//
//  Example:
//
//    May 31 2001 09:45:54 AM
//
//  Modified:
//
//    04 October 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    None
//
{
# define TIME_SIZE 40

  static char time_buffer[TIME_SIZE];
  const struct tm *tm;
  size_t len;
  time_t now;

  now = time ( NULL );
  tm = localtime ( &now );

  len = strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm );

  cout << time_buffer << "\n";

  return;
# undef TIME_SIZE
}
//****************************************************************************80

char *timestring ( void )

//****************************************************************************80
//
//  Purpose:
//
//    TIMESTRING returns the current YMDHMS date as a string.
//
//  Example:
//
//    May 31 2001 09:45:54 AM
//
//  Modified:
//
//    24 September 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, char *TIMESTRING, a string containing the current YMDHMS date.
//
{
# define TIME_SIZE 40

  const struct tm *tm;
  size_t len;
  time_t now;
  char *s;

  now = time ( NULL );
  tm = localtime ( &now );

  s = new char[TIME_SIZE];

  len = strftime ( s, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm );

  return s;
# undef TIME_SIZE
}
//****************************************************************************80

double van_der_corput ( void )

//****************************************************************************80
//
//  Purpose:
//
//    VAN_DER_CORPUT computes the next element in the van der Corput sequence.
//
//  Discussion:
//
//    The internal variables van_der_corput_SEED and van_der_corput_BASE 
//    control the van der Corput sequence.
//
//    The value of van_der_corput_SEED is incremented by one on each call.
//
//  Modified:
//
//    25 February 2003
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    J G van der Corput,
//    Verteilungsfunktionen I & II,
//    Nederl. Akad. Wetensch. Proc.,
//    Volume 38, 1935, pages 813-820, pages 1058-1066.
//
//  Parameters:
//
//    Output, double VAN_DER_CORPUT, the next element of the van der Corput sequence.
//
{
  double r;

  r = i4_to_van_der_corput ( van_der_corput_SEED, van_der_corput_BASE );

  van_der_corput_SEED = van_der_corput_SEED + 1;

  return r;
}
//****************************************************************************80

int van_der_corput_base_get ( void )

//****************************************************************************80
//
//  Purpose:
//
//    VAN_DER_CORPUT_BASE_GET gets the base for a van der Corput sequence.
//
//  Modified:
//
//    28 August 2002
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    J G van der Corput,
//    Verteilungsfunktionen I & II,
//    Nederl. Akad. Wetensch. Proc.,
//    Volume 38, 1935, pages 813-820, pages 1058-1066.
//
//  Parameters:
//
//    Output, int VAN_DER_CORPUT_BASE_GET, the base for the 
//    van der Corput sequence.
//
{
  return van_der_corput_BASE;
}
//****************************************************************************80

void van_der_corput_base_set ( int base )

//****************************************************************************80
//
//  Purpose:
//
//    VAN_DER_CORPUT_BASE_SET sets the base for a van der Corput sequence.
//
//
//  Modified:
//
//    28 August 2002
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    J G van der Corput,
//    Verteilungsfunktionen I & II,
//    Nederl. Akad. Wetensch. Proc.,
//    Volume 38, 1935, pages 813-820, pages 1058-1066.
//
//  Parameters:
//
//    Input, int BASE, the base for the van der Corput sequence.
//    BASE must be greater than 1.
//
{
  if ( base <= 1 )
  {
    cout << "\n";
    cout << "VAN_DER_CORPUT_BASE_SET - Fatal error!\n";
    cout << "  The input base BASE is <= 1!\n";
    cout << "  BASE = " << base << "\n";
    exit ( 1 );
  }

  van_der_corput_BASE = base;

  return;
}
//****************************************************************************80

int van_der_corput_seed_get ( void )

//****************************************************************************80
//
//  Purpose:
//
//    VAN_DER_CORPUT_SEED_GET gets the "seed" for the van der Corput sequence.
//
//  Modified:
//
//    25 February 2003
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    J G van der Corput,
//    Verteilungsfunktionen I & II,
//    Nederl. Akad. Wetensch. Proc.,
//    Volume 38, 1935, pages 813-820, pages 1058-1066.
//
//  Parameters:
//
//    Output, int VAN_DER_CORPUT_SEED_GET, the current seed for 
//    the van der Corput sequence.
//
{
  return van_der_corput_SEED;
}
//****************************************************************************80

void van_der_corput_seed_set ( int seed )

//****************************************************************************80
//
//  Purpose:
//
//    VAN_DER_CORPUT_SEED_SET sets the "seed" for the van der Corput sequence.
//
//  Discussion:
//
//    Calling VAN_DER_CORPUT repeatedly returns the elements of the 
//    van der Corput sequence in order, starting with element number 1.  
//    An internal counter, called SEED, keeps track of the next element 
//    to return.  Each time the routine is called, the SEED-th element
//    is computed, and then SEED is incremented by 1.
//
//    To restart the van der Corput sequence, it is only necessary to reset 
//    SEED to 1.  It might also be desirable to reset SEED to some other value.
//    This routine allows the user to specify any value of SEED.
//
//    The default value of SEED is 1, which restarts the van der Corput sequence.
//
//  Modified:
//
//    25 February 2003
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    J G van der Corput,
//    Verteilungsfunktionen I & II,
//    Nederl. Akad. Wetensch. Proc.,
//    Volume 38, 1935, pages 813-820, pages 1058-1066.
//
//  Parameters:
//
//    Input, int SEED, the seed for the van der Corput sequence.
//    SEED must be nonnegative.
//
{
  if ( seed < 0 ) 
  {
    cout << "\n";
    cout << "VAN_DER_CORPUT_SEED_SET - Fatal error!\n";
    cout << "  SEED < 0.";
    cout << "  SEED = " << seed << "\n";
    exit ( 1 );
  }

  van_der_corput_SEED = seed;

  return;
}
//****************************************************************************80

void van_der_corput_sequence ( int n, double r[] )

//****************************************************************************80
//
//  Purpose:
//
//    VAN_DER_CORPUT_SEQUENCE: next N elements in the van der Corput sequence.
//
//  Modified:
//
//    27 February 2003
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    J G van der Corput,
//    Verteilungsfunktionen I & II,
//    Nederl. Akad. Wetensch. Proc.,
//    Volume 38, 1935, pages 813-820, pages 1058-1066.
//
//  Parameters:
//
//    Input, int N, the number of elements desired.
//
//    Output, double R[N], the next N elements of the van 
//    der Corput sequence.
//
{
  int base;
  int seed;

  seed = van_der_corput_seed_get ( );

  base = van_der_corput_base_get ( );

  i4_to_van_der_corput_sequence ( seed, base, n, r );

  seed = seed + n;
  van_der_corput_seed_set ( seed );

  return;
}
//****************************************************************************80

void van_der_corput_write ( int n, int seed, int base, double r[], 
  char *file_out_name )

//****************************************************************************80
//
//  Purpose:
//
//    VAN_DER_CORPUT_WRITE writes a van der Corput dataset to a file.
//
//  Discussion:
//
//    The initial lines of the file are comments, which begin with a
//    "#" character.
//
//    Thereafter, each line of the file contains the NDIM-dimensional
//    components of the next entry of the sequence.
//
//  Modified:
//
//    30 April 2005
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the number of elements in the subsequence.
//
//    Input, int SEED, the initial seed.
//
//    Input, int BASE, the base.
//
//    Input, double R[M*N], the points.
//
//    Input, char *FILE_OUT_NAME, the name of the output file.
//
{
  ofstream file_out;
  int i;
  char *s;

  file_out.open ( file_out_name );

  if ( !file_out )
  {
    cout << "\n";
    cout << "VAN_DER_CORPUT_WRITE - Fatal error!\n";
    cout << "  Could not open the output file.\n";
    exit ( 1 );
  }

  s = timestring ( );

  file_out << "#  " << file_out_name << "\n";
  file_out << "#  created by routine VAN_DER_CORPUT.C" << "\n";
  file_out << "#  at " << s << "\n";
  file_out << "#\n";
  file_out << "#  N =    " << setw(12) << n    << "\n";
  file_out << "#  SEED = " << setw(12) << seed << "\n";
  file_out << "#  BASE = " << setw(12) << base    << "\n";
  file_out << "#  EPSILON (unit roundoff) = " << r8_epsilon ( ) << "\n";
  file_out << "#\n";

  for ( i = 0; i < n; i++ )
  {
    file_out << setw(10) << r[i] << "\n";
  }

  file_out.close ( );

  delete [] s;

  return;
}