/* nag_zeros_poly_real_fpml (c02abc) Example Program.
 *
 * Copyright 2022 Numerical Algorithms Group.
 *
 * Mark 28.3, 2022.
 */

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

/* Declarations */
int ex1_basic(void);
int ex2_polishing(void);

int main(void) {

  int exit_status = 0;
  int polish_example = 0;

  printf("nag_zeros_poly_real_fpml (c02abc) Example Program Results\n");

  exit_status += ex1_basic();
  if (polish_example)
    exit_status += ex2_polishing();

  return exit_status;
}

int ex1_basic() {

  Complex *z = 0;
  double *a = 0, *berr = 0, *cond = 0;
  Integer *conv = 0, exit_status = 0, i, itmax, n;
  /* Nag Types */
  Nag_Root_Polish polish;
  NagError fail;

  printf("\n Basic Problem\n\n");

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

  /* Read polynomial degree and allocate */
  scanf("%" NAG_IFMT " %*[^\n]", &n);
  if (n < 0) {
    printf(" Invalid n = %" NAG_IFMT "\n", n);
    return -1;
  }
  if (!(a = NAG_ALLOC(n + 1, double)) || !(berr = NAG_ALLOC(n, double)) ||
      !(cond = NAG_ALLOC(n, double)) || !(conv = NAG_ALLOC(n, Integer)) ||
      !(z = NAG_ALLOC(n, Complex))) {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }

  /* Read polynomial coefficients */
  for (i = 0; i <= n; i++) {
    scanf(" %lf%*[^\n]", &a[i]);
  }

  /* Find roots of the polynomial */
  itmax = 30;
  polish = Nag_Root_Polish_Simple;
  INIT_FAIL(fail);
  nag_zeros_poly_real_fpml(a, n, itmax, polish, z, berr, cond, conv, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_zeros_poly_real_fpml (c02abc).\n%s\n", fail.message);
    exit_status = -1;
    goto END;
  }

  /* Print output */
  printf("  i    z                      conv  berr      cond\n");
  printf(" -----------------------------------------------------\n");
  for (i = 0; i < n; i++) {
    printf(" %2" NAG_IFMT "  ( %9.2e, %9.2e)  %3" NAG_IFMT "  %9.2e %9.2e\n", i,
           z[i].re, z[i].im, conv[i], berr[i], cond[i]);
  }

END:

  NAG_FREE(a);
  NAG_FREE(berr);
  NAG_FREE(cond);
  NAG_FREE(conv);
  NAG_FREE(z);
  return exit_status;
}

int ex2_polishing() {

  const double eps = nag_machine_precision, rmax = nag_machine_real_largest;
  const double zero = 0.0;

  char polish_names[3][7] = {"None  ", "Simple", "Comp  "};
  Complex *ac, az, pz, *z = 0, *zact = 0;
  double *berr = 0, *cond = 0, *a = 0, delta, err, fwderr, maxfwderr, maxrelerr,
         relerr;
  Integer *conv = 0, exit_status = 0, i, j, k, n, p;
  /* Nag Types */
  Nag_Boolean *matched = 0;
  Nag_Root_Polish polish[3] = {Nag_Root_Polish_None, Nag_Root_Polish_Simple,
                               Nag_Root_Polish_Compensated};
  NagError fail;

  printf("\n Polishing Modes\n\n");

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

  /* Read polynomial degree and allocate */
  scanf("%" NAG_IFMT " %*[^\n]", &n);
  if (n < 0) {
    printf(" Invalid n = %" NAG_IFMT "\n", n);
    return -1;
  }
  if (!(a = NAG_ALLOC(n + 1, double)) || !(berr = NAG_ALLOC(n, double)) ||
      !(cond = NAG_ALLOC(n, double)) || !(conv = NAG_ALLOC(n, Integer)) ||
      !(matched = NAG_ALLOC(n, Nag_Boolean)) || !(z = NAG_ALLOC(n, Complex)) ||
      !(zact = NAG_ALLOC(n, Complex))) {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }

  /* Set expected roots */
  for (i = 0; i < n; i++)
    zact[i] = nag_complex_create(i + 1, 0.0);

  /* Multiply out (z-1)...(z-n) for coefficients */
  INIT_FAIL(fail);
  nag_blast_dload(n, zero, a, 1, &fail);
  a[n] = 1.0;

  /* ac used for temporary complex operations*/
  ac = NAG_ALLOC(n + 1, Complex);

  for (i = 0; i < n + 1; i++) {
    ac[i] = nag_complex_create(a[i], 0.0);
  }

  for (i = 0; i < n; i++) {
    for (j = n - i; j <= n; j++) {
      az = nag_complex_multiply(ac[j - 1], zact[i]);
      ac[j - 1] = nag_complex_subtract(ac[j], az);
    }
    az = nag_complex_multiply(ac[n], zact[i]);
    ac[n] = nag_complex_negate(az);
  }

  for (i = 0; i < n + 1; i++) {
    a[i] = ac[i].re;
  }

  printf("  polish  relerr    fwderr\n");
  printf(" ----------------------------\n");

  for (p = 0; p <= 2; p++) {
    /* Find roots */
    INIT_FAIL(fail);
    nag_zeros_poly_real_fpml(a, n, 30, polish[p], z, berr, cond, conv, &fail);
    if (fail.code != NE_NOERROR) {
      printf("Error from nag_zeros_poly_real_fpml (c02abc).\n%s\n",
             fail.message);
      exit_status = -1;
      goto END;
    }

    /* Calculate the maximum relative errors of the roots, and the maximum
     * forward error evaluating the polynomial at those roots. Errors are
     * capped at machine precision. */
    maxfwderr = maxrelerr = eps;
    for (i = 0; i < n; i++)
      matched[i] = Nag_FALSE;
    for (i = 0; i < n; i++) {

      /* Evaluate polynomial at this root */
      pz = nag_complex_create(a[0], 0.0);
      for (j = 1; j <= n; j++) {
        pz = nag_complex_multiply(pz, z[i]);
        pz = nag_complex_add(pz, nag_complex_create(a[j], 0.0));
      }

      /* Match to an expected root */
      k = 0;
      err = rmax;
      for (j = 0; j < n; j++) {
        if (!matched[j]) {
          delta = nag_complex_abs(nag_complex_subtract(z[i], zact[j]));
          if (delta <= err) {
            err = delta;
            k = j;
          }
        }
      }

      /* Mark as matched and update max errors */
      matched[k] = Nag_TRUE;
      relerr = err / nag_complex_abs(zact[k]);
      fwderr = nag_complex_abs(pz);
      if (relerr > maxrelerr)
        maxrelerr = relerr;
      if (fwderr > maxfwderr)
        maxfwderr = fwderr;
    }

    /* Print output */
    printf("  %s%10.2e%10.2e\n", polish_names[p], maxrelerr, maxfwderr);
  }

END:

  NAG_FREE(a);
  NAG_FREE(berr);
  NAG_FREE(cond);
  NAG_FREE(conv);
  NAG_FREE(matched);
  NAG_FREE(z);
  NAG_FREE(zact);
  return exit_status;
}