/* nag_quad_dim1_gauss_1 (d01tac) 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 fun1(double x, Nag_User *comm);
  static double NAG_CALL fun2(double x, Nag_User *comm);
  static double NAG_CALL fun3(double x, Nag_User *comm);
  static double NAG_CALL fun4(double x, Nag_User *comm);
#ifdef __cplusplus
}
#endif

int main(void)
{
  static Integer use_comm[4] = { 1, 1, 1, 1 };
  Integer exit_status = 0;
  static Integer nstor[3] = { 4, 8, 16 };
  double a, b;
  Integer i;
  double ans;
  Nag_GaussFormulae gaussformula;
  Nag_User comm;
  NagError fail;

  fail.print = Nag_TRUE;

  INIT_FAIL(fail);

  printf("nag_quad_dim1_gauss_1 (d01tac) Example Program Results\n");

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

  printf("\nGauss-Legendre example\n\n");
  for (i = 0; i < 3; ++i) {
    a = 0.0;
    b = 1.0;
    gaussformula = Nag_Legendre;
    /* nag_quad_dim1_gauss_1 (d01tac).
     * One-dimensional Gaussian quadrature rule evaluation,
     * thread-safe
     */
    ans = d01tac(gaussformula, fun1, a, b, nstor[i],
                                       &comm, &fail);
    if (fail.code == NE_NOERROR || fail.code == NE_QUAD_GAUSS_NPTS_RULE)
      printf("%" NAG_IFMT " Points    Answer = %10.5f\n\n", nstor[i], ans);
    else {
      printf("%s\n", fail.message);
      exit_status = 1;
      goto END;
    }
  }
  printf("\nGauss-Rational example\n\n");
  for (i = 0; i < 3; ++i) {
    a = 2.0;
    b = 0.0;
    gaussformula = Nag_Rational;
    /* nag_quad_dim1_gauss_1 (d01tac), see above. */
    ans = d01tac(gaussformula, fun2, a, b, nstor[i],
                                       &comm, &fail);
    if (fail.code == NE_NOERROR || fail.code == NE_QUAD_GAUSS_NPTS_RULE)
      printf("%" NAG_IFMT " Points    Answer = %10.5f\n\n", nstor[i], ans);
    else {
      printf("%s\n", fail.message);
      exit_status = 1;
      goto END;
    }
  }
  printf("\nGauss-Laguerre example\n\n");
  for (i = 0; i < 3; ++i) {
    a = 2.0;
    b = 1.0;
    gaussformula = Nag_Laguerre;
    /* nag_quad_dim1_gauss_1 (d01tac), see above. */
    ans = d01tac(gaussformula, fun3, a, b, nstor[i],
                                       &comm, &fail);
    if (fail.code == NE_NOERROR || fail.code == NE_QUAD_GAUSS_NPTS_RULE)
      printf("%" NAG_IFMT " Points    Answer = %10.5f\n\n", nstor[i], ans);
    else {
      printf("%s\n", fail.message);
      exit_status = 1;
      goto END;
    }
  }
  printf("\nGauss-Hermite  example\n\n");
  for (i = 0; i < 3; ++i) {
    a = -1.0;
    b = 3.0;
    gaussformula = Nag_Hermite;
    /* nag_quad_dim1_gauss_1 (d01tac), see above. */
    ans = d01tac(gaussformula, fun4, a, b, nstor[i],
                                       &comm, &fail);
    if (fail.code == NE_NOERROR || fail.code == NE_QUAD_GAUSS_NPTS_RULE)
      printf("%" NAG_IFMT " Points    Answer = %10.5f\n\n", nstor[i], ans);
    else {
      printf("%s\n", fail.message);
      exit_status = 1;
      goto END;
    }
  }

END:
  return exit_status;
}

static double NAG_CALL fun1(double x, Nag_User *comm)
{
  Integer *use_comm = (Integer *) comm->p;

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

  return 4.0 / (x * x + 1.0);
}

static double NAG_CALL fun2(double x, Nag_User *comm)
{
  Integer *use_comm = (Integer *) comm->p;

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

  return 1.0 / (x * x * log(x));
}

static double NAG_CALL fun3(double x, Nag_User *comm)
{
  Integer *use_comm = (Integer *) comm->p;

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

  return exp(-x) / x;
}

static double NAG_CALL fun4(double x, Nag_User *comm)
{
  Integer *use_comm = (Integer *) comm->p;

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

  return exp(x * (-3.0) * x - x * 4.0 - 1.0);
}