Integration

O2scl

Integration contents

One-dimensional integration based on GSL, CERNLIB, and Boost

Several classes integrate arbitrary one-dimensional functions:

\[W(x)=(x-a)^{\alpha}(b-x)^{\beta}\log^{\mu}(x-a)\log^{\nu}(b-x)\]

is performed by inte_qaws_gsl.

Note that some of the integrators support multiprecision, see Multiprecision Support.

There are two competing factors that can slow down an adaptive integration algorithm: (1) each evaluation of the integrand can be numerically expensive, depending on how the function is defined, and (2) the process of subdividing regions and recalculating values is almost always numerically expensive in its own right. For integrands that are very smooth (e.g., analytic functions), a high-order Gauss-Kronrod rule (e.g., 61-point) will achieve the desired error tolerance with relatively few subdivisions. For integrands with discontinuities or singular derivatives, a low-order rule (e.g., 15-point) is often more efficient.

GSL-based integration details

For the GSL-based integration routines, the variables o2scl::inte::tol_abs and o2scl::inte::tol_rel have the same role as the quantities usually denoted in the GSL integration routines by epsabs and epsrel. In particular, the integration classes attempt to ensure that

\[|\mathrm{result}-I| \leq \mathrm{Max}(\mathrm{tol\_abs}, \mathrm{tol\_rel}|I|)\]

and returns an error to attempt to ensure that

\[|\mathrm{result}-I| \leq \mathrm{abserr} \leq \mathrm{Max}(\mathrm{tol\_abs},\mathrm{tol\_rel}|I|)\]

where I is the integral to be evaluated. Even when the corresponding descendant of o2scl::inte::integ() returns success, these inequalities may fail for sufficiently difficult functions. All of the GSL integration routines except for inte_qng_gsl use a workspace given in inte_workspace_gsl which holds the results of the various subdivisions of the original interval.

The GSL routines were originally based on QUADPACK, which is available at http://www.netlib.org/quadpack .

For adaptive GSL integration classes, the type of Gauss-Kronrod quadrature rule that is used to approximate the integral and estimate the error of a subinterval is set by o2scl::inte_kronrod_gsl::set_rule().

The number of subdivisions of the integration region is limited by the size of the workspace, set in o2scl::inte_kronrod_gsl::set_limit(). The number of subdivisions required for the most recent call to o2scl::inte::integ() or o2scl::inte::integ_err() is given in o2scl::inte::last_iter. This number will always be less than or equal to the workspace size.

Note

The GSL integration routines can sometimes lose precision if the integrand is everywhere much smaller than unity. Some rescaling may be required in these cases.

GSL-based integration error messages

The error messages given by the adaptive GSL integration routines tend to follow a standard form and are documented here. There are several error messages which indicate improper usage and cause the error handler to be called regardless of the value of o2scl::inte::err_nonconv:

  • Iteration limit exceeds workspace in class::function(). [ exc_einval ]

  • Could not integrate function in class::function() (it may have returned a non-finite result). [ exc_efailed ] This often occurs when the user-specified function returns inf or nan.

  • Tolerance cannot be achieved with given value of tol_abs and tol_rel in class::function(). [ exc_ebadtol ] This occurs if the user supplies unreasonable values for o2scl::inte::tol_abs and o2scl::inte::tol_rel. All positive values for o2scl::inte::tol_abs are allowed. If zero-tolerance for o2scl::inte::tol_abs is desired, then o2scl::inte::tol_rel must be at least \(50 \cdot \epsilon_\mathrm{mach}\) (\(\approx 1.11 \times 10^{-14}\) ).

  • Cannot integrate with singularity on endpoint in inte_qawc_gsl::qawc(). [ exc_einval ] The class ref o2scl::inte_qawc_gsl cannot handle the case when a singularity is one of the endpoints of the integration.

