/* nag_quad_dim1_gauss_recm (d01tec) Example Program.
 *
 * Copyright 2023 Numerical Algorithms Group.
 *
 * Mark 29.2, 2023.
 */

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

int main(void) {
  /* Scalars */
  Integer exit_status = 0;
  Integer n, i;
  double ri, muzero;
  /* Arrays */
  double *a = 0, *b = 0, *c = 0, *abscissae = 0, *weights = 0, *mu = 0;
  /* Nag Types */
  NagError fail;

  INIT_FAIL(fail);

  printf("nag_quad_dim1_gauss_recm (d01tec) Example Program Results\n");
  /* Skip heading in data file */
  scanf("%*[^\n] ");
  /* Input number of abscissae required, n */
  scanf("%" NAG_IFMT "%*[^\n] ", &n);

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

  /* Set up moments w.r.t w(x)x^j in array mu */
  mu[0] = 2.0;
  for (i = 1; i < 2 * n; i = i + 2) {
    ri = (double)(i + 2);
    mu[i] = 0.0;
    mu[i + 1] = 2.0 / ri;
  }

  /* nag_quad_dim1_gauss_recm (d01tec).
   * Compute three term recurrence coefficients from moments in w(x)x^j.
   */
  nag_quad_dim1_gauss_recm(n, mu, a, b, c, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_quad_dim1_gauss_recm (d01tec).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }

  printf("\n       a         b         c\n");
  for (i = 0; i < n; i++) {
    printf("%10.5f%10.5f%10.5f\n", a[i], b[i], c[i]);
  }

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

  printf("\n   weights       abscissae\n");
  for (i = 0; i < n; i++) {
    printf("%10.5f     %10.5f\n", weights[i], abscissae[i]);
  }

END:
  NAG_FREE(a);
  NAG_FREE(b);
  NAG_FREE(c);
  NAG_FREE(mu);
  NAG_FREE(abscissae);
  NAG_FREE(weights);

  return exit_status;
}