#include <math.h>
#include <float.h>
#include <Rmath.h> /* for fmax2, fmin2, imin2 */
#include <RcppArmadillo.h>
#include <vector>

/*
  - Adapted from src/appl/integrate.c, which is a modified f2c translation of QUADPACK.
  - Principle: change as little of the code as possible
  - Re-factor to carefully use vec and mat -- but use C arrays if possible and loop when required
  - f2c uses 1-based indexes -- the re-factor uses some 0-based indexes (elist, rlist, elistSum)
  - As a reminder, function f: R -> R^ny for one integration point and ny function values
    (whereas integrate.c has f: R^n -> R^n for n integration points)
  - Potential issue: use of any() or all() in the conditions may be incorrect. 
    Principle: we want to be conservative.
  - Potential issue: limit*ny and 52*ny may be large:(
  - Limitation: in some cases, this is not much faster than using integrate() in an R for loop
  - TODO: can this be extended to allow for vectorised integration in C++?

  - GPL>=3
  - Mark Clements 2023-07-31
*/

using namespace arma;
using std::vector;

static void vrdqk21(const Rcpp::Function f,
		    const double, const double, vec *, vec *, vec *, vec *);
static void vrdqk15i(const Rcpp::Function f,
		     const double boun, // scalar or vec?
		     const int inf, double *a, double *b,
		     vec *result,
		     vec *abserr, vec *resabs, vec *resasc);

static void rdqpsrt(const int, int *, int *, double *, double *, int *, int *);

static void rdqelg(int *, double *, double *, double *, double *, int *);

static void vrdqagse(Rcpp::Function f, const double a, const double b, const double 
		     epsabs, const double epsrel, const int limit, const int ny,
		     double *result, double *abserr, int *neval, int *ier, double *alist,
		     double *blist, double *rlist, double *elist, int *
		     iord, int *last);

static void vrdqagie(Rcpp::Function f, const double bound, const int inf,
		     const double epsabs, const double epsrel, const int limit, const int ny,
		     double *resultp, double *abserrp, int *neval, int *ier, double *alist,
		     double *blist, double *rlistp, double *elistp, int *iord, int *last);

void vRdqagi(Rcpp::Function f, const double bound, const int inf,
	     const double epsabs, const double epsrel, const int limit, const int ny,
	     double *result, double *abserr, int *neval, int *ier,
	     int *lenw, int *last,
	     int *iwork, double *work)
{
  int l1, l2, l3;
  *ier = 6;
  *neval = 0;
  *last = 0;
  for (int i=0; i<ny; i++) {
    result[i] = 0.0;
    abserr[i] = 0.0;
  }
  if (limit < 1 || *lenw < limit*2 + limit*ny*2) return;

  /*         prepare call for dqagse. */

  // double work[2*limit + 2*limit*ny];
  l1 = limit;
  l2 = limit + l1;
  l3 = limit*ny + l2;

  vrdqagie(f, bound, inf, epsabs, epsrel, limit, ny,
	   result, abserr, neval, ier,
	   work, &work[l1], &work[l2], &work[l3], iwork, last);

  return;
} /* vRdqags_ */


static void vRdqags(Rcpp::Function f, const double a, const double b,
		    const double epsabs, const double epsrel, const int ny,
		    double *result, double *abserr, int *neval, int *ier,
		    const int limit, int *lenw, int *last, int *iwork, double *work)
{
  int l1, l2, l3;

  /*         check validity of limit and lenw. */

  *ier = 6;
  *neval = 0;
  *last = 0;
  for (int i=0; i<ny; i++) {
    result[i] = 0.0;
    abserr[i] = 0.0;
  }
  if (limit < 1 || *lenw < limit*2 + limit*ny*2) return;

  /*         prepare call for dqagse. */

  // double work[2*limit + 2*limit*ny];
  l1 = limit;
  l2 = limit + l1;
  l3 = limit*ny + l2;

  vrdqagse(f, a, b, epsabs, epsrel, limit, ny, result, abserr, neval, ier,
	   work, &work[l1], &work[l2], &work[l3], iwork, last);

  return;
} /* vRdqags_ */

