/* E02BC_P0W_F C++ Header Example Program.
 *
 * Copyright 2023 Numerical Algorithms Group.
 * Mark 29.0, 2023.
 */

#include <dco.hpp>
#include <iostream>
#include <math.h>
#include <nag.h>
#include <nagad.h>
#include <stdio.h>
#include <string>
using namespace std;

int main()
{
  Integer exit_status = 0;

  cout << "E02BC_P0W_F C++ Header Example Program\n";

  // Skip heading in data file
  string mystr;
  getline(cin, mystr);

  Integer ncap, m;
  cin >> ncap;
  cin >> m;

  Integer ncap7 = ncap + 7;

  if (m <= 0)
  {
    printf("Invalid m.\n");
    exit_status = 1;
    return exit_status;
  }
  if (ncap <= 0)
  {
    printf("Invalid ncap.\n");
    exit_status = 1;
    return exit_status;
  }

  // Allocate arrays
  double *c = 0, *lamda = 0;
  double *dsdl = 0, *dsdc = 0;
  c     = new double[ncap7];
  lamda = new double[ncap7];
  dsdc  = new double[2 * m * ncap7];
  dsdl  = new double[2 * m * ncap7];

  // Read knots and spline coefficients from file
  double tmp;
  for (int j = 0; j < ncap7; j++)
  {
    cin >> tmp;
    lamda[j] = tmp;
  }
  for (int j = 0; j < ncap + 3; j++)
  {
    cin >> tmp;
    c[j] = tmp;
  }

  cout << "\n        x             Spline    1st deriv";
  cout << "   2nd deriv   3rd deriv\n";
  cout.setf(ios::scientific, ios::floatfield);
  cout.setf(ios::right);
  cout.precision(4);

  Integer           ifail = 0;
  nag::ad::handle_t ad_handle;
  double            s[4], x;

  for (int i = 0; i < m; i++)
  {
    cin >> tmp;
    x = tmp;

    for (Integer left = 1; left <= 2; left++)
    {

      nag::ad::e02bc(ad_handle, ncap7, lamda, c, x, left, s, ifail);

      if (ifail != 0)
      {
        printf("\nError from nag::ad::e02bc.\n%" NAG_IFMT "\n", ifail);
        exit_status = 1;
        return exit_status;
      }

      double xv = x;
      cout.width(11);
      cout << xv;

      if (left == 1)
      {
        cout << "   Left";
      }
      else
      {
        cout << "  Right";
      }
      for (int l = 0; l < 4; l++)
      {
        tmp = s[l];
        cout.width(12);
        cout << tmp;
      }

      cout << endl;
    }
    cout << endl;
  }

  delete[] c;
  delete[] lamda;
  delete[] dsdc;
  delete[] dsdl;

  return exit_status;
}