There are also convergence errors which will call the error handler unless o2scl::inte::err_nonconv is false. See What is an error? for more discussion on convergence errors versus fatal errors:

  • Cannot reach tolerance because of roundoff error on first attempt in class::function(). [ exc_eround ] Each integration attempt tests for round-off errors by comparing the computed integral with that of the integrand’s absolute value (i.e., \(L^1\) -norm). A highly oscillatory integrand may cause this error.

  • A maximum of 1 iteration was insufficient in class::function(). [ exc_emaxiter ] This occurs if the workspace is allocated for one interval and a single Gauss-Kronrod integration does not yield the accuracy demanded by o2scl::inte::tol_abs and o2scl::inte::tol_rel.

  • Bad integrand behavior in class::function(). [ exc_esing ] This occurs if the integrand is (effectively) singular in a region, causing the subdivided intervals to become too small for floating-point precision.

  • Maximum number of subdivisions 'value' reached in class::function(). [ exc_emaxiter ] This occurs if the refinement algorithm runs out of allocated workspace. The number of iterations required for the most recent call to o2scl::inte::integ() or o2scl::inte::integ_err() is given in o2scl::inte::last_iter. This number will always be less than or equal to the workspace size.

  • Roundoff error prevents tolerance from being achieved in class::function(). [ exc_eround ] The refinement procedure counts round-off errors as they occur and terminates if too many such errors accumulate.

  • Roundoff error detected in extrapolation table in inte_singular_gsl::qags(). [ exc_eround ] This occurs when error-terms from the \(\varepsilon\) -algorithm are are monitored and compared with the error-terms from the refinement procedure. The algorithm terminates if these sequences differ by too many orders of magnitude. See o2scl::inte_singular_gsl::qelg().

  • Integral is divergent or slowly convergent in inte_singular_gsl::qags(). [ exc_ediverge ] This occurs if the approximations produced by the refinement algorithm and the extrapolation algorithm differ by too many orders of magnitude.

  • Exceeded limit of trigonometric table in inte_qawo_gsl_sin()::qawo(). [ exc_etable ] This occurs if the maximum level of the table of Chebyshev moments is reached.

One-dimensional integration example

This example computes the integral \(\int_{-\infty}^{\infty} e^{-x^2} ~dx\) with inte_qagi_gsl, the integral \(\int_0^{\infty} e^{-x^2} ~dx\) with inte_qagiu_gsl, the integral \(\int_{-\infty}^{0} e^{-x^2} ~dx\) with inte_qagil_gsl, and the integral \(\int_0^1 \left[ \sin (2 x) + \frac{1}{2} \right]~dx\) with both inte_qag_gsl and inte_adapt_cern_tl, and compares the computed results with the exact results.

/* Example: ex_inte.cpp
   -------------------------------------------------------------------
   An example to demonstrate numerical integration.
*/

#include <cmath>
#include <o2scl/test_mgr.h>
#include <o2scl/constants.h>
#include <o2scl/funct.h>
#include <o2scl/inte_qag_gsl.h>
#include <o2scl/inte_qagi_gsl.h>
#include <o2scl/inte_qagiu_gsl.h>
#include <o2scl/inte_qagil_gsl.h>
#include <o2scl/inte_adapt_cern.h>

using namespace std;
using namespace o2scl;
using namespace o2scl_const;

class cl {

public:

  // We'll use this to count the number of function
  // evaulations required by the integration routines
  int nf;

  // A function to be integrated
  double integrand(double x) {
    nf++;
    return exp(-x*x);
  }

  // Another function to be integrated
  double integrand2(double x) {
    nf++;
    return sin(2.0*x)+0.5;
  }
};

