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

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

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

int main(void)
{
  /* Scalars */
  Integer                  i, idim, ntimes, skip;
  Integer                  exit_status = 0;
  double                   sum, val1, val2;
  NagError                 fail;
  static Nag_QuasiRandom   GF;
  Nag_QuasiRandom_Sequence seq;
  Nag_QuasiRandom_State    state;
  Nag_Distributions        dist;

  /* Arrays */
  double                   *mean, *quasi, *std;

#define STD(I)   std[(I) -1]
#define QUASI(I) quasi[(I) -1]
#define MEAN(I)  mean[(I) -1]

  INIT_FAIL(fail);

  printf(
          "nag_quasi_random_normal (g05ybc) Example Program Results\n\n");

  idim = 20;
  /* Allocate memory */
  if (!(mean = NAG_ALLOC(idim, double)) ||
      !(quasi = NAG_ALLOC(idim, double)) ||
      !(std = NAG_ALLOC(idim, double)))
    {
      printf("Allocation error \n");
      exit_status = -1;
      goto END;
    }
  ntimes = 10000;
  dist = Nag_Normal;
  seq = Nag_QuasiRandom_Faure;
  if (seq == Nag_QuasiRandom_Nied)
    skip = 1000;
  else
    skip = 0;

  for (i = 1; i <= idim; ++i)
    {
      MEAN(i) = 2.0;
      STD(i) = 1.0;
    }

  /* Initialise quasi-random generator*/
  state = Nag_QuasiRandom_Init;
  /* nag_quasi_random_normal (g05ybc).
   * Multi-dimensional quasi-random number generator with a
   * Gaussian or log-normal probability distribution
   */
  nag_quasi_random_normal(state, seq, dist, &MEAN(1), &STD(1), skip, idim,
                          &QUASI(1), &GF, &fail);
  if (fail.code != NE_NOERROR)
    {
      printf(
              "Initialization Error from nag_quasi_random_normal (g05ybc)."
              "\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_normal (g05ybc), see above. */
      nag_quasi_random_normal(state, seq, dist, &MEAN(1), &STD(1), skip,
                              idim, &QUASI(1), &GF, &fail);
      sum += fun(idim, &MEAN(1), &STD(1), &QUASI(1));
    }
  if (fail.code != NE_NOERROR)
    {
      printf("Error from nag_quasi_random_normal (g05ybc).\n%s\n",
              fail.message);
      exit_status = 1;
      goto END;
    }
  val1 = sum / (double) ntimes;
  printf("Calculated value of the integral = %8.3f\n\n", val1);
  val2 = (double) idim * 3.0;
  printf("Exact value of the integral      = %8.3f\n", val2);

  /* Finish quasi-random generator */
  state = Nag_QuasiRandom_Finish;
  /* nag_quasi_random_normal (g05ybc), see above. */
  nag_quasi_random_normal(state, seq, dist, &MEAN(1), &STD(1), skip, idim,
                          &QUASI(1), &GF, &fail);
  if (fail.code != NE_NOERROR)
    {
      printf(
              "Finish Error from nag_quasi_random_normal (g05ybc).\n%s\n",
              fail.message);
      exit_status = 1;
      goto END;
    }
 END:
  if (mean) NAG_FREE(mean);
  if (quasi) NAG_FREE(quasi);
  if (std) NAG_FREE(std);
  return exit_status;
}

static double fun(Integer idim, double  mean[], double std[], double x[])
{
  Integer j;
  double  tmp1, tmp2;
#define X(I)    x[(I) -1]
#define STD(I)  std[(I) -1]
#define MEAN(I) mean[(I) -1]

  tmp1 = 0.0;
  for (j = 1; j <= idim; ++j)
    {
      tmp2 = (X(j) - MEAN(j)) / STD(j);
      tmp1 += tmp2 * tmp2 * tmp2 * tmp2;
    }
  return tmp1;
}