/* nag_quad_dim1_gauss_wrec (d01tdc) Example Program.
 *
 * Copyright 2023 Numerical Algorithms Group.
 *
 * Mark 29.0, 2023.
 */

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

int main(void) {
  /* Scalars */
  Integer exit_status = 0;
  Integer n, i;
  double muzero, alpha, beta, a1, b1, ab1, ab2, d, ri, fa1, fb1, fab2;
  /* Arrays */
  char nag_enum_arg[40];
  double *a = 0, *b = 0, *c = 0, *abscis = 0, *weight = 0;
  const char *str_type;
  /* Nag Types */
  Nag_QuadType quadtype;
  NagError fail, fail1, fail2;

  INIT_FAIL(fail);

  printf("nag_quad_dim1_gauss_wrec (d01tdc) Example Program Results\n");
  /* Skip heading in data file */
  scanf("%*[^\n] ");
  /* Input number of abscissae required, n */
  scanf("%" NAG_IFMT "%*[^\n] ", &n);
  /* Input quadrature rule to simulate, Nag_QuadType */
  scanf("%39s%*[^\n] ", nag_enum_arg);
  quadtype = (Nag_QuadType)nag_enum_name_to_value(nag_enum_arg);

  /* Allocate coefficient, weight and abscissae arrays */
  if (!(a = NAG_ALLOC(n, double)) || !(b = NAG_ALLOC(n, double)) ||
      !(c = NAG_ALLOC(n, double)) || !(weight = NAG_ALLOC(n, double)) ||
      !(abscis = NAG_ALLOC(n, double))) {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }

  /* Convert QuadType to string using
   * nag_enum_value_to_name (x04nbc).
   * Converts NAG enum member value to its name
   */
  str_type = nag_enum_value_to_name(quadtype);
  printf("\nQuadrature type = %s\n\n", str_type);

  /* Generate recursion coefficients a, b, c from quadtype */
  switch (quadtype) {
  case Nag_Quad_Gauss_Legendre: {
    muzero = 2.0;
    for (i = 0; i < n; ++i) {
      ri = (double)(i);
      a[i] = (2.0 * ri + 1.0) / (ri + 1.0);
      b[i] = 0.0;
      c[i] = ri / (ri + 1.0);
    }
  } break;
  case Nag_Quad_Gauss_Chebyshev_first: {
    muzero = X01AAC;
    for (i = 0; i < n; ++i) {
      a[i] = 2.0;
      b[i] = 0.0;
      c[i] = 1.0;
    }
    a[0] = 1.0;
  } break;
  case Nag_Quad_Gauss_Chebyshev_second: {
    muzero = 0.5 * X01AAC;
    for (i = 0; i < n; ++i) {
      a[i] = 2.0;
      b[i] = 0.0;
      c[i] = 1.0;
    }
  } break;
  case Nag_Quad_Gauss_Jacobi: {
    /* Input quadrature paramaters alpha and beta */
    scanf("%lf%lf\n", &alpha, &beta);
    printf("Using parameters alpha = %10.5f, beta = %10.5f\n\n", alpha, beta);
    a1 = alpha + 1.0;
    b1 = beta + 1.0;
    ab1 = alpha + beta + 1.0;
    ab2 = a1 + b1;
    INIT_FAIL(fail1);
    INIT_FAIL(fail2);
    /* nag_specfun_gamma (s14aac).
     * Gamma function Gamma(x)
     */
    fa1 = nag_specfun_gamma(a1, &fail);
    fb1 = nag_specfun_gamma(b1, &fail1);
    fab2 = nag_specfun_gamma(ab2, &fail2);
    if (fail.code != NE_NOERROR || fail1.code != NE_NOERROR ||
        fail2.code != NE_NOERROR) {
      printf("Error from nag_specfun_gamma (s14aac).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }
    muzero = pow(2.0, ab1) * fa1 * fb1 / fab2;
    a[0] = 0.5 * ab2;
    b[0] = 0.5 * (alpha - beta);
    c[0] = 0.0;
    for (i = 1; i < n; ++i) {
      ri = (double)i;
      d = (2.0 * ri + 2.0) * (ri + ab1);
      a[i] = (2.0 * ri + ab1) * (2.0 * ri + ab2) / d;
      d = (2.0 * ri + alpha + beta) * d;
      b[i] = (2.0 * ri + ab1) * (alpha * alpha - beta * beta) / d;
      c[i] = 2.0 * (ri + alpha) * (ri + beta) * (2.0 * ri + ab2) / d;
    }
  } break;
  case Nag_Quad_Gauss_Laguerre: {
    /* Input quadrature paramaters alpha */
    scanf("%lf\n", &alpha);
    printf("Using parameter alpha = %10.5f\n\n", alpha);
    a1 = alpha + 1.0;
    /* nag_specfun_gamma (s14aac).
     * gamma(x)
     */
    muzero = nag_specfun_gamma(a1, &fail);
    if (fail.code != NE_NOERROR) {
      printf("Error from nag_specfun_gamma (s14aac).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }
    for (i = 0; i < n; ++i) {
      ri = (double)i;
      a[i] = -1.0 / (ri + 1.0);
      b[i] = (2.0 * ri + a1) / (ri + 1.0);
      c[i] = (ri + alpha) / (ri + 1.0);
    }
  } break;
  case Nag_Quad_Gauss_Hermite: {
    muzero = sqrt(X01AAC);
    for (i = 0; i < n; ++i) {
      a[i] = 2.0;
      b[i] = 0.0;
      c[i] = (double)2 * i;
    }
  } break;
  default: {
    exit_status = 1;

    printf("The quadrature type %s is not handled here\n", str_type);
    goto END;
  }
  }

  /* nag_quad_dim1_gauss_wrec (d01tdc).
   * Compute weights and abscissae for a Gaussian quadrature rule
   * governed by a three-term recurrence relation.
   */
  nag_quad_dim1_gauss_wrec(n, a, b, c, muzero, weight, abscis, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_quad_dim1_gauss_wrec (d01tdc).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }

  printf(" %3" NAG_IFMT " points\n\n    Abscissae        Weights\n\n", n);
  for (i = 0; i < n; i++) {
    printf("%15.6e", abscis[i]);
    printf("%15.6e\n", weight[i]);
  }
  printf("\n");

END:
  NAG_FREE(a);
  NAG_FREE(b);
  NAG_FREE(c);
  NAG_FREE(abscis);
  NAG_FREE(weight);

  return exit_status;
}