/* nag_interp_dim1_cheb (e01aec) Example Program.
 *
 * Copyright 2023 Numerical Algorithms Group.
 *
 * Mark 29.2, 2023.
 */

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

int main(void) {
  /* Scalars */
  double xmax, xmin;
  Integer exit_status, i, pmax, ires, iy, j, k, m, n, itmin, itmax, num_iter;
  NagError fail;

  /* Arrays */
  double *a = 0, *perf = 0, *x = 0, *y = 0;
  Integer *p = 0;

  exit_status = 0;

  INIT_FAIL(fail);

  printf("nag_interp_dim1_cheb (e01aec) Example Program Results\n");

  /* Skip heading in data file */
  scanf("%*[^\n] ");
  while (scanf("%" NAG_IFMT "%lf%lf%*[^\n] ", &m, &xmin, &xmax) != EOF)
  {
    if (m > 0) {
      /* Allocate memory for p and x. */
      if (!(p = NAG_ALLOC(m, Integer)) || !(x = NAG_ALLOC(m, double))) {
        printf("Allocation failure\n");
        exit_status = -1;
        goto END;
      }

      /* Read p, x and y arrays */
      n = 0;
      pmax = 0;
      for (i = 1; i <= m; ++i) {
        scanf("%" NAG_IFMT "%lf", &p[i - 1], &x[i - 1]);
        k = n + p[i - 1] + 1;
        /* We need to extend array y as we go along */
        if (!(y = NAG_REALLOC(y, k, double))) {
          printf("Allocation failure\n");
          exit_status = -1;
          goto END;
        }
        for (j = n + 1; j <= k; ++j)
          scanf("%lf", &y[j - 1]);
        scanf("%*[^\n] ");
        if (p[i - 1] > pmax)
          pmax = p[i - 1];
        n = k;
      }

      /* Allocate array a */
      if (!(a = NAG_ALLOC(n, double)) ||
          !(perf = NAG_ALLOC(pmax + n + 1, double))) {
        printf("Allocation failure\n");
        exit_status = -1;
        goto END;
      }

      itmin = -1;
      itmax = -1;
      /* nag_interp_dim1_cheb (e01aec).
       * Interpolating functions, polynomial interpolant, data may
       * include derivative values, one variable
       */
      nag_interp_dim1_cheb(m, xmin, xmax, x, y, p, itmin, itmax, a, perf,
                           &num_iter, &fail);

      printf("\n");
      if (fail.code == NE_NOERROR) {
        printf("Total number of interpolating conditions ="
               " %" NAG_IFMT "\n",
               n);
        printf("\n");
        printf("Interpolating polynomial\n");
        printf("\n");
        printf("   i    Chebyshev Coefficient a(i+1)\n");

        for (i = 1; i <= n; ++i)
          printf("%4" NAG_IFMT "%20.4f\n", i - 1, a[i - 1]);

        printf("\n");

        printf("  x    R   Rth derivative    Residual\n");
        iy = 0;
        ires = pmax + 1;
        for (i = 1; i <= m; ++i) {
          for (j = 1; j <= p[i - 1] + 1; ++j) {
            ++iy;
            ++ires;
            if (j - 1 != 0)
              printf("    %4" NAG_IFMT "%12.1f%17.6f\n", j - 1, y[iy - 1],
                     perf[ires - 1]);
            else
              printf("%4.1f   0%12.1f%17.6f\n", x[i - 1], y[iy - 1],
                     perf[ires - 1]);
          }
        }
      } else {
        printf("Error from nag_interp_dim1_cheb (e01aec).\n%s\n", fail.message);
        exit_status = 1;
      }
    }
  }

END:
  NAG_FREE(a);
  NAG_FREE(x);
  NAG_FREE(y);
  NAG_FREE(p);
  NAG_FREE(perf);

  return exit_status;
}