static void vrdqagse(Rcpp::Function f, const double a, const double b, const double 
		     epsabs, const double epsrel, const int limit, const int ny,
		     double *resultp, double *abserrp, int *neval, int *ier, double *alist,
		     double *blist, double *rlistp, double *elistp, int *
		     iord, int *last)
{
  /* Local variables */
  Rboolean noext, extrap;
  int k,ksgn, nres;
  int ierro;
  int ktmin, nrmax;
  int iroff1, iroff2, iroff3;
  int id;
  int numrl2;
  int jupbnd;
  int maxerr;
  double res3la[3*ny];
  double rlist2[52*ny];
  vec abseps= zeros(ny), area1= zeros(ny), area2= zeros(ny), area12= zeros(ny);
  double epmach;
  double a1, a2, b1, b2, oflow, uflow;
  vec defab1= zeros(ny), defab2= zeros(ny), reseps= zeros(ny);
  vec area= zeros(ny), defabs= zeros(ny), resabs= zeros(ny), dres= zeros(ny), errbnd= zeros(ny);
  vec error1= zeros(ny), error2= zeros(ny), erro12= zeros(ny), correc= zeros(ny), erlarg= zeros(ny), errsum= zeros(ny);
  vec errmax= zeros(ny), erlast= zeros(ny), ertest= zeros(ny);
  double errmaxsum, small = 0.0;

  // vec result(resultp, limit, false); 
  // vec abserr(abserrp, limit, false); 
  mat elist(elistp, ny, limit, false);
  mat rlist(rlistp, ny, limit, false);
  double elistSum[limit];
    
  /* Parameter adjustments */
  --iord;
  // --elist;
  // --rlist;
  --blist;
  --alist;

  /* Function Body */
  epmach = DBL_EPSILON;

  /*            test on validity of parameters */
  /*            ------------------------------ */
  *ier = 0;
  *neval = 0;
  *last = 0;
  vec result = zeros(ny);
  vec abserr = zeros(ny);
  alist[1] = a;
  blist[1] = b;
  rlist.col(0) = zeros(ny);
  elist.col(0) = zeros(ny);
  elistSum[0] = 0.0;
  if (epsabs <= 0. && epsrel < R::fmax2(epmach * 50., 5e-29)) {
    *ier = 6;
    return;
  }

  /*           first approximation to the integral */
  /*           ----------------------------------- */

  uflow = DBL_MIN;
  oflow = DBL_MAX;
  ierro = 0;
  vrdqk21(f, a, b, &result, &abserr, &defabs, &resabs);

  /*           test on accuracy. */

  dres = abs(result);
  errbnd = max(ones(ny)*epsabs, epsrel * dres);
  *last = 1;
  rlist.col(0) = result;
  elist.col(0) = abserr;
  elistSum[0] = max(abserr);
  iord[1] = 1;
  if (all(abserr <= epmach * 100. * defabs && abserr > errbnd))
    *ier = 2;
  if (limit == 1)
    *ier = 1;
  if (*ier != 0 || all(abserr <= errbnd && abserr != resabs) // all or any??
      || all(abserr == 0.)) goto L140;

  /*           initialization */
  /*           -------------- */

  for (int i=0; i<ny; ++i)
    rlist2[i*52] = result[i];
  errmax = abserr;
  maxerr = 1; // base 1
  area = result;
  errsum = abserr;
  for (int i=0; i<ny; ++i)
    abserr[i] = oflow;
  nrmax = 1;
  nres = 0;
  numrl2 = 2;
  ktmin = 0;
  extrap = FALSE;
  noext = FALSE;
  iroff1 = 0;
  iroff2 = 0;
  iroff3 = 0;
  ksgn = -1;
  if (any(dres >= (1. - epmach * 50.) * defabs)) {
    ksgn = 1;
  }

  /*           main do-loop */
  /*           ------------ */

  for (*last = 2; *last <= limit; ++(*last)) {

    /*           bisect the subinterval with the nrmax-th largest error estimate. */

    a1 = alist[maxerr];
    b1 = (alist[maxerr] + blist[maxerr]) * .5;
    a2 = b1;
    b2 = blist[maxerr];
    erlast = errmax;
    vrdqk21(f, a1, b1, &area1, &error1, &resabs, &defab1);
    vrdqk21(f, a2, b2, &area2, &error2, &resabs, &defab2);

    /*           improve previous approximations to integral
		 and error and test for accuracy. */

    area12 = area1 + area2;
    erro12 = error1 + error2;
    errsum = errsum + erro12 - errmax;
    area = area + area12 - rlist.col(maxerr-1);
    if (all(defab1 != error1 && defab2 != error2)) {

      if (all(abs(rlist.col(maxerr-1) - area12) <= abs(area12) * 1e-5 &&
	      erro12 >= errmax * .99)) {
	if (extrap)
	  ++iroff2;
	else /* if(! extrap) */
	  ++iroff1;
      }
      if (*last > 10 && all(erro12 > errmax))
	++iroff3;
    }
    rlist.col(maxerr-1) = area1;
    rlist.col(*last-1) = area2;
    errbnd = max(ones(ny)*epsabs, epsrel * abs(area));

    /*           test for roundoff error and eventually set error flag. */
    if (iroff1 + iroff2 >= 10 || iroff3 >= 20)
      *ier = 2;
    if (iroff2 >= 5)
      ierro = 3;

    /* set error flag in the case that the number of subintervals equals limit. */
    if (*last == limit)
      *ier = 1;

    /*           set error flag in the case of bad integrand behaviour
		 at a point of the integration range. */

    if (R::fmax2(fabs(a1), fabs(b2)) <=
	(epmach * 100. + 1.) * (fabs(a2) + uflow * 1e3)) {
      *ier = 4;
    }

    /*           append the newly-created intervals to the list. */
    if (max(error2) > max(error1)) {
      alist[maxerr] = a2;
      alist[*last] = a1;
      blist[*last] = b1;
      rlist.col(maxerr-1) = area2;
      rlist.col(*last-1) = area1;
      elist.col(maxerr-1) = error2;
      elist.col(*last-1) = error1;
      elistSum[maxerr-1] = max(error2);
      elistSum[*last-1] = max(error1);
    } else {
      alist[*last] = a2;
      blist[maxerr] = b1;
      blist[*last] = b2;
      elist.col(maxerr-1) = error1;
      elist.col(*last-1) = error2;
      elistSum[maxerr-1] = max(error1);
      elistSum[*last-1] = max(error2);
    }

    /*           call subroutine dqpsrt to maintain the descending ordering
		 in the list of error estimates and select the subinterval
		 with nrmax-th largest error estimate (to be bisected next). */

    /*L30:*/
    // Constant values: limit, last, elist
    // Updated values: ermax, maxerr, iord, nrmax
    errmaxsum = max(errmax);
    rdqpsrt(limit, last, &maxerr, &errmaxsum, elistSum, &iord[1], &nrmax);
    errmax = elist.col(maxerr-1);

    if (all(errsum <= errbnd))   goto L115;/* ***jump out of do-loop */
    if (*ier != 0)		break;
    if (*last == 2)	{ /* L80: */
      small = fabs(b - a) * .375;
      erlarg = errsum;
      ertest = errbnd;
      for (int i=0; i<ny; ++i)
	rlist2[i*52+1] = area(i);
      continue;
    }
    if (noext)		continue;

    erlarg -= erlast;
    if (fabs(b1 - a1) > small) {
      erlarg += erro12;
    }
    if (!extrap) {

      /*          test whether the interval to be bisected next is the
		  smallest interval. */

      if (fabs(blist[maxerr] - alist[maxerr]) > small) {
	continue;
      }
      extrap = TRUE;
      nrmax = 2;
    }

    if (ierro != 3 && any(erlarg > ertest)) {

      /*           the smallest interval has the largest error.
		   before bisecting decrease the sum of the errors over the
		   larger intervals (erlarg) and perform extrapolation. */

      id = nrmax;
      jupbnd = *last;
      if (*last > limit / 2 + 2) {
	jupbnd = limit + 3 - *last;
      }
      for (k = id; k <= jupbnd; ++k) {
	maxerr = iord[nrmax];
	errmax = elist.col(maxerr-1);
	if (fabs(blist[maxerr] - alist[maxerr]) > small) {
	  goto L90;
	}
	++nrmax;
	/* L50: */
      }
    }
    /*           perform extrapolation.  L60: */

    ++numrl2;
    for (int i=0; i<ny; ++i)
      rlist2[i*52+(numrl2 - 1)] = area[i];
    for (int i=0; i<ny; ++i) {
      double resepsi = 0.0;
      double absepsi = 0.0;
      int numrl2i = numrl2;
      int nresi = nres;
      rdqelg(&numrl2i, &rlist2[i*52], &resepsi, &absepsi, &res3la[i*3], &nresi);
      reseps(i) = resepsi;
      abseps(i) = absepsi;
      if (i==ny-1) numrl2 = numrl2i;
    }
    ++ktmin;
    if (ktmin > 5 && all(abserr < errsum * .001)) {
      *ier = 5;
    }
    if (any(abseps < abserr)) {
      ktmin = 0;
      abserr = abseps;
      result = reseps;
      correc = erlarg;
      ertest = max(ones(ny)*epsabs, epsrel * abs(reseps));
      if (all(abserr <= ertest)) {
	break;
      }
    }

    /*           prepare bisection of the smallest interval.  L70: */

    if (numrl2 == 1) {
      noext = TRUE;
    }
    if (*ier == 5) {
      break;
    }
    maxerr = iord[1];
    errmax = elist.col(maxerr-1);
    nrmax = 1;
    extrap = FALSE;
    small *= .5;
    erlarg = errsum;
  L90:
    ;
  }


  /* L100:	set final result and error estimate. */
  /*		------------------------------------ */

  if (all(abserr == oflow)) 	goto L115; // ier should be ivec :(
  if (*ier + ierro != 0) {
    if (ierro == 3)
      abserr += correc;
    if (*ier == 0)
      *ier = 3;
    if (all(result == 0. || area == 0.)) {
      if (all(abserr > errsum)) 	goto L115;
      if (all(area == 0.)) 		goto L130;
    }
    else { /* L105:*/
      if (all(abserr / abs(result) > errsum / abs(area)))
	goto L115;
    }
  }

  /* L110: test on divergence. */
  if (!(ksgn == -1 && all(max(abs(result), abs(area)) <= defabs * .01)) &&
      all(.01 > result / area || result / area > 100. || errsum > abs(area))) {
    *ier = 5;
  }
  goto L130;

 L115:/*		compute global integral sum. */
  result = zeros(ny);
  for (k = 1; k <= *last; ++k)
    result += rlist.col(k-1);
  abserr = errsum;
 L130:
 L140:
  *neval = *last * 42 - 21;
  for (int i=0; i<ny; ++i) {
    resultp[i] = result[i];
    abserrp[i] = abserr[i];
  }
  return;
} /* vrdqagse_ */

