/* nag_quad_dim1_gauss_vec (d01uac) Example Program.
 *
 * Copyright 2022 Numerical Algorithms Group.
 *
 * Mark 28.7, 2022.
 */
#include <math.h>
#include <nag.h>

#ifdef __cplusplus
extern "C" {
#endif
static void NAG_CALL f(const double x[], const Integer nx, double fv[],
                       Integer *iflag, Nag_Comm *comm);
#ifdef __cplusplus
}
#endif

int main(void) {
  /* Scalars */
  Integer exit_status = 0;
  double a, b, dinest;
  Integer funid, i, nstor;
  /* Arrays */
  Integer iuser[1];
  /* Nag Types */
  Nag_Comm comm;
  NagError fail;
  Nag_QuadType quad_type;

  printf("nag_quad_dim1_gauss_vec (d01uac) Example Program Results\n");

  INIT_FAIL(fail);

  /* Use comm to pass information to f. */
  comm.iuser = iuser;

  for (funid = 1; funid <= 6; funid++) {
    switch (funid) {
    case 1: {
      printf("\nGauss-Legendre example\n");
      a = 0.0;
      b = 1.0;
      quad_type = Nag_Quad_Gauss_Legendre;
      break;
    }
    case 2: {
      printf("\nRational Gauss example\n");
      a = 2.0;
      b = 0.0;
      quad_type = Nag_Quad_Gauss_Rational_Adjusted;
      break;
    }
    case 3: {
      printf("\nGauss-Laguerre example (adjusted weights)\n");
      a = 2.0;
      b = 1.0;
      quad_type = Nag_Quad_Gauss_Laguerre_Adjusted;
      break;
    }
    case 4: {
      printf("\nGauss-Laguerre example (normal weights)\n");
      a = 2.0;
      b = 1.0;
      quad_type = Nag_Quad_Gauss_Laguerre;
      break;
    }
    case 5: {
      printf("\nGauss-Hermite example (adjusted weights)\n");
      a = -1.0;
      b = 3.0;
      quad_type = Nag_Quad_Gauss_Hermite_Adjusted;
      break;
    }
    case 6: {
      printf("\nGauss-Hermite example (normal weights)\n");
      a = -1.0;
      b = 3.0;
      quad_type = Nag_Quad_Gauss_Hermite;
      break;
    }
    }

    iuser[0] = funid;
    for (i = 0; i < 6; i++) {
      nstor = pow(2, i + 1);
      /* Compute the one-dimensional integral employing Gaussian quadrature,
       * with quadrature type and weights specified in quad_type, using
       * nag_quad_dim1_gauss_vec (d01uac).
       */
      nag_quad_dim1_gauss_vec(quad_type, a, b, nstor, f, &dinest, &comm, &fail);
      switch (fail.code) {
      case NE_NOERROR:
      case NE_QUAD_GAUSS_NPTS_RULE:
      case NE_UNDERFLOW:
      case NE_WEIGHT_ZERO: {
        /* The definite integral has been estimated. */
        printf("%5" NAG_IFMT "   Points  Answer = %10.5f\n", nstor, dinest);
        break;
      }
      default: {
        /* A solution could not be calculated due to an illegal parameter
         * or a requested exit.
         */
        printf("%s\n", fail.message);
        exit_status++;
      }
      }
    }
  }
  return exit_status;
}

static void NAG_CALL f(const double x[], const Integer nx, double fv[],
                       Integer *iflag, Nag_Comm *comm) {
  Integer i, funid;

  funid = comm->iuser[0];
  switch (funid) {
  case 1:
    for (i = 0; i < nx; i++)
      fv[i] = 4.0 / (1.0 + x[i] * x[i]);
    break;
  case 2:
    for (i = 0; i < nx; i++)
      fv[i] = 1.0 / (x[i] * x[i] * log(x[i]));
    break;
  case 3:
    for (i = 0; i < nx; i++)
      fv[i] = exp(-x[i]) / x[i];
    break;
  case 4:
    for (i = 0; i < nx; i++)
      fv[i] = 1.0 / x[i];
    break;
  case 5:
    for (i = 0; i < nx; i++)
      fv[i] = exp(-3.0 * x[i] * x[i] - 4.0 * x[i] - 1.0);
    break;
  case 6:
    for (i = 0; i < nx; i++)
      fv[i] = exp(2.0 * x[i] + 2.0);
    break;
  default:
    *iflag = -1;
  }
}