72#define GSL_FN_EVAL(F,x) (*((F)->function))(x,(F)->params)
99 double epsabs,
double epsrel,
size_t limit,
102 double *
result,
double *abserr);
107 double epsabs,
double epsrel,
size_t limit,
109 double *
result,
double * abserr) ;
113 double epsabs,
double epsrel,
size_t limit,
115 double *
result,
double *abserr) ;
120 double epsabs,
double epsrel,
size_t limit,
122 double *
result,
double *abserr) ;
127 double epsabs,
double epsrel,
size_t limit,
129 double *
result,
double *abserr) ;
158 RooRealVar maxSeg(
"maxSeg",
"maximum number of segments", 100);
159 RooCategory method(
"method",
"Integration method for each segment");
169 oocoutI(
nullptr,
Integration) <<
"RooAdaptiveGaussKronrodIntegrator1D has been registered " << std::endl;
189 _epsAbs(config.epsRel()),
190 _epsRel(config.epsAbs()),
211 _epsAbs(config.epsRel()),
212 _epsRel(config.epsAbs()),
275 coutE(
Integration) <<
"RooAdaptiveGaussKronrodIntegrator1D::setLimits: cannot override integrand's limits" << endl;
302 if (!infLo && !infHi) {
304 }
else if (infLo && infHi) {
306 }
else if (infLo && !infHi) {
401#define GSL_EBADTOL 13
403#define GSL_ERROR(a,b) oocoutE(nullptr,Integration) << "RooAdaptiveGaussKronrodIntegrator1D::integral() ERROR: " << a << endl ; return b ;
404#define GSL_DBL_MIN 2.2250738585072014e-308
405#define GSL_DBL_MAX 1.7976931348623157e+308
406#define GSL_DBL_EPSILON 2.2204460492503131e-16
409#define GSL_EMAXITER 3
412#define GSL_EDIVERGE 6
415#define GSL_ERROR_VAL(reason, gsl_errno, value) return value ;
417#define GSL_MAX(a,b) ((a) > (b) ? (a) : (b))
433#define GSL_COERCE_DBL(x) (gsl_coerce_double(x))
446 double *
result,
double *abserr,
447 double *defabs,
double *resabs);
450 double *
result,
double *abserr,
451 double *resabs,
double *resasc);
454 double *
result,
double *abserr,
455 double *resabs,
double *resasc);
458 double *
result,
double *abserr,
459 double *resabs,
double *resasc);
462 double *
result,
double *abserr,
463 double *resabs,
double *resasc);
466 double *
result,
double *abserr,
467 double *resabs,
double *resasc);
470 double *
result,
double *abserr,
471 double *resabs,
double *resasc);
474 double *cheb12,
double *cheb24);
492 const double wg[],
const double wgk[],
493 double fv1[],
double fv2[],
495 double *
result,
double * abserr,
496 double * resabs,
double * resasc);
501 double epsabs,
double epsrel,
size_t limit,
504 double *
result,
double *abserr);
515 workspace->
nrmax = 0;
519 workspace->
rlist[0] = 0.0;
520 workspace->
elist[0] = 0.0;
521 workspace->
order[0] = 0;
522 workspace->
level[0] = 0;
530 double result,
double error);
534 double result,
double error)
538 workspace->
elist[0] = error;
548 const size_t last = workspace->
size - 1;
549 const size_t limit = workspace->
limit;
551 double * elist = workspace->
elist;
552 size_t * order = workspace->
order;
558 size_t i_nrmax = workspace->
nrmax;
559 size_t i_maxerr = order[i_nrmax] ;
567 workspace->
i = i_maxerr ;
571 errmax = elist[i_maxerr] ;
578 while (i_nrmax > 0 && errmax > elist[order[i_nrmax - 1]])
580 order[i_nrmax] = order[i_nrmax - 1] ;
588 if(last < (limit/2 + 2))
594 top = limit - last + 1;
605 while (i < top && errmax < elist[order[i]])
607 order[i-1] = order[i] ;
611 order[i-1] = i_maxerr ;
615 errmin = elist[last] ;
619 while (k > i - 2 && errmin >= elist[order[k]])
621 order[k+1] = order[k] ;
629 i_maxerr = order[i_nrmax] ;
631 workspace->
i = i_maxerr ;
632 workspace->
nrmax = i_nrmax ;
640 double a1,
double b1,
double area1,
double error1,
641 double a2,
double b2,
double area2,
double error2);
645 double *
a,
double *
b,
double *
r,
double *
e);
651 double a1,
double b1,
double area1,
double error1,
652 double a2,
double b2,
double area2,
double error2)
654 double * alist = workspace->
alist ;
655 double * blist = workspace->
blist ;
656 double * rlist = workspace->
rlist ;
657 double * elist = workspace->
elist ;
658 size_t * level = workspace->
level ;
660 const size_t i_max = workspace->
i ;
661 const size_t i_new = workspace->
size ;
663 const size_t new_level = workspace->
level[i_max] + 1;
670 rlist[i_max] = area2;
671 elist[i_max] = error2;
672 level[i_max] = new_level;
676 rlist[i_new] = area1;
677 elist[i_new] = error1;
678 level[i_new] = new_level;
683 rlist[i_max] = area1;
684 elist[i_max] = error1;
685 level[i_max] = new_level;
689 rlist[i_new] = area2;
690 elist[i_new] = error2;
691 level[i_new] = new_level;
706 double *
a,
double *
b,
double *
r,
double *
e)
708 const size_t i = workspace->
i;
709 double * alist = workspace->
alist;
710 double * blist = workspace->
blist;
711 double * rlist = workspace->
rlist;
712 double * elist = workspace->
elist;
726 const double *
const rlist = workspace->
rlist ;
727 const size_t n = workspace->
size;
730 double result_sum = 0;
732 for (k = 0; k <
n; k++)
734 result_sum += rlist[k];
749 double tmp = (1 + 100 *
e) * (
fabs (a2) + 1000 * u);
751 int status =
fabs (a1) <= tmp &&
fabs (b2) <= tmp;
759 const double a,
const double b,
760 const double epsabs,
const double epsrel,
763 double *
result,
double * abserr,
769 double epsabs,
double epsrel,
size_t limit,
772 double *
result,
double * abserr)
808 status =
qag (
f,
a,
b, epsabs, epsrel, limit,
818 const double a,
const double b,
819 const double epsabs,
const double epsrel,
822 double *
result,
double *abserr,
826 double result0, abserr0, resabs0, resasc0;
828 size_t iteration = 0;
829 int roundoff_type1 = 0, roundoff_type2 = 0, error_type = 0;
840 if (limit > workspace->
limit)
845 if (epsabs <= 0 && (epsrel < 50 *
GSL_DBL_EPSILON || epsrel < 0.5e-28))
847 GSL_ERROR (
"tolerance cannot be acheived with given epsabs and epsrel",
853 q (
f,
a,
b, &result0, &abserr0, &resabs0, &resasc0);
865 if (abserr0 <= round_off && abserr0 > tolerance)
870 GSL_ERROR (
"cannot reach tolerance because of roundoff error "
873 else if ((abserr0 <= tolerance && abserr0 != resasc0) || abserr0 == 0.0)
895 double a1, b1, a2, b2;
896 double a_i, b_i, r_i, e_i;
897 double area1 = 0, area2 = 0, area12 = 0;
898 double error1 = 0, error2 = 0, error12 = 0;
899 double resasc1, resasc2;
900 double resabs1, resabs2;
904 retrieve (workspace, &a_i, &b_i, &r_i, &e_i);
907 b1 = 0.5 * (a_i + b_i);
911 q (
f, a1, b1, &area1, &error1, &resabs1, &resasc1);
912 q (
f, a2, b2, &area2, &error2, &resabs2, &resasc2);
914 area12 = area1 + area2;
915 error12 = error1 + error2;
917 errsum += (error12 - e_i);
918 area += area12 - r_i;
920 if (resasc1 != error1 && resasc2 != error2)
922 double delta = r_i - area12;
924 if (
fabs (delta) <= 1.0e-5 *
fabs (area12) && error12 >= 0.99 * e_i)
928 if (iteration >= 10 && error12 > e_i)
936 if (errsum > tolerance)
938 if (roundoff_type1 >= 6 || roundoff_type2 >= 20)
952 update (workspace, a1, b1, area1, error1, a2, b2, area2, error2);
954 retrieve (workspace, &a_i, &b_i, &r_i, &e_i);
959 while (iteration < limit && !error_type && errsum > tolerance);
964 if (errsum <= tolerance)
968 else if (error_type == 2)
970 GSL_ERROR (
"roundoff error prevents tolerance from being achieved",
973 else if (error_type == 3)
975 GSL_ERROR (
"bad integrand behavior found in the integration interval",
978 else if (iteration == limit)
988static double rescale_error (
double err,
const double result_abs,
const double result_asc) ;
995 if (result_asc != 0 && err != 0)
997 double scale =
TMath::Power((200 * err / result_asc), 1.5) ;
1001 err = result_asc * scale ;
1024 const double xgk[],
const double wg[],
const double wgk[],
1025 double fv1[],
double fv2[],
1027 double *
result,
double *abserr,
1028 double *resabs,
double *resasc)
1031 const double center = 0.5 * (
a +
b);
1032 const double half_length = 0.5 * (
b -
a);
1033 const double abs_half_length =
fabs (half_length);
1036 double result_gauss = 0;
1037 double result_kronrod = f_center * wgk[
n - 1];
1039 double result_abs =
fabs (result_kronrod);
1040 double result_asc = 0;
1041 double mean = 0, err = 0;
1047 result_gauss = f_center * wg[
n / 2 - 1];
1050 for (j = 0; j < (
n - 1) / 2; j++)
1052 const int jtw = j * 2 + 1;
1053 const double abscissa = half_length * xgk[jtw];
1054 const double fval1 =
GSL_FN_EVAL (
f, center - abscissa);
1055 const double fval2 =
GSL_FN_EVAL (
f, center + abscissa);
1056 const double fsum = fval1 + fval2;
1059 result_gauss += wg[j] * fsum;
1060 result_kronrod += wgk[jtw] * fsum;
1061 result_abs += wgk[jtw] * (
fabs (fval1) +
fabs (fval2));
1064 for (j = 0; j <
n / 2; j++)
1067 const double abscissa = half_length * xgk[jtwm1];
1068 const double fval1 =
GSL_FN_EVAL (
f, center - abscissa);
1069 const double fval2 =
GSL_FN_EVAL (
f, center + abscissa);
1072 result_kronrod += wgk[jtwm1] * (fval1 + fval2);
1073 result_abs += wgk[jtwm1] * (
fabs (fval1) +
fabs (fval2));
1076 mean = result_kronrod * 0.5;
1078 result_asc = wgk[
n - 1] *
fabs (f_center - mean);
1080 for (j = 0; j <
n - 1; j++)
1082 result_asc += wgk[j] * (
fabs (fv1[j] - mean) +
fabs (fv2[j] - mean));
1087 err = (result_kronrod - result_gauss) * half_length;
1089 result_kronrod *= half_length;
1090 result_abs *= abs_half_length;
1091 result_asc *= abs_half_length;
1093 *
result = result_kronrod;
1094 *resabs = result_abs;
1095 *resasc = result_asc;
1106 0.991455371120812639206854697526329,
1107 0.949107912342758524526189684047851,
1108 0.864864423359769072789712788640926,
1109 0.741531185599394439863864773280788,
1110 0.586087235467691130294144838258730,
1111 0.405845151377397166906606412076961,
1112 0.207784955007898467600689403773245,
1113 0.000000000000000000000000000000000
1121 0.129484966168869693270611432679082,
1122 0.279705391489276667901467771423780,
1123 0.381830050505118944950369775488975,
1124 0.417959183673469387755102040816327
1129 0.022935322010529224963732008058970,
1130 0.063092092629978553290700663189204,
1131 0.104790010322250183839876322541518,
1132 0.140653259715525918745189590510238,
1133 0.169004726639267902826583426598550,
1134 0.190350578064785409913256402421014,
1135 0.204432940075298892414161999234649,
1136 0.209482141084727828012999174891714
1141 double *
result,
double *abserr,
1142 double *resabs,
double *resasc)
1144 double fv1[8], fv2[8];
1146 gsl_integration_qk (8,
xgkA,
wgA,
wgkA, fv1, fv2,
f,
a,
b,
result, abserr, resabs, resasc);
1155 0.995657163025808080735527280689003,
1156 0.973906528517171720077964012084452,
1157 0.930157491355708226001207180059508,
1158 0.865063366688984510732096688423493,
1159 0.780817726586416897063717578345042,
1160 0.679409568299024406234327365114874,
1161 0.562757134668604683339000099272694,
1162 0.433395394129247190799265943165784,
1163 0.294392862701460198131126603103866,
1164 0.148874338981631210884826001129720,
1165 0.000000000000000000000000000000000
1173 0.066671344308688137593568809893332,
1174 0.149451349150580593145776339657697,
1175 0.219086362515982043995534934228163,
1176 0.269266719309996355091226921569469,
1177 0.295524224714752870173892994651338
1182 0.011694638867371874278064396062192,
1183 0.032558162307964727478818972459390,
1184 0.054755896574351996031381300244580,
1185 0.075039674810919952767043140916190,
1186 0.093125454583697605535065465083366,
1187 0.109387158802297641899210590325805,
1188 0.123491976262065851077958109831074,
1189 0.134709217311473325928054001771707,
1190 0.142775938577060080797094273138717,
1191 0.147739104901338491374841515972068,
1192 0.149445554002916905664936468389821
1198 double *
result,
double *abserr,
1199 double *resabs,
double *resasc)
1201 double fv1[11], fv2[11];
1203 gsl_integration_qk (11,
xgkB,
wgB,
wgkB, fv1, fv2,
f,
a,
b,
result, abserr, resabs, resasc);
1212 0.998002298693397060285172840152271,
1213 0.987992518020485428489565718586613,
1214 0.967739075679139134257347978784337,
1215 0.937273392400705904307758947710209,
1216 0.897264532344081900882509656454496,
1217 0.848206583410427216200648320774217,
1218 0.790418501442465932967649294817947,
1219 0.724417731360170047416186054613938,
1220 0.650996741297416970533735895313275,
1221 0.570972172608538847537226737253911,
1222 0.485081863640239680693655740232351,
1223 0.394151347077563369897207370981045,
1224 0.299180007153168812166780024266389,
1225 0.201194093997434522300628303394596,
1226 0.101142066918717499027074231447392,
1227 0.000000000000000000000000000000000
1235 0.030753241996117268354628393577204,
1236 0.070366047488108124709267416450667,
1237 0.107159220467171935011869546685869,
1238 0.139570677926154314447804794511028,
1239 0.166269205816993933553200860481209,
1240 0.186161000015562211026800561866423,
1241 0.198431485327111576456118326443839,
1242 0.202578241925561272880620199967519
1247 0.005377479872923348987792051430128,
1248 0.015007947329316122538374763075807,
1249 0.025460847326715320186874001019653,
1250 0.035346360791375846222037948478360,
1251 0.044589751324764876608227299373280,
1252 0.053481524690928087265343147239430,
1253 0.062009567800670640285139230960803,
1254 0.069854121318728258709520077099147,
1255 0.076849680757720378894432777482659,
1256 0.083080502823133021038289247286104,
1257 0.088564443056211770647275443693774,
1258 0.093126598170825321225486872747346,
1259 0.096642726983623678505179907627589,
1260 0.099173598721791959332393173484603,
1261 0.100769845523875595044946662617570,
1262 0.101330007014791549017374792767493
1267 double *
result,
double *abserr,
1268 double *resabs,
double *resasc)
1270 double fv1[16], fv2[16];
1272 gsl_integration_qk (16,
xgkC,
wgC,
wgkC, fv1, fv2,
f,
a,
b,
result, abserr, resabs, resasc);
1281 0.998859031588277663838315576545863,
1282 0.993128599185094924786122388471320,
1283 0.981507877450250259193342994720217,
1284 0.963971927277913791267666131197277,
1285 0.940822633831754753519982722212443,
1286 0.912234428251325905867752441203298,
1287 0.878276811252281976077442995113078,
1288 0.839116971822218823394529061701521,
1289 0.795041428837551198350638833272788,
1290 0.746331906460150792614305070355642,
1291 0.693237656334751384805490711845932,
1292 0.636053680726515025452836696226286,
1293 0.575140446819710315342946036586425,
1294 0.510867001950827098004364050955251,
1295 0.443593175238725103199992213492640,
1296 0.373706088715419560672548177024927,
1297 0.301627868114913004320555356858592,
1298 0.227785851141645078080496195368575,
1299 0.152605465240922675505220241022678,
1300 0.076526521133497333754640409398838,
1301 0.000000000000000000000000000000000
1309 0.017614007139152118311861962351853,
1310 0.040601429800386941331039952274932,
1311 0.062672048334109063569506535187042,
1312 0.083276741576704748724758143222046,
1313 0.101930119817240435036750135480350,
1314 0.118194531961518417312377377711382,
1315 0.131688638449176626898494499748163,
1316 0.142096109318382051329298325067165,
1317 0.149172986472603746787828737001969,
1318 0.152753387130725850698084331955098
1323 0.003073583718520531501218293246031,
1324 0.008600269855642942198661787950102,
1325 0.014626169256971252983787960308868,
1326 0.020388373461266523598010231432755,
1327 0.025882133604951158834505067096153,
1328 0.031287306777032798958543119323801,
1329 0.036600169758200798030557240707211,
1330 0.041668873327973686263788305936895,
1331 0.046434821867497674720231880926108,
1332 0.050944573923728691932707670050345,
1333 0.055195105348285994744832372419777,
1334 0.059111400880639572374967220648594,
1335 0.062653237554781168025870122174255,
1336 0.065834597133618422111563556969398,
1337 0.068648672928521619345623411885368,
1338 0.071054423553444068305790361723210,
1339 0.073030690332786667495189417658913,
1340 0.074582875400499188986581418362488,
1341 0.075704497684556674659542775376617,
1342 0.076377867672080736705502835038061,
1343 0.076600711917999656445049901530102
1348 double *
result,
double *abserr,
1349 double *resabs,
double *resasc)
1351 double fv1[21], fv2[21];
1353 gsl_integration_qk (21,
xgkD,
wgD,
wgkD, fv1, fv2,
f,
a,
b,
result, abserr, resabs, resasc);
1362 0.999262104992609834193457486540341,
1363 0.995556969790498097908784946893902,
1364 0.988035794534077247637331014577406,
1365 0.976663921459517511498315386479594,
1366 0.961614986425842512418130033660167,
1367 0.942974571228974339414011169658471,
1368 0.920747115281701561746346084546331,
1369 0.894991997878275368851042006782805,
1370 0.865847065293275595448996969588340,
1371 0.833442628760834001421021108693570,
1372 0.797873797998500059410410904994307,
1373 0.759259263037357630577282865204361,
1374 0.717766406813084388186654079773298,
1375 0.673566368473468364485120633247622,
1376 0.626810099010317412788122681624518,
1377 0.577662930241222967723689841612654,
1378 0.526325284334719182599623778158010,
1379 0.473002731445714960522182115009192,
1380 0.417885382193037748851814394594572,
1381 0.361172305809387837735821730127641,
1382 0.303089538931107830167478909980339,
1383 0.243866883720988432045190362797452,
1384 0.183718939421048892015969888759528,
1385 0.122864692610710396387359818808037,
1386 0.061544483005685078886546392366797,
1387 0.000000000000000000000000000000000
1395 0.011393798501026287947902964113235,
1396 0.026354986615032137261901815295299,
1397 0.040939156701306312655623487711646,
1398 0.054904695975835191925936891540473,
1399 0.068038333812356917207187185656708,
1400 0.080140700335001018013234959669111,
1401 0.091028261982963649811497220702892,
1402 0.100535949067050644202206890392686,
1403 0.108519624474263653116093957050117,
1404 0.114858259145711648339325545869556,
1405 0.119455763535784772228178126512901,
1406 0.122242442990310041688959518945852,
1407 0.123176053726715451203902873079050
1412 0.001987383892330315926507851882843,
1413 0.005561932135356713758040236901066,
1414 0.009473973386174151607207710523655,
1415 0.013236229195571674813656405846976,
1416 0.016847817709128298231516667536336,
1417 0.020435371145882835456568292235939,
1418 0.024009945606953216220092489164881,
1419 0.027475317587851737802948455517811,
1420 0.030792300167387488891109020215229,
1421 0.034002130274329337836748795229551,
1422 0.037116271483415543560330625367620,
1423 0.040083825504032382074839284467076,
1424 0.042872845020170049476895792439495,
1425 0.045502913049921788909870584752660,
1426 0.047982537138836713906392255756915,
1427 0.050277679080715671963325259433440,
1428 0.052362885806407475864366712137873,
1429 0.054251129888545490144543370459876,
1430 0.055950811220412317308240686382747,
1431 0.057437116361567832853582693939506,
1432 0.058689680022394207961974175856788,
1433 0.059720340324174059979099291932562,
1434 0.060539455376045862945360267517565,
1435 0.061128509717053048305859030416293,
1436 0.061471189871425316661544131965264,
1437 0.061580818067832935078759824240066
1444 double *
result,
double *abserr,
1445 double *resabs,
double *resasc)
1447 double fv1[26], fv2[26];
1449 gsl_integration_qk (26,
xgkE,
wgE,
wgkE, fv1, fv2,
f,
a,
b,
result, abserr, resabs, resasc);
1458 0.999484410050490637571325895705811,
1459 0.996893484074649540271630050918695,
1460 0.991630996870404594858628366109486,
1461 0.983668123279747209970032581605663,
1462 0.973116322501126268374693868423707,
1463 0.960021864968307512216871025581798,
1464 0.944374444748559979415831324037439,
1465 0.926200047429274325879324277080474,
1466 0.905573307699907798546522558925958,
1467 0.882560535792052681543116462530226,
1468 0.857205233546061098958658510658944,
1469 0.829565762382768397442898119732502,
1470 0.799727835821839083013668942322683,
1471 0.767777432104826194917977340974503,
1472 0.733790062453226804726171131369528,
1473 0.697850494793315796932292388026640,
1474 0.660061064126626961370053668149271,
1475 0.620526182989242861140477556431189,
1476 0.579345235826361691756024932172540,
1477 0.536624148142019899264169793311073,
1478 0.492480467861778574993693061207709,
1479 0.447033769538089176780609900322854,
1480 0.400401254830394392535476211542661,
1481 0.352704725530878113471037207089374,
1482 0.304073202273625077372677107199257,
1483 0.254636926167889846439805129817805,
1484 0.204525116682309891438957671002025,
1485 0.153869913608583546963794672743256,
1486 0.102806937966737030147096751318001,
1487 0.051471842555317695833025213166723,
1488 0.000000000000000000000000000000000
1496 0.007968192496166605615465883474674,
1497 0.018466468311090959142302131912047,
1498 0.028784707883323369349719179611292,
1499 0.038799192569627049596801936446348,
1500 0.048402672830594052902938140422808,
1501 0.057493156217619066481721689402056,
1502 0.065974229882180495128128515115962,
1503 0.073755974737705206268243850022191,
1504 0.080755895229420215354694938460530,
1505 0.086899787201082979802387530715126,
1506 0.092122522237786128717632707087619,
1507 0.096368737174644259639468626351810,
1508 0.099593420586795267062780282103569,
1509 0.101762389748405504596428952168554,
1510 0.102852652893558840341285636705415
1515 0.001389013698677007624551591226760,
1516 0.003890461127099884051267201844516,
1517 0.006630703915931292173319826369750,
1518 0.009273279659517763428441146892024,
1519 0.011823015253496341742232898853251,
1520 0.014369729507045804812451432443580,
1521 0.016920889189053272627572289420322,
1522 0.019414141193942381173408951050128,
1523 0.021828035821609192297167485738339,
1524 0.024191162078080601365686370725232,
1525 0.026509954882333101610601709335075,
1526 0.028754048765041292843978785354334,
1527 0.030907257562387762472884252943092,
1528 0.032981447057483726031814191016854,
1529 0.034979338028060024137499670731468,
1530 0.036882364651821229223911065617136,
1531 0.038678945624727592950348651532281,
1532 0.040374538951535959111995279752468,
1533 0.041969810215164246147147541285970,
1534 0.043452539701356069316831728117073,
1535 0.044814800133162663192355551616723,
1536 0.046059238271006988116271735559374,
1537 0.047185546569299153945261478181099,
1538 0.048185861757087129140779492298305,
1539 0.049055434555029778887528165367238,
1540 0.049795683427074206357811569379942,
1541 0.050405921402782346840893085653585,
1542 0.050881795898749606492297473049805,
1543 0.051221547849258772170656282604944,
1544 0.051426128537459025933862879215781,
1545 0.051494729429451567558340433647099
1550 double *
result,
double *abserr,
1551 double *resabs,
double *resasc)
1553 double fv1[31], fv2[31];
1555 gsl_integration_qk (31,
xgkF,
wgF,
wgkF, fv1, fv2,
f,
a,
b,
result, abserr, resabs, resasc);
1565 GSL_ERROR_VAL (
"workspace length n must be positive integer",
1574 GSL_ERROR_VAL (
"failed to allocate space for workspace struct",
1578 w->alist = (
double *)
malloc (
n *
sizeof (
double));
1588 w->blist = (
double *)
malloc (
n *
sizeof (
double));
1599 w->rlist = (
double *)
malloc (
n *
sizeof (
double));
1612 w->elist = (
double *)
malloc (
n *
sizeof (
double));
1625 w->order = (
size_t *)
malloc (
n *
sizeof (
size_t));
1639 w->level = (
size_t *)
malloc (
n *
sizeof (
size_t));
1656 w->maximum_level = 0 ;
1682 workspace->
nrmax = 0;
1683 workspace->
i = workspace->
order[0] ;
1699 int id = workspace->
nrmax;
1702 const size_t * level = workspace->
level;
1703 const size_t * order = workspace->
order;
1705 size_t limit = workspace->
limit ;
1706 size_t last = workspace->
size - 1 ;
1708 if (last > (1 + limit / 2))
1710 jupbnd = limit + 1 - last;
1717 for (k =
id; k <= jupbnd; k++)
1719 size_t i_max = order[workspace->
nrmax];
1721 workspace->
i = i_max ;
1737 size_t i = workspace->
i ;
1738 const size_t * level = workspace->
level;
1801 double *epstab = table->
rlist2;
1802 double *res3la = table->
res3la;
1803 const size_t n = table->
n - 1;
1805 const double current = epstab[
n];
1810 const size_t newelm =
n / 2;
1811 const size_t n_orig =
n;
1815 const size_t nres_orig = table->
nres;
1827 epstab[
n + 2] = epstab[
n];
1830 for (i = 0; i < newelm; i++)
1832 double res = epstab[
n - 2 * i + 2];
1833 double e0 = epstab[
n - 2 * i - 2];
1834 double e1 = epstab[
n - 2 * i - 1];
1837 double e1abs =
fabs (e1);
1838 double delta2 = e2 - e1;
1839 double err2 =
fabs (delta2);
1841 double delta3 = e1 - e0;
1842 double err3 =
fabs (delta3);
1845 double e3, delta1, err1, tol1, ss;
1847 if (err2 <= tol2 && err3 <= tol3)
1853 absolute = err2 + err3;
1859 e3 = epstab[
n - 2 * i];
1860 epstab[
n - 2 * i] = e1;
1862 err1 =
fabs (delta1);
1868 if (err1 <= tol1 || err2 <= tol2 || err3 <= tol3)
1874 ss = (1 / delta1 + 1 / delta2) - 1 / delta3;
1880 if (
fabs (ss * e1) <= 0.0001)
1890 epstab[
n - 2 * i] = res;
1893 const double error = err2 +
fabs (res - e2) + err3;
1895 if (error <= *abserr)
1906 const size_t limexp = 50 - 1;
1908 if (n_final == limexp)
1910 n_final = 2 * (limexp / 2);
1914 if (n_orig % 2 == 1)
1916 for (i = 0; i <= newelm; i++)
1918 epstab[1 + i * 2] = epstab[i * 2 + 3];
1923 for (i = 0; i <= newelm; i++)
1925 epstab[i * 2] = epstab[i * 2 + 2];
1929 if (n_orig != n_final)
1931 for (i = 0; i <= n_final; i++)
1933 epstab[i] = epstab[n_orig - n_final + i];
1937 table->
n = n_final + 1;
1941 res3la[nres_orig] = *
result;
1949 res3la[0] = res3la[1];
1950 res3la[1] = res3la[2];
1960 table->
nres = nres_orig + 1;
1984 b,
const double epsabs,
const double epsrel,
const size_t limit,
1991 double epsabs,
double epsrel,
size_t limit,
1993 double *
result,
double * abserr)
1995 int status =
qags (
f,
a,
b, epsabs, epsrel, limit,
2009static double i_transform (
double t,
void *params);
2013 double epsabs,
double epsrel,
size_t limit,
2015 double *
result,
double *abserr)
2024 status =
qags (&f_transform, 0.0, 1.0,
2025 epsabs, epsrel, limit,
2037 double x = (1 - t) / t;
2057 double epsabs,
double epsrel,
size_t limit,
2059 double *
result,
double *abserr)
2066 transform_params.
b =
b ;
2067 transform_params.
f =
f ;
2070 f_transform.
params = &transform_params;
2072 status =
qags (&f_transform, 0.0, 1.0,
2073 epsabs, epsrel, limit,
2087 double x =
b - (1 - t) / t;
2106 double epsabs,
double epsrel,
size_t limit,
2108 double *
result,
double *abserr)
2115 transform_params.
a =
a ;
2116 transform_params.
f =
f ;
2119 f_transform.
params = &transform_params;
2121 status =
qags (&f_transform, 0.0, 1.0,
2122 epsabs, epsrel, limit,
2136 double x =
a + (1 - t) / t;
2145 const double a,
const double b,
2146 const double epsabs,
const double epsrel,
2149 double *
result,
double *abserr,
2152 double area, errsum;
2153 double res_ext, err_ext;
2154 double result0, abserr0, resabs0, resasc0;
2158 double error_over_large_intervals = 0;
2159 double reseps = 0, abseps = 0, correc = 0;
2161 int roundoff_type1 = 0, roundoff_type2 = 0, roundoff_type3 = 0;
2162 int error_type = 0, error_type2 = 0;
2164 size_t iteration = 0;
2166 int positive_integrand = 0;
2167 int extrapolate = 0;
2168 int disallow_extrapolation = 0;
2179 if (limit > workspace->
limit)
2186 if (epsabs <= 0 && (epsrel < 50 *
GSL_DBL_EPSILON || epsrel < 0.5e-28))
2188 GSL_ERROR (
"tolerance cannot be acheived with given epsabs and epsrel",
2194 q (
f,
a,
b, &result0, &abserr0, &resabs0, &resasc0);
2200 if (abserr0 <= 100 * GSL_DBL_EPSILON * resabs0 && abserr0 > tolerance)
2205 GSL_ERROR (
"cannot reach tolerance because of roundoff error"
2208 else if ((abserr0 <= tolerance && abserr0 != resasc0) || abserr0 == 0.0)
2215 else if (limit == 1)
2240 size_t current_level;
2241 double a1, b1, a2, b2;
2242 double a_i, b_i, r_i, e_i;
2243 double area1 = 0, area2 = 0, area12 = 0;
2244 double error1 = 0, error2 = 0, error12 = 0;
2245 double resasc1, resasc2;
2246 double resabs1, resabs2;
2251 retrieve (workspace, &a_i, &b_i, &r_i, &e_i);
2253 current_level = workspace->
level[workspace->
i] + 1;
2256 b1 = 0.5 * (a_i + b_i);
2262 q (
f, a1, b1, &area1, &error1, &resabs1, &resasc1);
2263 q (
f, a2, b2, &area2, &error2, &resabs2, &resasc2);
2265 area12 = area1 + area2;
2266 error12 = error1 + error2;
2276 errsum = errsum + error12 - e_i;
2277 area = area + area12 - r_i;
2281 if (resasc1 != error1 && resasc2 != error2)
2283 double delta = r_i - area12;
2285 if (
fabs (delta) <= 1.0e-5 *
fabs (area12) && error12 >= 0.99 * e_i)
2296 if (iteration > 10 && error12 > e_i)
2304 if (roundoff_type1 + roundoff_type2 >= 10 || roundoff_type3 >= 20)
2309 if (roundoff_type2 >= 5)
2324 update (workspace, a1, b1, area1, error1, a2, b2, area2, error2);
2326 if (errsum <= tolerance)
2328 goto compute_result;
2336 if (iteration >= limit - 1)
2344 error_over_large_intervals = errsum;
2350 if (disallow_extrapolation)
2355 error_over_large_intervals += -last_e_i;
2357 if (current_level < workspace->maximum_level)
2359 error_over_large_intervals += error12;
2371 workspace->
nrmax = 1;
2374 if (!error_type2 && error_over_large_intervals > ertest)
2384 qelg (&table, &reseps, &abseps);
2388 if (ktmin > 5 && err_ext < 0.001 * errsum)
2393 if (abseps < err_ext)
2398 correc = error_over_large_intervals;
2400 if (err_ext <= ertest)
2408 disallow_extrapolation = 1;
2411 if (error_type == 5)
2420 error_over_large_intervals = errsum;
2423 while (iteration < limit);
2429 goto compute_result;
2431 if (error_type || error_type2)
2441 if (res_ext != 0.0 && area != 0.0)
2443 if (err_ext /
fabs (res_ext) > errsum /
fabs (area))
2444 goto compute_result;
2446 else if (err_ext > errsum)
2448 goto compute_result;
2450 else if (area == 0.0)
2461 if (!positive_integrand && max_area < 0.01 * resabs0)
2466 double ratio = res_ext / area;
2468 if (ratio < 0.01 || ratio > 100.0 || errsum >
fabs (area))
2486 if (error_type == 0)
2490 else if (error_type == 1)
2494 else if (error_type == 2)
2496 GSL_ERROR (
"cannot reach tolerance because of roundoff error",
2499 else if (error_type == 3)
2501 GSL_ERROR (
"bad integrand behavior found in the integration interval",
2504 else if (error_type == 4)
2506 GSL_ERROR (
"roundoff error detected in the extrapolation table",
2509 else if (error_type == 5)
2511 GSL_ERROR (
"integral is divergent, or slowly convergent",
static void update(gsl_integration_workspace *workspace, double a1, double b1, double area1, double error1, double a2, double b2, double area2, double error2)
void gsl_integration_qk61(const gsl_function *f, double a, double b, double *result, double *abserr, double *resabs, double *resasc)
gsl_integration_workspace * gsl_integration_workspace_alloc(const size_t n)
static int increase_nrmax(gsl_integration_workspace *workspace)
static double rescale_error(double err, const double result_abs, const double result_asc)
void gsl_integration_qk51(const gsl_function *f, double a, double b, double *result, double *abserr, double *resabs, double *resasc)
static void retrieve(const gsl_integration_workspace *workspace, double *a, double *b, double *r, double *e)
static void set_initial_result(gsl_integration_workspace *workspace, double result, double error)
static const double wgD[11]
static void qelg(struct extrapolation_table *table, double *result, double *abserr)
static double il_transform(double t, void *params)
void gsl_integration_qk(const int n, const double xgk[], const double wg[], const double wgk[], double fv1[], double fv2[], const gsl_function *f, double a, double b, double *result, double *abserr, double *resabs, double *resasc)
static void initialise_table(struct extrapolation_table *table)
static const double wgC[8]
static int test_positivity(double result, double resabs)
static const double wgkE[26]
static const double wgkB[11]
static int large_interval(gsl_integration_workspace *workspace)
static const double wgkF[31]
static void qpsrt(gsl_integration_workspace *workspace)
#define GSL_ERROR_VAL(reason, gsl_errno, value)
double GSL_MAX_DBL(double a, double b)
static const double wgB[5]
double RooAdaptiveGaussKronrodIntegrator1D_GSL_GlueFunction(double x, void *data)
Glue function interacing to GSL code.
static int qag(const gsl_function *f, const double a, const double b, const double epsabs, const double epsrel, const size_t limit, gsl_integration_workspace *workspace, double *result, double *abserr, gsl_integration_rule *q)
static const double xgkE[26]
int gsl_integration_qagil(gsl_function *f, double b, double epsabs, double epsrel, size_t limit, gsl_integration_workspace *workspace, double *result, double *abserr)
int gsl_integration_qag(const gsl_function *f, double a, double b, double epsabs, double epsrel, size_t limit, int key, gsl_integration_workspace *workspace, double *result, double *abserr)
static void reset_nrmax(gsl_integration_workspace *workspace)
#define GSL_COERCE_DBL(x)
static const double wgkC[16]
static const double xgkA[8]
static const double wgF[15]
void gsl_integration_qk21(const gsl_function *f, double a, double b, double *result, double *abserr, double *resabs, double *resasc)
int gsl_integration_qagi(gsl_function *f, double epsabs, double epsrel, size_t limit, gsl_integration_workspace *workspace, double *result, double *abserr)
void gsl_integration_rule(const gsl_function *f, double a, double b, double *result, double *abserr, double *defabs, double *resabs)
static int qags(const gsl_function *f, const double a, const double b, const double epsabs, const double epsrel, const size_t limit, gsl_integration_workspace *workspace, double *result, double *abserr, gsl_integration_rule *q)
static double sum_results(const gsl_integration_workspace *workspace)
void gsl_integration_qk15(const gsl_function *f, double a, double b, double *result, double *abserr, double *resabs, double *resasc)
static int subinterval_too_small(double a1, double a2, double b2)
int gsl_integration_qagiu(gsl_function *f, double a, double epsabs, double epsrel, size_t limit, gsl_integration_workspace *workspace, double *result, double *abserr)
static const double xgkD[21]
static const double wgA[4]
void gsl_integration_qk41(const gsl_function *f, double a, double b, double *result, double *abserr, double *resabs, double *resasc)
void gsl_integration_qcheb(gsl_function *f, double a, double b, double *cheb12, double *cheb24)
static const double wgkD[21]
void gsl_integration_workspace_free(gsl_integration_workspace *w)
int gsl_integration_qags(const gsl_function *f, double a, double b, double epsabs, double epsrel, size_t limit, gsl_integration_workspace *workspace, double *result, double *abserr)
#define GSL_FN_EVAL(F, x)
static const double xgkF[31]
static const double wgE[13]
static double i_transform(double t, void *params)
double gsl_coerce_double(const double x)
void gsl_integration_qk31(const gsl_function *f, double a, double b, double *result, double *abserr, double *resabs, double *resasc)
static double iu_transform(double t, void *params)
static const double xgkC[16]
static void append_table(struct extrapolation_table *table, double y)
static const double xgkB[11]
static const double wgkA[8]
static void initialise(gsl_integration_workspace *workspace, double a, double b)
winID h TVirtualViewer3D TVirtualGLPainter p
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t b
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t r
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t result
Int_t getCatIndex(const char *name, Int_t defVal=0, bool verbose=false) const
Get index value of a RooAbsCategory stored in set with given name.
double getRealValue(const char *name, double defVal=0, bool verbose=false) const
Get value of a RooAbsReal stored in set with given name.
Abstract interface for evaluating a real-valued function of one real variable and performing numerica...
virtual double getMaxLimit(UInt_t dimension) const =0
virtual double getMinLimit(UInt_t dimension) const =0
UInt_t getDimension() const
RooAbsIntegrator is the abstract interface for integrators of real-valued functions that implement th...
bool isValid() const
Is integrator in valid state.
const RooAbsFunc * _function
Pointer to function binding of integrand.
const RooAbsFunc * integrand() const
Return integrand function binding.
bool _valid
Is integrator in valid state?
RooAdaptiveGaussKronrodIntegrator1D implements the Gauss-Kronrod integration algorithm.
double integral(const double *yvec=0) override
Calculate and return integral at at given parameter values.
RooAbsIntegrator * clone(const RooAbsFunc &function, const RooNumIntConfig &config) const override
Virtual constructor.
bool initialize()
Initialize integrator allocate buffers and setup GSL workspace.
double _epsAbs
Current coordinate.
TClass * IsA() const override
bool setLimits(double *xmin, double *xmax) override
Change our integration limits.
RooAdaptiveGaussKronrodIntegrator1D()
coverity[UNINIT_CTOR] Default constructor
bool checkLimits() const override
Check that our integration range is finite and otherwise return false.
friend double RooAdaptiveGaussKronrodIntegrator1D_GSL_GlueFunction(double x, void *data)
Glue function interacing to GSL code.
~RooAdaptiveGaussKronrodIntegrator1D() override
Destructor.
double _xmax
Lower integration bound.
static void registerIntegrator(RooNumIntFactory &fact)
Register this class with RooNumIntConfig as a possible choice of numeric integrator for one-dimension...
RooArgSet is a container object that can hold multiple RooAbsArg objects.
RooCategory is an object to represent discrete states.
bool setIndex(Int_t index, bool printError=true) override
Set value by specifying the index code of the desired state.
bool defineType(const std::string &label)
Define a state with given name.
RooNumIntConfig holds the configuration parameters of the various numeric integrators used by RooReal...
const RooArgSet & getConfigSection(const char *name) const
Retrieve configuration information specific to integrator with given name.
RooNumIntFactory is a factory to instantiate numeric integrators from a given function binding and a ...
bool storeProtoIntegrator(RooAbsIntegrator *proto, const RooArgSet &defConfig, const char *depName="")
Method accepting registration of a prototype numeric integrator along with a RooArgSet of its default...
static RooNumIntFactory & instance()
Static method returning reference to singleton instance of factory.
static Int_t isInfinite(double x)
Return true if x is infinite by RooNumBer internal specification.
RooRealVar represents a variable that can be changed from the outside.
virtual const char * GetName() const
Returns name of object.
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
void function(const Char_t *name_, T fun, const Char_t *docstring=0)
static Roo_reg_AGKInteg1D instance
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
static void registerIntegrator()
double(* function)(double x, void *params)