// example of using:
//   nagcpp::quad::md_gauss (d01fb)
//   nagcpp::quad::dim1_gauss_wres (d01tb)
#include <iostream>
#include <math.h>
#include <vector>

#include "d01/nagcpp_d01fb.hpp"
#include "d01/nagcpp_d01tb.hpp"

// function of interest
double f(const std::vector<double> &x) {
  double p1 = 6;
  double p2 = 8;
  return pow(x[0] * x[1] * x[2], p1) / pow(x[3] + 2, p2) *
         exp(-2 * x[1] - 0.5 * x[2] * x[2]);
}

int main(void) {
  std::cout << "nagcpp::quad::md_gauss Example" << std::endl;

  // number of points per dimension
  std::vector<nagcpp::types::f77_integer> nptvec = {4, 4, 4, 4};

  // preallocate the arrays for the combined weights and abscis
  size_t lwa = 0;
  for (auto n : nptvec)
    lwa += n;
  std::vector<double> weight(lwa);
  std::vector<double> abscis(lwa);

  // calculate some weights and abscis ...
  std::vector<nagcpp::types::f77_integer> key = {0, -3, -4, -5};
  std::vector<double> a = {1.0, 0.0, 0.0, 1.0};
  std::vector<double> b = {2.0, 2.0, 0.5, 2.0};

  // loop over dimension
  size_t k = 0;
  for (auto i = 0u; i < nptvec.size(); ++i) {
    std::vector<double> this_weight;
    std::vector<double> this_abscis;

    try {
      nagcpp::quad::dim1_gauss_wres(key[i], a[i], b[i], nptvec[i], this_weight,
                                    this_abscis);

    } catch (nagcpp::error_handler::Exception &e) {
      std::cout << e.msg << std::endl;
      return 1;
    }

    // concatenate the weights and abscis
    for (auto j = 0; j < nptvec[i]; ++j, ++k) {
      weight[k] = this_weight[j];
      abscis[k] = this_abscis[j];
    }
  }
  // ... calculate some weights and abscis

  // approximate the integral of f
  double mdint;
  try {
    mdint = nagcpp::quad::md_gauss(nptvec, weight, abscis, f);

  } catch (nagcpp::error_handler::Exception &e) {
    std::cout << e.msg << std::endl;
    return 1;
  }

  std::cout.precision(5);
  std::cout << "Answer = " << mdint << std::endl;

  return 0;
}