int main(void) {
  cl acl;
  test_mgr t;

  t.set_output_level(1);

  funct f1=std::bind(std::mem_fn<double(double)>
		       (&cl::integrand),&acl,std::placeholders::_1);
  funct f2=std::bind(std::mem_fn<double(double)>
		       (&cl::integrand2),&acl,std::placeholders::_1);

  // We don't need to specify the function type in the integration
  // objects, because we're using the default function type (type
  // funct).
  inte_qag_gsl<> g;
  inte_qagi_gsl<> gi;
  inte_qagiu_gsl<> gu;
  inte_qagil_gsl<> gl;
  inte_adapt_cern ca;

  // The result and the uncertainty
  double res, err;
  
  // An integral from -infinity to +infinity (the limits are ignored)
  acl.nf=0;
  int ret1=gi.integ_err(f1,0.0,0.0,res,err);
  cout << "inte_qagi_gsl: " << endl;
  cout << "Return value: " << ret1 << endl;
  cout << "Result: " << res << " Uncertainty: " << err << endl;
  cout << "Number of iterations: " << gi.last_iter << endl;
  cout << "Number of function evaluations: " << acl.nf << endl;
  cout << endl;
  t.test_rel(res,sqrt(pi),1.0e-8,"inte 1");
  
  // An integral from 0 to +infinity (the second limit argument is
  // ignored in the line below)
  acl.nf=0;
  gu.integ_err(f1,0.0,0.0,res,err);
  cout << "inte_qagiu_gsl: " << endl;
  cout << "Return value: " << ret1 << endl;
  cout << "Result: " << res << " Uncertainty: " << err << endl;
  cout << "Number of iterations: " << gu.last_iter << endl;
  cout << "Number of function evaluations: " << acl.nf << endl;
  cout << endl;
  t.test_rel(res,sqrt(pi)/2.0,1.0e-8,"inte 2");

  // An integral from -infinity to zero (the first limit argument is
  // ignored in the line below)
  acl.nf=0;
  gl.integ_err(f1,0.0,0.0,res,err);
  cout << "inte_qagil_gsl: " << endl;
  cout << "Return value: " << ret1 << endl;
  cout << "Result: " << res << " Uncertainty: " << err << endl;
  cout << "Number of iterations: " << gl.last_iter << endl;
  cout << "Number of function evaluations: " << acl.nf << endl;
  cout << endl;
  t.test_rel(res,sqrt(pi)/2.0,1.0e-8,"inte 3");

  // An integral from 0 to 1 with the GSL integrator
  acl.nf=0;
  g.integ_err(f2,0.0,1.0,res,err);
  cout << "inte_qag_gsl: " << endl;
  cout << "Return value: " << ret1 << endl;
  cout << "Result: " << res << " Uncertainty: " << err << endl;
  cout << "Number of iterations: " << g.last_iter << endl;
  cout << "Number of function evaluations: " << acl.nf << endl;
  cout << endl;
  t.test_rel(res,0.5+sin(1.0)*sin(1.0),1.0e-8,"inte 4");

  // The same integral with the CERNLIB integrator
  acl.nf=0;
  ca.integ_err(f2,0.0,1.0,res,err);
  cout << "inte_adapt_cern: " << endl;
  cout << "Return value: " << ret1 << endl;
  cout << "Result: " << res << " Uncertainty: " << err << endl;
  cout << "Number of iterations: " << ca.last_iter << endl;
  cout << "Number of function evaluations: " << acl.nf << endl;
  cout << endl;
  t.test_rel(res,0.5+sin(1.0)*sin(1.0),1.0e-8,"inte 5");

  t.report();
  return 0;
}

Gauss-Kronrod integration coefficients

Top

namespace o2scl_inte_gk_coeffs

A namespace for the GSL adaptive Gauss-Kronrod integration coefficients.

The storage scheme follows that of QUADPACK: symmetry about the

origin permits storing only those abscissae in the interval

\( [0, 1). \) Similiarly for the weights. The odd-indexed abscissae xgk[1], xgk[3],… belong to the low-order Gauss rule. For the full Gauss-Kronrod rule, the array wgk[] holds the weights corresponding to the abscissae xgk[]. For the low-order rule, wg[i] is the weight correpsonding to xgk[2*i+1].

Documentation from GSL

:

Gauss quadrature weights and kronrod quadrature abscissae and weights as evaluated with 80 decimal digit arithmetic by L. W. Fullerton, Bell Labs, Nov. 1981.

