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

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

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

int main(void)
{
  /* 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);
  printf("nag_quasi_random_uniform (g05yac) Example Program Results\n\n");

  idim = 22;

  /* Allocate memory */
  if ( !(quasi = NAG_ALLOC(idim, double)))
    {
      Vprintf("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)
    {
      Vprintf("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)
    {
      Vprintf("Error from nag_quasi_random_uniform (g05yac).\n%s\n",
              fail.message);
      exit_status = 1;
      goto END;
    }
  vsbl = sum / (double) ntimes;
  Vprintf("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)
    {
      Vprintf("Finish Error from nag_quasi_random_uniform (g05yac).\n%s\n",
              fail.message);
      exit_status = 1;
      goto END;
    }

 END:
  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;
}