integrate.cpp

Go to the documentation of this file.

/*
 *  R : A Computer Language for Statistical Data Analysis
 *  Copyright (C) 2001-2014  The R Core Team
 *
 *  This program is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation; either version 2 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program; if not, a copy is available at
 *  https://www.R-project.org/Licenses/
 *
 *  Most of this file is C translations of Fortran routines in
 *  QUADPACK: the latter is part of SLATEC 'and therefore in the public
 *  domain' (https://en.wikipedia.org/wiki/QUADPACK).
 *
 *
 */

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#ifdef _MSC_VER
  #ifdef small
    #undef small
  #endif
#endif

#include <math.h>
#include <float.h>
//#include <Rmath.h> /* for fmax2, fmin2, imin2 */
//#include <R_ext/Applic.h> /* exporting the API , particularly */
/*--- typedef void integr_fn(double *x, int n, void *ex) ---
 * vectorizing function   f(x[1:n], ...) -> x[]  {overwriting x[]}.
 * Vectorization can be used to speed up the integrand
 * instead of calling it  n  times.
*/

/* f2c-ed translations + modifications of QUADPACK functions */

template<class Float, class integr_fn> static void rdqagie(integr_fn f, void *ex,
        Float *, int *, Float * , Float *, int *,
        Float *, Float *, int *,
        int *, Float *, Float *, Float *, Float *,
        int *, int *);

template<class Float, class integr_fn> static void rdqk15i(integr_fn f, void *ex,
        Float *, int *, Float * , Float *,
        Float *, Float *, Float *, Float *);

template<class Float, class integr_fn> static void rdqagse(integr_fn f, void *ex, Float *, Float *,
        Float *, Float *, int *, Float *, Float *,
        int *, int *, Float *, Float *, Float *,
        Float *, int *, int *);

template<class Float, class integr_fn> static void rdqk21(integr_fn f, void *ex,
       Float *, Float *, Float *, Float *, Float *, Float *);

template<class Float> static void rdqpsrt(int *, int *, int *, Float *, Float *, int *, int *);

template<class Float> static void rdqelg(int *, Float *, Float *, Float *, Float *, int *);

template<class Float, class integr_fn>
void Rdqagi(integr_fn f, void *ex, Float *bound, int *inf,
      Float *epsabs, Float *epsrel,
      Float *result, Float *abserr, int *neval, int *ier,
      int *limit, int *lenw, int *last,
      int *iwork, Float *work)
{
    int l1, l2, l3;

/*
***begin prologue  dqagi
***date written   800101   (yymmdd)
***revision date  830518   (yymmdd)
***category no.  h2a3a1,h2a4a1
***keywords  automatic integrator, infinite intervals,
            general-purpose, transformation, extrapolation,
            globally adaptive
***author  piessens,robert,appl. math. & progr. div. - k.u.leuven
          de doncker,elise,appl. math. & progr. div. -k.u.leuven
***purpose  the routine calculates an approximation result to a given
           integral   i = integral of f over (bound,+infinity)
           or i = integral of f over (-infinity,bound)
           or i = integral of f over (-infinity,+infinity)
           hopefully satisfying following claim for accuracy
           abs(i-result) <= max(epsabs,epsrel*abs(i)).
***description

       integration over infinite intervals
       standard fortran subroutine

       parameters
        on entry
           f      - double precision
                    function subprogram defining the integrand
                    function f(x). the actual name for f needs to be
                    declared e x t e r n a l in the driver program.

           bound  - double precision
                    finite bound of integration range
                    (has no meaning if interval is doubly-infinite)

           inf    - int
                    indicating the kind of integration range involved
                    inf = 1 corresponds to  (bound,+infinity),
                    inf = -1            to  (-infinity,bound),
                    inf = 2             to (-infinity,+infinity).

           epsabs - double precision
                    absolute accuracy requested
           epsrel - double precision
                    relative accuracy requested
                    if  epsabs <= 0
                    and epsrel < max(50*rel.mach.acc.,0.5d-28),
                    the routine will end with ier = 6.


        on return
           result - double precision
                    approximation to the integral

           abserr - double precision
                    estimate of the modulus of the absolute error,
                    which should equal or exceed abs(i-result)

           neval  - int
                    number of integrand evaluations

           ier    - int
                    ier = 0 normal and reliable termination of the
                            routine. it is assumed that the requested
                            accuracy has been achieved.
                  - ier > 0 abnormal termination of the routine. the
                            estimates for result and error are less
                            reliable. it is assumed that the requested
                            accuracy has not been achieved.
           error messages
                    ier = 1 maximum number of subdivisions allowed
                            has been achieved. one can allow more
                            subdivisions by increasing the value of
                            limit (and taking the according dimension
                            adjustments into account). however, if
                            this yields no improvement it is advised
                            to analyze the integrand in order to
                            determine the integration difficulties. if
                            the position of a local difficulty can be
                            determined (e.g. singularity,
                            discontinuity within the interval) one
                            will probably gain from splitting up the
                            interval at this point and calling the
                            integrator on the subranges. if possible,
                            an appropriate special-purpose integrator
                            should be used, which is designed for
                            handling the type of difficulty involved.
                        = 2 the occurrence of roundoff error is
                            detected, which prevents the requested
                            tolerance from being achieved.
                            the error may be under-estimated.
                        = 3 extremely bad integrand behaviour occurs
                            at some points of the integration
                            interval.
                        = 4 the algorithm does not converge.
                            roundoff error is detected in the
                            extrapolation table.
                            it is assumed that the requested tolerance
                            cannot be achieved, and that the returned
                            result is the best which can be obtained.
                        = 5 the integral is probably divergent, or
                            slowly convergent. it must be noted that
                            divergence can occur with any other value
                            of ier.
                        = 6 the input is invalid, because
                            (epsabs <= 0 and
                             epsrel < max(50*rel.mach.acc.,0.5d-28))
                             or limit < 1 or leniw < limit*4.
                            result, abserr, neval, last are set to
                            zero. exept when limit or leniw is
                            invalid, iwork(1), work(limit*2+1) and
                            work(limit*3+1) are set to zero, work(1)
                            is set to a and work(limit+1) to b.

        dimensioning parameters
           limit - int
                   dimensioning parameter for iwork
                   limit determines the maximum number of subintervals
                   in the partition of the given integration interval
                   (a,b), limit >= 1.
                   if limit < 1, the routine will end with ier = 6.

           lenw  - int
                   dimensioning parameter for work
                   lenw must be at least limit*4.
                   if lenw < limit*4, the routine will end
                   with ier = 6.

           last  - int
                   on return, last equals the number of subintervals
                   produced in the subdivision process, which
                   determines the number of significant elements
                   actually in the work arrays.

        work arrays
           iwork - int
                   vector of dimension at least limit, the first
                   k elements of which contain pointers
                   to the error estimates over the subintervals,
                   such that work(limit*3+iwork(1)),... ,
                   work(limit*3+iwork(k)) form a decreasing
                   sequence, with k = last if last <= (limit/2+2), and
                   k = limit+1-last otherwise

           work  - double precision
                   vector of dimension at least lenw
                   on return
                   work(1), ..., work(last) contain the left
                    end points of the subintervals in the
                    partition of (a,b),
                   work(limit+1), ..., work(limit+last) contain
                    the right end points,
                   work(limit*2+1), ...,work(limit*2+last) contain the
                    integral approximations over the subintervals,
                   work(limit*3+1), ..., work(limit*3)
                    contain the error estimates.

***routines called  dqagie
***end prologue  dqagi */

    *ier = 6;
    *neval = 0;
    *last = 0;
    *result = 0.;
    *abserr = 0.;
    if (*limit < 1 || *lenw < *limit << 2) return;

    l1 = *limit;
    l2 = *limit + l1;
    l3 = *limit + l2;

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

    return;
} /* Rdqagi */

