# ifdef ANSI_HEADERS
#   include <cstdlib>
#   include <cmath>
#   include <ctime>
# else
#   include <stdlib.h>
#   include <math.h>
#   include <time.h>
# endif

# include <iostream>
# include <fstream>
# include <iomanip>

using namespace std;

# include "fsu.hpp"

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

void fsu_hammersley  ( int dim_num, int n, int step, int seed[], int leap[],
  int base[], double r[] )

//****************************************************************************80
//
//  Purpose:
//
//    FSU_HAMMERSLEY computes N elements of a leaped Hammersley subsequence.
//
//  License:
//
//    Copyright (C) 2004  John Burkardt and Max Gunzburger
//
//    This library is free software; you can redistribute it and/or
//    modify it under the terms of the GNU Lesser General Public
//    License as published by the Free Software Foundation; either
//    version 2.1 of the License, or (at your option) any later version.
//
//    This library is distributed in the hope that it will be useful,
//    but WITHOUT ANY WARRANTY; without even the implied warranty of
//    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
//    Lesser General Public License for more details.
//
//    You should have received a copy of the GNU Lesser General Public
//    License along with this library; if not, write to the Free Software
//    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
//
//  Discussion:
//
//    The DIM_NUM-dimensional Hammersley sequence is really DIM_NUM separate
//    sequences, each generated by a particular base.  If the base is 
//    greater than 1, a standard 1-dimensional
//    van der Corput sequence is generated.  But if the base is 
//    negative, this is a signal that the much simpler sequence J/(-BASE) 
//    is to be generated.  For the standard Hammersley sequence, the
//    first spatial coordinate uses a base of (-N), and subsequent
//    coordinates use bases of successive primes (2, 3, 5, 7, 11, ...).
//    This program allows the user to specify any combination of bases,
//    included nonprimes and repeated values.
//
//    This routine selects elements of a "leaped" subsequence of the
//    Hammersley sequence.  The subsequence elements are indexed by a
//    quantity called STEP, which starts at 0.  The STEP-th subsequence
//    element is simply element
//
//      SEED(1:DIM_NUM) + STEP * LEAP(1:DIM_NUM)
//
//    of the original Hammersley sequence.
//
//
//    The data to be computed has two dimensions.
//
//    The number of data items is DIM_NUM * N, where DIM_NUM is the spatial dimension
//    of each element of the sequence, and N is the number of elements of the sequence.
//
//    The data is stored in a one dimensional array R.  The first element of the
//    sequence is stored in the first DIM_NUM entries of R, followed by the DIM_NUM entries
//    of the second element, and so on.
//
//    In particular, the J-th element of the sequence is stored in entries
//    0+(J-1)*DIM_NUM through (DIM_NUM-1) + (J-1)*DIM_NUM.
//
//  Modified:
//
//    10 November 2006
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    J M Hammersley,
//    Monte Carlo methods for solving multivariable problems,
//    Proceedings of the New York Academy of Science,
//    Volume 86, 1960, pages 844-874.
//
//    Ladislav Kocis and William Whiten,
//    Computational Investigations of Low-Discrepancy Sequences,
//    ACM Transactions on Mathematical Software,
//    Volume 23, Number 2, 1997, pages 266-294.
//
//  Parameters:
//
//    Input, int DIM_NUM, the spatial dimension.
//
//    Input, int N, the number of elements of the sequence.
//
//    Input, int STEP, the index of the subsequence element.
//    0 <= STEP is required
//
//    Input, int SEED[DIM_NUM], the Hammersley sequence index corresponding
//    to STEP = 0. 
//
//    Input, int LEAP[DIM_NUM], the succesive jumps in the Hammersley sequence.
//
//    Input, int BASE[DIM_NUM], the Hammersley bases.
//
//    Output, double R[DIM_NUM*N], the next N elements of the
//    leaped Hammersley subsequence, beginning with element STEP.
//
{
# define FIDDLE 1.0

  double base_inv;
  int digit;
  int i;
  int j;
  int *seed2;
  int temp;
//
//  Check the input.
//
  if ( !halham_dim_num_check ( dim_num ) )
  {
    exit ( 1 );
  }

  if ( !halham_n_check ( n ) )
  {
    exit ( 1 );
  }

  if ( !halham_step_check ( step ) )
  {
    exit ( 1 );
  }

  if ( !halham_seed_check ( dim_num, seed ) )
  {
    exit ( 1 );
  }

  if ( !halham_leap_check ( dim_num, leap ) )
  {
    exit ( 1 );
  }

  if ( !hammersley_base_check ( dim_num, base ) )
  {
    exit ( 1 );
  }
//
//  Calculate the data.
//
  seed2 = new int[n];

  for ( i = 0; i < dim_num; i++ )
  {
    if ( 1 < base[i] )
    {
      for ( j = 0; j < n; j++ )
      {
        seed2[j] = seed[i] + ( step + j ) * leap[i];
      }

      for ( j = 0; j < n; j++ )
      {
        r[i+j*dim_num] = 0.0E+00;
      }

      for ( j = 0; j < n; j++ )
      {
        base_inv = 1.0E+00 / ( ( double ) base[i] );
  
        while ( seed2[j] != 0 )
        {
          digit = seed2[j] % base[i];
          r[i+j*dim_num] = r[i+j*dim_num] + ( ( double ) digit ) * base_inv;
          base_inv = base_inv / ( ( double ) base[i] );
          seed2[j] = seed2[j] / base[i];
        }
      }
    }
//
//  In the following computation, the value of FIDDLE can be:
//
//    0,   for the sequence 0/N, 1/N, ..., N-1/N
//    1,   for the sequence 1/N, 2/N, ..., N/N
//    1/2, for the sequence 1/(2N), 3/(2N), ..., (2*N-1)/(2N)
//
    else
    {
      for ( j = 0; j < n; j++ )
      {
        temp = ( seed[i] + ( step + j - 1 ) * leap[i] ) % ( -base[i] );

        r[i+j*dim_num] = ( ( double ) ( temp ) + FIDDLE )
                    / ( double ) ( -base[i] );
      }
    }
  }

  delete [] seed2;

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

bool hammersley_base_check ( int dim_num, int base[] )

//****************************************************************************80
//
//  Purpose:
//
//    HAMMERSLEY_BASE_CHECK checks BASE for a Hammersley sequence.
//
//  License:
//
//    Copyright (C) 2004  John Burkardt and Max Gunzburger
//
//    This library is free software; you can redistribute it and/or
//    modify it under the terms of the GNU Lesser General Public
//    License as published by the Free Software Foundation; either
//    version 2.1 of the License, or (at your option) any later version.
//
//    This library is distributed in the hope that it will be useful,
//    but WITHOUT ANY WARRANTY; without even the implied warranty of
//    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
//    Lesser General Public License for more details.
//
//    You should have received a copy of the GNU Lesser General Public
//    License along with this library; if not, write to the Free Software
//    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
//
//  Modified:
//
//    10 November 2006
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int DIM_NUM, the spatial dimension.
//
//    Input, int BASE[DIM_NUM], the bases.
//
//    Output, bool HAMMERSLEY_BASE_CHECK, is true if BASE is legal.
{
  int i;
  bool value;

  value = true;

  for ( i = 0; i < dim_num; i++ )
  {
    if ( base[i] == 0 || base[i] == 1 ) 
    {
      cout << "\n";
      cout << "HAMMERSLEY_BASE_CHECK - Fatal error!\n";
      cout << "  Bases may not be 0 or 1.\n";
      i4vec_transpose_print ( dim_num, base, "BASE:  " );
      value = false;
      break;
    }
  }

  return value;
}