/* nag_fit_dim1_cheb_glp (e02afc) Example Program.
 *
 * Copyright 2023 Numerical Algorithms Group.
 *
 * Mark 29.2, 2023.
 *
 */

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

int main(void) {
  Integer exit_status = 0;
  double *an = 0, d1, *f = 0, fit, pi, piby2n, *xcap = 0;

  Integer i, j, n;
  Integer r;
  NagError fail;

  INIT_FAIL(fail);

  printf("nag_fit_dim1_cheb_glp (e02afc) Example Program Results \n");

  /* Skip heading in data file */
  scanf("%*[^\n]");

  /* nag_math_pi (x01aac).
   * pi
   */
  pi = nag_math_pi;
  while ((scanf("%" NAG_IFMT "", &n)) != EOF)
  {
    if (n > 0) {
      if (!(an = NAG_ALLOC(n + 1, double)) || !(f = NAG_ALLOC(n + 1, double)) ||
          !(xcap = NAG_ALLOC(n + 1, double))) {
        printf("Allocation failure\n");
        exit_status = -1;
        goto END;
      }
    } else {
      printf("Invalid n.\n");
      exit_status = 1;
      return exit_status;
    }

    piby2n = pi * 0.5 / (double)n;
    for (r = 0; r < n + 1; ++r)
      scanf("%lf", &f[r]);

    for (r = 0; r < n + 1; ++r) {
      i = r;
      /* The following method of evaluating  xcap = cos(pi*i/n)
       *  ensures that the computed value has a small relative error
       *  and, moreover, is bounded in modulus by unity for all
       *  i = 0, 1, ..., n.  (It is assumed that the sine routine
       *  produces a result with a small relative error for values
       *  of the argument between  -PI/4  and  PI/4).
       */
      if (2 * i <= n) {
        d1 = sin(piby2n * i);
        xcap[i] = 1.0 - d1 * d1 * 2.0;
      } else if (2 * i > n * 3) {
        d1 = sin(piby2n * (n - i));
        xcap[i] = d1 * d1 * 2.0 - 1.0;
      } else {
        xcap[i] = sin(piby2n * (n - 2 * i));
      }
    }
    /* nag_fit_dim1_cheb_glp (e02afc).
     * Computes the coefficients of a Chebyshev series
     * polynomial for interpolated data
     */
    nag_fit_dim1_cheb_glp(n + 1, f, an, &fail);
    if (fail.code != NE_NOERROR) {
      printf("Error from nag_fit_dim1_cheb_glp (e02afc).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }

    printf("\n          Chebyshev \n");
    printf("   J   coefficient A(J) \n");
    for (j = 0; j < n + 1; ++j)
      printf(" %3" NAG_IFMT "%14.7f\n", j + 1, an[j]);
    printf("\n   R    Abscissa   Ordinate      Fit \n");
    for (r = 0; r < n + 1; ++r) {
      /* nag_fit_dim1_cheb_eval (e02aec).
       * Evaluates the coefficients of a Chebyshev series
       * polynomial
       */
      nag_fit_dim1_cheb_eval(n + 1, an, xcap[r], &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(" %3" NAG_IFMT "%11.4f%11.4f%11.4f\n", r + 1, xcap[r], f[r], fit);
    }
  END:
    NAG_FREE(an);
    NAG_FREE(f);
    NAG_FREE(xcap);
  }
  return exit_status;
}