template<class Float, class integr_fn> static
void rdqagie(integr_fn f, void *ex, Float *bound, int *inf, Float *
       epsabs, Float *epsrel, int *limit, Float *result,
       Float *abserr, int *neval, int *ier, Float *alist,
       Float *blist, Float *rlist, Float *elist, int *
       iord, int *last)
{
    /* Local variables */
    Float area, dres;
    int ksgn;
    Float boun;
    int nres;
    Float area1, area2, area12;
    int k;
    Float small, erro12;
    small = 0.0;
    int ierro;
    Float a1, a2, b1, b2, defab1, defab2, oflow;
    int ktmin, nrmax;
    Float uflow;
    bool noext;
    int iroff1, iroff2, iroff3;
    Float res3la[3], error1, error2;
    int id;
    Float rlist2[52];
    int numrl2;
    Float defabs, epmach, erlarg, abseps, correc, errbnd, resabs;
    erlarg = 0.0;
    correc = 0.0;
    int jupbnd;
    Float erlast, errmax;
    int maxerr;
    Float reseps;
    bool extrap;
    Float ertest, errsum;
    ertest = 0.0;

/* ***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;
    *result = 0.;
    *abserr = 0.;
    alist[1] = 0.;
    blist[1] = 1.;
    rlist[1] = 0.;
    elist[1] = 0.;
    iord[1] = 0;
    if (*epsabs <= 0. && (*epsrel < 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.;
    }

    /* Table of constant values */
    static Float c_b6 = 0.;
    static Float c_b7 = 1.;

    rdqk15i(f, ex, &boun, inf, &c_b6, &c_b7, result, abserr, &defabs, &resabs);

/*           test on accuracy */

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

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

    uflow = DBL_MIN;
    oflow = DBL_MAX;
    rlist2[0] = *result;
    errmax = *abserr;
    maxerr = 1;
    area = *result;
    errsum = *abserr;
    *abserr = oflow;
    nrmax = 1;
    nres = 0;
    ktmin = 0;
    numrl2 = 2;
    extrap = FALSE;
    noext = FALSE;
    ierro = 0;
    iroff1 = 0;
    iroff2 = 0;
    iroff3 = 0;
    ksgn = -1;
    if (dres >= (1. - epmach * 50.) * defabs) {
  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;
  rdqk15i(f, ex, &boun, inf, &a1, &b1, &area1, &error1, &resabs, &defab1);
  rdqk15i(f, ex, &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[maxerr];
  if (!(defab1 == error1 || defab2 == error2)) {
      if (fabs(rlist[maxerr] - area12) <= fabs(area12) * 1e-5 &&
    erro12 >= errmax * .99) {
    if (extrap)
        ++iroff2;
    else /* if (! extrap) */
        ++iroff1;
      }
      if (*last > 10 && erro12 > errmax)
    ++iroff3;
  }

  rlist[maxerr] = area1;
  rlist[*last] = area2;
  errbnd = fmax2(*epsabs, *epsrel * fabs(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 (fmax2(fabs(a1), fabs(b2)) <=
      (epmach * 100. + 1.) * (fabs(a2) + uflow * 1e3))
    {
      *ier = 4;
  }

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

  if (error2 <= error1) {
      alist[*last] = a2;
      blist[maxerr] = b1;
      blist[*last] = b2;
      elist[maxerr] = error1;
      elist[*last] = error2;
  }
  else {
      alist[maxerr] = a2;
      alist[*last] = a1;
      blist[*last] = b1;
      rlist[maxerr] = area2;
      rlist[*last] = area1;
      elist[maxerr] = error2;
      elist[*last] = 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). */

  rdqpsrt(limit, last, &maxerr, &errmax, &elist[1], &iord[1], &nrmax);
  if (errsum <= errbnd) {
      goto L115;
  }
  if (*ier != 0)      break;
  if (*last == 2) { /* L80: */
      small = .375;
      erlarg = errsum;
      ertest = errbnd;
      rlist2[1] = area; 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 && 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[maxerr];
    if (fabs(blist[maxerr] - alist[maxerr]) > small) {
        goto L90;
    }
    ++nrmax;
    /* L50: */
      }
  }
/*           perform extrapolation.  L60: */
  ++numrl2;
  rlist2[numrl2 - 1] = area;
  rdqelg(&numrl2, rlist2, &reseps, &abseps, res3la, &nres);
  ++ktmin;
  if (ktmin > 5 && *abserr < errsum * .001) {
      *ier = 5;
  }
  if (abseps >= *abserr) {
      goto L70;
  }
  ktmin = 0;
  *abserr = abseps;
  *result = reseps;
  correc = erlarg;
  ertest = fmax2(*epsabs, *epsrel * fabs(reseps));
  if (*abserr <= ertest) {
      break;
  }

/*            prepare bisection of the smallest interval. */

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

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

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

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

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

/*           compute global integral sum. */

L115:
    *result = 0.;
    for (k = 1; k <= *last; ++k)
  *result += rlist[k];

    *abserr = errsum;
L130:
    *neval = *last * 30 - 15;
    if (*inf == 2) {
  *neval <<= 1;
    }
    if (*ier > 2) {
  --(*ier);
    }
    return;
} /* rdqagie_ */

template<class Float, class integr_fn>
void Rdqags(integr_fn f, void *ex, Float *a, Float *b,
      Float *epsabs, Float *epsrel,
      Float *result, Float *abserr, int *neval, int *ier,
      int *limit, int *lenw, int *last, int *iwork, Float *work)
{
    int l1, l2, l3;

/*
***begin prologue  dqags
***date written   800101   (yymmdd)
***revision date  830518   (yymmdd)
***category no.  h2a1a1
***keywords  automatic integrator, general-purpose,
            (end-point) singularities, extrapolation,
            globally adaptive
***author  piessens,robert,appl. math. & progr. div. - k.u.leuven
          de doncker,elise,appl. math. & prog. div. - k.u.leuven
***purpose  the routine calculates an approximation result to a given
           definite integral  i = integral of f over (a,b),
           hopefully satisfying following claim for accuracy
           abs(i-result) <= max(epsabs,epsrel*abs(i)).
***description

       computation of a definite integral
       standard fortran subroutine
       double precision version


       parameters
        on entry
           f      - double precision
                    function subprogram defining the integrand
                    function f(x). the actual name for f needs to be
                    declared e x t e r n a l in the driver program.

           a      - double precision
                    lower limit of integration

           b      - double precision
                    upper limit of integration

           epsabs - double precision
                    absolute accuracy requested
           epsrel - double precision
                    relative accuracy requested
                    if  epsabs <= 0
                    and epsrel < max(50*rel.mach.acc.,0.5d-28),
                    the routine will end with ier = 6.

        on return
           result - double precision
                    approximation to the integral

           abserr - double precision
                    estimate of the modulus of the absolute error,
                    which should equal or exceed abs(i-result)

           neval  - int
                    number of integrand evaluations

           ier    - int
                    ier = 0 normal and reliable termination of the
                            routine. it is assumed that the requested
                            accuracy has been achieved.
                    ier > 0 abnormal termination of the routine
                            the estimates for integral and error are
                            less reliable. it is assumed that the
                            requested accuracy has not been achieved.
           error messages
                    ier = 1 maximum number of subdivisions allowed
                            has been achieved. one can allow more sub-
                            divisions by increasing the value of limit
                            (and taking the according dimension
                            adjustments into account. however, if
                            this yields no improvement it is advised
                            to analyze the integrand in order to
                            determine the integration difficulties. if
                            the position of a local difficulty can be
                            determined (e.g. singularity,
                            discontinuity within the interval) one
                            will probably gain from splitting up the
                            interval at this point and calling the
                            integrator on the subranges. if possible,
                            an appropriate special-purpose integrator
                            should be used, which is designed for
                            handling the type of difficulty involved.
                        = 2 the occurrence of roundoff error is detec-
                            ted, which prevents the requested
                            tolerance from being achieved.
                            the error may be under-estimated.
                        = 3 extremely bad integrand behaviour
                            occurs at some points of the integration
                            interval.
                        = 4 the algorithm does not converge.
                            roundoff error is detected in the
                            extrapolation table. it is presumed that
                            the requested tolerance cannot be
                            achieved, and that the returned result is
                            the best which can be obtained.
                        = 5 the integral is probably divergent, or
                            slowly convergent. it must be noted that
                            divergence can occur with any other value
                            of ier.
                        = 6 the input is invalid, because
                            (epsabs <= 0 and
                             epsrel < max(50*rel.mach.acc.,0.5d-28)
                            or limit < 1 or lenw < limit*4.
                            result, abserr, neval, last are set to
                            zero.except when limit or lenw is invalid,
                            iwork(1), work(limit*2+1) and
                            work(limit*3+1) are set to zero, work(1)
                            is set to a and work(limit+1) to b.

        dimensioning parameters
           limit - int
                   dimensioning parameter for iwork
                   limit determines the maximum number of subintervals
                   in the partition of the given integration interval
                   (a,b), limit >= 1.
                   if limit < 1, the routine will end with ier = 6.

           lenw  - int
                   dimensioning parameter for work
                   lenw must be at least limit*4.
                   if lenw < limit*4, the routine will end
                   with ier = 6.

           last  - int
                   on return, last equals the number of subintervals
                   produced in the subdivision process, detemines the
                   number of significant elements actually in the work
                   arrays.

        work arrays
           iwork - int
                   vector of dimension at least limit, the first k
                   elements of which contain pointers
                   to the error estimates over the subintervals
                   such that work(limit*3+iwork(1)),... ,
                   work(limit*3+iwork(k)) form a decreasing
                   sequence, with k = last if last <= (limit/2+2),
                   and k = limit+1-last otherwise

           work  - double precision
                   vector of dimension at least lenw
                   on return
                   work(1), ..., work(last) contain the left
                    end-points of the subintervals in the
                    partition of (a,b),
                   work(limit+1), ..., work(limit+last) contain
                    the right end-points,
                   work(limit*2+1), ..., work(limit*2+last) contain
                    the integral approximations over the subintervals,
                   work(limit*3+1), ..., work(limit*3+last)
                    contain the error estimates.

***routines called  dqagse
***end prologue  dqags */

/*         check validity of limit and lenw. */

    *ier = 6;
    *neval = 0;
    *last = 0;
    *result = 0.;
    *abserr = 0.;
    if (*limit < 1 || *lenw < *limit *4) return;

/*         prepare call for dqagse. */

    l1 = *limit;
    l2 = *limit + l1;
    l3 = *limit + l2;

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

    return;
} /* rdqags_ */

template<class Float, class integr_fn> static
void rdqagse(integr_fn f, void *ex, Float *a, Float *b, Float *
       epsabs, Float *epsrel, int *limit, Float *result,
       Float *abserr, int *neval, int *ier, Float *alist,
       Float *blist, Float *rlist, Float *elist, int *
       iord, int *last)
{
    /* Local variables */
    bool noext, extrap;
    int k,ksgn, nres;
    int ierro;
    int ktmin, nrmax;
    int iroff1, iroff2, iroff3;
    int id;
    int numrl2;
    int jupbnd;
    int maxerr;
    Float res3la[3];
    Float rlist2[52];
    Float abseps, area, area1, area2, area12, dres, epmach;
    Float a1, a2, b1, b2, defabs, defab1, defab2, oflow, uflow, resabs, reseps;
    Float error1, error2, erro12, errbnd, erlast, errmax, errsum;

    Float correc = 0.0, erlarg = 0.0, ertest = 0.0, small = 0.0;
/*
***begin prologue  dqagse
***date written   800101   (yymmdd)
***revision date  830518   (yymmdd)
***category no.  h2a1a1
***keywords  automatic integrator, general-purpose,
            (end point) singularities, extrapolation,
            globally adaptive
***author  piessens,robert,appl. math. & progr. div. - k.u.leuven
          de doncker,elise,appl. math. & progr. div. - k.u.leuven
***purpose  the routine calculates an approximation result to a given
           definite integral i = integral of f over (a,b),
           hopefully satisfying following claim for accuracy
           abs(i-result) <= max(epsabs,epsrel*abs(i)).
***description

       computation of a definite integral
       standard fortran subroutine
       double precision version

       parameters
        on entry
           f      - double precision
                    function subprogram defining the integrand
                    function f(x). the actual name for f needs to be
                    declared e x t e r n a l in the driver program.

           a      - double precision
                    lower limit of integration

           b      - double precision
                    upper limit of integration

           epsabs - double precision
                    absolute accuracy requested
           epsrel - double precision
                    relative accuracy requested
                    if  epsabs <= 0
                    and epsrel < max(50*rel.mach.acc.,0.5d-28),
                    the routine will end with ier = 6.

           limit  - int
                    gives an upperbound on the number of subintervals
                    in the partition of (a,b)

        on return
           result - double precision
                    approximation to the integral

           abserr - double precision
                    estimate of the modulus of the absolute error,
                    which should equal or exceed abs(i-result)

           neval  - int
                    number of integrand evaluations

           ier    - int
                    ier = 0 normal and reliable termination of the
                            routine. it is assumed that the requested
                            accuracy has been achieved.
                    ier > 0 abnormal termination of the routine
                            the estimates for integral and error are
                            less reliable. it is assumed that the
                            requested accuracy has not been achieved.
           error messages
                        = 1 maximum number of subdivisions allowed
                            has been achieved. one can allow more sub-
                            divisions by increasing the value of limit
                            (and taking the according dimension
                            adjustments into account). however, if
                            this yields no improvement it is advised
                            to analyze the integrand in order to
                            determine the integration difficulties. if
                            the position of a local difficulty can be
                            determined (e.g. singularity,
                            discontinuity within the interval) one
                            will probably gain from splitting up the
                            interval at this point and calling the
                            integrator on the subranges. if possible,
                            an appropriate special-purpose integrator
                            should be used, which is designed for
                            handling the type of difficulty involved.
                        = 2 the occurrence of roundoff error is detec-
                            ted, which prevents the requested
                            tolerance from being achieved.
                            the error may be under-estimated.
                        = 3 extremely bad integrand behaviour
                            occurs at some points of the integration
                            interval.
                        = 4 the algorithm does not converge.
                            roundoff error is detected in the
                            extrapolation table.
                            it is presumed that the requested
                            tolerance cannot be achieved, and that the
                            returned result is the best which can be
                            obtained.
                        = 5 the integral is probably divergent, or
                            slowly convergent. it must be noted that
                            divergence can occur with any other value
                            of ier.
                        = 6 the input is invalid, because
                            epsabs <= 0 and
                            epsrel < max(50*rel.mach.acc.,0.5d-28).
                            result, abserr, neval, last, rlist(1),
                            iord(1) and elist(1) are set to zero.
                            alist(1) and blist(1) are set to a and b
                            respectively.

           alist  - double precision
                    vector of dimension at least limit, the first
                     last  elements of which are the left end points
                    of the subintervals in the partition of the
                    given integration range (a,b)

           blist  - double precision
                    vector of dimension at least limit, the first
                     last  elements of which are the right end points
                    of the subintervals in the partition of the given
                    integration range (a,b)

           rlist  - double precision
                    vector of dimension at least limit, the first
                     last  elements of which are the integral
                    approximations on the subintervals

           elist  - double precision
                    vector of dimension at least limit, the first
                     last  elements of which are the moduli of the
                    absolute error estimates on the subintervals

           iord   - int
                    vector of dimension at least limit, the first k
                    elements of which are pointers to the
                    error estimates over the subintervals,
                    such that elist(iord(1)), ..., elist(iord(k))
                    form a decreasing sequence, with k = last
                    if last <= (limit/2+2), and k = limit+1-last
                    otherwise

           last   - int
                    number of subintervals actually produced in the
                    subdivision process

***references  (none)
***routines called  dqelg,dqk21,dqpsrt
***end prologue  dqagse



           the dimension of rlist2 is determined by the value of
           limexp in subroutine dqelg (rlist2 should be of dimension
           (limexp+2) at least).

           list of major variables
           -----------------------

          alist     - list of left end points of all subintervals
                      considered up to now
          blist     - list of right end points of all subintervals
                      considered up to now
          rlist(i)  - approximation to the integral over
                      (alist(i),blist(i))
          rlist2    - array of dimension at least limexp+2 containing
                      the part of the epsilon table which is still
                      needed for further computations
          elist(i)  - error estimate applying to rlist(i)
          maxerr    - pointer to the interval with largest error
                      estimate
          errmax    - elist(maxerr)
          erlast    - error on the interval currently subdivided
                      (before that subdivision has taken place)
          area      - sum of the integrals over the subintervals
          errsum    - sum of the errors over the subintervals
          errbnd    - requested accuracy max(epsabs,epsrel*
                      abs(result))
          *****1    - variable for the left interval
          *****2    - variable for the right interval
          last      - index for subdivision
          nres      - number of calls to the extrapolation routine
          numrl2    - number of elements currently in rlist2. if an
                      appropriate approximation to the compounded
                      integral has been obtained it is put in
                      rlist2(numrl2) after numrl2 has been increased
                      by one.
          small     - length of the smallest interval considered up
                      to now, multiplied by 1.5
          erlarg    - sum of the errors over the intervals larger
                      than the smallest interval considered up to now
          extrap    - logical variable denoting that the routine is
                      attempting to perform extrapolation i.e. before
                      subdividing the smallest interval we try to
                      decrease the value of erlarg.
          noext     - logical variable denoting that extrapolation
                      is no longer allowed (true value)

           machine dependent constants
           ---------------------------

          epmach is the largest relative spacing.
          uflow is the smallest positive magnitude.
          oflow is the largest positive magnitude. */

/* ***first executable statement  dqagse */
    /* Parameter adjustments */
    --iord;
    --elist;
    --rlist;
    --blist;
    --alist;

    /* Function Body */
    epmach = DBL_EPSILON;

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

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

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

/*           test on accuracy. */

    dres = fabs(*result);
    errbnd = fmax2(*epsabs, *epsrel * dres);
    *last = 1;
    rlist[1] = *result;
    elist[1] = *abserr;
    iord[1] = 1;
    if (*abserr <= epmach * 100. * defabs && *abserr > errbnd)
  *ier = 2;
    if (*limit == 1)
  *ier = 1;
    if (*ier != 0 || (*abserr <= errbnd && *abserr != resabs)
  || *abserr == 0.) goto L140;

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

    rlist2[0] = *result;
    errmax = *abserr;
    maxerr = 1;
    area = *result;
    errsum = *abserr;
    *abserr = oflow;
    nrmax = 1;
    nres = 0;
    numrl2 = 2;
    ktmin = 0;
    extrap = FALSE;
    noext = FALSE;
    iroff1 = 0;
    iroff2 = 0;
    iroff3 = 0;
    ksgn = -1;
    if (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;
  rdqk21(f, ex, &a1, &b1, &area1, &error1, &resabs, &defab1);
  rdqk21(f, ex, &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[maxerr];
  if (!(defab1 == error1 || defab2 == error2)) {

      if (fabs(rlist[maxerr] - area12) <= fabs(area12) * 1e-5 &&
    erro12 >= errmax * .99) {
    if (extrap)
        ++iroff2;
    else /* if(! extrap) */
        ++iroff1;
      }
      if (*last > 10 && erro12 > errmax)
    ++iroff3;
  }
  rlist[maxerr] = area1;
  rlist[*last] = area2;
  errbnd = fmax2(*epsabs, *epsrel * fabs(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 (fmax2(fabs(a1), fabs(b2)) <=
      (epmach * 100. + 1.) * (fabs(a2) + uflow * 1e3)) {
      *ier = 4;
  }

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

  if (error2 > error1) {
      alist[maxerr] = a2;
      alist[*last] = a1;
      blist[*last] = b1;
      rlist[maxerr] = area2;
      rlist[*last] = area1;
      elist[maxerr] = error2;
      elist[*last] = error1;
  } else {
      alist[*last] = a2;
      blist[maxerr] = b1;
      blist[*last] = b2;
      elist[maxerr] = error1;
      elist[*last] = 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:*/
  rdqpsrt(limit, last, &maxerr, &errmax, &elist[1], &iord[1], &nrmax);

  if (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;
      rlist2[1] = area; 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 && 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[maxerr];
    if (fabs(blist[maxerr] - alist[maxerr]) > small) {
        goto L90;
    }
    ++nrmax;
    /* L50: */
      }
  }
/*           perform extrapolation.  L60: */

  ++numrl2;
  rlist2[numrl2 - 1] = area;
  rdqelg(&numrl2, rlist2, &reseps, &abseps, res3la, &nres);
  ++ktmin;
  if (ktmin > 5 && *abserr < errsum * .001) {
      *ier = 5;
  }
  if (abseps < *abserr) {
      ktmin = 0;
      *abserr = abseps;
      *result = reseps;
      correc = erlarg;
      ertest = fmax2(*epsabs, *epsrel * fabs(reseps));
      if (*abserr <= ertest) {
    break;
      }
  }

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

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


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

    if (*abserr == oflow)   goto L115;
    if (*ier + ierro == 0)  goto L110;
    if (ierro == 3)
  *abserr += correc;
    if (*ier == 0)
  *ier = 3;
    if (*result == 0. || area == 0.) {
  if (*abserr > errsum)   goto L115;
  if (area == 0.)   goto L130;
    }
    else { /* L105:*/
  if (*abserr / fabs(*result) > errsum / fabs(area))
      goto L115;
    }

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

L115:/*   compute global integral sum. */
    *result = 0.;
    for (k = 1; k <= *last; ++k)
  *result += rlist[k];
    *abserr = errsum;
L130:
    if (*ier > 2)
L140:
    *neval = *last * 42 - 21;
    return;
} /* rdqagse_ */


template<class Float, class integr_fn> static void rdqk15i(integr_fn f, void *ex,
        Float *boun, int *inf, Float *a, Float *b,
        Float *result,
        Float *abserr, Float *resabs, Float *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 */
    Float absc, dinf, resg, resk, fsum, absc1, absc2, fval1, fval2;
    int j;
    Float hlgth, centr, reskh, uflow;
    Float tabsc1, tabsc2, fc, epmach;
    Float fv1[7], fv2[7], vec[15], vec2[15];

/*
***begin prologue  dqk15i
***date written   800101   (yymmdd)
***revision date  830518   (yymmdd)
***category no.  h2a3a2,h2a4a2
***keywords  15-point transformed gauss-kronrod rules
***author  piessens,robert,appl. math. & progr. div. - k.u.leuven
          de doncker,elise,appl. math. & progr. div. - k.u.leuven
***purpose  the original (infinite integration range is mapped
           onto the interval (0,1) and (a,b) is a part of (0,1).
           it is the purpose to compute
           i = integral of transformed integrand over (a,b),
           j = integral of abs(transformed integrand) over (a,b).
***description

          integration rule
          standard fortran subroutine
          double precision version

          parameters
           on entry
             f      - double precision
                      fuction subprogram defining the integrand
                      function f(x). the actual name for f needs to be
                      declared e x t e r n a l in the calling program.

             boun   - double precision
                      finite bound of original integration
                      range (set to zero if inf = +2)

             inf    - int
                      if inf = -1, the original interval is
                                  (-infinity,bound),
                      if inf = +1, the original interval is
                                  (bound,+infinity),
                      if inf = +2, the original interval is
                                  (-infinity,+infinity) and
                      the integral is computed as the sum of two
                      integrals, one over (-infinity,0) and one over
                      (0,+infinity).

             a      - double precision
                      lower limit for integration over subrange
                      of (0,1)

             b      - double precision
                      upper limit for integration over subrange
                      of (0,1)

           on return
             result - double precision
                      approximation to the integral i
                      result is computed by applying the 15-point
                      kronrod rule(resk) obtained by optimal addition
                      of abscissae to the 7-point gauss rule(resg).

             abserr - double precision
                      estimate of the modulus of the absolute error,
                      which should equal or exceed abs(i-result)

             resabs - double precision
                      approximation to the integral j

             resasc - double precision
                      approximation to the integral of
                      abs((transformed integrand)-i/(b-a)) over (a,b)

***references  (none)
***end prologue  dqk15i


          the abscissae and weights are supplied for the interval
          (-1,1).  because of symmetry only the positive abscissae and
          their corresponding weights are given.

          xgk    - abscissae of the 15-point kronrod rule
                   xgk(2), xgk(4), ... abscissae of the 7-point
                   gauss rule
                   xgk(1), xgk(3), ...  abscissae which are optimally
                   added to the 7-point gauss rule

          wgk    - weights of the 15-point kronrod rule

          wg     - weights of the 7-point gauss rule, corresponding
                   to the abscissae xgk(2), xgk(4), ...
                   wg(1), wg(3), ... are set to zero.





          list of major variables
          -----------------------

          centr  - mid point of the interval
          hlgth  - half-length of the interval
          absc*  - abscissa
          tabsc* - transformed abscissa
          fval*  - function value
          resg   - result of the 7-point gauss formula
          resk   - result of the 15-point kronrod formula
          reskh  - approximation to the mean value of the transformed
                   integrand over (a,b), i.e. to i/(b-a)

          machine dependent constants
          ---------------------------

          epmach is the largest relative spacing.
          uflow is the smallest positive magnitude.
*/

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

    centr = (*a + *b) * .5;
    hlgth = (*b - *a) * .5;
    tabsc1 = *boun + dinf * (1. - centr) / centr;
    vec[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;
  vec[(j << 1) - 1] = tabsc1;
  vec[j * 2] = tabsc2;
  if (*inf == 2) {
      vec2[(j << 1) - 1] = -tabsc1;
      vec2[j * 2] = -tabsc2;
  }
/* L5: */
    }
    f(vec, 15, ex); /* -> new vec[] overwriting old vec[] */
    if (*inf == 2) f(vec2, 15, ex);
    fval1 = vec[0];
    if (*inf == 2) fval1 += vec2[0];
    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 = fabs(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 = vec[(j << 1) - 1];
  fval2 = vec[j * 2];
  if (*inf == 2) {
      fval1 += vec2[(j << 1) - 1];
  }
  if (*inf == 2) {
      fval2 += 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] * (fabs(fval1) + fabs(fval2));
/* L10: */
    }
    reskh = resk * .5;
    *resasc = wgk[7] * fabs(fc - reskh);
    for (j = 1; j <= 7; ++j) {
  *resasc += wgk[j - 1] * (fabs(fv1[j - 1] - reskh) +
         fabs(fv2[j - 1] - reskh));
/* L20: */
    }
    *result = resk * hlgth;
    *resasc *= hlgth;
    *resabs *= hlgth;
    *abserr = fabs((resk - resg) * hlgth);
    if (*resasc != 0. && *abserr != 0.) {
  *abserr = *resasc * fmin2(1., pow(*abserr * 200. / *resasc, 1.5));
    }
    if (*resabs > uflow / (epmach * 50.)) {
  *abserr = fmax2(epmach * 50. * *resabs, *abserr);
    }
    return;
} /* rdqk15i_ */

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

/* ***begin prologue  dqelg
***refer to  dqagie,dqagoe,dqagpe,dqagse
***revision date  830518   (yymmdd)
***keywords  epsilon algorithm, convergence acceleration,
            extrapolation
***author  piessens,robert,appl. math. & progr. div. - k.u.leuven
          de doncker,elise,appl. math & progr. div. - k.u.leuven
***purpose  the routine determines the limit of a given sequence of
           approximations, by means of the epsilon algorithm of
           p.wynn. an estimate of the absolute error is also given.
           the condensed epsilon table is computed. only those
           elements needed for the computation of the next diagonal
           are preserved.
***description

          epsilon algorithm
          standard fortran subroutine
          double precision version

          parameters
             n      - int
                      epstab(n) contains the new element in the
                      first column of the epsilon table.

             epstab - double precision
                      vector of dimension 52 containing the elements
                      of the two lower diagonals of the triangular
                      epsilon table. the elements are numbered
                      starting at the right-hand corner of the
                      triangle.

             result - double precision
                      resulting approximation to the integral

             abserr - double precision
                      estimate of the absolute error computed from
                      result and the 3 previous results

             res3la - double precision
                      vector of dimension 3 containing the last 3
                      results

             nres   - int
                      number of calls to the routine
                      (should be zero at first call)

***end prologue  dqelg


          list of major variables
          -----------------------

          e0     - the 4 elements on which the computation of a new
          e1       element in the epsilon table is based
          e2
          e3                 e0
                       e3    e1    new
                             e2

          newelm - number of elements to be computed in the new diagonal
          errA   - errA = abs(e1-e0)+abs(e2-e1)+abs(new-e2)
          result - the element in the new diagonal with least value of errA

          machine dependent constants
          ---------------------------

          epmach is the largest relative spacing.
          oflow is the largest positive magnitude.
          limexp is the maximum number of elements the epsilon
          table can contain. if this number is reached, the upper
          diagonal of the epsilon table is deleted. */

/* ***first executable statement  dqelg */
    /* 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 = fmax2(fabs(e2), e1abs) * epmach;
  delta3 = e1 - e0;
  err3 = fabs(delta3);
  tol3 = 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 = 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 = fmax2(*abserr, epmach * 5. * fabs(*result));
    return;
} /* rdqelg_ */

template<class Float, class integr_fn> static void  rdqk21(integr_fn f, void *ex, Float *a, Float *b, Float *result,
        Float *abserr, Float *resabs, Float *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 */
    Float fv1[10], fv2[10], vec[21];
    Float absc, resg, resk, fsum, fval1, fval2;
    Float hlgth, centr, reskh, uflow;
    Float fc, epmach, dhlgth;
    int j, jtw, jtwm1;

/* ***begin prologue  dqk21
***date written   800101   (yymmdd)
***revision date  830518   (yymmdd)
***category no.  h2a1a2
***keywords  21-point gauss-kronrod rules
***author  piessens,robert,appl. math. & progr. div. - k.u.leuven
          de doncker,elise,appl. math. & progr. div. - k.u.leuven
***purpose  to compute i = integral of f over (a,b), with error
                          estimate
                      j = integral of abs(f) over (a,b)
***description

          integration rules
          standard fortran subroutine
          double precision version

          parameters
           on entry
             f      - double precision
                      function subprogram defining the integrand
                      function f(x). the actual name for f needs to be
                      declared e x t e r n a l in the driver program.

             a      - double precision
                      lower limit of integration

             b      - double precision
                      upper limit of integration

           on return
             result - double precision
                      approximation to the integral i
                      result is computed by applying the 21-point
                      kronrod rule (resk) obtained by optimal addition
                      of abscissae to the 10-point gauss rule (resg).

             abserr - double precision
                      estimate of the modulus of the absolute error,
                      which should not exceed abs(i-result)

             resabs - double precision
                      approximation to the integral j

             resasc - double precision
                      approximation to the integral of abs(f-i/(b-a))
                      over (a,b)

***references  (none)
***end prologue  dqk21



          the abscissae and weights are given for the interval (-1,1).
          because of symmetry only the positive abscissae and their
          corresponding weights are given.

          xgk    - abscissae of the 21-point kronrod rule
                   xgk(2), xgk(4), ...  abscissae of the 10-point
                   gauss rule
                   xgk(1), xgk(3), ...  abscissae which are optimally
                   added to the 10-point gauss rule

          wgk    - weights of the 21-point kronrod rule

          wg     - weights of the 10-point gauss rule


gauss quadrature weights and kronron quadrature abscissae and weights
as evaluated with 80 decimal digit arithmetic by l. w. fullerton,
bell labs, nov. 1981.





          list of major variables
          -----------------------

          centr  - mid point of the interval
          hlgth  - half-length of the interval
          absc   - abscissa
          fval*  - function value
          resg   - result of the 10-point gauss formula
          resk   - result of the 21-point kronrod formula
          reskh  - approximation to the mean value of f over (a,b),
                   i.e. to i/(b-a)


          machine dependent constants
          ---------------------------

          epmach is the largest relative spacing.
          uflow is the smallest positive magnitude. */

/* ***first executable statement  dqk21 */
    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. */

    resg = 0.;
    vec[0] = centr;
    for (j = 1; j <= 5; ++j) {
  jtw = j << 1;
  absc = hlgth * xgk[jtw - 1];
  vec[(j << 1) - 1] = centr - absc;
/* L5: */
  vec[j * 2] = centr + absc;
    }
    for (j = 1; j <= 5; ++j) {
  jtwm1 = (j << 1) - 1;
  absc = hlgth * xgk[jtwm1 - 1];
  vec[(j << 1) + 9] = centr - absc;
  vec[(j << 1) + 10] = centr + absc;
    }
    f(vec, 21, ex);
    fc = vec[0];
    resk = wgk[10] * fc;
    *resabs = fabs(resk);
    for (j = 1; j <= 5; ++j) {
  jtw = j << 1;
  absc = hlgth * xgk[jtw - 1];
  fval1 = vec[(j << 1) - 1];
  fval2 = vec[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] * (fabs(fval1) + fabs(fval2));
/* L10: */
    }
    for (j = 1; j <= 5; ++j) {
  jtwm1 = (j << 1) - 1;
  absc = hlgth * xgk[jtwm1 - 1];
  fval1 = vec[(j << 1) + 9];
  fval2 = vec[(j << 1) + 10];
  fv1[jtwm1 - 1] = fval1;
  fv2[jtwm1 - 1] = fval2;
  fsum = fval1 + fval2;
  resk += wgk[jtwm1 - 1] * fsum;
  *resabs += wgk[jtwm1 - 1] * (fabs(fval1) + fabs(fval2));
/* L15: */
    }
    reskh = resk * .5;
    *resasc = wgk[10] * fabs(fc - reskh);
    for (j = 1; j <= 10; ++j) {
  *resasc += wgk[j - 1] * (fabs(fv1[j - 1] - reskh) +
         fabs(fv2[j - 1] - reskh));
/* L20: */
    }
    *result = resk * hlgth;
    *resabs *= dhlgth;
    *resasc *= dhlgth;
    *abserr = fabs((resk - resg) * hlgth);
    if (*resasc != 0. && *abserr != 0.) {
  *abserr = *resasc * fmin2(1., pow(*abserr * 200. / *resasc, 1.5));
    }
    if (*resabs > uflow / (epmach * 50.)) {
  *abserr = fmax2(epmach * 50. * *resabs, *abserr);
    }
    return;
} /* rdqk21_ */

template<class Float> static void rdqpsrt(int *limit, int *last, int *maxerr,
        Float *ermax, Float *elist, int *iord, int *nrmax)
{
    /* Local variables */
    int i, j, k, ido, jbnd, isucc, jupbn;
    Float errmin, errmax;

/* ***begin prologue  dqpsrt
 ***refer to  dqage,dqagie,dqagpe,dqawse
 ***routines called  (none)
 ***revision date  810101   (yymmdd)
 ***keywords  sequential sorting
 ***author  piessens,robert,appl. math. & progr. div. - k.u.leuven
           de doncker,elise,appl. math. & progr. div. - k.u.leuven
 ***purpose  this routine maintains the descending ordering in the
            list of the local error estimated resulting from the
            interval subdivision process. at each call two error
            estimates are inserted using the sequential search
            method, top-down for the largest error estimate and
            bottom-up for the smallest error estimate.
 ***description

           ordering routine
           standard fortran subroutine
           double precision version

           parameters (meaning at output)
              limit  - int
                       maximum number of error estimates the list
                       can contain

              last   - int
                       number of error estimates currently in the list

              maxerr - int
                       maxerr points to the nrmax-th largest error
                       estimate currently in the list

              ermax  - double precision
                       nrmax-th largest error estimate
                       ermax = elist(maxerr)

              elist  - double precision
                       vector of dimension last containing
                       the error estimates

              iord   - int
                       vector of dimension last, the first k elements
                       of which contain pointers to the error
                       estimates, such that
                       elist(iord(1)),...,  elist(iord(k))
                       form a decreasing sequence, with
                       k = last if last <= (limit/2+2), and
                       k = limit+1-last otherwise

              nrmax  - int
                       maxerr = iord(nrmax)

***end prologue  dqpsrt
*/


    /* 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_ */