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

#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);
  Vprintf("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)) )
    {
      Vprintf("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;
  g05ybc(state, seq, dist, &MEAN(1), &STD(1), skip, idim,
         &QUASI(1), &GF, &fail);
  if (fail.code != NE_NOERROR)
    {
      Vprintf("Initialization Error from 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)
    {
      g05ybc(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)
    {
      Vprintf("Error from g05ybc.\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }
  val1 = sum / (double) ntimes;
  Vprintf("Calculated value of the integral = %8.3f\n\n",val1);
  val2 = (double) idim * 3.0;
  Vprintf("Exact value of the integral      = %8.3f\n",val2);

  /* Finish quasi-random generator */
  state = Nag_QuasiRandom_Finish;
  g05ybc(state, seq, dist, &MEAN(1), &STD(1), skip, idim,
         &QUASI(1), &GF, &fail);
  if (fail.code != NE_NOERROR)
    {
      Vprintf("Finish Error from 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;
}