rule15gauss_evalError()

static int rule15gauss_evalError ( rule *  r, unsigned  fdim, integrand_v  f, void *  fdata, unsigned  nR, region *  R  )

static

{
     /* Gauss quadrature weights and kronrod quadrature abscissae and
        weights as evaluated with 80 decimal digit arithmetic by
        L. W. Fullerton, Bell Labs, Nov. 1981. */
     const unsigned n = 8;
     const double xgk[8] = {  /* abscissae of the 15-point kronrod rule */
          0.991455371120812639206854697526329,
          0.949107912342758524526189684047851,
          0.864864423359769072789712788640926,
          0.741531185599394439863864773280788,
          0.586087235467691130294144838258730,
          0.405845151377397166906606412076961,
          0.207784955007898467600689403773245,
          0.000000000000000000000000000000000
          /* xgk[1], xgk[3], ... abscissae of the 7-point gauss rule. 
             xgk[0], xgk[2], ... to optimally extend the 7-point gauss rule */
     };
     static const double wg[4] = {  /* weights of the 7-point gauss rule */
          0.129484966168869693270611432679082,
          0.279705391489276667901467771423780,
          0.381830050505118944950369775488975,
          0.417959183673469387755102040816327
     };
     static const double wgk[8] = { /* weights of the 15-point kronrod rule */
          0.022935322010529224963732008058970,
          0.063092092629978553290700663189204,
          0.104790010322250183839876322541518,
          0.140653259715525918745189590510238,
          0.169004726639267902826583426598550,
          0.190350578064785409913256402421014,
          0.204432940075298892414161999234649,
          0.209482141084727828012999174891714
     };
     unsigned j, k, iR, npts = 0;
     double *pts, *vals;

     if (alloc_rule_pts(r, nR)) return FAILURE;
     pts = r->pts; vals = r->vals;

     for (iR = 0; iR < nR; ++iR) {
          const double center = R[iR].h.data[0];
          const double halfwidth = R[iR].h.data[1];

          pts[npts++] = center;

          for (j = 0; j < (n - 1) / 2; ++j) {
               int j2 = 2*j + 1;
               double w = halfwidth * xgk[j2];
               pts[npts++] = center - w;
               pts[npts++] = center + w;
          }
          for (j = 0; j < n/2; ++j) {
               int j2 = 2*j;
               double w = halfwidth * xgk[j2];
               pts[npts++] = center - w;
               pts[npts++] = center + w;
          }

          R[iR].splitDim = 0; /* no choice but to divide 0th dimension */
     }

     f(1, npts, pts, fdata, fdim, vals);
     
     for (k = 0; k < fdim; ++k) {
          for (iR = 0; iR < nR; ++iR) {
               const double halfwidth = R[iR].h.data[1];
               double result_gauss = vals[0] * wg[n/2 - 1];
               double result_kronrod = vals[0] * wgk[n - 1];
               double result_abs = fabs(result_kronrod);
               double result_asc, mean, err;

               /* accumulate integrals */
               npts = 1;
               for (j = 0; j < (n - 1) / 2; ++j) {
                    int j2 = 2*j + 1;
                    double v = vals[npts] + vals[npts+1];
                    result_gauss += wg[j] * v;
                    result_kronrod += wgk[j2] * v;
                    result_abs += wgk[j2] * (fabs(vals[npts]) 
                                             + fabs(vals[npts+1]));
                    npts += 2;
               }
               for (j = 0; j < n/2; ++j) {
                    int j2 = 2*j;
                    result_kronrod += wgk[j2] * (vals[npts] + vals[npts+1]);
                    result_abs += wgk[j2] * (fabs(vals[npts]) 
                                             + fabs(vals[npts+1]));
                    npts += 2;
               }
               
               /* integration result */
               R[iR].ee[k].val = result_kronrod * halfwidth;

               /* error estimate 
                  (from GSL, probably dates back to QUADPACK
                  ... not completely clear to me why we don't just use
                  fabs(result_kronrod - result_gauss) * halfwidth */
               mean = result_kronrod * 0.5;
               result_asc = wgk[n - 1] * fabs(vals[0] - mean);
               npts = 1;
               for (j = 0; j < (n - 1) / 2; ++j) {
                    int j2 = 2*j + 1;
                    result_asc += wgk[j2] * (fabs(vals[npts]-mean)
                                             + fabs(vals[npts+1]-mean));
                    npts += 2;
               }
               for (j = 0; j < n/2; ++j) {
                    int j2 = 2*j;
                    result_asc += wgk[j2] * (fabs(vals[npts]-mean)
                                             + fabs(vals[npts+1]-mean));
                    npts += 2;
               }
               err = fabs(result_kronrod - result_gauss) * halfwidth;
               result_abs *= halfwidth;
               result_asc *= halfwidth;
               if (result_asc != 0 && err != 0) {
                    double scale = pow((200 * err / result_asc), 1.5);
                    err = (scale < 1) ? result_asc * scale : result_asc;
               }
               if (result_abs > DBL_MIN / (50 * DBL_EPSILON)) {
                    double min_err = 50 * DBL_EPSILON * result_abs;
                    if (min_err > err) err = min_err;
               }
               R[iR].ee[k].err = err;
               
               /* increment vals to point to next batch of results */
               vals += 15;
          }
     }
     return SUCCESS;
}