/* nag_quad_multid_quad_monte_carlo_1 (d01xbc) Example Program.
 *
 * Copyright 2019 Numerical Algorithms Group.
 *
 * Mark 27.0, 2019.
 *
 */

#include <nag.h>
#include <stdio.h>
#include <math.h>

#ifdef __cplusplus
extern "C"
{
#endif
  static double NAG_CALL f(Integer ndim, const double x[], Nag_User *comm);
#ifdef __cplusplus
}
#endif

#define MAXCLS 20000

int main(void)
{

  static Integer use_comm[1] = { 1 };
  Integer exit_status = 0, k, maxcls = MAXCLS, mincls, ndim = 4;
  NagError fail;
  Nag_MCMethod method;
  Nag_Start cont;
  Nag_User comm;
  double *a = 0, acc, *b = 0, *comm_arr = 0, eps, finest;

  INIT_FAIL(fail);

  printf("nag_quad_multid_quad_monte_carlo_1 (d01xbc) Example Program Results\n");

  /* For communication with user-supplied functions: */
  comm.p = (Pointer) &use_comm;

  if (ndim >= 1) {
    if (!(a = NAG_ALLOC(ndim, double)) || !(b = NAG_ALLOC(ndim, double)))
    {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }
  }
  else {
    printf("Invalid ndim.\n");
    exit_status = 1;
    return exit_status;
  }
  for (k = 0; k < ndim; ++k) {
    a[k] = 0.0;
    b[k] = 1.0;
  }
  eps = 0.01;
  mincls = 1000;
  method = Nag_ManyIterations;
  cont = Nag_Cold;

  /* nag_quad_multid_quad_monte_carlo_1 (d01xbc).
   * Multi-dimensional quadrature, using Monte Carlo method,
   * thread-safe
   */
  nag_quad_multid_quad_monte_carlo_1(ndim, f, method, cont, a, b, &mincls, maxcls,
                                eps, &finest, &acc, &comm_arr, &comm, &fail);
  if (fail.code == NE_NOERROR || fail.code == NE_QUAD_MAX_INTEGRAND_EVAL) {
    if (fail.code == NE_QUAD_MAX_INTEGRAND_EVAL) {
      printf("Error from nag_quad_multid_quad_monte_carlo_1 (d01xbc).\n%s\n",
             fail.message);
      exit_status = 2;
    }
    printf("Requested accuracy    = %11.2e\n", eps);
    printf("Estimated value       = %10.5f\n", finest);
    printf("Estimated accuracy    = %11.2e\n", acc);
    printf("Number of evaluations = %5" NAG_IFMT "\n", mincls);
  }
  else {
    printf("Error from nag_quad_multid_quad_monte_carlo_1 (d01xbc).\n%s\n",
           fail.message);
    printf("%s\n", fail.message);
    exit_status = 1;
  }
END:
  NAG_FREE(a);
  NAG_FREE(b);
  /* Free memory allocated internally */
  NAG_FREE(comm_arr);
  return exit_status;
}

static double NAG_CALL f(Integer ndim, const double x[], Nag_User *comm)
{
  Integer *use_comm = (Integer *) comm->p;

  if (use_comm[0]) {
    printf("(User-supplied callback f, first invocation.)\n");
    use_comm[0] = 0;
  }

  return x[0] * 4.0 * (x[2] * x[2]) * exp(x[0] * 2.0 * x[2]) /
         ((x[1] + 1.0 + x[ndim - 1]) * (x[1] + 1.0 + x[ndim - 1]));
}