/* nag_numdiff_rcomm (d04bac) Example Program.
 *
 * Copyright 2023 Numerical Algorithms Group.
 *
 * Mark 29.0, 2023.
 */

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

int main(void) {
/* DER_COMP is used to store results for printing */
#define DER_COMP(I, J, K) der_comp[I + J * n_hbase + K * n_hbase * n_der_comp]

  Integer exit_status = 0;
  double x_0;
  Integer n_der_comp, n_display, n_hbase;
  double hbase;
  Integer i, j, k;
  double der[14], erest[14], fval[21], xval[21];
  double *actder = 0, *der_comp = 0;
  char *clabs = 0, **clabsc = 0;
  NagError fail;

  INIT_FAIL(fail);

  printf("nag_numdiff_rcomm (d04bac) Example Program Results\n");

  printf(
      "\n"
      "Find the derivatives of the polygamma (psi) function\n"
      "using function values generated by nag_specfun_psi_deriv_real (s14aec).\n\n"
      "Demonstrate the effect of successively reducing hbase.\n\n");
  fflush(stdout);

  n_der_comp = 3;
  n_display = 3;
  n_hbase = 4;

  if (!(clabs = NAG_ALLOC(n_der_comp * 10, char)) ||
      !(clabsc = NAG_ALLOC(n_der_comp, char *)) ||
      !(actder = NAG_ALLOC(n_display, double)) ||
      !(der_comp = NAG_ALLOC(n_hbase * n_der_comp * 14, double))) {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }

  x_0 = 0.05;

  /* Select an initial separation distance hbase. */
  hbase = 0.0025;

  /* Compute the actual derivatives at target location x_0
   * using nag_specfun_psi_deriv_real (s14aec), for comparison.
   */
  for (j = 0; j < n_display; j++) {
    /* nag_specfun_psi_deriv_real (s14aec).
     * Derivative of the psi function psi(x).
     */
    actder[j] = nag_specfun_psi_deriv_real(x_0, j + 1, &fail);
    if (fail.code != NE_NOERROR) {
      printf("Error from nag_specfun_psi_deriv_real (s14aec).\n%s\n",
             fail.message);
      exit_status = 1;
      goto END;
    }
  }

  /* Attempt n_hbase approximations, reducing hbase by factor 0.1 each time. */
  for (j = 0; j < n_hbase; j++) {
    /* nag_numdiff_sample (d04bbc).
     * Generates sample points for function evaluations
     * by nag_numdiff_rcomm (d04bac).
     */
    nag_numdiff_sample(x_0, hbase, xval);

    /* Calculate the corresponding objective function values. */
    for (k = 0; k < 21; k++) {
      fval[k] = nag_specfun_psi_deriv_real(xval[k], 0, &fail);
      if (fail.code != NE_NOERROR) {
        printf("Error from nag_specfun_psi_deriv_real (s14aec).\n%s\n",
               fail.message);
        exit_status = 1;
        goto END;
      }
    }

    /* Calculate the derivative estimates. */

    /* nag_numdiff_rcomm (d04bac).
     * Numerical differentiation, user-supplied function values,
     * derivatives up to order 14,
     * derivatives with respect to one real variable.
     */
    nag_numdiff_rcomm(xval, fval, der, erest, &fail);
    if (fail.code != NE_NOERROR) {
      printf("Error from nag_numdiff_rcomm (d04bac).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }

    /* Store results in DER_COMP. */
    for (i = 0; i < 14; i++) {
      DER_COMP(j, 0, i) = hbase;
      DER_COMP(j, 1, i) = der[i];
      DER_COMP(j, 2, i) = erest[i];
    }

    /* Decrease hbase for next loop. */
    hbase = hbase * 0.1;
  }

  /* Display results for first n_display derivatives. */
  for (j = 0; j < n_display; j++) {
    printf(" Derivative (%" NAG_IFMT ") calculated using "
           "nag_specfun_psi_deriv_real (s14aec): %13.4e\n",
           j + 1, actder[j]);
    printf(" Derivative and error estimates for derivative (%" NAG_IFMT ")\n",
           j + 1);
    printf("          hbase        der[%" NAG_IFMT "]      erest[%" NAG_IFMT
           "]\n",
           j, j);
    for (k = 0; k < n_hbase; k++)
      printf("%15.4e %13.4e %13.1e\n", DER_COMP(k, 0, j), DER_COMP(k, 1, j),
             DER_COMP(k, 2, j));
    printf("\n");
  }

END:
  NAG_FREE(actder);
  NAG_FREE(der_comp);
  NAG_FREE(clabs);
  NAG_FREE(clabsc);

  return exit_status;
}