static double c_b6 = 0.;
static double c_b7 = 1.;

static
void vrdqagie(Rcpp::Function f, const double bound, const int inf,
	      const double epsabs, const double epsrel, const int limit, const int ny,
	      double *resultp, double *abserrp, int *neval, int *ier, double *alist,
	      double *blist, double *rlistp, double *elistp, int *iord, int *last) {
  /* Local variables */
  vec area=zeros(ny), dres=zeros(ny);
  int ksgn;
  double boun;
  int nres;
  vec area1=zeros(ny), area2=zeros(ny), area12=zeros(ny), erro12=zeros(ny);
  int k;
  double small = 0.0;
  int ierro;
  double a1, a2, b1, b2, oflow;
  vec defab1=zeros(ny), defab2=zeros(ny);
  int ktmin, nrmax;
  double uflow;
  Rboolean noext;
  int iroff1, iroff2, iroff3;
  double res3la[3*ny];
  vec error1=zeros(ny), error2=zeros(ny);
  int id;
  double rlist2[52*ny];
  int numrl2;
  vec correc=zeros(ny), abseps=zeros(ny), errbnd=zeros(ny), resabs=zeros(ny), erlarg=zeros(ny), defabs=zeros(ny);
  double epmach;
  int jupbnd;
  vec erlast=zeros(ny), errmax=zeros(ny), reseps=zeros(ny);
  int maxerr;
  double errmaxsum;
  Rboolean extrap;
  vec ertest=zeros(ny), errsum=zeros(ny);
  vec result = zeros(ny), abserr = zeros(ny);

  mat elist(elistp, ny, limit, false);
  mat rlist(rlistp, ny, limit, false);
  double elistSum[limit];
    
  /* ***first executable statement  dqagie */
  /* Parameter adjustments */
  --iord;
  // --elist;
  // --rlist;
  --blist;
  --alist;

  /* Function Body */
  epmach = DBL_EPSILON;

  /*           test on validity of parameters */
  /*           ----------------------------- */

  *ier = 0;
  *neval = 0;
  *last = 0;
  alist[1] = 0.;
  blist[1] = 1.;
  rlist.col(0) = zeros(ny);
  elist.col(0) = zeros(ny);
  elistSum[0] = 0.0;
  iord[1] = 0;
  if (epsabs <= 0. && (epsrel < R::fmax2(epmach * 50., 5e-29))) *ier = 6;
  if (*ier == 6) return;

  /*           first approximation to the integral */
  /*           ----------------------------------- */

  /*         determine the interval to be mapped onto (0,1).
	     if inf = 2 the integral is computed as i = i1+i2, where
	     i1 = integral of f over (-infinity,0),
	     i2 = integral of f over (0,+infinity). */

  boun = bound;
  if (inf == 2) {
    boun = 0.;
  }
  vrdqk15i(f, boun, inf, &c_b6, &c_b7, &result, &abserr, &defabs, &resabs);

  /*           test on accuracy */

  *last = 1;
  rlist.col(0) = result;
  elist.col(0) = abserr;
  elistSum[0] = max(abserr);
  iord[1] = 1;
  dres = abs(result);
  errbnd = max(ones(ny)*epsabs, epsrel * dres);
  if (all(abserr <= epmach * 100. * defabs && abserr > errbnd)) *ier = 2;
  if (limit == 1) *ier = 1;
  if (*ier != 0 || all(abserr <= errbnd && abserr != resabs)
      || all(abserr == 0.)) goto L130;

  /*           initialization */
  /*           -------------- */

  uflow = DBL_MIN;
  oflow = DBL_MAX;
  for (int i=0; i<ny; ++i)
    rlist2[i*52] = result[i];
  errmax = abserr;
  maxerr = 1;
  area = result;
  errsum = abserr;
  for (int i=0; i<ny; ++i)
    abserr[i] = oflow;
  nrmax = 1;
  nres = 0;
  ktmin = 0;
  numrl2 = 2;
  extrap = FALSE;
  noext = FALSE;
  ierro = 0;
  iroff1 = 0;
  iroff2 = 0;
  iroff3 = 0;
  ksgn = -1;
  if (any(dres >= (1. - epmach * 50.) * defabs)) { // any or all or do separately?
    ksgn = 1;
  }

  /*           main do-loop */
  /*           ------------ */

  for (*last = 2; *last <= limit; ++(*last)) {

    /*           bisect the subinterval with nrmax-th largest error estimate. */

    a1 = alist[maxerr];
    b1 = (alist[maxerr] + blist[maxerr]) * .5;
    a2 = b1;
    b2 = blist[maxerr];
    erlast = errmax;
    vrdqk15i(f, boun, inf, &a1, &b1, &area1, &error1, &resabs, &defab1);
    vrdqk15i(f, boun, inf, &a2, &b2, &area2, &error2, &resabs, &defab2);

    /*           improve previous approximations to integral
		 and error and test for accuracy. */

    area12 = area1 + area2;
    erro12 = error1 + error2;
    errsum = errsum + erro12 - errmax;
    area = area + area12 - rlist.col(maxerr-1);
    if (all(defab1 != error1 && defab2 != error2)) {
      if (all(abs(rlist.col(maxerr-1) - area12) <= abs(area12) * 1e-5 &&
	      erro12 >= errmax * .99)) {
	if (extrap)
	  ++iroff2;
	else /* if (! extrap) */
	  ++iroff1;
      }
      if (*last > 10 && all(erro12 > errmax))
	++iroff3;
    }

    rlist.col(maxerr-1) = area1;
    rlist.col(*last-1) = area2;
    errbnd = max(ones(ny)*epsabs, epsrel * abs(area));

    /*           test for roundoff error and eventually set error flag. */

    if (iroff1 + iroff2 >= 10 || iroff3 >= 20)
      *ier = 2;
    if (iroff2 >= 5)
      ierro = 3;

    /*           set error flag in the case that the number of
		 subintervals equals limit. */

    if (*last == limit)
      *ier = 1;

    /*           set error flag in the case of bad integrand behaviour
		 at some points of the integration range. */

    if (R::fmax2(fabs(a1), fabs(b2)) <=
	(epmach * 100. + 1.) * (fabs(a2) + uflow * 1e3))
      {
	*ier = 4;
      }

    /*           append the newly-created intervals to the list. */

    if (max(error2) <= max(error1)) {
      alist[*last] = a2;
      blist[maxerr] = b1;
      blist[*last] = b2;
      elist.col(maxerr-1) = error1;
      elist.col(*last-1) = error2;
      elistSum[maxerr-1] = max(error1);
      elistSum[*last-1] = max(error2);
    }
    else {
      alist[maxerr] = a2;
      alist[*last] = a1;
      blist[*last] = b1;
      rlist.col(maxerr-1) = area2;
      rlist.col(*last-1) = area1;
      elist.col(maxerr-1) = error2;
      elist.col(*last-1) = error1;
      elistSum[maxerr-1] = max(error2);
      elistSum[*last-1] = max(error1);
    }

    /*           call subroutine dqpsrt to maintain the descending ordering
		 in the list of error estimates and select the subinterval
		 with nrmax-th largest error estimate (to be bisected next). */

    errmaxsum = max(errmax);
    rdqpsrt(limit, last, &maxerr, &errmaxsum, elistSum, &iord[1], &nrmax);
    errmax = elist.col(maxerr-1);
    if (all(errsum <= errbnd)) {
      goto L115;
    }
    if (*ier != 0)	    break;
    if (*last == 2) { /* L80: */
      small = .375;
      erlarg = errsum;
      ertest = errbnd;
      for (int i=0; i<ny; ++i)
	rlist2[i*52+1] = area(i);
      continue;
    }
    if (noext) 	    continue;

    erlarg -= erlast;
    if (fabs(b1 - a1) > small) {
      erlarg += erro12;
    }
    if (!extrap) {

      /*           test whether the interval to be bisected next is the
		   smallest interval. */

      if (fabs(blist[maxerr] - alist[maxerr]) > small) {
	continue;
      }
      extrap = TRUE;
      nrmax = 2;
    }

    if (ierro != 3 && any(erlarg > ertest)) {

      /*	    the smallest interval has the largest error.
		    before bisecting decrease the sum of the errors over the
		    larger intervals (erlarg) and perform extrapolation. */

      id = nrmax;
      jupbnd = *last;
      if (*last > limit / 2 + 2) {
	jupbnd = limit + 3 - *last;
      }
      for (k = id; k <= jupbnd; ++k) {
	maxerr = iord[nrmax];
	errmax = elist.col(maxerr-1);
	if (fabs(blist[maxerr] - alist[maxerr]) > small) {
	  goto L90;
	}
	++nrmax;
	/* L50: */
      }
    }
    /*           perform extrapolation.  L60: */
    ++numrl2;
    for (int i=0; i<ny; ++i)
      rlist2[i*52+(numrl2 - 1)] = area[i];
    // rdqelg(&numrl2, rlist2, &reseps, &abseps, res3la, &nres);
    for (int i=0; i<ny; ++i) {
      double resepsi = 0.0;
      double absepsi = 0.0;
      int numrl2i = numrl2;
      int nresi = nres;
      rdqelg(&numrl2i, &rlist2[i*52], &resepsi, &absepsi, &res3la[i*3], &nresi);
      reseps(i) = resepsi;
      abseps(i) = absepsi;
      if (i==ny-1) numrl2 = numrl2i;
    }

    ++ktmin;
    if (ktmin > 5 && all(abserr < errsum * .001)) {
      *ier = 5;
    }
    if (all(abseps >= abserr)) {
      goto L70;
    }
    ktmin = 0;
    abserr = abseps;
    result = reseps;
    correc = erlarg;
    ertest = max(ones(ny)*epsabs, epsrel * abs(reseps));
    if (all(abserr <= ertest)) {
      break;
    }

    /*            prepare bisection of the smallest interval. */

  L70:
    if (numrl2 == 1) {
      noext = TRUE;
    }
    if (*ier == 5) {
      break;
    }
    maxerr = iord[1];
    errmax = elist.col(maxerr-1);
    nrmax = 1;
    extrap = FALSE;
    small *= .5;
    erlarg = errsum;
  L90:
    ;
  }

  /* L100:     set final result and error estimate. */
  /*	     ------------------------------------ */

  if (all(abserr == oflow)) {
    goto L115;
  }
  if (*ier + ierro == 0) {
    goto L110;
  }
  if (ierro == 3) {
    abserr += correc;
  }
  if (*ier == 0) {
    *ier = 3;
  }
  if (all(result == 0. || area == 0.)) {
    if (all(abserr > errsum))
      goto L115;

    if (all(area == 0.))
      goto L130;
  }
  else { /* L105: */
    if (all(abserr / abs(result) > errsum / abs(area))) {
      goto L115;
    }
  }

  /*           test on divergence */
 L110:
  if (ksgn == -1 && all(max(abs(result), abs(area)) <= defabs * .01)) {
    goto L130;
  }
  if (all(.01 > result / area || result / area > 100. || errsum > abs(area))) {
    *ier = 6;
  }
  goto L130;

  /*           compute global integral sum. */

 L115:
  result = zeros(ny);
  for (k = 1; k <= *last; ++k)
    result += rlist.col(k-1);

  abserr = errsum;
 L130:
  *neval = *last * 30 - 15;
  if (inf == 2) {
    *neval <<= 1;
  }
  if (*ier > 2) {
    --(*ier);
  }
  for (int i=0; i<ny; ++i) {
    resultp[i] = result[i];
    abserrp[i] = abserr[i];
  }
  return;
} /* vrdqagie_ */


