/* nag_fit_dim1_cheb_arb (e02adc) Example Program.
 *
 * Copyright 2019 Numerical Algorithms Group.
 *
 * Mark 27.0, 2019.
 *
 */

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

int main(void)
{
#define A(I, J) a[(I) *tda + J]

  Integer exit_status = 0, i, iwght, j, k, m, r, tda;
  NagError fail;
  double *a = 0, *ak = 0, d1, fit, *s = 0, *w = 0, *x = 0, x1, xarg, xcapr,
         xm, *y = 0;

  INIT_FAIL(fail);

  printf("nag_fit_dim1_cheb_arb (e02adc) Example Program Results \n");

  /* Skip heading in data file */
  scanf("%*[^\n]");
  while ((scanf("%" NAG_IFMT "", &m)) != EOF)
  {
    if (m >= 2) {
      if (!(x = NAG_ALLOC(m, double)) ||
          !(y = NAG_ALLOC(m, double)) || !(w = NAG_ALLOC(m, double)))
      {
        printf("Allocation failure\n");
        exit_status = -1;
        goto END;
      }
    }
    else {
      printf("Invalid m.\n");
      exit_status = 1;
      return exit_status;
    }
    scanf("%" NAG_IFMT "", &k);
    if (k >= 1) {
      if (!(a = NAG_ALLOC((k + 1) * (k + 1), double)) ||
          !(s = NAG_ALLOC(k + 1, double)) || !(ak = NAG_ALLOC(k + 1, double)))
      {
        printf("Allocation failure\n");
        exit_status = -1;
        goto END;
      }
      tda = k + 1;
    }
    else {
      printf("Invalid k.\n");
      exit_status = 1;
      return exit_status;
    }
    scanf("%" NAG_IFMT "", &iwght);
    for (r = 0; r < m; ++r) {
      if (iwght != 1) {
        scanf("%lf", &x[r]);
        scanf("%lf", &y[r]);
        scanf("%lf", &w[r]);
      }
      else {
        scanf("%lf", &x[r]);
        scanf("%lf", &y[r]);
        w[r] = 1.0;
      }
    }
    /* nag_fit_dim1_cheb_arb (e02adc).
     * Computes the coefficients of a Chebyshev series
     * polynomial for arbitrary data
     */
    nag_fit_dim1_cheb_arb(m, k + 1, tda, x, y, w, a, s, &fail);
    if (fail.code != NE_NOERROR) {
      printf("Error from nag_fit_dim1_cheb_arb (e02adc).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }

    for (i = 0; i <= k; ++i) {
      printf("\n");
      printf(" %s%4" NAG_IFMT "%s%12.2e\n", "Degree", i,
             "   R.M.S. residual =", s[i]);
      printf("\n   J  Chebyshev coeff A(J) \n");
      for (j = 0; j < i + 1; ++j)
        printf(" %3" NAG_IFMT "%15.4f\n", j + 1, A(i, j));
    }
    for (j = 0; j < k + 1; ++j)
      ak[j] = A(k, j);
    x1 = x[0];
    xm = x[m - 1];
    printf("\n %s%4" NAG_IFMT "\n\n",
           "Polynomial approximation and residuals for degree", k);
    printf("   R   Abscissa     Weight   Ordinate  Polynomial  Residual \n");
    for (r = 1; r <= m; ++r) {
      xcapr = (x[r - 1] - x1 - (xm - x[r - 1])) / (xm - x1);
      /* nag_fit_dim1_cheb_eval (e02aec).
       * Evaluates the coefficients of a Chebyshev series
       * polynomial
       */
      nag_fit_dim1_cheb_eval(k + 1, ak, xcapr, &fit, &fail);
      if (fail.code != NE_NOERROR) {
        printf("Error from nag_fit_dim1_cheb_eval (e02aec).\n%s\n", fail.message);
        exit_status = 1;
        goto END;
      }

      d1 = fit - y[r - 1];
      printf(" %3" NAG_IFMT "%11.4f%11.4f%11.4f%11.4f%11.2e\n", r, x[r - 1],
             w[r - 1], y[r - 1], fit, d1);
      if (r < m) {
        xarg = (x[r - 1] + x[r]) * 0.5;
        xcapr = (xarg - x1 - (xm - xarg)) / (xm - x1);
        /* nag_fit_dim1_cheb_eval (e02aec), see above. */
        nag_fit_dim1_cheb_eval(k + 1, ak, xcapr, &fit, &fail);
        if (fail.code != NE_NOERROR) {
          printf("Error from nag_fit_dim1_cheb_eval (e02aec).\n%s\n", fail.message);
          exit_status = 1;
          goto END;
        }
        printf("    %11.4f                      %11.4f\n", xarg, fit);
      }
    }
  END:
    NAG_FREE(a);
    NAG_FREE(x);
    NAG_FREE(y);
    NAG_FREE(w);
    NAG_FREE(s);
    NAG_FREE(ak);
  }
  return exit_status;
}