Variables

static const double qk15_xgk[8] = {0.991455371120812639206854697526329, 0.949107912342758524526189684047851, 0.864864423359769072789712788640926, 0.741531185599394439863864773280788, 0.586087235467691130294144838258730, 0.405845151377397166906606412076961, 0.207784955007898467600689403773245, 0.000000000000000000000000000000000}

Abscissae of the 15-point Kronrod rule.

static const double qk15_wg[4] = {0.129484966168869693270611432679082, 0.279705391489276667901467771423780, 0.381830050505118944950369775488975, 0.417959183673469387755102040816327}

Weights of the 7-point Gauss rule.

static const double qk15_wgk[8] = {0.022935322010529224963732008058970, 0.063092092629978553290700663189204, 0.104790010322250183839876322541518, 0.140653259715525918745189590510238, 0.169004726639267902826583426598550, 0.190350578064785409913256402421014, 0.204432940075298892414161999234649, 0.209482141084727828012999174891714}

Weights of the 15-point Kronrod rule.

static const double qk21_xgk[11] = {0.995657163025808080735527280689003, 0.973906528517171720077964012084452, 0.930157491355708226001207180059508, 0.865063366688984510732096688423493, 0.780817726586416897063717578345042, 0.679409568299024406234327365114874, 0.562757134668604683339000099272694, 0.433395394129247190799265943165784, 0.294392862701460198131126603103866, 0.148874338981631210884826001129720, 0.000000000000000000000000000000000}

Abscissae of the 21-point Kronrod rule.

static const double qk21_wg[5] = {0.066671344308688137593568809893332, 0.149451349150580593145776339657697, 0.219086362515982043995534934228163, 0.269266719309996355091226921569469, 0.295524224714752870173892994651338}

Weights of the 10-point Gauss rule.

static const double qk21_wgk[11] = {0.011694638867371874278064396062192, 0.032558162307964727478818972459390, 0.054755896574351996031381300244580, 0.075039674810919952767043140916190, 0.093125454583697605535065465083366, 0.109387158802297641899210590325805, 0.123491976262065851077958109831074, 0.134709217311473325928054001771707, 0.142775938577060080797094273138717, 0.147739104901338491374841515972068, 0.149445554002916905664936468389821}

Weights of the 21-point Kronrod rule.

static const double qk31_xgk[16] = {0.998002298693397060285172840152271, 0.987992518020485428489565718586613, 0.967739075679139134257347978784337, 0.937273392400705904307758947710209, 0.897264532344081900882509656454496, 0.848206583410427216200648320774217, 0.790418501442465932967649294817947, 0.724417731360170047416186054613938, 0.650996741297416970533735895313275, 0.570972172608538847537226737253911, 0.485081863640239680693655740232351, 0.394151347077563369897207370981045, 0.299180007153168812166780024266389, 0.201194093997434522300628303394596, 0.101142066918717499027074231447392, 0.000000000000000000000000000000000}

Abscissae of the 31-point Kronrod rule.

static const double qk31_wg[8] = {0.030753241996117268354628393577204, 0.070366047488108124709267416450667, 0.107159220467171935011869546685869, 0.139570677926154314447804794511028, 0.166269205816993933553200860481209, 0.186161000015562211026800561866423, 0.198431485327111576456118326443839, 0.202578241925561272880620199967519}

Weights of the 15-point Gauss rule.

static const double qk31_wgk[16] = {0.005377479872923348987792051430128, 0.015007947329316122538374763075807, 0.025460847326715320186874001019653, 0.035346360791375846222037948478360, 0.044589751324764876608227299373280, 0.053481524690928087265343147239430, 0.062009567800670640285139230960803, 0.069854121318728258709520077099147, 0.076849680757720378894432777482659, 0.083080502823133021038289247286104, 0.088564443056211770647275443693774, 0.093126598170825321225486872747346, 0.096642726983623678505179907627589, 0.099173598721791959332393173484603, 0.100769845523875595044946662617570, 0.101330007014791549017374792767493}

