EIC code displayed by LXR master/​include/​NumA/​integration/​one_dimension/​gauss_kronrod_adaptive/​GaussKronrodAdaptive.h	Source navigation Diff markup Identifier search General search
	Version: master test Revert   Change

File indexing completed on 2026-06-02 08:48:22

 #ifndef GAUSS_KRONROD_ADAPTIVE_H
 #define GAUSS_KRONROD_ADAPTIVE_H
 
 /**
  * @file GaussKronrodAdaptive.h
  * @author Bryan BERTHOU (SPhN / CEA Saclay)
  * @date 27 January 2016
  * @version 1.0
  */
 
 #include <cmath>
 #include <iostream>
 #include <limits>
 #include <vector>
 #include <stddef.h>
 #include <stdlib.h>
 
 #include "../../../functor/one_dimension/FunctionType1D.h"
 //#include "../../../utils/Errors.h"
 //#include "../../../utils/Tolerances.h"
 #include "../Integrator1D.h"
 
 #define GSL_DBL_EPSILON        2.2204460492503131e-16
 #define GSL_DBL_MIN        2.2250738585072014e-308
 #define GSL_DBL_MAX        1.7976931348623157e+308
 
 #define GSL_ERROR(reason, gsl_errno) \
        do { \
        gsl_error (reason, __FILE__, __LINE__, gsl_errno) ; \
        return gsl_errno ; \
        } while (0)
 
 #define GSL_ERROR_VAL(reason, gsl_errno, value) \
        do { \
        gsl_error (reason, __FILE__, __LINE__, gsl_errno) ; \
        return value ; \
        } while (0)
 
 namespace NumA {
 
 /**
  * @class GaussKronrodAdaptive
  *
  * @brief
  */
 
 class GaussKronrodAdaptive: public Integrator1D {
 public:
     void gsl_error(const char * reason, const char * file, int line,
             int gsl_errno) {
         std::cerr << "GSL_ERROR = " << gsl_errno << std::endl;
     }
 
     enum {
         GSL_SUCCESS = 0, GSL_FAILURE = -1, GSL_CONTINUE = -2, /* iteration has not converged */
         GSL_EDOM = 1, /* input domain error, e.g sqrt(-1) */
         GSL_ERANGE = 2, /* output range error, e.g. exp(1e100) */
         GSL_EFAULT = 3, /* invalid pointer */
         GSL_EINVAL = 4, /* invalid argument supplied by user */
         GSL_EFAILED = 5, /* generic failure */
         GSL_EFACTOR = 6, /* factorization failed */
         GSL_ESANITY = 7, /* sanity check failed - shouldn't happen */
         GSL_ENOMEM = 8, /* malloc failed */
         GSL_EBADFUNC = 9, /* problem with user-supplied function */
         GSL_ERUNAWAY = 10, /* iterative process is out of control */
         GSL_EMAXITER = 11, /* exceeded max number of iterations */
         GSL_EZERODIV = 12, /* tried to divide by zero */
         GSL_EBADTOL = 13, /* user specified an invalid tolerance */
         GSL_ETOL = 14, /* failed to reach the specified tolerance */
         GSL_EUNDRFLW = 15, /* underflow */
         GSL_EOVRFLW = 16, /* overflow  */
         GSL_ELOSS = 17, /* loss of accuracy */
         GSL_EROUND = 18, /* failed because of roundoff error */
         GSL_EBADLEN = 19, /* matrix, vector lengths are not conformant */
         GSL_ENOTSQR = 20, /* matrix not square */
         GSL_ESING = 21, /* apparent singularity detected */
         GSL_EDIVERGE = 22, /* integral or series is divergent */
         GSL_EUNSUP = 23, /* requested feature is not supported by the hardware */
         GSL_EUNIMPL = 24, /* requested feature not (yet) implemented */
         GSL_ECACHE = 25, /* cache limit exceeded */
         GSL_ETABLE = 26, /* table limit exceeded */
         GSL_ENOPROG = 27, /* iteration is not making progress towards solution */
         GSL_ENOPROGJ = 28, /* jacobian evaluations are not improving the solution */
         GSL_ETOLF = 29, /* cannot reach the specified tolerance in F */
         GSL_ETOLX = 30, /* cannot reach the specified tolerance in X */
         GSL_ETOLG = 31, /* cannot reach the specified tolerance in gradient */
         GSL_EOF = 32 /* end of file */
     };
 
     GaussKronrodAdaptive() :
             Integrator1D() {
 
     }
 
     virtual ~GaussKronrodAdaptive() {
 
     }
 
     virtual void configure() {
 
     }
 
     typedef struct {
         size_t limit;
         size_t size;
         size_t nrmax;
         size_t i;
         size_t maximum_level;
         double *alist;
         double *blist;
         double *rlist;
         double *elist;
         size_t *order;
         size_t *level;
     } gsl_integration_workspace;
 
     struct extrapolation_table {
         size_t n;
         double rlist2[52];
         size_t nres;
         double res3la[3];
     };
 
     inline double integrate(FunctionType1D* pFunction, double a, double b,
             std::vector<double> &parameters) {
         double result = 0.;
         const size_t limit = 1000;
 
         gsl_integration_workspace * workspace = gsl_integration_workspace_alloc(
                 limit);
 
         compute(pFunction, a, b, parameters, result, workspace, limit);
 
         gsl_integration_workspace_free(workspace);
 
         return result;
     }
 
     inline int compute(FunctionType1D* pFunction, const double a,
             const double b, std::vector<double> &parameters, double &result,
             gsl_integration_workspace * workspace, const double limit) {
 
         double epsabs = m_tolerances.getAbsoluteTolerance();
         double epsrel = m_tolerances.getRelativeTolerance();
 
         //TODO when use those calls results are differents !!!
 //        double epsabs = m_tolerances.getAbsoluteTolerance();
 //        double epsrel = m_tolerances.getRelativeTolerance();
 
 // ##### ROOT MathMore section #####
         if (epsabs <= 0.) {
             // from static variable in ROOT::GSLIntegrator
             epsabs = 1e-06;
         }
         if (epsrel <= 0.) {
             // from static variable in ROOT::GSLIntegrator
             epsrel = 1e-09;
         }
         // ##################################
 
         double area, errsum;
         double res_ext, err_ext;
         double result0, abserr0, resabs0, resasc0;
         double tolerance;
 
         double ertest = 0;
         double error_over_large_intervals = 0;
         double reseps = 0, abseps = 0, correc = 0;
         size_t ktmin = 0;
         int roundoff_type1 = 0, roundoff_type2 = 0, roundoff_type3 = 0;
         int error_type = 0, error_type2 = 0;
 
         size_t iteration = 0;
 
         int positive_integrand = 0;
         int extrapolate = 0;
         int disallow_extrapolation = 0;
 
         struct extrapolation_table table;
 
         /* Initialize results */
 
         initialise(workspace, a, b);
 
         result = 0;
         double abserr = 0;
 
         if (limit > workspace->limit) {
             GSL_ERROR ("iteration limit exceeds available workspace", GSL_EINVAL)
             ;
         }
 
         /* Test on accuracy */
 
         if (epsabs <= 0
                 && (epsrel < 50 * GSL_DBL_EPSILON || epsrel < 0.5e-28)) {
             GSL_ERROR(
                     "tolerance cannot be acheived with given epsabs and epsrel",
                     GSL_EBADTOL);
         }
 
         /* Perform the first integration */
 
         gsl_integration_21qk(pFunction, a, b, result0, abserr0, resabs0,
                 resasc0, parameters);
 
         set_initial_result(workspace, result0, abserr0);
 
         tolerance = GSL_MAX_DBL(epsabs, epsrel * fabs(result0));
 
         if (abserr0 <= 100 * GSL_DBL_EPSILON * resabs0 && abserr0 > tolerance) {
             result = result0;
             abserr = abserr0;
             m_errors.setAbsolute(abserr);
 
             GSL_ERROR("cannot reach tolerance because of roundoff error"
                     "on first attempt", GSL_EROUND);
         } else if ((abserr0 <= tolerance && abserr0 != resasc0)
                 || abserr0 == 0.0) {
             result = result0;
             abserr = abserr0;
             m_errors.setAbsolute(abserr);
 
             return GSL_SUCCESS;
         } else if (limit == 1) {
             result = result0;
             abserr = abserr0;
             m_errors.setAbsolute(abserr);
 
             GSL_ERROR("a maximum of one iteration was insufficient",
                     GSL_EMAXITER);
         }
 
         /* Initialization */
 
         initialise_table(&table);
         append_table(&table, result0);
 
         area = result0;
         errsum = abserr0;
 
         res_ext = result0;
         err_ext = GSL_DBL_MAX;
 
         positive_integrand = test_positivity(result0, resabs0);
 
         iteration = 1;
 
         do {
             size_t current_level;
             double a1, b1, a2, b2;
             double a_i, b_i, r_i, e_i;
             double area1 = 0, area2 = 0, area12 = 0;
             double error1 = 0, error2 = 0, error12 = 0;
             double resasc1, resasc2;
             double resabs1, resabs2;
             double last_e_i;
 
             /* Bisect the subinterval with the largest error estimate */
 
             retrieve(workspace, &a_i, &b_i, &r_i, &e_i);
 
             current_level = workspace->level[workspace->i] + 1;
 
             a1 = a_i;
             b1 = 0.5 * (a_i + b_i);
             a2 = b1;
             b2 = b_i;
 
             iteration++;
 
             gsl_integration_21qk(pFunction, a1, b1, area1, error1, resabs1,
                     resasc1, parameters);
             gsl_integration_21qk(pFunction, a2, b2, area2, error2, resabs2,
                     resasc2, parameters);
 
             area12 = area1 + area2;
             error12 = error1 + error2;
             last_e_i = e_i;
 
             /* Improve previous approximations to the integral and test for
              accuracy.
 
              We write these expressions in the same way as the original
              QUADPACK code so that the rounding errors are the same, which
              makes testing easier. */
 
             errsum = errsum + error12 - e_i;
             area = area + area12 - r_i;
 
             tolerance = GSL_MAX_DBL(epsabs, epsrel * fabs(area));
 
             if (resasc1 != error1 && resasc2 != error2) {
                 double delta = r_i - area12;
 
                 if (fabs(delta) <= 1.0e-5 * fabs(area12)
                         && error12 >= 0.99 * e_i) {
                     if (!extrapolate) {
                         roundoff_type1++;
                     } else {
                         roundoff_type2++;
                     }
                 }
                 if (iteration > 10 && error12 > e_i) {
                     roundoff_type3++;
                 }
             }
 
             /* Test for roundoff and eventually set error flag */
 
             if (roundoff_type1 + roundoff_type2 >= 10 || roundoff_type3 >= 20) {
                 error_type = 2; /* round off error */
             }
 
             if (roundoff_type2 >= 5) {
                 error_type2 = 1;
             }
 
             /* set error flag in the case of bad integrand behaviour at
              a point of the integration range */
 
             if (subinterval_too_small(a1, a2, b2)) {
                 error_type = 4;
             }
 
             /* append the newly-created intervals to the list */
 
             update(workspace, a1, b1, area1, error1, a2, b2, area2, error2);
 
             if (errsum <= tolerance) {
                 goto compute_result;
             }
 
             if (error_type) {
                 break;
             }
 
             if (iteration >= limit - 1) {
                 error_type = 1;
                 break;
             }
 
             if (iteration == 2) /* set up variables on first iteration */
             {
                 error_over_large_intervals = errsum;
                 ertest = tolerance;
                 append_table(&table, area);
                 continue;
             }
 
             if (disallow_extrapolation) {
                 continue;
             }
 
             error_over_large_intervals += -last_e_i;
 
             if (current_level < workspace->maximum_level) {
                 error_over_large_intervals += error12;
             }
 
             if (!extrapolate) {
                 /* test whether the interval to be bisected next is the
                  smallest interval. */
 
                 if (large_interval(workspace))
                     continue;
 
                 extrapolate = 1;
                 workspace->nrmax = 1;
             }
 
             if (!error_type2 && error_over_large_intervals > ertest) {
                 if (increase_nrmax(workspace))
                     continue;
             }
 
             /* Perform extrapolation */
 
             append_table(&table, area);
 
             qelg(&table, reseps, abseps);
 
             ktmin++;
 
             if (ktmin > 5 && err_ext < 0.001 * errsum) {
                 error_type = 5;
             }
 
             if (abseps < err_ext) {
                 ktmin = 0;
                 err_ext = abseps;
                 res_ext = reseps;
                 correc = error_over_large_intervals;
                 ertest = GSL_MAX_DBL(epsabs, epsrel * fabs(reseps));
                 if (err_ext <= ertest)
                     break;
             }
 
             /* Prepare bisection of the smallest interval. */
 
             if (table.n == 1) {
                 disallow_extrapolation = 1;
             }
 
             if (error_type == 5) {
                 break;
             }
 
             /* work on interval with largest error */
 
             reset_nrmax(workspace);
             extrapolate = 0;
             error_over_large_intervals = errsum;
 
         } while (iteration < limit);
 
         result = res_ext;
         abserr = err_ext;
         m_errors.setAbsolute(abserr);
 
         if (err_ext == GSL_DBL_MAX)
             goto compute_result;
 
         if (error_type || error_type2) {
             if (error_type2) {
                 err_ext += correc;
             }
 
             if (error_type == 0)
                 error_type = 3;
 
             if (res_ext != 0.0 && area != 0.0) {
                 if (err_ext / fabs(res_ext) > errsum / fabs(area))
                     goto compute_result;
             } else if (err_ext > errsum) {
                 goto compute_result;
             } else if (area == 0.0) {
                 goto return_error;
             }
         }
 
         /*  Test on divergence. */
 
         {
             double max_area = GSL_MAX_DBL(fabs(res_ext), fabs(area));
 
             if (!positive_integrand && max_area < 0.01 * resabs0)
                 goto return_error;
         }
 
         {
             double ratio = res_ext / area;
 
             if (ratio < 0.01 || ratio > 100.0 || errsum > fabs(area))
                 error_type = 6;
         }
 
         goto return_error;
 
         compute_result:
 
         result = sum_results(workspace);
         abserr = errsum;
         m_errors.setAbsolute(abserr);
 
         return_error:
 
         if (error_type > 2)
             error_type--;
 
         if (error_type == 0) {
             return GSL_SUCCESS;
         } else if (error_type == 1) {
             GSL_ERROR("number of iterations was insufficient", GSL_EMAXITER);
         } else if (error_type == 2) {
             GSL_ERROR("cannot reach tolerance because of roundoff error",
                     GSL_EROUND);
         } else if (error_type == 3) {
             GSL_ERROR(
                     "bad integrand behavior found in the integration interval",
                     GSL_ESING);
         } else if (error_type == 4) {
             GSL_ERROR("roundoff error detected in the extrapolation table",
                     GSL_EROUND);
         } else if (error_type == 5) {
             GSL_ERROR("integral is divergent, or slowly convergent",
                     GSL_EDIVERGE);
         } else {
             GSL_ERROR("could not integrate function", GSL_EFAILED);
         }
     }
 
     gsl_integration_workspace *
     gsl_integration_workspace_alloc(const size_t n) {
         gsl_integration_workspace * w;
 
         if (n == 0) {
             GSL_ERROR_VAL("workspace length n must be positive integer",
                     GSL_EDOM, 0);
         }
 
         w = (gsl_integration_workspace *) malloc(
                 sizeof(gsl_integration_workspace));
 
         if (w == 0) {
             GSL_ERROR_VAL("failed to allocate space for workspace struct",
                     GSL_ENOMEM, 0);
         }
 
         w->alist = (double *) malloc(n * sizeof(double));
 
         if (w->alist == 0) {
             free(w); /* exception in constructor, avoid memory leak */
 
             GSL_ERROR_VAL("failed to allocate space for alist ranges",
                     GSL_ENOMEM, 0);
         }
 
         w->blist = (double *) malloc(n * sizeof(double));
 
         if (w->blist == 0) {
             free(w->alist);
             free(w); /* exception in constructor, avoid memory leak */
 
             GSL_ERROR_VAL("failed to allocate space for blist ranges",
                     GSL_ENOMEM, 0);
         }
 
         w->rlist = (double *) malloc(n * sizeof(double));
 
         if (w->rlist == 0) {
             free(w->blist);
             free(w->alist);
             free(w); /* exception in constructor, avoid memory leak */
 
             GSL_ERROR_VAL("failed to allocate space for rlist ranges",
                     GSL_ENOMEM, 0);
         }
 
         w->elist = (double *) malloc(n * sizeof(double));
 
         if (w->elist == 0) {
             free(w->rlist);
             free(w->blist);
             free(w->alist);
             free(w); /* exception in constructor, avoid memory leak */
 
             GSL_ERROR_VAL("failed to allocate space for elist ranges",
                     GSL_ENOMEM, 0);
         }
 
         w->order = (size_t *) malloc(n * sizeof(size_t));
 
         if (w->order == 0) {
             free(w->elist);
             free(w->rlist);
             free(w->blist);
             free(w->alist);
             free(w); /* exception in constructor, avoid memory leak */
 
             GSL_ERROR_VAL("failed to allocate space for order ranges",
                     GSL_ENOMEM, 0);
         }
 
         w->level = (size_t *) malloc(n * sizeof(size_t));
 
         if (w->level == 0) {
             free(w->order);
             free(w->elist);
             free(w->rlist);
             free(w->blist);
             free(w->alist);
             free(w); /* exception in constructor, avoid memory leak */
 
             GSL_ERROR_VAL("failed to allocate space for order ranges",
                     GSL_ENOMEM, 0);
         }
 
         w->size = 0;
         w->limit = n;
         w->maximum_level = 0;
 
         return w;
     }
 
     void gsl_integration_workspace_free(gsl_integration_workspace * w) {
         if (w != 0) {
             free(w->level);
             free(w->order);
             free(w->elist);
             free(w->rlist);
             free(w->blist);
             free(w->alist);
             free(w);
         }
     }
 
     inline void initialise(gsl_integration_workspace * workspace, double a,
             double b) {
         workspace->size = 0;
         workspace->nrmax = 0;
         workspace->i = 0;
         workspace->alist[0] = a;
         workspace->blist[0] = b;
         workspace->rlist[0] = 0.0;
         workspace->elist[0] = 0.0;
         workspace->order[0] = 0;
         workspace->level[0] = 0;
 
         workspace->maximum_level = 0;
     }
 
     inline void retrieve(const gsl_integration_workspace * workspace,
             double * a, double * b, double * r, double * e) {
         const size_t i = workspace->i;
         double * alist = workspace->alist;
         double * blist = workspace->blist;
         double * rlist = workspace->rlist;
         double * elist = workspace->elist;
 
         *a = alist[i];
         *b = blist[i];
         *r = rlist[i];
         *e = elist[i];
     }
 
     inline
     void update(gsl_integration_workspace * workspace, double a1, double b1,
             double area1, double error1, double a2, double b2, double area2,
             double error2) {
         double * alist = workspace->alist;
         double * blist = workspace->blist;
         double * rlist = workspace->rlist;
         double * elist = workspace->elist;
         size_t * level = workspace->level;
 
         const size_t i_max = workspace->i;
         const size_t i_new = workspace->size;
 
         const size_t new_level = workspace->level[i_max] + 1;
 
         /* append the newly-created intervals to the list */
 
         if (error2 > error1) {
             alist[i_max] = a2; /* blist[maxerr] is already == b2 */
             rlist[i_max] = area2;
             elist[i_max] = error2;
             level[i_max] = new_level;
 
             alist[i_new] = a1;
             blist[i_new] = b1;
             rlist[i_new] = area1;
             elist[i_new] = error1;
             level[i_new] = new_level;
         } else {
             blist[i_max] = b1; /* alist[maxerr] is already == a1 */
             rlist[i_max] = area1;
             elist[i_max] = error1;
             level[i_max] = new_level;
 
             alist[i_new] = a2;
             blist[i_new] = b2;
             rlist[i_new] = area2;
             elist[i_new] = error2;
             level[i_new] = new_level;
         }
 
         workspace->size++;
 
         if (new_level > workspace->maximum_level) {
             workspace->maximum_level = new_level;
         }
 
         qpsrt(workspace);
     }
 
     inline
     void qpsrt(gsl_integration_workspace * workspace) {
         const size_t last = workspace->size - 1;
         const size_t limit = workspace->limit;
 
         double * elist = workspace->elist;
         size_t * order = workspace->order;
 
         double errmax;
         double errmin;
         int i, k, top;
 
         size_t i_nrmax = workspace->nrmax;
         size_t i_maxerr = order[i_nrmax];
 
         /* Check whether the list contains more than two error estimates */
 
         if (last < 2) {
             order[0] = 0;
             order[1] = 1;
             workspace->i = i_maxerr;
             return;
         }
 
         errmax = elist[i_maxerr];
 
         /* 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. */
 
         while (i_nrmax > 0 && errmax > elist[order[i_nrmax - 1]]) {
             order[i_nrmax] = order[i_nrmax - 1];
             i_nrmax--;
         }
 
         /* 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)) {
             top = last;
         } else {
             top = limit - last + 1;
         }
 
         /* Insert errmax by traversing the list top-down, starting
          comparison from the element elist(order(i_nrmax+1)). */
 
         i = i_nrmax + 1;
 
         /* The order of the tests in the following line is important to
          prevent a segmentation fault */
 
         while (i < top && errmax < elist[order[i]]) {
             order[i - 1] = order[i];
             i++;
         }
 
         order[i - 1] = i_maxerr;
 
         /* Insert errmin by traversing the list bottom-up */
 
         errmin = elist[last];
 
         k = top - 1;
 
         while (k > i - 2 && errmin >= elist[order[k]]) {
             order[k + 1] = order[k];
             k--;
         }
 
         order[k + 1] = last;
 
         /* Set i_max and e_max */
 
         i_maxerr = order[i_nrmax];
 
         workspace->i = i_maxerr;
         workspace->nrmax = i_nrmax;
     }
 
     int increase_nrmax(gsl_integration_workspace * workspace) {
         int k;
         int id = workspace->nrmax;
         int jupbnd;
 
         const size_t * level = workspace->level;
         const size_t * order = workspace->order;
 
         size_t limit = workspace->limit;
         size_t last = workspace->size - 1;
 
         if (last > (1 + limit / 2)) {
             jupbnd = limit + 1 - last;
         } else {
             jupbnd = last;
         }
 
         for (k = id; k <= jupbnd; k++) {
             size_t i_max = order[workspace->nrmax];
 
             workspace->i = i_max;
 
             if (level[i_max] < workspace->maximum_level) {
                 return 1;
             }
 
             workspace->nrmax++;
 
         }
 
         return 0;
     }
 
     inline void reset_nrmax(gsl_integration_workspace * workspace) {
         workspace->nrmax = 0;
         workspace->i = workspace->order[0];
     }
 
     inline double sum_results(const gsl_integration_workspace * workspace) {
         const double * const rlist = workspace->rlist;
         const size_t n = workspace->size;
 
         size_t k;
         double result_sum = 0;
 
         for (k = 0; k < n; k++) {
             result_sum += rlist[k];
         }
 
         return result_sum;
     }
 
     int large_interval(gsl_integration_workspace * workspace) {
         size_t i = workspace->i;
         const size_t * level = workspace->level;
 
         if (level[i] < workspace->maximum_level) {
             return 1;
         } else {
             return 0;
         }
     }
 
     inline void set_initial_result(gsl_integration_workspace * workspace,
             double result, double error) {
         workspace->size = 1;
         workspace->rlist[0] = result;
         workspace->elist[0] = error;
     }
 
     void initialise_table(struct extrapolation_table *table) {
         table->n = 0;
         table->nres = 0;
     }
 
     void append_table(struct extrapolation_table *table, double y) {
         size_t n;
         n = table->n;
         table->rlist2[n] = y;
         table->n++;
     }
 
     inline void qelg(struct extrapolation_table *table, double &result,
             double &abserr) {
         double *epstab = table->rlist2;
         double *res3la = table->res3la;
         const size_t n = table->n - 1;
 
         const double current = epstab[n];
 
         double absolute = std::numeric_limits<double>::max();
         double relative = 5 * std::numeric_limits<double>::epsilon()
                 * fabs(current);
 
         const size_t newelm = n / 2;
         const size_t n_orig = n;
         size_t n_final = n;
         size_t i;
 
         const size_t nres_orig = table->nres;
 
         result = current;
         abserr = std::numeric_limits<double>::max();
 
         if (n < 2) {
             result = current;
             abserr = GSL_MAX_DBL(absolute, relative);
             return;
         }
 
         epstab[n + 2] = epstab[n];
         epstab[n] = std::numeric_limits<double>::max();
 
         for (i = 0; i < newelm; i++) {
             double res = epstab[n - 2 * i + 2];
             double e0 = epstab[n - 2 * i - 2];
             double e1 = epstab[n - 2 * i - 1];
             double e2 = res;
 
             double e1abs = fabs(e1);
             double delta2 = e2 - e1;
             double err2 = fabs(delta2);
             double tol2 = GSL_MAX_DBL(fabs(e2), e1abs)
                     * std::numeric_limits<double>::epsilon();
             double delta3 = e1 - e0;
             double err3 = fabs(delta3);
             double tol3 = GSL_MAX_DBL(e1abs, fabs(e0))
                     * std::numeric_limits<double>::epsilon();
 
             double e3, delta1, err1, tol1, ss;
 
             if (err2 <= tol2 && err3 <= tol3) {
                 /* If e0, e1 and e2 are equal to within machine accuracy,
                  convergence is assumed.  */
 
                 result = res;
                 absolute = err2 + err3;
                 relative = 5 * std::numeric_limits<double>::epsilon()
                         * fabs(res);
                 abserr = GSL_MAX_DBL(absolute, relative);
                 return;
             }
 
             e3 = epstab[n - 2 * i];
             epstab[n - 2 * i] = e1;
             delta1 = e1 - e3;
             err1 = fabs(delta1);
             tol1 = GSL_MAX_DBL(e1abs, fabs(e3))
                     * std::numeric_limits<double>::epsilon();
 
             /* 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) {
                 n_final = 2 * i;
                 break;
             }
 
             ss = (1 / delta1 + 1 / delta2) - 1 / delta3;
 
             /* Test to detect irregular behaviour in the table, and
              eventually omit a part of the table by adjusting the value of
              n. */
 
             if (fabs(ss * e1) <= 0.0001) {
                 n_final = 2 * i;
                 break;
             }
 
             /* Compute a new element and eventually adjust the value of
              result. */
 
             res = e1 + 1 / ss;
             epstab[n - 2 * i] = res;
 
             {
                 const double error = err2 + fabs(res - e2) + err3;
 
                 if (error <= abserr) {
                     abserr = error;
                     result = res;
                 }
             }
         }
 
         /* Shift the table */
 
         {
             const size_t limexp = 50 - 1;
 
             if (n_final == limexp) {
                 n_final = 2 * (limexp / 2);
             }
         }
 
         if (n_orig % 2 == 1) {
             for (i = 0; i <= newelm; i++) {
                 epstab[1 + i * 2] = epstab[i * 2 + 3];
             }
         } else {
             for (i = 0; i <= newelm; i++) {
                 epstab[i * 2] = epstab[i * 2 + 2];
             }
         }
 
         if (n_orig != n_final) {
             for (i = 0; i <= n_final; i++) {
                 epstab[i] = epstab[n_orig - n_final + i];
             }
         }
 
         table->n = n_final + 1;
 
         if (nres_orig < 3) {
             res3la[nres_orig] = result;
             abserr = std::numeric_limits<double>::max();
         } else { /* Compute error estimate */
             abserr = (fabs(result - res3la[2]) + fabs(result - res3la[1])
                     + fabs(result - res3la[0]));
 
             res3la[0] = res3la[1];
             res3la[1] = res3la[2];
             res3la[2] = result;
         }
 
         /* In QUADPACK the variable table->nres is incremented at the top of
          qelg, so it increases on every call. This leads to the array
          res3la being accessed when its elements are still undefined, so I
          have moved the update to this point so that its value more
          useful. */
 
         table->nres = nres_orig + 1;
 
         abserr = GSL_MAX_DBL(abserr,
                 5 * std::numeric_limits<double>::epsilon() * fabs(result));
 
         return;
     }
 
     inline int test_positivity(double result, double resabs) {
         int status = (fabs(result)
                 >= (1 - 50 * std::numeric_limits<double>::epsilon()) * resabs);
 
         return status;
     }
 
     inline int subinterval_too_small(double a1, double a2, double b2) {
         const double e = std::numeric_limits<double>::epsilon();
         const double u = std::numeric_limits<double>::min();
 
         double tmp = (1 + 100 * e) * (fabs(a2) + 1000 * u);
 
         int status = fabs(a1) <= tmp && fabs(b2) <= tmp;
 
         return status;
     }
 
     inline double GSL_MAX_DBL(double a, double b) {
         return GSL_MAX(a, b);
     }
 
     void gsl_integration_21qk(FunctionType1D* pFunction, double a, double b,
             double &result, double &abserr, double &resabs, double &resasc,
             std::vector<double> &parameters) {
 
         double x = 0.;
         const int n = 11;
         double fv1[11], fv2[11];
 
         const double center = 0.5 * (a + b);
         const double half_length = 0.5 * (b - a);
         const double abs_half_length = fabs(half_length);
 
         const double f_center = (*pFunction)(center, parameters);
 
         double result_gauss = 0;
         double result_kronrod = f_center * wgk[n - 1];
 
         double result_abs = fabs(result_kronrod);
         double result_asc = 0;
         double mean = 0, err = 0;
 
         int j;
 
         if (n % 2 == 0) {
             result_gauss = f_center * wg[n / 2 - 1];
         }
 
         for (j = 0; j < (n - 1) / 2; j++) {
             const int jtw = j * 2 + 1; /* in original fortran j=1,2,3 jtw=2,4,6 */
             const double abscissa = half_length * xgk[jtw];
 
             x = center - abscissa;
             const double fval1 = (*pFunction)(x, parameters);
 
             x = center + abscissa;
             const double fval2 = (*pFunction)(x, parameters);
 
             const double fsum = fval1 + fval2;
             fv1[jtw] = fval1;
             fv2[jtw] = fval2;
             result_gauss += wg[j] * fsum;
             result_kronrod += wgk[jtw] * fsum;
             result_abs += wgk[jtw] * (fabs(fval1) + fabs(fval2));
         }
 
         for (j = 0; j < n / 2; j++) {
             int jtwm1 = j * 2;
             const double abscissa = half_length * xgk[jtwm1];
 
             x = center - abscissa;
             const double fval1 = (*pFunction)(x, parameters);
 
             x = center + abscissa;
             const double fval2 = (*pFunction)(x, parameters);
             fv1[jtwm1] = fval1;
             fv2[jtwm1] = fval2;
             result_kronrod += wgk[jtwm1] * (fval1 + fval2);
             result_abs += wgk[jtwm1] * (fabs(fval1) + fabs(fval2));
         };
 
         mean = result_kronrod * 0.5;
 
         result_asc = wgk[n - 1] * fabs(f_center - mean);
 
         for (j = 0; j < n - 1; j++) {
             result_asc += wgk[j] * (fabs(fv1[j] - mean) + fabs(fv2[j] - mean));
         }
 
         /* scale by the width of the integration region */
 
         err = (result_kronrod - result_gauss) * half_length;
 
         result_kronrod *= half_length;
         result_abs *= abs_half_length;
         result_asc *= abs_half_length;
 
         result = result_kronrod;
         resabs = result_abs;
         resasc = result_asc;
         abserr = rescale_error(err, result_abs, result_asc);
     }
 
     double rescale_error(double err, const double result_abs,
             const double result_asc) {
         err = fabs(err);
 
         if (result_asc != 0 && err != 0) {
             double scale = pow((200 * err / result_asc), 1.5);
 
             if (scale < 1) {
                 err = result_asc * scale;
             } else {
                 err = result_asc;
             }
         }
         if (result_abs
                 > std::numeric_limits<double>::min()
                         / (50 * std::numeric_limits<double>::epsilon())) {
             double min_err = 50 * std::numeric_limits<double>::epsilon()
                     * result_abs;
 
             if (min_err > err) {
                 err = min_err;
             }
         }
 
         return err;
     }
 
     double GSL_MAX(double a, double b) {
         return (a > b) ? a : b;
     }
 
     static const std::vector<double> xgk;
     static const std::vector<double> wg;
     static const std::vector<double> wgk;
 
     static const std::vector<double> make_xgk_21();
     static const std::vector<double> make_wg_21();
     static const std::vector<double> make_wgk_21();
 
     virtual GaussKronrodAdaptive* clone() const {
         return new GaussKronrodAdaptive(*this);
     }
 
 protected:
     GaussKronrodAdaptive(const GaussKronrodAdaptive &other) :
             Integrator1D(other) {
     }
 };
 
 } // namespace NumA
 
 #endif /* GAUSS_KRONROD_ADAPTIVE_H */

This page was automatically generated by the 2.3.7 LXR engine. The LXR team