static vec eval_f(const Rcpp::Function f, const double x) {
  return Rcpp::as<vec>(Rcpp::wrap(f(x)));
}

static void vrdqk15i(const Rcpp::Function f,
		     const double boun, // scalar or vec?
		     const int inf, double *a, double *b,
		     vec *result,
		     vec *abserr, vec *resabs, vec *resasc)
{
  /* Initialized data */

  static double wg[8] = {
    0., .129484966168869693270611432679082,
    0., .27970539148927666790146777142378,
    0., .381830050505118944950369775488975,
    0., .417959183673469387755102040816327 };
  static double xgk[8] = {
    .991455371120812639206854697526329,
    .949107912342758524526189684047851,
    .864864423359769072789712788640926,
    .741531185599394439863864773280788,
    .58608723546769113029414483825873,
    .405845151377397166906606412076961,
    .207784955007898467600689403773245, 0. };
  static double wgk[8] = {
    .02293532201052922496373200805897,
    .063092092629978553290700663189204,
    .104790010322250183839876322541518,
    .140653259715525918745189590510238,
    .16900472663926790282658342659855,
    .190350578064785409913256402421014,
    .204432940075298892414161999234649,
    .209482141084727828012999174891714 };

  /* Local variables */
  double epmach, uflow, dinf, centr, hlgth;
  std::vector<vec> fv1(7), fv2(7);
  double vect[15], vec2[15];
  double tabsc1, tabsc2, absc, absc1, absc2;
  vec fval1, fval2, fc, resg, resk, fsum, reskh;
  int j, ny;

  /* ***first executable statement  dqk15i */
  epmach = DBL_EPSILON;
  uflow = DBL_MIN;
  dinf = (double) R::imin2(1, inf);

  centr = (*a + *b) * .5;
  hlgth = (*b - *a) * .5;
  tabsc1 = boun + dinf * (1. - centr) / centr;
  vect[0] = tabsc1;
  if (inf == 2) {
    vec2[0] = -tabsc1;
  }
  for (j = 1; j <= 7; ++j) {
    absc = hlgth * xgk[j - 1];
    absc1 = centr - absc;
    absc2 = centr + absc;
    tabsc1 = boun + dinf * (1. - absc1) / absc1;
    tabsc2 = boun + dinf * (1. - absc2) / absc2;
    vect[(j << 1) - 1] = tabsc1;
    vect[j * 2] = tabsc2;
    if (inf == 2) {
      vec2[(j << 1) - 1] = -tabsc1;
      vec2[j * 2] = -tabsc2;
    }
    /* L5: */
  }
  fval1 = eval_f(f, vect[0]);
  ny = fval1.size();
  if (inf == 2) fval2 = eval_f(f,vec2[0]);
  if (inf == 2) fval1 += fval2;
  fc = fval1 / centr / centr;

  /*           compute the 15-point kronrod approximation to
	       the integral, and estimate the error. */

  resg = wg[7] * fc;
  resk = wgk[7] * fc;
  *resabs = abs(resk);
  for (j = 1; j <= 7; ++j) {
    absc = hlgth * xgk[j - 1];
    absc1 = centr - absc;
    absc2 = centr + absc;
    tabsc1 = boun + dinf * (1. - absc1) / absc1;
    tabsc2 = boun + dinf * (1. - absc2) / absc2;
    fval1 = eval_f(f,vect[(j << 1) - 1]);
    fval2 = eval_f(f,vect[j * 2]);
    if (inf == 2) {
      fval1 += eval_f(f,vec2[(j << 1) - 1]);
      fval2 += eval_f(f,vec2[j * 2]);
    }
    fval1 = fval1 / absc1 / absc1;
    fval2 = fval2 / absc2 / absc2;
    fv1[j - 1] = fval1;
    fv2[j - 1] = fval2;
    fsum = fval1 + fval2;
    resg += wg[j - 1] * fsum;
    resk += wgk[j - 1] * fsum;
    *resabs += wgk[j - 1] * (abs(fval1) + abs(fval2));
    /* L10: */
  }
  reskh = resk * .5;
  *resasc = wgk[7] * abs(fc - reskh);
  for (j = 1; j <= 7; ++j) {
    *resasc += wgk[j - 1] * (abs(fv1[j - 1] - reskh) +
			     abs(fv2[j - 1] - reskh));
    /* L20: */
  }
  *result = resk * hlgth;
  *resasc *= hlgth;
  *resabs *= hlgth;
  *abserr = abs((resk - resg) * hlgth);
  for (int i=0; i<ny; ++i) {
    if ((*resasc)[i] != 0. && (*abserr)[i] != 0.) {
      (*abserr)[i] = (*resasc)[i] * R::fmin2(1.0, pow((*abserr)[i] * 200. / (*resasc)[i], 1.5));
    }
    if ((*resabs)[i] > uflow / (epmach * 50.)) {
      (*abserr)[i] = R::fmax2(epmach * 50. * (*resabs)[i], (*abserr)[i]);
    }
  }
  return;
} /* vrdqk15i_ */