Weights of the 31-point Kronrod rule.

static const double qk41_xgk[21] = {0.998859031588277663838315576545863, 0.993128599185094924786122388471320, 0.981507877450250259193342994720217, 0.963971927277913791267666131197277, 0.940822633831754753519982722212443, 0.912234428251325905867752441203298, 0.878276811252281976077442995113078, 0.839116971822218823394529061701521, 0.795041428837551198350638833272788, 0.746331906460150792614305070355642, 0.693237656334751384805490711845932, 0.636053680726515025452836696226286, 0.575140446819710315342946036586425, 0.510867001950827098004364050955251, 0.443593175238725103199992213492640, 0.373706088715419560672548177024927, 0.301627868114913004320555356858592, 0.227785851141645078080496195368575, 0.152605465240922675505220241022678, 0.076526521133497333754640409398838, 0.000000000000000000000000000000000}

Abscissae of the 41-point Kronrod rule.

static const double qk41_wg[11] = {0.017614007139152118311861962351853, 0.040601429800386941331039952274932, 0.062672048334109063569506535187042, 0.083276741576704748724758143222046, 0.101930119817240435036750135480350, 0.118194531961518417312377377711382, 0.131688638449176626898494499748163, 0.142096109318382051329298325067165, 0.149172986472603746787828737001969, 0.152753387130725850698084331955098}

Weights of the 20-point Gauss rule.

static const double qk41_wgk[21] = {0.003073583718520531501218293246031, 0.008600269855642942198661787950102, 0.014626169256971252983787960308868, 0.020388373461266523598010231432755, 0.025882133604951158834505067096153, 0.031287306777032798958543119323801, 0.036600169758200798030557240707211, 0.041668873327973686263788305936895, 0.046434821867497674720231880926108, 0.050944573923728691932707670050345, 0.055195105348285994744832372419777, 0.059111400880639572374967220648594, 0.062653237554781168025870122174255, 0.065834597133618422111563556969398, 0.068648672928521619345623411885368, 0.071054423553444068305790361723210, 0.073030690332786667495189417658913, 0.074582875400499188986581418362488, 0.075704497684556674659542775376617, 0.076377867672080736705502835038061, 0.076600711917999656445049901530102}

Weights of the 41-point Kronrod rule.

static const double qk51_xgk[26] = {0.999262104992609834193457486540341, 0.995556969790498097908784946893902, 0.988035794534077247637331014577406, 0.976663921459517511498315386479594, 0.961614986425842512418130033660167, 0.942974571228974339414011169658471, 0.920747115281701561746346084546331, 0.894991997878275368851042006782805, 0.865847065293275595448996969588340, 0.833442628760834001421021108693570, 0.797873797998500059410410904994307, 0.759259263037357630577282865204361, 0.717766406813084388186654079773298, 0.673566368473468364485120633247622, 0.626810099010317412788122681624518, 0.577662930241222967723689841612654, 0.526325284334719182599623778158010, 0.473002731445714960522182115009192, 0.417885382193037748851814394594572, 0.361172305809387837735821730127641, 0.303089538931107830167478909980339, 0.243866883720988432045190362797452, 0.183718939421048892015969888759528, 0.122864692610710396387359818808037, 0.061544483005685078886546392366797, 0.000000000000000000000000000000000}

Abscissae of the 51-point Kronrod rule.

static const double qk51_wg[13] = {0.011393798501026287947902964113235, 0.026354986615032137261901815295299, 0.040939156701306312655623487711646, 0.054904695975835191925936891540473, 0.068038333812356917207187185656708, 0.080140700335001018013234959669111, 0.091028261982963649811497220702892, 0.100535949067050644202206890392686, 0.108519624474263653116093957050117, 0.114858259145711648339325545869556, 0.119455763535784772228178126512901, 0.122242442990310041688959518945852, 0.123176053726715451203902873079050}

