/* nag_quasi_random_uniform (g05yac) Example Program.
 *
 * Copyright 2001 Numerical Algorithms Group.
 *
 * Mark 7, 2001.
 * Mark 7b revised, 2004.
 */

#include <nag.h>
#include <nagx04.h>
#include <nag_types.h>
#include <nag_stdlib.h>
#include <nagg05.h>
#include <stdio.h>

static double fun(Integer idim, double  x[]);

int main(int argc, char *argv[])
{
  FILE                     *fpout;
  /* Scalars */
  Integer                  i, idim, ntimes, skip;
  Integer                  exit_status = 0;
  double                   sum, vsbl;
  NagError                 fail;
  Nag_QuasiRandom          GF;
  Nag_QuasiRandom_Sequence seq;
  Nag_QuasiRandom_State    state;

  /* Arrays */
  double                   *quasi;

#define QUASI(I) quasi[(I) -1]

  INIT_FAIL(fail);

  /* Check for command-line IO options */
  fpout = nag_example_file_io(argc, argv, "-results", NULL);
  fprintf(fpout,
          "nag_quasi_random_uniform (g05yac) Example Program Results\n\n");

  idim = 22;

  /* Allocate memory */
  if (!(quasi = NAG_ALLOC(idim, double)))
    {
      fprintf(fpout, "Allocation error \n");
      exit_status = -1;
      goto END;
    }

  ntimes = 50000;
  seq = Nag_QuasiRandom_Faure;
  if (seq == Nag_QuasiRandom_Nied)
    skip = 1000;
  else
    skip = 0;
  /* Initialise quasi-random generator*/
  state = Nag_QuasiRandom_Init;
  /* nag_quasi_random_uniform (g05yac).
   * Multi-dimensional quasi-random number generator with a
   * uniform probability distribution
   */
  nag_quasi_random_uniform(state, seq, skip, idim, &QUASI(1), &GF, &fail);
  if (fail.code != NE_NOERROR)
    {
      fprintf(fpout,
              "Initialization Error from nag_quasi_random_uniform (g05yac)."
              "\n%s\n",
              fail.message);
      exit_status = 1;
      goto END;
    }

  /* Evaluate integrand at quasi-random locations and sum */
  sum = 0.0;
  state = Nag_QuasiRandom_Cont;
  for (i = 1; i <= ntimes; ++i)
    {
      /* nag_quasi_random_uniform (g05yac), see above. */
      nag_quasi_random_uniform(state, seq, skip, idim, &QUASI(1), &GF, &fail);
      sum += fun(idim, &QUASI(1));
    }
  if (fail.code != NE_NOERROR)
    {
      fprintf(fpout, "Error from nag_quasi_random_uniform (g05yac).\n%s\n",
              fail.message);
      exit_status = 1;
      goto END;
    }
  vsbl = sum / (double) ntimes;
  fprintf(fpout, "Value of integral = %8.3f\n", vsbl);

  /* Finish quasi-random generator */
  state = Nag_QuasiRandom_Finish;
  /* nag_quasi_random_uniform (g05yac), see above. */
  nag_quasi_random_uniform(state, seq, skip, idim, &QUASI(1), &GF, &fail);
  if (fail.code != NE_NOERROR)
    {
      fprintf(fpout,
              "Finish Error from nag_quasi_random_uniform (g05yac).\n%s\n",
              fail.message);
      exit_status = 1;
      goto END;
    }

 END:
  if (fpout != stdout) fclose(fpout);
  if (quasi) NAG_FREE(quasi);
  return exit_status;
}

static double fun(Integer idim, double  x[])
{
  Integer j;
  double  tmp;

#define X(I) x[(I) -1]

  tmp = 1.0;
  for (j = 1; j <= idim; ++j)
    {
      tmp = tmp * ABS(X(j) * 4.0 - 2.0);
    }
  return tmp;
}