static void vrdqk21(const Rcpp::Function f, const double a, const double b, vec *result,
		    vec *abserr, vec *resabs, vec *resasc)
{
  /* Initialized data */

  static double wg[5] = { .066671344308688137593568809893332,
    .149451349150580593145776339657697,
    .219086362515982043995534934228163,
    .269266719309996355091226921569469,
    .295524224714752870173892994651338 };
  static double xgk[11] = { .995657163025808080735527280689003,
    .973906528517171720077964012084452,
    .930157491355708226001207180059508,
    .865063366688984510732096688423493,
    .780817726586416897063717578345042,
    .679409568299024406234327365114874,
    .562757134668604683339000099272694,
    .433395394129247190799265943165784,
    .294392862701460198131126603103866,
    .14887433898163121088482600112972,0. };
  static double wgk[11] = { .011694638867371874278064396062192,
    .03255816230796472747881897245939,
    .05475589657435199603138130024458,
    .07503967481091995276704314091619,
    .093125454583697605535065465083366,
    .109387158802297641899210590325805,
    .123491976262065851077958109831074,
    .134709217311473325928054001771707,
    .142775938577060080797094273138717,
    .147739104901338491374841515972068,
    .149445554002916905664936468389821 };

  /* Local variables */
  std::vector<vec> fv1(10), fv2(10);
  double absc, vect[21];
  vec resg, resk, fsum, fval1, fval2;
  vec fc, reskh;
  double hlgth, centr, uflow;
  double epmach, dhlgth;
  int j, jtw, jtwm1, ny;

  epmach = DBL_EPSILON;
  uflow = DBL_MIN;

  centr = (a + b) * .5;
  hlgth = (b - a) * .5;
  dhlgth = fabs(hlgth);

  /*           compute the 21-point kronrod approximation to
	       the integral, and estimate the absolute error. */

  vect[0] = centr;
  for (j = 1; j <= 5; ++j) {
    jtw = j << 1;
    absc = hlgth * xgk[jtw - 1];
    vect[(j << 1) - 1] = centr - absc;
    /* L5: */
    vect[j * 2] = centr + absc;
  }
  for (j = 1; j <= 5; ++j) {
    jtwm1 = (j << 1) - 1;
    absc = hlgth * xgk[jtwm1 - 1];
    vect[(j << 1) + 9] = centr - absc;
    vect[(j << 1) + 10] = centr + absc;
  }
  fc = eval_f(f,vect[0]);
  ny = fc.size();
  resg = zeros(ny);
  resk = wgk[10] * fc;
  *resabs = abs(resk);
  for (j = 1; j <= 5; ++j) {
    jtw = j << 1;
    absc = hlgth * xgk[jtw - 1];
    fval1 = eval_f(f,vect[(j << 1) - 1]);
    fval2 = eval_f(f,vect[j * 2]);
    fv1[jtw - 1] = fval1;
    fv2[jtw - 1] = fval2;
    fsum = fval1 + fval2;
    resg += wg[j - 1] * fsum;
    resk += wgk[jtw - 1] * fsum;
    *resabs += wgk[jtw - 1] * (abs(fval1) + abs(fval2));
    /* L10: */
  }
  for (j = 1; j <= 5; ++j) {
    jtwm1 = (j << 1) - 1;
    absc = hlgth * xgk[jtwm1 - 1];
    fval1 = eval_f(f,vect[(j << 1) + 9]);
    fval2 = eval_f(f,vect[(j << 1) + 10]);
    fv1[jtwm1 - 1] = fval1;
    fv2[jtwm1 - 1] = fval2;
    fsum = fval1 + fval2;
    resk += wgk[jtwm1 - 1] * fsum;
    *resabs += wgk[jtwm1 - 1] * (abs(fval1) + abs(fval2));
    /* L15: */
  }
  reskh = resk * .5;
  *resasc = wgk[10] * abs(fc - reskh);
  for (j = 1; j <= 10; ++j) {
    *resasc += wgk[j - 1] * (abs(fv1[j - 1] - reskh) +
			     abs(fv2[j - 1] - reskh));
    /* L20: */
  }
  *result = resk * hlgth;
  *resabs *= dhlgth;
  *resasc *= dhlgth;
  *abserr = abs((resk - resg) * hlgth);
  for (int i=0; i<ny; ++i) {
    if ((*resasc)[i] != 0. && (*abserr)[i] != 0.) {
      (*abserr)[i] = (*resasc)[i] * R::fmin2(1.0, pow((*abserr)[i] * 200. / (*resasc)[i], 1.5));
    }
    if ((*resabs)[i] > uflow / (epmach * 50.)) {
      (*abserr)[i] = R::fmax2(epmach * 50. * (*resabs)[i], (*abserr)[i]);
    }
  }
  return;
} /* vrdqk21_ */