Weights of the 25-point Gauss rule.

static const double qk51_wgk[26] = {0.001987383892330315926507851882843, 0.005561932135356713758040236901066, 0.009473973386174151607207710523655, 0.013236229195571674813656405846976, 0.016847817709128298231516667536336, 0.020435371145882835456568292235939, 0.024009945606953216220092489164881, 0.027475317587851737802948455517811, 0.030792300167387488891109020215229, 0.034002130274329337836748795229551, 0.037116271483415543560330625367620, 0.040083825504032382074839284467076, 0.042872845020170049476895792439495, 0.045502913049921788909870584752660, 0.047982537138836713906392255756915, 0.050277679080715671963325259433440, 0.052362885806407475864366712137873, 0.054251129888545490144543370459876, 0.055950811220412317308240686382747, 0.057437116361567832853582693939506, 0.058689680022394207961974175856788, 0.059720340324174059979099291932562, 0.060539455376045862945360267517565, 0.061128509717053048305859030416293, 0.061471189871425316661544131965264, 0.061580818067832935078759824240066}

Weights of the 51-point Kronrod rule.

static const double qk61_xgk[31]

Abscissae of the 61-point Kronrod rule.

static const double qk61_wg[15] = {0.007968192496166605615465883474674, 0.018466468311090959142302131912047, 0.028784707883323369349719179611292, 0.038799192569627049596801936446348, 0.048402672830594052902938140422808, 0.057493156217619066481721689402056, 0.065974229882180495128128515115962, 0.073755974737705206268243850022191, 0.080755895229420215354694938460530, 0.086899787201082979802387530715126, 0.092122522237786128717632707087619, 0.096368737174644259639468626351810, 0.099593420586795267062780282103569, 0.101762389748405504596428952168554, 0.102852652893558840341285636705415}

Weights of the 30-point Gauss rule.

static const double qk61_wgk[31]

Weights of the 61-point Kronrod rule.

Non-adaptive quadrature integration coefficients

Top

namespace o2scl_inte_qng_coeffs

A namespace for the quadrature coefficients for non-adaptive integration.

Documentation from GSL

:

Gauss-Kronrod-Patterson quadrature coefficients for use in quadpack routine qng. These coefficients were calculated with 101 decimal digit arithmetic by L. W. Fullerton, Bell Labs, Nov

Variables

static const double x1[5] = {0.973906528517171720077964012084452, 0.865063366688984510732096688423493, 0.679409568299024406234327365114874, 0.433395394129247190799265943165784, 0.148874338981631210884826001129720}

Abcissae common to the 10-, 21-, 43- and 87-point rule

static const double w10[5] = {0.066671344308688137593568809893332, 0.149451349150580593145776339657697, 0.219086362515982043995534934228163, 0.269266719309996355091226921569469, 0.295524224714752870173892994651338}

Weights of the 10-point formula

static const double x2[5] = {0.995657163025808080735527280689003, 0.930157491355708226001207180059508, 0.780817726586416897063717578345042, 0.562757134668604683339000099272694, 0.294392862701460198131126603103866}

Abcissae common to the 21-, 43- and 87-point rule

static const double w21a[5] = {0.032558162307964727478818972459390, 0.075039674810919952767043140916190, 0.109387158802297641899210590325805, 0.134709217311473325928054001771707, 0.147739104901338491374841515972068}

Weights of the 21-point formula for abcissae x1

static const double w21b[6] = {0.011694638867371874278064396062192, 0.054755896574351996031381300244580, 0.093125454583697605535065465083366, 0.123491976262065851077958109831074, 0.142775938577060080797094273138717, 0.149445554002916905664936468389821}

Weights of the 21-point formula for abcissae x2

static const double x3[11] = {0.999333360901932081394099323919911, 0.987433402908088869795961478381209, 0.954807934814266299257919200290473, 0.900148695748328293625099494069092, 0.825198314983114150847066732588520, 0.732148388989304982612354848755461, 0.622847970537725238641159120344323, 0.499479574071056499952214885499755, 0.364901661346580768043989548502644, 0.222254919776601296498260928066212, 0.074650617461383322043914435796506}