// Constant values: limit, last, elist
// Updated values: ermax, maxerr, iord, nrmax
static void rdqpsrt(const int limit, int *last, int *maxerr,
		    double *ermax, double *elist, int *iord, int *nrmax)
{
  /* Local variables */
  int i, j, k, ido, jbnd, isucc, jupbn;
  double errmin, errmax;

  /* Parameter adjustments */
  --iord;
  --elist;

  /* Function Body */

  /*           check whether the list contains more than
	       two error estimates. */
  if (*last <= 2) {
    iord[1] = 1;
    iord[2] = 2;
    goto Last;
  }
  /*           this part of the routine is only executed if, due to a
	       difficult integrand, subdivision increased the error
	       estimate. in the normal case the insert procedure should
	       start after the nrmax-th largest error estimate. */

  errmax = elist[*maxerr];
  if (*nrmax > 1) {
    ido = *nrmax - 1;
    for (i = 1; i <= ido; ++i) {
      isucc = iord[*nrmax - 1];
      if (errmax <= elist[isucc])
	break; /* out of for-loop */
      iord[*nrmax] = isucc;
      --(*nrmax);
      /* L20: */
    }
  }

  /*L30:       compute the number of elements in the list to be maintained
    in descending order. this number depends on the number of
    subdivisions still allowed. */
  if (*last > limit / 2 + 2)
    jupbn = limit + 3 - *last;
  else
    jupbn = *last;

  errmin = elist[*last];

  /*           insert errmax by traversing the list top-down,
	       starting comparison from the element elist(iord(nrmax+1)). */

  jbnd = jupbn - 1;
  for (i = *nrmax + 1; i <= jbnd; ++i) {
    isucc = iord[i];
    if (errmax >= elist[isucc]) {/* ***jump out of do-loop */
      /* L60: insert errmin by traversing the list bottom-up. */
      iord[i - 1] = *maxerr;
      for (j = i, k = jbnd; j <= jbnd; j++, k--) {
	isucc = iord[k];
	if (errmin < elist[isucc]) {
	  /* goto L80; ***jump out of do-loop */
	  iord[k + 1] = *last;
	  goto Last;
	}
	iord[k + 1] = isucc;
      }
      iord[i] = *last;
      goto Last;
    }
    iord[i - 1] = isucc;
  }

  iord[jbnd] = *maxerr;
  iord[jupbn] = *last;

 Last:/* set maxerr and ermax. */

  *maxerr = iord[*nrmax];
  *ermax = elist[*maxerr];
  return;
} /* rdqpsrt_ */

// Constant values:
// Updated values: n, epstab, result, abserr, res3la, nres
static void rdqelg(int *n, double *epstab, double *
		   result, double *abserr, double *res3la, int *nres)
{
  /* Local variables */
  int i__, indx, ib, ib2, ie, k1, k2, k3, num, newelm, limexp;
  double delta1, delta2, delta3, e0, e1, e1abs, e2, e3, epmach, epsinf;
  double oflow, ss, res;
  double errA, err1, err2, err3, tol1, tol2, tol3;

  /* Parameter adjustments */
  --res3la;
  --epstab;

  /* Function Body */
  epmach = DBL_EPSILON;
  oflow = DBL_MAX;
  ++(*nres);
  *abserr = oflow;
  *result = epstab[*n];
  if (*n < 3) {
    goto L100;
  }
  limexp = 50;
  epstab[*n + 2] = epstab[*n];
  newelm = (*n - 1) / 2;
  epstab[*n] = oflow;
  num = *n;
  k1 = *n;
  for (i__ = 1; i__ <= newelm; ++i__) {
    k2 = k1 - 1;
    k3 = k1 - 2;
    res = epstab[k1 + 2];
    e0 = epstab[k3];
    e1 = epstab[k2];
    e2 = res;
    e1abs = fabs(e1);
    delta2 = e2 - e1;
    err2 = fabs(delta2);
    tol2 = R::fmax2(fabs(e2), e1abs) * epmach;
    delta3 = e1 - e0;
    err3 = fabs(delta3);
    tol3 = R::fmax2(e1abs, fabs(e0)) * epmach;
    if (err2 <= tol2 && err3 <= tol3) {
      /*           if e0, e1 and e2 are equal to within machine
		   accuracy, convergence is assumed. */
      *result = res;/*		result = e2 */
      *abserr = err2 + err3;/*	abserr = fabs(e1-e0)+fabs(e2-e1) */

      goto L100;	/* ***jump out of do-loop */
    }

    e3 = epstab[k1];
    epstab[k1] = e1;
    delta1 = e1 - e3;
    err1 = fabs(delta1);
    tol1 = R::fmax2(e1abs, fabs(e3)) * epmach;

    /*           if two elements are very close to each other, omit
		 a part of the table by adjusting the value of n */

    if (err1 > tol1 && err2 > tol2 && err3 > tol3) {
      ss = 1. / delta1 + 1. / delta2 - 1. / delta3;
      epsinf = fabs(ss * e1);

      /*           test to detect irregular behaviour in the table, and
		   eventually omit a part of the table adjusting the value of n. */

      if (epsinf > 1e-4) {
	goto L30;
      }
    }

    *n = i__ + i__ - 1;
    goto L50;/* ***jump out of do-loop */


  L30:/* compute a new element and eventually adjust the value of result. */

    res = e1 + 1. / ss;
    epstab[k1] = res;
    k1 += -2;
    errA = err2 + fabs(res - e2) + err3;
    if (errA <= *abserr) {
      *abserr = errA;
      *result = res;
    }
  }

  /*           shift the table. */

 L50:
  if (*n == limexp) {
    *n = (limexp / 2 << 1) - 1;
  }

  if (num / 2 << 1 == num) ib = 2; else ib = 1;
  ie = newelm + 1;
  for (i__ = 1; i__ <= ie; ++i__) {
    ib2 = ib + 2;
    epstab[ib] = epstab[ib2];
    ib = ib2;
  }
  if (num != *n) {
    indx = num - *n + 1;
    for (i__ = 1; i__ <= *n; ++i__) {
      epstab[i__] = epstab[indx];
      ++indx;
    }
  }
  /*L80:*/
  if (*nres >= 4) {
    /* L90: */
    *abserr = fabs(*result - res3la[3]) +
      fabs(*result - res3la[2]) +
      fabs(*result - res3la[1]);
    res3la[1] = res3la[2];
    res3la[2] = res3la[3];
    res3la[3] = *result;
  } else {
    res3la[*nres] = *result;
    *abserr = oflow;
  }

 L100:/* compute error estimate */
  *abserr = R::fmax2(*abserr, epmach * 5. * fabs(*result));
  return;
} /* rdqelg_ */