Abscissae common to the 43- and 87-point rule

static const double w43a[10] = {0.016296734289666564924281974617663, 0.037522876120869501461613795898115, 0.054694902058255442147212685465005, 0.067355414609478086075553166302174, 0.073870199632393953432140695251367, 0.005768556059769796184184327908655, 0.027371890593248842081276069289151, 0.046560826910428830743339154433824, 0.061744995201442564496240336030883, 0.071387267268693397768559114425516}

Weights of the 43-point formula for abscissae x1, x3

static const double w43b[12] = {0.001844477640212414100389106552965, 0.010798689585891651740465406741293, 0.021895363867795428102523123075149, 0.032597463975345689443882222526137, 0.042163137935191811847627924327955, 0.050741939600184577780189020092084, 0.058379395542619248375475369330206, 0.064746404951445885544689259517511, 0.069566197912356484528633315038405, 0.072824441471833208150939535192842, 0.074507751014175118273571813842889, 0.074722147517403005594425168280423}

Weights of the 43-point formula for abscissae x3

static const double x4[22] = {0.999902977262729234490529830591582, 0.997989895986678745427496322365960, 0.992175497860687222808523352251425, 0.981358163572712773571916941623894, 0.965057623858384619128284110607926, 0.943167613133670596816416634507426, 0.915806414685507209591826430720050, 0.883221657771316501372117548744163, 0.845710748462415666605902011504855, 0.803557658035230982788739474980964, 0.757005730685495558328942793432020, 0.706273209787321819824094274740840, 0.651589466501177922534422205016736, 0.593223374057961088875273770349144, 0.531493605970831932285268948562671, 0.466763623042022844871966781659270, 0.399424847859218804732101665817923, 0.329874877106188288265053371824597, 0.258503559202161551802280975429025, 0.185695396568346652015917141167606, 0.111842213179907468172398359241362, 0.037352123394619870814998165437704}

Abscissae of the 87-point rule

static const double w87a[21] = {0.008148377384149172900002878448190, 0.018761438201562822243935059003794, 0.027347451050052286161582829741283, 0.033677707311637930046581056957588, 0.036935099820427907614589586742499, 0.002884872430211530501334156248695, 0.013685946022712701888950035273128, 0.023280413502888311123409291030404, 0.030872497611713358675466394126442, 0.035693633639418770719351355457044, 0.000915283345202241360843392549948, 0.005399280219300471367738743391053, 0.010947679601118931134327826856808, 0.016298731696787335262665703223280, 0.021081568889203835112433060188190, 0.025370969769253827243467999831710, 0.029189697756475752501446154084920, 0.032373202467202789685788194889595, 0.034783098950365142750781997949596, 0.036412220731351787562801163687577, 0.037253875503047708539592001191226}

Weights of the 87-point formula for abscissae x1, x2, x3

static const double w87b[23] = {0.000274145563762072350016527092881, 0.001807124155057942948341311753254, 0.004096869282759164864458070683480, 0.006758290051847378699816577897424, 0.009549957672201646536053581325377, 0.012329447652244853694626639963780, 0.015010447346388952376697286041943, 0.017548967986243191099665352925900, 0.019938037786440888202278192730714, 0.022194935961012286796332102959499, 0.024339147126000805470360647041454, 0.026374505414839207241503786552615, 0.028286910788771200659968002987960, 0.030052581128092695322521110347341, 0.031646751371439929404586051078883, 0.033050413419978503290785944862689, 0.034255099704226061787082821046821, 0.035262412660156681033782717998428, 0.036076989622888701185500318003895, 0.036698604498456094498018047441094, 0.037120549269832576114119958413599, 0.037334228751935040321235449094698, 0.037361073762679023410321241766599}

Weights of the 87-point formula for abscissae x4