// [[Rcpp::export]]
Rcpp::List vdqagsRcpp(const Rcpp::Function f, const double a, const double b,
		      const double epsrel, const double epsabs, const int limit,
		      const int ny) {
// RcppExport SEXP vdqagsRcpp(SEXP _f, SEXP _a, SEXP _b,
// 			   SEXP _epsrel, SEXP _epsabs, SEXP _limit,
// 			   SEXP _ny) {
  using namespace Rcpp;
  // Function f = as<Function>(_f);
  // double a = as<double>(_a);
  // double b = as<double>(_b);
  // double epsrel = as<double>(_epsrel);
  // double epsabs = as<double>(_epsabs);
  // int limit = as<int>(_limit);
  // int ny = as<int>(_ny);
  double result[ny], abserr[ny];
  int lenw, ier, neval, last;
  lenw = 2*limit*ny + 2*limit;
  int iwork[limit];
  double work[lenw];
  vRdqags(f, a, b, epsabs, epsrel, ny,
	  result, abserr, &neval, &ier, limit, &lenw, &last, iwork, work);
  NumericVector resultnv(result, result+ny);
  NumericVector abserrnv(abserr, abserr+ny);
  return wrap(List::create(_("result") = resultnv,
			   _("abserr") = abserrnv,
			   _("subdivisions") = last,
			   _("ierr") = ier
			   ));
}

// void vRdqagi(Rcpp::Function f, const double bound, const int inf,
// 	     const double epsabs, const double epsrel, const int limit, const int ny,
// 	    double *result, double *abserr, int *neval, int *ier,
// 	    int *lenw, int *last,
// 	    int *iwork, double *work)

// [[Rcpp::export]]
Rcpp::List vdqagiRcpp(const Rcpp::Function f, const double bound, const int inf,
		      const double epsrel, const double epsabs, const int limit, const int ny) {
// RcppExport SEXP vdqagiRcpp(SEXP _f, SEXP _bound, SEXP _inf,
// 			   SEXP _epsrel, SEXP _epsabs, SEXP _limit, SEXP _ny) {
  using namespace Rcpp;
  // Function f = as<Function>(_f);
  // double bound = as<double>(_bound);
  // int inf = as<int>(_inf);
  // double epsrel = as<double>(_epsrel);
  // double epsabs = as<double>(_epsabs);
  // int limit = as<int>(_limit);
  // int ny = as<int>(_ny);
  double result[ny], abserr[ny];
  int lenw, ier, neval, last;
  lenw = 2*limit*ny + 2*limit;
  int iwork[limit];
  double work[lenw];
  vRdqagi(f, bound, inf, epsabs, epsrel, limit, ny,
	  result, abserr, &neval, &ier, &lenw, &last, iwork, work);
  NumericVector resultnv(result, result+ny);
  NumericVector abserrnv(abserr, abserr+ny);
  return wrap(List::create(_("result") = resultnv,
			   _("abserr") = abserrnv,
			   _("subdivisions") = last,
			   _("ierr") = ier
			   ));
}

// [[Rcpp::export]]
Rcpp::List vrdqk21Rcpp(const Rcpp::Function f, const double a, const double b) {
// RcppExport SEXP vrdqk21Rcpp(SEXP _f, SEXP _a, SEXP _b) {
  using namespace Rcpp;
  // Function f = as<Function>(_f);
  // double a = as<double>(_a);
  // double b = as<double>(_b);
  vec result, abserr, resabs, resasc;
  vrdqk21(f, a, b, &result, &abserr, &resabs, &resasc);
  return wrap(List::create(_("result") = result,
			   _("abserr") = abserr,
			   _("resabs") = resabs,
			   _("resasc") = resasc));
}

// [[Rcpp::export]]
Rcpp::List vrdqk15Rcpp(const Rcpp::Function f, const double boun, const int inf,
		       double a, double b) {
// RcppExport SEXP vrdqk15Rcpp(SEXP _f, SEXP _boun, SEXP _inf, SEXP _a, SEXP _b) {
  using namespace Rcpp;
  // Function f = as<Function>(_f);
  // double boun = as<double>(_boun);
  // int inf = as<int>(_inf);
  // double a = as<double>(_a);
  // double b = as<double>(_b);
  vec result, abserr, resabs, resasc;
  vrdqk15i(f, boun, inf, &a, &b, &result, &abserr, &resabs, &resasc);
  return wrap(List::create(_("result") = result,
			   _("abserr") = abserr,
			   _("resabs") = resabs,
			   _("resasc") = resasc));
}