7 const long *nerr,
const long *level,
long)
12 fprintf(
ioQQQ,
"Error [%ld]: %*s\n",*nerr,
int(*nmesg),mesg);
17 fprintf(
ioQQQ,
"Warning [%ld]: %*s\n",*nerr,
int(*nmesg),mesg);
37 return FLT_EPSILON/FLT_RADIX;
41 return log10(
double(FLT_RADIX));
43 fprintf(stderr,
"Error in input to r1mach");
63 return DBL_EPSILON/FLT_RADIX;
67 return log10(
double(FLT_RADIX));
75 const static long c__4 = 4;
76 const static long c__1 = 1;
80 const static long c__2 = 2;
81 const static long c__0 = 0;
89 static void dgtsl_(
const long *,
double*,
double*,
double*,
92 void dqage_(
const D_fp& f,
const double *a,
const double *b,
const double *epsabs,
93 const double *epsrel,
const long *key,
const long *limit,
94 double *result,
double *abserr,
long *neval,
long *ier,
double *
95 alist__,
double *blist,
double *rlist,
double *elist,
96 long *iord,
long *last)
106 double area1, area2, area12, erro12, defab1, defab2;
110 double error1, error2, defabs, epmach, errbnd, resabs, errmax;
317 if (*epsabs <= 0. && *epsrel <
max(d__1,5e-29)) {
336 dqk15_(f, a, b, result, abserr, &defabs, &resabs);
339 dqk21_(f, a, b, result, abserr, &defabs, &resabs);
342 dqk31_(f, a, b, result, abserr, &defabs, &resabs);
345 dqk41_(f, a, b, result, abserr, &defabs, &resabs);
348 dqk51_(f, a, b, result, abserr, &defabs, &resabs);
351 dqk61_(f, a, b, result, abserr, &defabs, &resabs);
361 d__1 = *epsabs, d__2 = *epsrel * fabs(*result);
362 errbnd =
max(d__1,d__2);
363 if (*abserr <= epmach * 50. * defabs && *abserr > errbnd) {
369 if (*ier != 0 || ( *abserr <= errbnd && *abserr != resabs ) || *abserr == 0.)
390 for (*last = 2; *last <= i__1; ++(*last)) {
394 a1 = alist__[maxerr];
395 b1 = (alist__[maxerr] + blist[maxerr]) * .5;
399 dqk15_(f, &a1, &b1, &area1, &error1, &resabs, &defab1);
402 dqk21_(f, &a1, &b1, &area1, &error1, &resabs, &defab1);
405 dqk31_(f, &a1, &b1, &area1, &error1, &resabs, &defab1);
408 dqk41_(f, &a1, &b1, &area1, &error1, &resabs, &defab1);
411 dqk51_(f, &a1, &b1, &area1, &error1, &resabs, &defab1);
414 dqk61_(f, &a1, &b1, &area1, &error1, &resabs, &defab1);
417 dqk15_(f, &a2, &b2, &area2, &error2, &resabs, &defab2);
420 dqk21_(f, &a2, &b2, &area2, &error2, &resabs, &defab2);
423 dqk31_(f, &a2, &b2, &area2, &error2, &resabs, &defab2);
426 dqk41_(f, &a2, &b2, &area2, &error2, &resabs, &defab2);
429 dqk51_(f, &a2, &b2, &area2, &error2, &resabs, &defab2);
432 dqk61_(f, &a2, &b2, &area2, &error2, &resabs, &defab2);
439 area12 = area1 + area2;
440 erro12 = error1 + error2;
441 errsum = errsum + erro12 - errmax;
442 area = area + area12 - rlist[maxerr];
443 if (defab1 == error1 || defab2 == error2) {
446 if ((d__1 = rlist[maxerr] - area12, fabs(d__1)) <= fabs(area12) * 1e-5
447 && erro12 >= errmax * .99) {
450 if (*last > 10 && erro12 > errmax) {
454 rlist[maxerr] = area1;
455 rlist[*last] = area2;
457 d__1 = *epsabs, d__2 = *epsrel * fabs(area);
458 errbnd =
max(d__1,d__2);
459 if (errsum <= errbnd) {
465 if (iroff1 >= 6 || iroff2 >= 20) {
472 if (*last == *limit) {
480 d__1 = fabs(a1), d__2 = fabs(b2);
481 if (
max(d__1,d__2) <= (epmach * 100. + 1.) * (fabs(a2) + uflow * 1e3))
489 if (error2 > error1) {
495 elist[maxerr] = error1;
496 elist[*last] = error2;
499 alist__[maxerr] =
a2;
502 rlist[maxerr] = area2;
503 rlist[*last] = area1;
504 elist[maxerr] = error2;
505 elist[*last] = error1;
512 dqpsrt_(limit, last, &maxerr, &errmax, &elist[1], &iord[1], &nrmax);
514 if (*ier != 0 || errsum <= errbnd) {
526 for (k = 1; k <= i__1; ++k) {
533 *neval = (keyf * 10 + 1) * ((*neval << 1) + 1);
536 *neval = *neval * 30 + 15;
542 void dqag_(
const D_fp& f,
const double *a,
const double *b,
const double *epsabs,
543 const double *epsrel,
const long *key,
double *result,
544 double *abserr,
long *neval,
long *ier,
long *limit,
545 const long *lenw,
long *last,
long *iwork,
double *work)
547 long l1, l2, l3, lvl;
712 if (*limit < 1 || *lenw < *limit << 2) {
722 dqage_(f, a, b, epsabs, epsrel, key, limit, result, abserr, neval,
723 ier, &work[1], &work[l1], &work[l2], &work[l3], &iwork[1], last);
733 xerror_(
"abnormal return from dqag ", &
c__26, ier, &lvl, 26);
739 const double *epsabs,
const double *epsrel,
const long *limit,
740 double *result,
double *abserr,
long *neval,
long *ier,
741 double *alist__,
double *blist,
double *rlist,
double *elist,
742 long *iord,
long *last)
756 double area1, area2, area12;
757 double small=0., erro12;
759 double defab1, defab2;
763 long iroff1, iroff2, iroff3;
764 double res3la[3], error1, error2, rlist2[52];
766 double defabs, epmach, erlarg=0., abseps, correc=0., errbnd, resabs;
768 double erlast, errmax;
772 double ertest=0., errsum;
1002 d__1 = epmach * 50.;
1003 if (*epsabs <= 0. && *epsrel <
max(d__1,5e-29)) {
1032 dres = fabs(*result);
1034 d__1 = *epsabs, d__2 = *epsrel * dres;
1035 errbnd =
max(d__1,d__2);
1036 if (*abserr <= epmach * 100. * defabs && *abserr > errbnd) {
1042 if (*ier != 0 || ( *abserr <= errbnd && *abserr != resabs ) || *abserr == 0.)
1052 rlist2[0] = *result;
1069 if (dres >= (1. - epmach * 50.) * defabs) {
1077 for (*last = 2; *last <= i__1; ++(*last)) {
1081 a1 = alist__[maxerr];
1082 b1 = (alist__[maxerr] + blist[maxerr]) * .5;
1086 dqk15i_(f, &boun, inf, &a1, &b1, &area1, &error1, &resabs, &
1088 dqk15i_(f, &boun, inf, &a2, &b2, &area2, &error2, &resabs, &
1094 area12 = area1 + area2;
1095 erro12 = error1 + error2;
1096 errsum = errsum + erro12 - errmax;
1097 area = area + area12 - rlist[maxerr];
1098 if (defab1 == error1 || defab2 == error2) {
1101 if ((d__1 = rlist[maxerr] - area12, fabs(d__1)) > fabs(area12) * 1e-5 ||
1102 erro12 < errmax * .99) {
1112 if (*last > 10 && erro12 > errmax) {
1116 rlist[maxerr] = area1;
1117 rlist[*last] = area2;
1119 d__1 = *epsabs, d__2 = *epsrel * fabs(area);
1120 errbnd =
max(d__1,d__2);
1124 if (iroff1 + iroff2 >= 10 || iroff3 >= 20) {
1134 if (*last == *limit) {
1142 d__1 = fabs(a1), d__2 = fabs(b2);
1143 if (
max(d__1,d__2) <= (epmach * 100. + 1.) * (fabs(a2) + uflow * 1e3))
1150 if (error2 > error1) {
1153 alist__[*last] =
a2;
1156 elist[maxerr] = error1;
1157 elist[*last] = error2;
1160 alist__[maxerr] =
a2;
1161 alist__[*last] =
a1;
1163 rlist[maxerr] = area2;
1164 rlist[*last] = area1;
1165 elist[maxerr] = error2;
1166 elist[*last] = error1;
1173 dqpsrt_(limit, last, &maxerr, &errmax, &elist[1], &iord[1], &nrmax);
1174 if (errsum <= errbnd) {
1187 if ((d__1 = b1 - a1, fabs(d__1)) > small) {
1197 if ((d__1 = blist[maxerr] - alist__[maxerr], fabs(d__1)) > small) {
1203 if (ierro == 3 || erlarg <= ertest) {
1213 if (*last > *limit / 2 + 2) {
1214 jupbnd = *limit + 3 - *last;
1217 for (k =
id; k <= i__2; ++k) {
1218 maxerr = iord[nrmax];
1219 errmax = elist[maxerr];
1220 if ((d__1 = blist[maxerr] - alist__[maxerr], fabs(d__1)) > small) {
1231 rlist2[numrl2 - 1] = area;
1232 dqelg_(&numrl2, rlist2, &reseps, &abseps, res3la, &nres);
1234 if (ktmin > 5 && *abserr < errsum * .001) {
1237 if (abseps >= *abserr) {
1245 d__1 = *epsabs, d__2 = *epsrel * fabs(reseps);
1246 ertest =
max(d__1,d__2);
1247 if (*abserr <= ertest) {
1261 errmax = elist[maxerr];
1280 if (*abserr == oflow) {
1283 if (*ier + ierro == 0) {
1292 if (*result != 0. && area != 0.) {
1295 if (*abserr > errsum) {
1303 if (*abserr / fabs(*result) > errsum / fabs(area)) {
1311 d__1 = fabs(*result), d__2 = fabs(area);
1312 if (ksgn == -1 &&
max(d__1,d__2) <= defabs * .01) {
1315 if (.01 > *result / area || *result / area > 100. || errsum > fabs(area)) {
1325 for (k = 1; k <= i__1; ++k) {
1326 *result += rlist[k];
1331 *neval = *last * 30 - 15;
1343 const double *epsabs,
const double *epsrel,
double *result,
1344 double *abserr,
long *neval,
long *ier,
long *limit,
1345 const long *lenw,
long *last,
long *iwork,
double *work)
1347 long l1, l2, l3, lvl;
1521 if (*limit < 1 || *lenw < *limit << 2) {
1531 dqagie_(f, bound, inf, epsabs, epsrel, limit, result, abserr, neval,
1532 ier, &work[1], &work[l1], &work[l2], &work[l3], &iwork[1], last);
1542 xerror_(
"abnormal return from dqagi", &
c__26, ier, &lvl, 26);
1547 void dqagpe_(
const D_fp& f,
const double *a,
const double *b,
const long *npts2,
1548 const double *points,
const double *epsabs,
const double *epsrel,
1549 const long *limit,
double *result,
double *abserr,
long *
1550 neval,
long *ier,
double *alist__,
double *blist,
1551 double *rlist,
double *elist,
double *pts,
long *iord,
1552 long *level,
long *ndin,
long *last)
1561 long id, ip1, ind1, ind2;
1563 double resa, dres,
sign;
1566 long nres,
nint, jlow, npts;
1567 double area1, area2, area12;
1570 double defab1, defab2;
1572 double oflow, uflow;
1574 long iroff1, iroff2, iroff3;
1577 double error1, error2, rlist2[52];
1579 double defabs, epmach, erlarg, abseps, correc=0., errbnd, resabs;
1584 long maxerr, levcur;
1587 double ertest, errsum;
1861 d__1 = epmach * 50.;
1862 if (*npts2 < 2 || *limit <= npts ||
1863 ( *epsabs <= 0. && *epsrel <
max(d__1,5e-29 ))) {
1877 pts[1] =
min(*a,*b);
1882 for (i__ = 1; i__ <= i__1; ++i__) {
1883 pts[i__ + 1] = points[i__];
1887 pts[npts + 2] =
max(*a,*b);
1895 for (i__ = 1; i__ <= i__1; ++i__) {
1898 for (j = ip1; j <= i__2; ++j) {
1899 if (pts[i__] <= pts[j]) {
1909 if (pts[1] !=
min(*a,*b) || pts[nintp1] !=
max(*a,*b)) {
1922 for (i__ = 1; i__ <= i__2; ++i__) {
1924 dqk21_(f, &a1, &b1, &area1, &error1, &defabs, &resa);
1928 if (error1 == resa && error1 != 0.) {
1933 elist[i__] = error1;
1943 for (i__ = 1; i__ <= i__2; ++i__) {
1944 if (ndin[i__] == 1) {
1945 elist[i__] = *abserr;
1947 errsum += elist[i__];
1955 dres = fabs(*result);
1957 d__1 = *epsabs, d__2 = *epsrel * dres;
1958 errbnd =
max(d__1,d__2);
1959 if (*abserr <= epmach * 100. * resabs && *abserr > errbnd) {
1966 for (i__ = 1; i__ <= i__2; ++i__) {
1970 for (j = jlow; j <= i__1; ++j) {
1972 if (elist[ind1] > elist[ind2]) {
1980 if (ind1 == iord[i__]) {
1983 iord[k] = iord[i__];
1988 if (*limit < *npts2) {
1992 if (*ier != 0 || *abserr <= errbnd) {
1999 rlist2[0] = *result;
2001 errmax = elist[maxerr];
2020 if (dres >= (1. - epmach * 50.) * resabs) {
2028 for (*last = *npts2; *last <= i__2; ++(*last)) {
2033 levcur = level[maxerr] + 1;
2034 a1 = alist__[maxerr];
2035 b1 = (alist__[maxerr] + blist[maxerr]) * .5;
2039 dqk21_(f, &a1, &b1, &area1, &error1, &resa, &defab1);
2040 dqk21_(f, &a2, &b2, &area2, &error2, &resa, &defab2);
2046 area12 = area1 + area2;
2047 erro12 = error1 + error2;
2048 errsum = errsum + erro12 - errmax;
2049 area = area + area12 - rlist[maxerr];
2050 if (defab1 == error1 || defab2 == error2) {
2053 if ((d__1 = rlist[maxerr] - area12, fabs(d__1)) > fabs(area12) * 1e-5 ||
2054 erro12 < errmax * .99) {
2064 if (*last > 10 && erro12 > errmax) {
2068 level[maxerr] = levcur;
2069 level[*last] = levcur;
2070 rlist[maxerr] = area1;
2071 rlist[*last] = area2;
2073 d__1 = *epsabs, d__2 = *epsrel * fabs(area);
2074 errbnd =
max(d__1,d__2);
2078 if (iroff1 + iroff2 >= 10 || iroff3 >= 20) {
2088 if (*last == *limit) {
2096 d__1 = fabs(a1), d__2 = fabs(b2);
2097 if (
max(d__1,d__2) <= (epmach * 100. + 1.) * (fabs(a2) + uflow * 1e3))
2104 if (error2 > error1) {
2107 alist__[*last] =
a2;
2110 elist[maxerr] = error1;
2111 elist[*last] = error2;
2114 alist__[maxerr] =
a2;
2115 alist__[*last] =
a1;
2117 rlist[maxerr] = area2;
2118 rlist[*last] = area1;
2119 elist[maxerr] = error2;
2120 elist[*last] = error1;
2127 dqpsrt_(limit, last, &maxerr, &errmax, &elist[1], &iord[1], &nrmax);
2129 if (errsum <= errbnd) {
2140 if (levcur + 1 <= levmax) {
2150 if (level[maxerr] + 1 <= levmax) {
2156 if (ierro == 3 || erlarg <= ertest) {
2166 if (*last > *limit / 2 + 2) {
2167 jupbnd = *limit + 3 - *last;
2170 for (k =
id; k <= i__1; ++k) {
2171 maxerr = iord[nrmax];
2172 errmax = elist[maxerr];
2174 if (level[maxerr] + 1 <= levmax) {
2185 rlist2[numrl2 - 1] = area;
2189 dqelg_(&numrl2, rlist2, &reseps, &abseps, res3la, &nres);
2191 if (ktmin > 5 && *abserr < errsum * .001) {
2194 if (abseps >= *abserr) {
2202 d__1 = *epsabs, d__2 = *epsrel * fabs(reseps);
2203 ertest =
max(d__1,d__2);
2205 if (*abserr < ertest) {
2220 errmax = elist[maxerr];
2234 if (*abserr == oflow) {
2237 if (*ier + ierro == 0) {
2246 if (*result != 0. && area != 0.) {
2249 if (*abserr > errsum) {
2257 if (*abserr / fabs(*result) > errsum / fabs(area)) {
2265 d__1 = fabs(*result), d__2 = fabs(area);
2266 if (ksgn == -1 &&
max(d__1,d__2) <= resabs * .01) {
2269 if (.01 > *result / area || *result / area > 100. || errsum > fabs(area)) {
2279 for (k = 1; k <= i__2; ++k) {
2280 *result += rlist[k];
2293 void dqagp_(
const D_fp& f,
const double *a,
const double *b,
const long *npts2,
2294 const double *points,
const double *epsabs,
const double *epsrel,
2295 double *result,
double *abserr,
long *neval,
long *ier,
2296 const long *leniw,
const long *lenw,
long *last,
long *iwork,
2299 long l1, l2, l3, l4, lvl, limit;
2503 if (*leniw < *npts2 * 3 - 2 || *lenw < (*leniw << 1) - *npts2 || *npts2 <
2510 limit = (*leniw - *npts2) / 2;
2516 dqagpe_(f, a, b, npts2, &points[1], epsabs, epsrel, &limit, result,
2517 abserr, neval, ier, &work[1], &work[l1], &work[l2], &work[l3], &
2518 work[l4], &iwork[1], &iwork[l1], &iwork[l2], last);
2528 xerror_(
"abnormal return from dqagp", &
c__26, ier, &lvl, 26);
2534 const double *epsabs,
const double *epsrel,
const long *limit,
2536 double *abserr,
long *neval,
long *ier,
double *alist__,
2537 double *blist,
double *rlist,
double *elist,
long *
2551 double area1, area2, area12;
2552 double small=0., erro12;
2554 double defab1, defab2;
2556 double oflow, uflow;
2558 long iroff1, iroff2, iroff3;
2559 double res3la[3], error1, error2, rlist2[52];
2561 double defabs, epmach, erlarg=0., abseps, correc=0., errbnd, resabs;
2563 double erlast, errmax;
2567 double ertest=0., errsum;
2794 d__1 = epmach * 50.;
2795 if (*epsabs <= 0. && *epsrel <
max(d__1,5e-29)) {
2808 dqk21_(f, a, b, result, abserr, &defabs, &resabs);
2812 dres = fabs(*result);
2814 d__1 = *epsabs, d__2 = *epsrel * dres;
2815 errbnd =
max(d__1,d__2);
2820 if (*abserr <= epmach * 100. * defabs && *abserr > errbnd) {
2826 if (*ier != 0 || ( *abserr <= errbnd && *abserr != resabs ) || *abserr == 0.)
2834 rlist2[0] = *result;
2850 if (dres >= (1. - epmach * 50.) * defabs) {
2858 for (*last = 2; *last <= i__1; ++(*last)) {
2863 a1 = alist__[maxerr];
2864 b1 = (alist__[maxerr] + blist[maxerr]) * .5;
2868 dqk21_(f, &a1, &b1, &area1, &error1, &resabs, &defab1);
2869 dqk21_(f, &a2, &b2, &area2, &error2, &resabs, &defab2);
2874 area12 = area1 + area2;
2875 erro12 = error1 + error2;
2876 errsum = errsum + erro12 - errmax;
2877 area = area + area12 - rlist[maxerr];
2878 if (defab1 == error1 || defab2 == error2) {
2881 if ((d__1 = rlist[maxerr] - area12, fabs(d__1)) > fabs(area12) * 1e-5 ||
2882 erro12 < errmax * .99) {
2892 if (*last > 10 && erro12 > errmax) {
2896 rlist[maxerr] = area1;
2897 rlist[*last] = area2;
2899 d__1 = *epsabs, d__2 = *epsrel * fabs(area);
2900 errbnd =
max(d__1,d__2);
2904 if (iroff1 + iroff2 >= 10 || iroff3 >= 20) {
2914 if (*last == *limit) {
2922 d__1 = fabs(a1), d__2 = fabs(b2);
2923 if (
max(d__1,d__2) <= (epmach * 100. + 1.) * (fabs(a2) + uflow * 1e3))
2930 if (error2 > error1) {
2933 alist__[*last] =
a2;
2936 elist[maxerr] = error1;
2937 elist[*last] = error2;
2940 alist__[maxerr] =
a2;
2941 alist__[*last] =
a1;
2943 rlist[maxerr] = area2;
2944 rlist[*last] = area1;
2945 elist[maxerr] = error2;
2946 elist[*last] = error1;
2953 dqpsrt_(limit, last, &maxerr, &errmax, &elist[1], &iord[1], &nrmax);
2955 if (errsum <= errbnd) {
2969 if ((d__1 = b1 - a1, fabs(d__1)) > small) {
2979 if ((d__1 = blist[maxerr] - alist__[maxerr], fabs(d__1)) > small) {
2985 if (ierro == 3 || erlarg <= ertest) {
2995 if (*last > *limit / 2 + 2) {
2996 jupbnd = *limit + 3 - *last;
2999 for (k =
id; k <= i__2; ++k) {
3000 maxerr = iord[nrmax];
3001 errmax = elist[maxerr];
3003 if ((d__1 = blist[maxerr] - alist__[maxerr], fabs(d__1)) > small) {
3014 rlist2[numrl2 - 1] = area;
3015 dqelg_(&numrl2, rlist2, &reseps, &abseps, res3la, &nres);
3017 if (ktmin > 5 && *abserr < errsum * .001) {
3020 if (abseps >= *abserr) {
3028 d__1 = *epsabs, d__2 = *epsrel * fabs(reseps);
3029 ertest =
max(d__1,d__2);
3031 if (*abserr <= ertest) {
3045 errmax = elist[maxerr];
3052 small = (d__1 = *b - *a, fabs(d__1)) * .375;
3064 if (*abserr == oflow) {
3067 if (*ier + ierro == 0) {
3076 if (*result != 0. && area != 0.) {
3079 if (*abserr > errsum) {
3087 if (*abserr / fabs(*result) > errsum / fabs(area)) {
3095 d__1 = fabs(*result), d__2 = fabs(area);
3096 if (ksgn == -1 &&
max(d__1,d__2) <= defabs * .01) {
3099 if (.01 > *result / area || *result / area > 100. || errsum > fabs(area)) {
3109 for (k = 1; k <= i__1; ++k) {
3110 *result += rlist[k];
3119 *neval = *last * 42 - 21;
3124 void dqags_(
const D_fp& f,
const double *a,
const double *b,
const double *epsabs,
3125 const double *epsrel,
double *result,
double *abserr,
3126 long *neval,
long *ier,
const long *limit,
const long *lenw,
long *
3127 last,
long *iwork,
double *work)
3129 long l1, l2, l3, lvl;
3300 if (*limit < 1 || *lenw < *limit << 2) {
3310 dqagse_(f, a, b, epsabs, epsrel, limit, result, abserr, neval, ier,
3311 &work[1], &work[l1], &work[l2], &work[l3], &iwork[1], last);
3321 xerror_(
"abnormal return from dqags", &
c__26, ier, &lvl, 26);
3326 void dqawce_(
const D_fp& f,
const double *a,
const double *b,
const double *c__,
3327 const double *epsabs,
const double *epsrel,
const long *limit,
3328 double *result,
double *abserr,
long *neval,
long *ier,
3329 double *alist__,
double *blist,
double *rlist,
double
3330 *elist,
long *iord,
long *last)
3340 double area, area1, area2, area12;
3344 long iroff1, iroff2;
3345 double error1, error2, epmach, errbnd, errmax;
3543 d__1 = epmach * 50.;
3544 if (*c__ == *a || *c__ == *b || ( *epsabs <= 0. && *epsrel <
max(d__1,5e-29) )
3562 dqc25c_(f, &aa, &bb, c__, result, abserr, &krule, neval);
3573 d__1 = *epsabs, d__2 = *epsrel * fabs(*result);
3574 errbnd =
max(d__1,d__2);
3579 d__1 = fabs(*result) * .01;
3580 if (*abserr <
min(d__1,errbnd) || *ier == 1) {
3602 for (*last = 2; *last <= i__1; ++(*last)) {
3607 a1 = alist__[maxerr];
3608 b1 = (alist__[maxerr] + blist[maxerr]) * .5;
3610 if (*c__ <= b1 && *c__ > a1) {
3611 b1 = (*c__ +
b2) * .5;
3613 if (*c__ > b1 && *c__ < b2) {
3614 b1 = (a1 + *c__) * .5;
3618 dqc25c_(f, &a1, &b1, c__, &area1, &error1, &krule, &nev);
3620 dqc25c_(f, &a2, &b2, c__, &area2, &error2, &krule, &nev);
3626 area12 = area1 + area2;
3627 erro12 = error1 + error2;
3628 errsum = errsum + erro12 - errmax;
3629 area = area + area12 - rlist[maxerr];
3630 if ((d__1 = rlist[maxerr] - area12, fabs(d__1)) < fabs(area12) * 1e-5 &&
3631 erro12 >= errmax * .99 && krule == 0) {
3634 if (*last > 10 && erro12 > errmax && krule == 0) {
3637 rlist[maxerr] = area1;
3638 rlist[*last] = area2;
3640 d__1 = *epsabs, d__2 = *epsrel * fabs(area);
3641 errbnd =
max(d__1,d__2);
3642 if (errsum <= errbnd) {
3648 if (iroff1 >= 6 && iroff2 > 20) {
3655 if (*last == *limit) {
3663 d__1 = fabs(a1), d__2 = fabs(b2);
3664 if (
max(d__1,d__2) <= (epmach * 100. + 1.) * (fabs(a2) + uflow * 1e3))
3672 if (error2 > error1) {
3675 alist__[*last] =
a2;
3678 elist[maxerr] = error1;
3679 elist[*last] = error2;
3682 alist__[maxerr] =
a2;
3683 alist__[*last] =
a1;
3685 rlist[maxerr] = area2;
3686 rlist[*last] = area1;
3687 elist[maxerr] = error2;
3688 elist[*last] = error1;
3695 dqpsrt_(limit, last, &maxerr, &errmax, &elist[1], &iord[1], &nrmax);
3697 if (*ier != 0 || errsum <= errbnd) {
3709 for (k = 1; k <= i__1; ++k) {
3710 *result += rlist[k];
3716 *result = -(*result);
3722 void dqawc_(
const D_fp& f,
const double *a,
const double *b,
const double *c__,
3723 const double *epsabs,
const double *epsrel,
double *result,
3724 double *abserr,
long *neval,
long *ier,
long *limit,
3725 const long *lenw,
long *last,
long *iwork,
double *work)
3727 long l1, l2, l3, lvl;
3889 if (*limit < 1 || *lenw < *limit << 2) {
3898 dqawce_(f, a, b, c__, epsabs, epsrel, limit, result, abserr, neval,
3899 ier, &work[1], &work[l1], &work[l2], &work[l3], &iwork[1], last);
3909 xerror_(
"abnormal return from dqawc", &
c__26, ier, &lvl, 26);
3915 const long *integr,
const double *epsabs,
const long *limlst,
3916 const long *limit,
const long *maxp1,
3917 double *result,
double *abserr,
long *neval,
3918 long *ier,
double *rslst,
double *erlst,
long *ierlst,
3919 long *lst,
double *alist__,
double *blist,
3920 double *rlist,
double *elist,
long *iord,
long *nnlog,
3926 double pi = 3.1415926535897932384626433832795;
3929 int chebmo_dim1, chebmo_offset, i__1;
3934 double c1, c2, p1, dl, ep;
3946 double abseps, correc;
3948 double reseps, errsum;
4182 chebmo_dim1 = *maxp1;
4183 chebmo_offset = 1 + chebmo_dim1;
4184 chebmo -= chebmo_offset;
4197 if ( (*integr != 1 && *integr != 2 ) || *epsabs <= 0. || *limlst < 3) {
4212 neval, ier, &alist__[1], &blist[1], &rlist[1], &elist[1], &
4225 l = (int) fabs(*omega);
4226 dl = (double) ((l << 1) + 1);
4227 cycle = dl * pi / fabs(*omega);
4238 if (*epsabs > uflow / p1) {
4251 for (*lst = 1; *lst <= i__1; ++(*lst)) {
4256 dqawoe_(f, &c1, &c2, omega, integr, &epsa, &
c_b20, limit, lst,
4257 maxp1, &rslst[*lst], &erlst[*lst], &nev, &ierlst[*lst], &last,
4258 &alist__[1], &blist[1], &rlist[1], &elist[1], &iord[1], &
4259 nnlog[1], &momcom, &chebmo[chebmo_offset]);
4262 errsum += erlst[*lst];
4263 drl = (d__1 = rslst[*lst], fabs(d__1)) * 50.;
4267 if (errsum + drl <= *epsabs && *lst >= 6) {
4271 d__1 = correc, d__2 = erlst[*lst];
4272 correc =
max(d__1,d__2);
4273 if (ierlst[*lst] != 0) {
4275 d__1 = ep, d__2 = correc * p1;
4276 eps =
max(d__1,d__2);
4278 if (ierlst[*lst] != 0) {
4281 if (*ier == 7 && errsum + drl <= correc * 10. && *lst > 5) {
4291 psum[numrl2 - 1] = psum[ll - 1] + rslst[*lst];
4298 if (*lst == *limlst) {
4304 dqelg_(&numrl2, psum, &reseps, &abseps, res3la, &nres);
4309 if (ktmin >= 15 && *abserr <= (errsum + drl) * .001) {
4312 if (abseps > *abserr && *lst != 3) {
4323 if (*abserr + correc * 10. <= *epsabs ||
4324 ( *abserr <= *epsabs && correc * 10. >= *epsabs ) ) {
4328 if (*ier != 0 && *ier != 7) {
4342 *abserr += correc * 10.;
4346 if (*result != 0. && psum[numrl2 - 1] != 0.) {
4349 if (*abserr > errsum) {
4352 if (psum[numrl2 - 1] == 0.) {
4356 if (*abserr / fabs(*result) > (errsum + drl) / (d__1 = psum[numrl2 - 1],
4360 if (*ier >= 1 && *ier != 7) {
4365 *result = psum[numrl2 - 1];
4366 *abserr = errsum + drl;
4371 void dqawf_(
const D_fp& f,
const double *a,
const double *omega,
const long *integr,
4372 const double *epsabs,
double *result,
double *abserr,
4373 long *neval,
long *ier,
long *limlst,
long *lst,
const long *
4374 leniw,
const long *maxp1,
const long *lenw,
long *iwork,
double *
4377 long l1, l2, l3, l4, l5, l6, ll2, lvl, limit;
4583 if (*limlst < 3 || *leniw < *limlst + 2 || *maxp1 < 1 || *lenw < (*leniw
4584 << 1) + *maxp1 * 25) {
4590 limit = (*leniw - *limlst) / 2;
4598 dqawfe_(f, a, omega, integr, epsabs, limlst, &limit, maxp1, result,
4599 abserr, neval, ier, &work[1], &work[l1], &iwork[1], lst, &work[l2]
4600 , &work[l3], &work[l4], &work[l5], &iwork[l1], &iwork[ll2], &work[
4611 xerror_(
"abnormal return from dqawf", &
c__26, ier, &lvl, 26);
4616 void dqawoe_(
const D_fp& f,
const double *a,
const double *b,
const double
4617 *omega,
const long *integr,
const double *epsabs,
4618 const double *epsrel,
4619 const long *limit,
const long *icall,
const long *maxp1,
4621 double *abserr,
long *neval,
long *ier,
long *last,
4622 double *alist__,
double *blist,
double *rlist,
double
4623 *elist,
long *iord,
long *nnlog,
long *momcom,
double *
4627 int chebmo_dim1, chebmo_offset, i__1, i__2;
4636 double area1, area2, area12;
4637 double small, erro12, width, defab1, defab2;
4643 long iroff1, iroff2, iroff3;
4644 double res3la[3], error1, error2, rlist2[52];
4646 double defabs, domega, epmach, erlarg=0., abseps, correc=0., errbnd,
4650 double erlast, errmax;
4654 double ertest=0., errsum;
4919 chebmo_dim1 = *maxp1;
4920 chebmo_offset = 1 + chebmo_dim1;
4921 chebmo -= chebmo_offset;
4941 d__1 = epmach * 50.;
4942 if ( ( *integr != 1 && *integr != 2 ) ||
4943 ( *epsabs <= 0. && *epsrel <
max(d__1, 5e-29) ) ||
4944 *icall < 1 || *maxp1 < 1) {
4954 domega = fabs(*omega);
4961 dqc25f_(f, a, b, &domega, integr, &nrmom, maxp1, &
c__0, result,
4962 abserr, neval, &defabs, &resabs, momcom, &chebmo[chebmo_offset]);
4966 dres = fabs(*result);
4968 d__1 = *epsabs, d__2 = *epsrel * dres;
4969 errbnd =
max(d__1,d__2);
4973 if (*abserr <= epmach * 100. * defabs && *abserr > errbnd) {
4979 if (*ier != 0 || *abserr <= errbnd) {
5001 small = (d__1 = *b - *a, fabs(d__1)) * .75;
5005 if ((d__1 = *b - *a, fabs(d__1)) * .5 * domega > 2.) {
5010 rlist2[0] = *result;
5012 if ((d__1 = *b - *a, fabs(d__1)) * .25 * domega <= 2.) {
5016 if (dres >= (1. - epmach * 50.) * defabs) {
5024 for (*last = 2; *last <= i__1; ++(*last)) {
5029 nrmom = nnlog[maxerr] + 1;
5030 a1 = alist__[maxerr];
5031 b1 = (alist__[maxerr] + blist[maxerr]) * .5;
5035 dqc25f_(f, &a1, &b1, &domega, integr, &nrmom, maxp1, &
c__0, &
5036 area1, &error1, &nev, &resabs, &defab1, momcom, &chebmo[
5039 dqc25f_(f, &a2, &b2, &domega, integr, &nrmom, maxp1, &
c__1, &
5040 area2, &error2, &nev, &resabs, &defab2, momcom, &chebmo[
5047 area12 = area1 + area2;
5048 erro12 = error1 + error2;
5049 errsum = errsum + erro12 - errmax;
5050 area = area + area12 - rlist[maxerr];
5051 if (defab1 == error1 || defab2 == error2) {
5054 if ((d__1 = rlist[maxerr] - area12, fabs(d__1)) > fabs(area12) * 1e-5 ||
5055 erro12 < errmax * .99) {
5065 if (*last > 10 && erro12 > errmax) {
5069 rlist[maxerr] = area1;
5070 rlist[*last] = area2;
5071 nnlog[maxerr] = nrmom;
5072 nnlog[*last] = nrmom;
5074 d__1 = *epsabs, d__2 = *epsrel * fabs(area);
5075 errbnd =
max(d__1,d__2);
5079 if (iroff1 + iroff2 >= 10 || iroff3 >= 20) {
5089 if (*last == *limit) {
5097 d__1 = fabs(a1), d__2 = fabs(b2);
5098 if (
max(d__1,d__2) <= (epmach * 100. + 1.) * (fabs(a2) + uflow * 1e3))
5105 if (error2 > error1) {
5108 alist__[*last] =
a2;
5111 elist[maxerr] = error1;
5112 elist[*last] = error2;
5115 alist__[maxerr] =
a2;
5116 alist__[*last] =
a1;
5118 rlist[maxerr] = area2;
5119 rlist[*last] = area1;
5120 elist[maxerr] = error2;
5121 elist[*last] = error1;
5128 dqpsrt_(limit, last, &maxerr, &errmax, &elist[1], &iord[1], &nrmax);
5130 if (errsum <= errbnd) {
5136 if (*last == 2 && extall) {
5146 if ((d__1 = b1 - a1, fabs(d__1)) > small) {
5157 width = (d__1 = blist[maxerr] - alist__[maxerr], fabs(d__1));
5158 if (width > small) {
5170 if (width * .25 * domega > 2.) {
5179 if (ierro == 3 || erlarg <= ertest) {
5188 if (*last > *limit / 2 + 2) {
5189 jupbnd = *limit + 3 - *last;
5193 for (k =
id; k <= i__2; ++k) {
5194 maxerr = iord[nrmax];
5195 errmax = elist[maxerr];
5196 if ((d__1 = blist[maxerr] - alist__[maxerr], fabs(d__1)) > small) {
5207 rlist2[numrl2 - 1] = area;
5211 dqelg_(&numrl2, rlist2, &reseps, &abseps, res3la, &nres);
5213 if (ktmin > 5 && *abserr < errsum * .001) {
5216 if (abseps >= *abserr) {
5224 d__1 = *epsabs, d__2 = *epsrel * fabs(reseps);
5225 ertest =
max(d__1,d__2);
5227 if (*abserr <= ertest) {
5242 errmax = elist[maxerr];
5251 rlist2[numrl2 - 1] = area;
5263 if (*abserr == oflow || nres == 0) {
5266 if (*ier + ierro == 0) {
5275 if (*result != 0. && area != 0.) {
5278 if (*abserr > errsum) {
5286 if (*abserr / fabs(*result) > errsum / fabs(area)) {
5294 d__1 = fabs(*result), d__2 = fabs(area);
5295 if (ksgn == -1 &&
max(d__1,d__2) <= defabs * .01) {
5298 if (.01 > *result / area || *result / area > 100. || errsum >= fabs(area))
5309 for (k = 1; k <= i__1; ++k) {
5310 *result += rlist[k];
5319 if (*integr == 2 && *omega < 0.) {
5320 *result = -(*result);
5326 void dqawo_(
const D_fp& f,
const double *a,
const double *b,
const double *omega,
5327 const long *integr,
const double *epsabs,
const double *epsrel,
5328 double *result,
double *abserr,
long *neval,
long *ier,
5329 const long *leniw,
long *maxp1,
const long *lenw,
long *last,
5330 long *iwork,
double *work)
5332 long l1, l2, l3, l4, lvl, limit;
5537 if (*leniw < 2 || *maxp1 < 1 || *lenw < (*leniw << 1) + *maxp1 * 25) {
5548 dqawoe_(f, a, b, omega, integr, epsabs, epsrel, &limit, &
c__1,
5549 maxp1, result, abserr, neval, ier, last, &work[1], &work[l1], &
5550 work[l2], &work[l3], &iwork[1], &iwork[l1], &momcom, &work[l4]);
5560 xerror_(
"abnormal return from dqawo", &
c__26, ier, &lvl, 26);
5565 void dqawse_(
const D_fp& f,
const double *a,
const double *b,
const double *alfa,
5566 const double *beta,
const long *integr,
const double *epsabs,
5567 const double *epsrel,
const long *limit,
double *result,
5568 double *abserr,
long *neval,
long *ier,
double *alist__,
5569 double *blist,
double *rlist,
double *elist,
long *iord,
5578 double a1,
a2,
b1,
b2, rg[25], rh[25], ri[25], rj[25];
5580 double area, area1, area2, area12;
5584 long iroff1, iroff2;
5585 double resas1, resas2, error1, error2, epmach, errbnd, centre;
5803 d__1 = epmach * 50.;
5804 if (*b <= *a || ( *epsabs == 0. && *epsrel <
max(d__1,5e-29) ) || *alfa <=
5805 -1. || *beta <= -1. || *integr < 1 || *integr > 4 || *limit < 2) {
5812 dqmomo_(alfa, beta, ri, rj, rg, rh, integr);
5816 centre = (*b + *a) * .5;
5817 dqc25s_(f, a, b, a, ¢re, alfa, beta, ri, rj, rg, rh, &area1, &
5818 error1, &resas1, integr, &nev);
5820 dqc25s_(f, a, b, ¢re, b, alfa, beta, ri, rj, rg, rh, &area2, &
5821 error2, &resas2, integr, &nev);
5824 *result = area1 + area2;
5825 *abserr = error1 + error2;
5830 d__1 = *epsabs, d__2 = *epsrel * fabs(*result);
5831 errbnd =
max(d__1,d__2);
5836 if (error2 > error1) {
5840 alist__[2] = centre;
5849 alist__[1] = centre;
5863 if (*abserr <= errbnd || *ier == 1) {
5878 for (*last = 3; *last <= i__1; ++(*last)) {
5882 a1 = alist__[maxerr];
5883 b1 = (alist__[maxerr] + blist[maxerr]) * .5;
5887 dqc25s_(f, a, b, &a1, &b1, alfa, beta, ri, rj, rg, rh, &area1, &
5888 error1, &resas1, integr, &nev);
5890 dqc25s_(f, a, b, &a2, &b2, alfa, beta, ri, rj, rg, rh, &area2, &
5891 error2, &resas2, integr, &nev);
5897 area12 = area1 + area2;
5898 erro12 = error1 + error2;
5899 errsum = errsum + erro12 - errmax;
5900 area = area + area12 - rlist[maxerr];
5901 if (*a == a1 || *b == b2) {
5904 if (resas1 == error1 || resas2 == error2) {
5910 if ((d__1 = rlist[maxerr] - area12, fabs(d__1)) < fabs(area12) * 1e-5 &&
5911 erro12 >= errmax * .99) {
5914 if (*last > 10 && erro12 > errmax) {
5918 rlist[maxerr] = area1;
5919 rlist[*last] = area2;
5924 d__1 = *epsabs, d__2 = *epsrel * fabs(area);
5925 errbnd =
max(d__1,d__2);
5926 if (errsum <= errbnd) {
5933 if (*last == *limit) {
5940 if (iroff1 >= 6 || iroff2 >= 20) {
5948 d__1 = fabs(a1), d__2 = fabs(b2);
5949 if (
max(d__1,d__2) <= (epmach * 100. + 1.) * (fabs(a2) + uflow * 1e3))
5957 if (error2 > error1) {
5960 alist__[*last] =
a2;
5963 elist[maxerr] = error1;
5964 elist[*last] = error2;
5967 alist__[maxerr] =
a2;
5968 alist__[*last] =
a1;
5970 rlist[maxerr] = area2;
5971 rlist[*last] = area1;
5972 elist[maxerr] = error2;
5973 elist[*last] = error1;
5980 dqpsrt_(limit, last, &maxerr, &errmax, &elist[1], &iord[1], &nrmax);
5982 if (*ier != 0 || errsum <= errbnd) {
5994 for (k = 1; k <= i__1; ++k) {
5995 *result += rlist[k];
6003 void dqaws_(
const D_fp& f,
const double *a,
const double *b,
const double *alfa,
6004 const double *beta,
const long *integr,
const double *epsabs,
6005 const double *epsrel,
double *result,
double *abserr,
6006 long *neval,
long *ier,
long *limit,
const long *lenw,
long *last,
6007 long *iwork,
double *work)
6009 long l1, l2, l3, lvl;
6192 if (*limit < 2 || *lenw < *limit << 2) {
6202 dqawse_(f, a, b, alfa, beta, integr, epsabs, epsrel, limit, result,
6203 abserr, neval, ier, &work[1], &work[l1], &work[l2], &work[l3], &
6214 xerror_(
"abnormal return from dqaws", &
c__26, ier, &lvl, 26);
6219 void dqc25c_(
const D_fp& f,
const double *a,
const double *b,
const double *c__,
6220 double *result,
double *abserr,
long *krul,
long *neval)
6224 double x[11] = { .991444861373810411144557526928563,
6225 .965925826289068286749743199728897,
6226 .923879532511286756128183189396788,
6227 .866025403784438646763723170752936,
6228 .793353340291235164579776961501299,
6229 .707106781186547524400844362104849,
6230 .608761429008720639416097542898164,.5,
6231 .382683432365089771728459984030399,
6232 .258819045102520762348898837624048,
6233 .130526192220051591548406227895489 };
6240 double u, p2, p3, p4, cc;
6242 double ak22, fval[25], res12, res24;
6244 double amom0, amom1, amom2, cheb12[13], cheb24[25], hlgth,
6246 double resabs, resasc;
6332 cc = (2. * *c__ - *b - *a) / (*b - *a);
6333 if (fabs(cc) < 1.1) {
6340 dqk15w_(f,
dqwgtc_, c__, &p2, &p3, &p4, &kp, a, b, result,
6341 abserr, &resabs, &resasc);
6343 if (resasc == *abserr) {
6351 hlgth = (*b - *a) * .5;
6352 centr = (*b + *a) * .5;
6354 d__1 = hlgth + centr;
6355 fval[0] = f(d__1) * .5;
6356 fval[12] = f(centr);
6357 d__1 = centr - hlgth;
6358 fval[24] = f(d__1) * .5;
6359 for (i__ = 2; i__ <= 12; ++i__) {
6360 u = hlgth * x[i__ - 2];
6363 fval[i__ - 1] = f(d__1);
6365 fval[isym - 1] = f(d__1);
6371 dqcheb_(x, fval, cheb12, cheb24);
6376 amom0 = log((d__1 = (1. - cc) / (cc + 1.), fabs(d__1)));
6377 amom1 = cc * amom0 + 2.;
6378 res12 = cheb12[0] * amom0 + cheb12[1] * amom1;
6379 res24 = cheb24[0] * amom0 + cheb24[1] * amom1;
6380 for (k = 3; k <= 13; ++k) {
6381 amom2 = cc * 2. * amom1 - amom0;
6382 ak22 = (double) ((k - 2) * (k - 2));
6383 if (k / 2 << 1 == k) {
6384 amom2 -= 4. / (ak22 - 1.);
6386 res12 += cheb12[k - 1] * amom2;
6387 res24 += cheb24[k - 1] * amom2;
6392 for (k = 14; k <= 25; ++k) {
6393 amom2 = cc * 2. * amom1 - amom0;
6394 ak22 = (double) ((k - 2) * (k - 2));
6395 if (k / 2 << 1 == k) {
6396 amom2 -= 4. / (ak22 - 1.);
6398 res24 += cheb24[k - 1] * amom2;
6404 *abserr = (d__1 = res24 - res12, fabs(d__1));
6409 void dqc25f_(
const D_fp& f,
const double *a,
const double *b,
const double
6410 *omega,
const long *integr,
const long *nrmom,
const long *maxp1,
6411 const long *ksave,
double *result,
double *abserr,
long *neval,
6412 double *resabs,
double *resasc,
long *momcom,
double *
6417 double x[11] = { .991444861373810411144557526928563,
6418 .965925826289068286749743199728897,
6419 .923879532511286756128183189396788,
6420 .866025403784438646763723170752936,
6421 .793353340291235164579776961501299,
6422 .707106781186547524400844362104849,
6423 .608761429008720639416097542898164,.5,
6424 .382683432365089771728459984030399,
6425 .258819045102520762348898837624048,
6426 .130526192220051591548406227895489 };
6429 int chebmo_dim1, chebmo_offset, i__1;
6434 long i__, j, k, m=0;
6435 double v[28], d1[25], d2[25], p2, p3, p4, ac, an, as, an2, ass,
6436 par2, conc, asap, par22, fval[25], estc, cons;
6440 double cheb12[13], cheb24[25], resc12, resc24, hlgth, centr;
6441 double ress12, ress24, oflow;
6548 chebmo_dim1 = *maxp1;
6549 chebmo_offset = 1 + chebmo_dim1;
6550 chebmo -= chebmo_offset;
6585 centr = (*b + *a) * .5;
6586 hlgth = (*b - *a) * .5;
6587 parint = *omega * hlgth;
6593 if (fabs(parint) > 2.) {
6597 result, abserr, resabs, resasc);
6605 conc = hlgth * cos(centr * *omega);
6606 cons = hlgth * sin(centr * *omega);
6613 if (*nrmom < *momcom || *ksave == 1) {
6620 par2 = parint * parint;
6622 sinpar = sin(parint);
6623 cospar = cos(parint);
6627 v[0] = sinpar * 2. / parint;
6628 v[1] = (cospar * 8. + (par2 + par2 - 8.) * sinpar / parint) / par2;
6629 v[2] = ((par2 - 12.) * 32. * cospar + ((par2 - 80.) * par2 + 192.) * 2. *
6630 sinpar / parint) / (par2 * par2);
6632 as = parint * 24. *
sinpar;
6633 if (fabs(parint) > 24.) {
6645 for (k = 1; k <= i__1; ++k) {
6647 d__[k - 1] = (an2 - 4.) * -2. * (par22 - an2 - an2);
6648 d2[k - 1] = (an - 1.) * (an - 2.) * par2;
6649 d1[k] = (an + 3.) * (an + 4.) * par2;
6650 v[k + 2] = as - (an2 - 4.) * ac;
6655 d__[noequ - 1] = (an2 - 4.) * -2. * (par22 - an2 - an2);
6656 v[noequ + 2] = as - (an2 - 4.) * ac;
6657 v[3] -= par2 * 56. * v[2];
6659 asap = (((((par2 * 210. - 1.) * cospar - (par2 * 105. - 63.) * ass) / an2
6660 - (1. - par2 * 15.) * cospar + ass * 15.) / an2 - cospar + ass *
6661 3.) / an2 - cospar) / an2;
6662 v[noequ + 2] -= asap * 2. * par2 * (an - 1.) * (an - 2.);
6670 dgtsl_(&noequ, d1, d__, d2, &v[3], &iers);
6678 for (i__ = 4; i__ <= 13; ++i__) {
6680 v[i__ - 1] = ((an2 - 4.) * ((par22 - an2 - an2) * 2. * v[i__ - 2] -
6681 ac) + as - par2 * (an + 1.) * (an + 2.) * v[i__ - 3]) / (par2
6682 * (an - 1.) * (an - 2.));
6687 for (j = 1; j <= 13; ++j) {
6688 chebmo[m + ((j << 1) - 1) * chebmo_dim1] = v[j - 1];
6694 v[0] = (sinpar - parint * cospar) * 2. / par2;
6695 v[1] = (18. - 48. / par2) * sinpar / par2 + (48. / par2 - 2.) * cospar /
6697 ac = parint * -24. * cospar;
6699 if (fabs(parint) > 24.) {
6709 for (k = 1; k <= i__1; ++k) {
6711 d__[k - 1] = (an2 - 4.) * -2. * (par22 - an2 - an2);
6712 d2[k - 1] = (an - 1.) * (an - 2.) * par2;
6713 d1[k] = (an + 3.) * (an + 4.) * par2;
6714 v[k + 1] = ac + (an2 - 4.) * as;
6719 d__[noequ - 1] = (an2 - 4.) * -2. * (par22 - an2 - an2);
6720 v[noequ + 1] = ac + (an2 - 4.) * as;
6721 v[2] -= par2 * 42. * v[1];
6722 ass = parint * cospar;
6723 asap = (((((par2 * 105. - 63.) * ass + (par2 * 210. - 1.) *
sinpar) / an2
6724 + (par2 * 15. - 1.) * sinpar - ass * 15.) / an2 - ass * 3. -
6725 sinpar) / an2 - sinpar) / an2;
6726 v[noequ + 1] -= asap * 2. * par2 * (an - 1.) * (an - 2.);
6734 dgtsl_(&noequ, d1, d__, d2, &v[2], &iers);
6741 for (i__ = 3; i__ <= 12; ++i__) {
6743 v[i__ - 1] = ((an2 - 4.) * ((par22 - an2 - an2) * 2. * v[i__ - 2] +
6744 as) + ac - par2 * (an + 1.) * (an + 2.) * v[i__ - 3]) / (par2
6745 * (an - 1.) * (an - 2.));
6750 for (j = 1; j <= 12; ++j) {
6751 chebmo[m + (j << 1) * chebmo_dim1] = v[j - 1];
6755 if (*nrmom < *momcom) {
6758 if (*momcom < *maxp1 - 1 && *nrmom >= *momcom) {
6765 d__1 = centr + hlgth;
6766 fval[0] = f(d__1) * .5;
6767 fval[12] = f(centr);
6768 d__1 = centr - hlgth;
6769 fval[24] = f(d__1) * .5;
6770 for (i__ = 2; i__ <= 12; ++i__) {
6772 d__1 = hlgth * x[i__ - 2] + centr;
6773 fval[i__ - 1] = f(d__1);
6774 d__1 = centr - hlgth * x[i__ - 2];
6775 fval[isym - 1] = f(d__1);
6778 dqcheb_(x, fval, cheb12, cheb24);
6782 resc12 = cheb12[12] * chebmo[m + chebmo_dim1 * 13];
6785 for (j = 1; j <= 6; ++j) {
6786 resc12 += cheb12[k - 1] * chebmo[m + k * chebmo_dim1];
6787 ress12 += cheb12[k] * chebmo[m + (k + 1) * chebmo_dim1];
6791 resc24 = cheb24[24] * chebmo[m + chebmo_dim1 * 25];
6793 *resabs = fabs(cheb24[24]);
6795 for (j = 1; j <= 12; ++j) {
6796 resc24 += cheb24[k - 1] * chebmo[m + k * chebmo_dim1];
6797 ress24 += cheb24[k] * chebmo[m + (k + 1) * chebmo_dim1];
6798 *resabs = (d__1 = cheb24[k - 1], fabs(d__1)) + (d__2 = cheb24[k], fabs(
6803 estc = (d__1 = resc24 - resc12, fabs(d__1));
6804 ests = (d__1 = ress24 - ress12, fabs(d__1));
6805 *resabs *= fabs(hlgth);
6809 *result = conc * resc24 - cons * ress24;
6810 *abserr = (d__1 = conc * estc, fabs(d__1)) + (d__2 = cons * ests, fabs(d__2)
6814 *result = conc * ress24 + cons * resc24;
6815 *abserr = (d__1 = conc * ests, fabs(d__1)) + (d__2 = cons * estc, fabs(d__2)
6822 const double *bl,
const double *br,
const double *alfa,
6823 const double *beta,
const double *ri,
const double *rj,
6824 const double *rg,
const double *rh,
6825 double *result,
double *abserr,
double *resasc,
6826 const long *integr,
long *nev)
6830 double x[11] = { .991444861373810411144557526928563,
6831 .965925826289068286749743199728897,
6832 .923879532511286756128183189396788,
6833 .866025403784438646763723170752936,
6834 .793353340291235164579776961501299,
6835 .707106781186547524400844362104849,
6836 .608761429008720639416097542898164,.5,
6837 .382683432365089771728459984030399,
6838 .258819045102520762348898837624048,
6839 .130526192220051591548406227895489 };
6846 double u, dc, fix, fval[25], res12, res24;
6848 double cheb12[13], cheb24[25], hlgth, centr;
6849 double factor, resabs;
6963 if (*bl == *a && (*alfa != 0. || *integr == 2 || *integr == 4)) {
6966 if (*br == *b && (*beta != 0. || *integr == 3 || *integr == 4)) {
6974 dqk15w_(f,
dqwgts_, a, b, alfa, beta, integr, bl, br, result,
6975 abserr, &resabs, resasc);
6988 hlgth = (*br - *bl) * .5;
6989 centr = (*br + *bl) * .5;
6991 d__1 = hlgth + centr;
6993 fval[0] = f(d__1) * .5 * pow(d__2, *beta);
6994 fval[12] = f(centr) * pow(fix, *beta);
6995 d__1 = centr - hlgth;
6997 fval[24] = f(d__1) * .5 * pow(d__2, *beta);
6998 for (i__ = 2; i__ <= 12; ++i__) {
6999 u = hlgth * x[i__ - 2];
7003 fval[i__ - 1] = f(d__1) * pow(d__2, *beta);
7006 fval[isym - 1] = f(d__1) * pow(d__2, *beta);
7010 factor = pow(hlgth, d__1);
7018 dqcheb_(x, fval, cheb12, cheb24);
7022 for (i__ = 1; i__ <= 13; ++i__) {
7023 res12 += cheb12[i__ - 1] * ri[i__];
7024 res24 += cheb24[i__ - 1] * ri[i__];
7027 for (i__ = 14; i__ <= 25; ++i__) {
7028 res24 += cheb24[i__ - 1] * ri[i__];
7037 dc = log(*br - *bl);
7038 *result = res24 * dc;
7039 *abserr = (d__1 = (res24 - res12) * dc, fabs(d__1));
7042 for (i__ = 1; i__ <= 13; ++i__) {
7043 res12 += cheb12[i__ - 1] * rg[i__];
7044 res24 = res12 + cheb24[i__ - 1] * rg[i__];
7047 for (i__ = 14; i__ <= 25; ++i__) {
7048 res24 += cheb24[i__ - 1] * rg[i__];
7058 fval[0] *= log(fix - hlgth);
7059 fval[12] *= log(fix);
7060 fval[24] *= log(fix + hlgth);
7061 for (i__ = 2; i__ <= 12; ++i__) {
7062 u = hlgth * x[i__ - 2];
7064 fval[i__ - 1] *= log(fix - u);
7065 fval[isym - 1] *= log(fix + u);
7068 dqcheb_(x, fval, cheb12, cheb24);
7072 for (i__ = 1; i__ <= 13; ++i__) {
7073 res12 += cheb12[i__ - 1] * ri[i__];
7074 res24 += cheb24[i__ - 1] * ri[i__];
7077 for (i__ = 14; i__ <= 25; ++i__) {
7078 res24 += cheb24[i__ - 1] * ri[i__];
7087 dc = log(*br - *bl);
7088 *result = res24 * dc;
7089 *abserr = (d__1 = (res24 - res12) * dc, fabs(d__1));
7092 for (i__ = 1; i__ <= 13; ++i__) {
7093 res12 += cheb12[i__ - 1] * rg[i__];
7094 res24 += cheb24[i__ - 1] * rg[i__];
7097 for (i__ = 14; i__ <= 25; ++i__) {
7098 res24 += cheb24[i__ - 1] * rg[i__];
7102 *result = (*result + res24) * factor;
7103 *abserr = (*abserr + (d__1 = res24 - res12, fabs(d__1))) * factor;
7115 hlgth = (*br - *bl) * .5;
7116 centr = (*br + *bl) * .5;
7118 d__1 = hlgth + centr;
7120 fval[0] = f(d__1) * .5 * pow(d__2, *alfa);
7121 fval[12] = f(centr) * pow(fix, *alfa);
7122 d__1 = centr - hlgth;
7124 fval[24] = f(d__1) * .5 * pow(d__2, *alfa);
7125 for (i__ = 2; i__ <= 12; ++i__) {
7126 u = hlgth * x[i__ - 2];
7130 fval[i__ - 1] = f(d__1) * pow(d__2, *alfa);
7133 fval[isym - 1] = f(d__1) * pow(d__2, *alfa);
7137 factor = pow(hlgth, d__1);
7142 if (*integr == 2 || *integr == 4) {
7148 dqcheb_(x, fval, cheb12, cheb24);
7149 for (i__ = 1; i__ <= 13; ++i__) {
7150 res12 += cheb12[i__ - 1] * rj[i__];
7151 res24 += cheb24[i__ - 1] * rj[i__];
7154 for (i__ = 14; i__ <= 25; ++i__) {
7155 res24 += cheb24[i__ - 1] * rj[i__];
7164 dc = log(*br - *bl);
7165 *result = res24 * dc;
7166 *abserr = (d__1 = (res24 - res12) * dc, fabs(d__1));
7169 for (i__ = 1; i__ <= 13; ++i__) {
7170 res12 += cheb12[i__ - 1] * rh[i__];
7171 res24 += cheb24[i__ - 1] * rh[i__];
7174 for (i__ = 14; i__ <= 25; ++i__) {
7175 res24 += cheb24[i__ - 1] * rh[i__];
7185 fval[0] *= log(hlgth + fix);
7186 fval[12] *= log(fix);
7187 fval[24] *= log(fix - hlgth);
7188 for (i__ = 2; i__ <= 12; ++i__) {
7189 u = hlgth * x[i__ - 2];
7191 fval[i__ - 1] *= log(u + fix);
7192 fval[isym - 1] *= log(fix - u);
7195 dqcheb_(x, fval, cheb12, cheb24);
7199 for (i__ = 1; i__ <= 13; ++i__) {
7200 res12 += cheb12[i__ - 1] * rj[i__];
7201 res24 += cheb24[i__ - 1] * rj[i__];
7204 for (i__ = 14; i__ <= 25; ++i__) {
7205 res24 += cheb24[i__ - 1] * rj[i__];
7211 dc = log(*br - *bl);
7212 *result = res24 * dc;
7213 *abserr = (d__1 = (res24 - res12) * dc, fabs(d__1));
7219 for (i__ = 1; i__ <= 13; ++i__) {
7220 res12 += cheb12[i__ - 1] * rh[i__];
7221 res24 += cheb24[i__ - 1] * rh[i__];
7224 for (i__ = 14; i__ <= 25; ++i__) {
7225 res24 += cheb24[i__ - 1] * rh[i__];
7229 *result = (*result + res24) * factor;
7230 *abserr = (*abserr + (d__1 = res24 - res12, fabs(d__1))) * factor;
7235 void dqcheb_(
const double *x,
double *fval,
double *cheb12,
double *cheb24)
7238 double v[12], alam, alam1, alam2, part1, part2, part3;
7294 for (i__ = 1; i__ <= 12; ++i__) {
7296 v[i__ - 1] = fval[i__] - fval[j];
7297 fval[i__] += fval[j];
7300 alam1 = v[0] - v[8];
7301 alam2 = x[6] * (v[2] - v[6] - v[10]);
7302 cheb12[4] = alam1 + alam2;
7303 cheb12[10] = alam1 - alam2;
7304 alam1 = v[1] - v[7] - v[9];
7305 alam2 = v[3] - v[5] - v[11];
7306 alam = x[3] * alam1 + x[9] * alam2;
7307 cheb24[4] = cheb12[4] + alam;
7308 cheb24[22] = cheb12[4] - alam;
7309 alam = x[9] * alam1 - x[3] * alam2;
7310 cheb24[10] = cheb12[10] + alam;
7311 cheb24[16] = cheb12[10] - alam;
7312 part1 = x[4] * v[4];
7313 part2 = x[8] * v[8];
7314 part3 = x[6] * v[6];
7315 alam1 = v[0] + part1 + part2;
7316 alam2 = x[2] * v[2] + part3 + x[10] * v[10];
7317 cheb12[2] = alam1 + alam2;
7318 cheb12[12] = alam1 - alam2;
7319 alam = x[1] * v[1] + x[3] * v[3] + x[5] * v[5] + x[7] * v[7] + x[9] * v[9]
7321 cheb24[2] = cheb12[2] + alam;
7322 cheb24[24] = cheb12[2] - alam;
7323 alam = x[11] * v[1] - x[9] * v[3] + x[7] * v[5] - x[5] * v[7] + x[3] * v[
7325 cheb24[12] = cheb12[12] + alam;
7326 cheb24[14] = cheb12[12] - alam;
7327 alam1 = v[0] - part1 + part2;
7328 alam2 = x[10] * v[2] - part3 + x[2] * v[10];
7329 cheb12[6] = alam1 + alam2;
7330 cheb12[8] = alam1 - alam2;
7331 alam = x[5] * v[1] - x[9] * v[3] - x[1] * v[5] - x[11] * v[7] + x[3] * v[
7333 cheb24[6] = cheb12[6] + alam;
7334 cheb24[20] = cheb12[6] - alam;
7335 alam = x[7] * v[1] - x[3] * v[3] - x[11] * v[5] + x[1] * v[7] - x[9] * v[
7337 cheb24[8] = cheb12[8] + alam;
7338 cheb24[18] = cheb12[8] - alam;
7339 for (i__ = 1; i__ <= 6; ++i__) {
7341 v[i__ - 1] = fval[i__] - fval[j];
7342 fval[i__] += fval[j];
7345 alam1 = v[0] + x[8] * v[4];
7346 alam2 = x[4] * v[2];
7347 cheb12[3] = alam1 + alam2;
7348 cheb12[11] = alam1 - alam2;
7349 cheb12[7] = v[0] - v[4];
7350 alam = x[2] * v[1] + x[6] * v[3] + x[10] * v[5];
7351 cheb24[3] = cheb12[3] + alam;
7352 cheb24[23] = cheb12[3] - alam;
7353 alam = x[6] * (v[1] - v[3] - v[5]);
7354 cheb24[7] = cheb12[7] + alam;
7355 cheb24[19] = cheb12[7] - alam;
7356 alam = x[10] * v[1] - x[6] * v[3] + x[2] * v[5];
7357 cheb24[11] = cheb12[11] + alam;
7358 cheb24[15] = cheb12[11] - alam;
7359 for (i__ = 1; i__ <= 3; ++i__) {
7361 v[i__ - 1] = fval[i__] - fval[j];
7362 fval[i__] += fval[j];
7365 cheb12[5] = v[0] + x[8] * v[2];
7366 cheb12[9] = fval[1] - x[8] * fval[3];
7368 cheb24[5] = cheb12[5] + alam;
7369 cheb24[21] = cheb12[5] - alam;
7370 alam = x[8] * fval[2] - fval[4];
7371 cheb24[9] = cheb12[9] + alam;
7372 cheb24[17] = cheb12[9] - alam;
7373 cheb12[1] = fval[1] + fval[3];
7374 alam = fval[2] + fval[4];
7375 cheb24[1] = cheb12[1] + alam;
7376 cheb24[25] = cheb12[1] - alam;
7377 cheb12[13] = v[0] - v[2];
7378 cheb24[13] = cheb12[13];
7379 alam = .16666666666666666;
7380 for (i__ = 2; i__ <= 12; ++i__) {
7381 cheb12[i__] *= alam;
7387 for (i__ = 2; i__ <= 24; ++i__) {
7388 cheb24[i__] *= alam;
7391 cheb24[1] = alam * .5 * cheb24[1];
7392 cheb24[25] = alam * .5 * cheb24[25];
7396 void dqelg_(
long *n,
double *epstab,
double *result,
7397 double *abserr,
double *res3la,
long *nres)
7401 double d__1, d__2, d__3;
7405 double e0,
e1,
e2, e3;
7406 long k1, k2, k3, ib, ie;
7411 double err1, err2, err3, tol1, tol2, tol3;
7413 double e1abs, oflow, error;
7414 double delta1, delta2, delta3, epmach, epsinf;
7415 long newelm, limexp;
7501 *result = epstab[*n];
7506 epstab[*n + 2] = epstab[*n];
7507 newelm = (*n - 1) / 2;
7512 for (i__ = 1; i__ <= i__1; ++i__) {
7515 res = epstab[k1 + 2];
7521 err2 = fabs(delta2);
7524 tol2 =
max(d__1,e1abs) * epmach;
7526 err3 = fabs(delta3);
7528 d__1 = e1abs, d__2 = fabs(e0);
7529 tol3 =
max(d__1,d__2) * epmach;
7530 if (err2 > tol2 || err3 > tol3) {
7540 *abserr = err2 + err3;
7547 err1 = fabs(delta1);
7549 d__1 = e1abs, d__2 = fabs(e3);
7550 tol1 =
max(d__1,d__2) * epmach;
7555 if (err1 <= tol1 || err2 <= tol2 || err3 <= tol3) {
7558 ss = 1. / delta1 + 1. / delta2 - 1. / delta3;
7559 epsinf = (d__1 = ss *
e1, fabs(d__1));
7565 if (epsinf > 1e-4) {
7580 error = err2 + (d__1 = res -
e2, fabs(d__1)) + err3;
7581 if (error > *abserr) {
7594 *n = (limexp / 2 << 1) - 1;
7597 if (num / 2 << 1 == num) {
7602 for (i__ = 1; i__ <= i__1; ++i__) {
7604 epstab[ib] = epstab[ib2];
7611 indx = num - *n + 1;
7613 for (i__ = 1; i__ <= i__1; ++i__) {
7614 epstab[i__] = epstab[indx];
7622 res3la[*nres] = *result;
7629 *abserr = (d__1 = *result - res3la[3], fabs(d__1)) + (d__2 = *result -
7630 res3la[2], fabs(d__2)) + (d__3 = *result - res3la[1], fabs(d__3));
7631 res3la[1] = res3la[2];
7632 res3la[2] = res3la[3];
7633 res3la[3] = *result;
7636 d__1 = *abserr, d__2 = epmach * 5. * fabs(*result);
7637 *abserr =
max(d__1,d__2);
7641 void dqk15_(
const D_fp& f,
const double *a,
const double *b,
double *
7642 result,
double *abserr,
double *resabs,
double *resasc)
7646 double wg[4] = { .129484966168869693270611432679082,
7647 .27970539148927666790146777142378,
7648 .381830050505118944950369775488975,
7649 .417959183673469387755102040816327 };
7650 double xgk[8] = { .991455371120812639206854697526329,
7651 .949107912342758524526189684047851,
7652 .864864423359769072789712788640926,
7653 .741531185599394439863864773280788,
7654 .58608723546769113029414483825873,
7655 .405845151377397166906606412076961,
7656 .207784955007898467600689403773245,0. };
7657 double wgk[8] = { .02293532201052922496373200805897,
7658 .063092092629978553290700663189204,
7659 .104790010322250183839876322541518,
7660 .140653259715525918745189590510238,
7661 .16900472663926790282658342659855,
7662 .190350578064785409913256402421014,
7663 .204432940075298892414161999234649,
7664 .209482141084727828012999174891714 };
7667 double d__1, d__2, d__3;
7671 double fc, fv1[7], fv2[7];
7673 double absc, resg, resk, fsum, fval1, fval2;
7675 double hlgth, centr, reskh, uflow;
7676 double epmach, dhlgth;
7776 centr = (*a + *b) * .5;
7777 hlgth = (*b - *a) * .5;
7778 dhlgth = fabs(hlgth);
7786 *resabs = fabs(resk);
7787 for (j = 1; j <= 3; ++j) {
7789 absc = hlgth * xgk[jtw - 1];
7790 d__1 = centr - absc;
7792 d__1 = centr + absc;
7794 fv1[jtw - 1] = fval1;
7795 fv2[jtw - 1] = fval2;
7796 fsum = fval1 + fval2;
7797 resg += wg[j - 1] * fsum;
7798 resk += wgk[jtw - 1] * fsum;
7799 *resabs += wgk[jtw - 1] * (fabs(fval1) + fabs(fval2));
7802 for (j = 1; j <= 4; ++j) {
7803 jtwm1 = (j << 1) - 1;
7804 absc = hlgth * xgk[jtwm1 - 1];
7805 d__1 = centr - absc;
7807 d__1 = centr + absc;
7809 fv1[jtwm1 - 1] = fval1;
7810 fv2[jtwm1 - 1] = fval2;
7811 fsum = fval1 + fval2;
7812 resk += wgk[jtwm1 - 1] * fsum;
7813 *resabs += wgk[jtwm1 - 1] * (fabs(fval1) + fabs(fval2));
7817 *resasc = wgk[7] * (d__1 = fc - reskh, fabs(d__1));
7818 for (j = 1; j <= 7; ++j) {
7819 *resasc += wgk[j - 1] * ((d__1 = fv1[j - 1] - reskh, fabs(d__1)) + (
7820 d__2 = fv2[j - 1] - reskh, fabs(d__2)));
7823 *result = resk * hlgth;
7826 *abserr = (d__1 = (resk - resg) * hlgth, fabs(d__1));
7827 if (*resasc != 0. && *abserr != 0.) {
7829 d__3 = *abserr * 200. / *resasc;
7830 d__1 = 1., d__2 = pow(d__3,
c_b270);
7831 *abserr = *resasc *
min(d__1,d__2);
7833 if (*resabs > uflow / (epmach * 50.)) {
7835 d__1 = epmach * 50. * *resabs;
7836 *abserr =
max(d__1,*abserr);
7842 const double *a,
const double *b,
double *result,
double *abserr,
7843 double *resabs,
double *resasc)
7847 double wg[8] = { 0.,.129484966168869693270611432679082,0.,
7848 .27970539148927666790146777142378,0.,
7849 .381830050505118944950369775488975,0.,
7850 .417959183673469387755102040816327 };
7851 double xgk[8] = { .991455371120812639206854697526329,
7852 .949107912342758524526189684047851,
7853 .864864423359769072789712788640926,
7854 .741531185599394439863864773280788,
7855 .58608723546769113029414483825873,
7856 .405845151377397166906606412076961,
7857 .207784955007898467600689403773245,0. };
7858 double wgk[8] = { .02293532201052922496373200805897,
7859 .063092092629978553290700663189204,
7860 .104790010322250183839876322541518,
7861 .140653259715525918745189590510238,
7862 .16900472663926790282658342659855,
7863 .190350578064785409913256402421014,
7864 .204432940075298892414161999234649,
7865 .209482141084727828012999174891714 };
7868 double d__1, d__2, d__3;
7872 double fc, fv1[7], fv2[7], absc, dinf, resg, resk, fsum, absc1,
7873 absc2, fval1, fval2, hlgth, centr, reskh, uflow;
7874 double tabsc1, tabsc2, epmach;
7990 dinf = (double)
min(1,*inf);
7992 centr = (*a + *b) * .5;
7993 hlgth = (*b - *a) * .5;
7994 tabsc1 = *boun + dinf * (1. - centr) / centr;
8000 fc = fval1 / centr / centr;
8007 *resabs = fabs(resk);
8008 for (j = 1; j <= 7; ++j) {
8009 absc = hlgth * xgk[j - 1];
8010 absc1 = centr - absc;
8011 absc2 = centr + absc;
8012 tabsc1 = *boun + dinf * (1. - absc1) / absc1;
8013 tabsc2 = *boun + dinf * (1. - absc2) / absc2;
8024 fval1 = fval1 / absc1 / absc1;
8025 fval2 = fval2 / absc2 / absc2;
8028 fsum = fval1 + fval2;
8029 resg += wg[j - 1] * fsum;
8030 resk += wgk[j - 1] * fsum;
8031 *resabs += wgk[j - 1] * (fabs(fval1) + fabs(fval2));
8035 *resasc = wgk[7] * (d__1 = fc - reskh, fabs(d__1));
8036 for (j = 1; j <= 7; ++j) {
8037 *resasc += wgk[j - 1] * ((d__1 = fv1[j - 1] - reskh, fabs(d__1)) + (
8038 d__2 = fv2[j - 1] - reskh, fabs(d__2)));
8041 *result = resk * hlgth;
8044 *abserr = (d__1 = (resk - resg) * hlgth, fabs(d__1));
8045 if (*resasc != 0. && *abserr != 0.) {
8047 d__3 = *abserr * 200. / *resasc;
8048 d__1 = 1., d__2 = pow(d__3,
c_b270);
8049 *abserr = *resasc *
min(d__1,d__2);
8051 if (*resabs > uflow / (epmach * 50.)) {
8053 d__1 = epmach * 50. * *resabs;
8054 *abserr =
max(d__1,*abserr);
8060 const double *p3,
const double *p4,
const long *kp,
8061 const double *a,
const double *b,
8062 double *result,
double *abserr,
double *resabs,
double *resasc)
8066 double xgk[8] = { .9914553711208126,.9491079123427585,
8067 .8648644233597691,.7415311855993944,.5860872354676911,
8068 .4058451513773972,.2077849550078985,0. };
8069 double wgk[8] = { .02293532201052922,.06309209262997855,
8070 .1047900103222502,.1406532597155259,.1690047266392679,
8071 .1903505780647854,.2044329400752989,.2094821410847278 };
8072 double wg[4] = { .1294849661688697,.2797053914892767,
8073 .3818300505051889,.4179591836734694 };
8076 double d__1, d__2, d__3;
8080 double fc, fv1[7], fv2[7];
8082 double absc, resg, resk, fsum, absc1, absc2, fval1, fval2;
8084 double hlgth, centr, reskh, uflow;
8085 double epmach, dhlgth;
8192 centr = (*a + *b) * .5;
8193 hlgth = (*b - *a) * .5;
8194 dhlgth = fabs(hlgth);
8199 fc = f(centr) * w(¢r, p1, p2, p3, p4, kp);
8202 *resabs = fabs(resk);
8203 for (j = 1; j <= 3; ++j) {
8205 absc = hlgth * xgk[jtw - 1];
8206 absc1 = centr - absc;
8207 absc2 = centr + absc;
8208 fval1 = f(absc1) * w(&absc1, p1, p2, p3, p4, kp);
8209 fval2 = f(absc2) * w(&absc2, p1, p2, p3, p4, kp);
8210 fv1[jtw - 1] = fval1;
8211 fv2[jtw - 1] = fval2;
8212 fsum = fval1 + fval2;
8213 resg += wg[j - 1] * fsum;
8214 resk += wgk[jtw - 1] * fsum;
8215 *resabs += wgk[jtw - 1] * (fabs(fval1) + fabs(fval2));
8218 for (j = 1; j <= 4; ++j) {
8219 jtwm1 = (j << 1) - 1;
8220 absc = hlgth * xgk[jtwm1 - 1];
8221 absc1 = centr - absc;
8222 absc2 = centr + absc;
8223 fval1 = f(absc1) * w(&absc1, p1, p2, p3, p4, kp);
8224 fval2 = f(absc2) * w(&absc2, p1, p2, p3, p4, kp);
8225 fv1[jtwm1 - 1] = fval1;
8226 fv2[jtwm1 - 1] = fval2;
8227 fsum = fval1 + fval2;
8228 resk += wgk[jtwm1 - 1] * fsum;
8229 *resabs += wgk[jtwm1 - 1] * (fabs(fval1) + fabs(fval2));
8233 *resasc = wgk[7] * (d__1 = fc - reskh, fabs(d__1));
8234 for (j = 1; j <= 7; ++j) {
8235 *resasc += wgk[j - 1] * ((d__1 = fv1[j - 1] - reskh, fabs(d__1)) + (
8236 d__2 = fv2[j - 1] - reskh, fabs(d__2)));
8239 *result = resk * hlgth;
8242 *abserr = (d__1 = (resk - resg) * hlgth, fabs(d__1));
8243 if (*resasc != 0. && *abserr != 0.) {
8245 d__3 = *abserr * 200. / *resasc;
8246 d__1 = 1., d__2 = pow(d__3,
c_b270);
8247 *abserr = *resasc *
min(d__1,d__2);
8249 if (*resabs > uflow / (epmach * 50.)) {
8251 d__1 = epmach * 50. * *resabs;
8252 *abserr =
max(d__1,*abserr);
8257 void dqk21_(
const D_fp& f,
const double *a,
const double *b,
double *
8258 result,
double *abserr,
double *resabs,
double *resasc)
8262 double wg[5] = { .066671344308688137593568809893332,
8263 .149451349150580593145776339657697,
8264 .219086362515982043995534934228163,
8265 .269266719309996355091226921569469,
8266 .295524224714752870173892994651338 };
8267 double xgk[11] = { .995657163025808080735527280689003,
8268 .973906528517171720077964012084452,
8269 .930157491355708226001207180059508,
8270 .865063366688984510732096688423493,
8271 .780817726586416897063717578345042,
8272 .679409568299024406234327365114874,
8273 .562757134668604683339000099272694,
8274 .433395394129247190799265943165784,
8275 .294392862701460198131126603103866,
8276 .14887433898163121088482600112972,0. };
8277 double wgk[11] = { .011694638867371874278064396062192,
8278 .03255816230796472747881897245939,
8279 .05475589657435199603138130024458,
8280 .07503967481091995276704314091619,
8281 .093125454583697605535065465083366,
8282 .109387158802297641899210590325805,
8283 .123491976262065851077958109831074,
8284 .134709217311473325928054001771707,
8285 .142775938577060080797094273138717,
8286 .147739104901338491374841515972068,
8287 .149445554002916905664936468389821 };
8290 double d__1, d__2, d__3;
8294 double fc, fv1[10], fv2[10];
8296 double absc, resg, resk, fsum, fval1, fval2;
8298 double hlgth, centr, reskh, uflow;
8299 double epmach, dhlgth;
8400 centr = (*a + *b) * .5;
8401 hlgth = (*b - *a) * .5;
8402 dhlgth = fabs(hlgth);
8409 resk = wgk[10] * fc;
8410 *resabs = fabs(resk);
8411 for (j = 1; j <= 5; ++j) {
8413 absc = hlgth * xgk[jtw - 1];
8414 d__1 = centr - absc;
8416 d__1 = centr + absc;
8418 fv1[jtw - 1] = fval1;
8419 fv2[jtw - 1] = fval2;
8420 fsum = fval1 + fval2;
8421 resg += wg[j - 1] * fsum;
8422 resk += wgk[jtw - 1] * fsum;
8423 *resabs += wgk[jtw - 1] * (fabs(fval1) + fabs(fval2));
8426 for (j = 1; j <= 5; ++j) {
8427 jtwm1 = (j << 1) - 1;
8428 absc = hlgth * xgk[jtwm1 - 1];
8429 d__1 = centr - absc;
8431 d__1 = centr + absc;
8433 fv1[jtwm1 - 1] = fval1;
8434 fv2[jtwm1 - 1] = fval2;
8435 fsum = fval1 + fval2;
8436 resk += wgk[jtwm1 - 1] * fsum;
8437 *resabs += wgk[jtwm1 - 1] * (fabs(fval1) + fabs(fval2));
8441 *resasc = wgk[10] * (d__1 = fc - reskh, fabs(d__1));
8442 for (j = 1; j <= 10; ++j) {
8443 *resasc += wgk[j - 1] * ((d__1 = fv1[j - 1] - reskh, fabs(d__1)) + (
8444 d__2 = fv2[j - 1] - reskh, fabs(d__2)));
8447 *result = resk * hlgth;
8450 *abserr = (d__1 = (resk - resg) * hlgth, fabs(d__1));
8451 if (*resasc != 0. && *abserr != 0.) {
8453 d__3 = *abserr * 200. / *resasc;
8454 d__1 = 1., d__2 = pow(d__3,
c_b270);
8455 *abserr = *resasc *
min(d__1,d__2);
8457 if (*resabs > uflow / (epmach * 50.)) {
8459 d__1 = epmach * 50. * *resabs;
8460 *abserr =
max(d__1,*abserr);
8465 void dqk31_(
const D_fp& f,
const double *a,
const double *b,
double *
8466 result,
double *abserr,
double *resabs,
double *resasc)
8470 double wg[8] = { .030753241996117268354628393577204,
8471 .070366047488108124709267416450667,
8472 .107159220467171935011869546685869,
8473 .139570677926154314447804794511028,
8474 .166269205816993933553200860481209,
8475 .186161000015562211026800561866423,
8476 .198431485327111576456118326443839,
8477 .202578241925561272880620199967519 };
8478 double xgk[16] = { .998002298693397060285172840152271,
8479 .987992518020485428489565718586613,
8480 .967739075679139134257347978784337,
8481 .937273392400705904307758947710209,
8482 .897264532344081900882509656454496,
8483 .848206583410427216200648320774217,
8484 .790418501442465932967649294817947,
8485 .724417731360170047416186054613938,
8486 .650996741297416970533735895313275,
8487 .570972172608538847537226737253911,
8488 .485081863640239680693655740232351,
8489 .394151347077563369897207370981045,
8490 .299180007153168812166780024266389,
8491 .201194093997434522300628303394596,
8492 .101142066918717499027074231447392,0. };
8493 double wgk[16] = { .005377479872923348987792051430128,
8494 .015007947329316122538374763075807,
8495 .025460847326715320186874001019653,
8496 .03534636079137584622203794847836,
8497 .04458975132476487660822729937328,
8498 .05348152469092808726534314723943,
8499 .062009567800670640285139230960803,
8500 .069854121318728258709520077099147,
8501 .076849680757720378894432777482659,
8502 .083080502823133021038289247286104,
8503 .088564443056211770647275443693774,
8504 .093126598170825321225486872747346,
8505 .096642726983623678505179907627589,
8506 .099173598721791959332393173484603,
8507 .10076984552387559504494666261757,
8508 .101330007014791549017374792767493 };
8511 double d__1, d__2, d__3;
8515 double fc, fv1[15], fv2[15];
8517 double absc, resg, resk, fsum, fval1, fval2;
8519 double hlgth, centr, reskh, uflow;
8520 double epmach, dhlgth;
8617 centr = (*a + *b) * .5;
8618 hlgth = (*b - *a) * .5;
8619 dhlgth = fabs(hlgth);
8626 resk = wgk[15] * fc;
8627 *resabs = fabs(resk);
8628 for (j = 1; j <= 7; ++j) {
8630 absc = hlgth * xgk[jtw - 1];
8631 d__1 = centr - absc;
8633 d__1 = centr + absc;
8635 fv1[jtw - 1] = fval1;
8636 fv2[jtw - 1] = fval2;
8637 fsum = fval1 + fval2;
8638 resg += wg[j - 1] * fsum;
8639 resk += wgk[jtw - 1] * fsum;
8640 *resabs += wgk[jtw - 1] * (fabs(fval1) + fabs(fval2));
8643 for (j = 1; j <= 8; ++j) {
8644 jtwm1 = (j << 1) - 1;
8645 absc = hlgth * xgk[jtwm1 - 1];
8646 d__1 = centr - absc;
8648 d__1 = centr + absc;
8650 fv1[jtwm1 - 1] = fval1;
8651 fv2[jtwm1 - 1] = fval2;
8652 fsum = fval1 + fval2;
8653 resk += wgk[jtwm1 - 1] * fsum;
8654 *resabs += wgk[jtwm1 - 1] * (fabs(fval1) + fabs(fval2));
8658 *resasc = wgk[15] * (d__1 = fc - reskh, fabs(d__1));
8659 for (j = 1; j <= 15; ++j) {
8660 *resasc += wgk[j - 1] * ((d__1 = fv1[j - 1] - reskh, fabs(d__1)) + (
8661 d__2 = fv2[j - 1] - reskh, fabs(d__2)));
8664 *result = resk * hlgth;
8667 *abserr = (d__1 = (resk - resg) * hlgth, fabs(d__1));
8668 if (*resasc != 0. && *abserr != 0.) {
8670 d__3 = *abserr * 200. / *resasc;
8671 d__1 = 1., d__2 = pow(d__3,
c_b270);
8672 *abserr = *resasc *
min(d__1,d__2);
8674 if (*resabs > uflow / (epmach * 50.)) {
8676 d__1 = epmach * 50. * *resabs;
8677 *abserr =
max(d__1,*abserr);
8682 void dqk41_(
const D_fp& f,
const double *a,
const double *b,
double *
8683 result,
double *abserr,
double *resabs,
double *resasc)
8687 double wg[10] = { .017614007139152118311861962351853,
8688 .040601429800386941331039952274932,
8689 .062672048334109063569506535187042,
8690 .083276741576704748724758143222046,
8691 .10193011981724043503675013548035,
8692 .118194531961518417312377377711382,
8693 .131688638449176626898494499748163,
8694 .142096109318382051329298325067165,
8695 .149172986472603746787828737001969,
8696 .152753387130725850698084331955098 };
8697 double xgk[21] = { .998859031588277663838315576545863,
8698 .99312859918509492478612238847132,
8699 .981507877450250259193342994720217,
8700 .963971927277913791267666131197277,
8701 .940822633831754753519982722212443,
8702 .912234428251325905867752441203298,
8703 .878276811252281976077442995113078,
8704 .839116971822218823394529061701521,
8705 .795041428837551198350638833272788,
8706 .746331906460150792614305070355642,
8707 .693237656334751384805490711845932,
8708 .636053680726515025452836696226286,
8709 .575140446819710315342946036586425,
8710 .510867001950827098004364050955251,
8711 .44359317523872510319999221349264,
8712 .373706088715419560672548177024927,
8713 .301627868114913004320555356858592,
8714 .227785851141645078080496195368575,
8715 .152605465240922675505220241022678,
8716 .076526521133497333754640409398838,0. };
8717 double wgk[21] = { .003073583718520531501218293246031,
8718 .008600269855642942198661787950102,
8719 .014626169256971252983787960308868,
8720 .020388373461266523598010231432755,
8721 .025882133604951158834505067096153,
8722 .031287306777032798958543119323801,
8723 .036600169758200798030557240707211,
8724 .041668873327973686263788305936895,
8725 .046434821867497674720231880926108,
8726 .050944573923728691932707670050345,
8727 .055195105348285994744832372419777,
8728 .059111400880639572374967220648594,
8729 .062653237554781168025870122174255,
8730 .065834597133618422111563556969398,
8731 .068648672928521619345623411885368,
8732 .07105442355344406830579036172321,
8733 .073030690332786667495189417658913,
8734 .074582875400499188986581418362488,
8735 .075704497684556674659542775376617,
8736 .076377867672080736705502835038061,
8737 .076600711917999656445049901530102 };
8740 double d__1, d__2, d__3;
8744 double fc, fv1[20], fv2[20];
8746 double absc, resg, resk, fsum, fval1, fval2;
8748 double hlgth, centr, reskh, uflow;
8749 double epmach, dhlgth;
8850 centr = (*a + *b) * .5;
8851 hlgth = (*b - *a) * .5;
8852 dhlgth = fabs(hlgth);
8859 resk = wgk[20] * fc;
8860 *resabs = fabs(resk);
8861 for (j = 1; j <= 10; ++j) {
8863 absc = hlgth * xgk[jtw - 1];
8864 d__1 = centr - absc;
8866 d__1 = centr + absc;
8868 fv1[jtw - 1] = fval1;
8869 fv2[jtw - 1] = fval2;
8870 fsum = fval1 + fval2;
8871 resg += wg[j - 1] * fsum;
8872 resk += wgk[jtw - 1] * fsum;
8873 *resabs += wgk[jtw - 1] * (fabs(fval1) + fabs(fval2));
8876 for (j = 1; j <= 10; ++j) {
8877 jtwm1 = (j << 1) - 1;
8878 absc = hlgth * xgk[jtwm1 - 1];
8879 d__1 = centr - absc;
8881 d__1 = centr + absc;
8883 fv1[jtwm1 - 1] = fval1;
8884 fv2[jtwm1 - 1] = fval2;
8885 fsum = fval1 + fval2;
8886 resk += wgk[jtwm1 - 1] * fsum;
8887 *resabs += wgk[jtwm1 - 1] * (fabs(fval1) + fabs(fval2));
8891 *resasc = wgk[20] * (d__1 = fc - reskh, fabs(d__1));
8892 for (j = 1; j <= 20; ++j) {
8893 *resasc += wgk[j - 1] * ((d__1 = fv1[j - 1] - reskh, fabs(d__1)) + (
8894 d__2 = fv2[j - 1] - reskh, fabs(d__2)));
8897 *result = resk * hlgth;
8900 *abserr = (d__1 = (resk - resg) * hlgth, fabs(d__1));
8901 if (*resasc != 0. && *abserr != 0.) {
8903 d__3 = *abserr * 200. / *resasc;
8904 d__1 = 1., d__2 = pow(d__3,
c_b270);
8905 *abserr = *resasc *
min(d__1,d__2);
8907 if (*resabs > uflow / (epmach * 50.)) {
8909 d__1 = epmach * 50. * *resabs;
8910 *abserr =
max(d__1,*abserr);
8915 void dqk51_(
const D_fp& f,
const double *a,
const double *b,
double *
8916 result,
double *abserr,
double *resabs,
double *resasc)
8920 double wg[13] = { .011393798501026287947902964113235,
8921 .026354986615032137261901815295299,
8922 .040939156701306312655623487711646,
8923 .054904695975835191925936891540473,
8924 .068038333812356917207187185656708,
8925 .080140700335001018013234959669111,
8926 .091028261982963649811497220702892,
8927 .100535949067050644202206890392686,
8928 .108519624474263653116093957050117,
8929 .114858259145711648339325545869556,
8930 .119455763535784772228178126512901,
8931 .122242442990310041688959518945852,
8932 .12317605372671545120390287307905 };
8933 double xgk[26] = { .999262104992609834193457486540341,
8934 .995556969790498097908784946893902,
8935 .988035794534077247637331014577406,
8936 .976663921459517511498315386479594,
8937 .961614986425842512418130033660167,
8938 .942974571228974339414011169658471,
8939 .920747115281701561746346084546331,
8940 .894991997878275368851042006782805,
8941 .86584706529327559544899696958834,
8942 .83344262876083400142102110869357,
8943 .797873797998500059410410904994307,
8944 .759259263037357630577282865204361,
8945 .717766406813084388186654079773298,
8946 .673566368473468364485120633247622,
8947 .626810099010317412788122681624518,
8948 .577662930241222967723689841612654,
8949 .52632528433471918259962377815801,
8950 .473002731445714960522182115009192,
8951 .417885382193037748851814394594572,
8952 .361172305809387837735821730127641,
8953 .303089538931107830167478909980339,
8954 .243866883720988432045190362797452,
8955 .183718939421048892015969888759528,
8956 .122864692610710396387359818808037,
8957 .061544483005685078886546392366797,0. };
8958 double wgk[26] = { .001987383892330315926507851882843,
8959 .005561932135356713758040236901066,
8960 .009473973386174151607207710523655,
8961 .013236229195571674813656405846976,
8962 .016847817709128298231516667536336,
8963 .020435371145882835456568292235939,
8964 .024009945606953216220092489164881,
8965 .027475317587851737802948455517811,
8966 .030792300167387488891109020215229,
8967 .034002130274329337836748795229551,
8968 .03711627148341554356033062536762,
8969 .040083825504032382074839284467076,
8970 .042872845020170049476895792439495,
8971 .04550291304992178890987058475266,
8972 .047982537138836713906392255756915,
8973 .05027767908071567196332525943344,
8974 .052362885806407475864366712137873,
8975 .054251129888545490144543370459876,
8976 .055950811220412317308240686382747,
8977 .057437116361567832853582693939506,
8978 .058689680022394207961974175856788,
8979 .059720340324174059979099291932562,
8980 .060539455376045862945360267517565,
8981 .061128509717053048305859030416293,
8982 .061471189871425316661544131965264,
8983 .061580818067832935078759824240066 };
8986 double d__1, d__2, d__3;
8990 double fc, fv1[25], fv2[25];
8992 double absc, resg, resk, fsum, fval1, fval2;
8994 double hlgth, centr, reskh, uflow;
8995 double epmach, dhlgth;
9096 centr = (*a + *b) * .5;
9097 hlgth = (*b - *a) * .5;
9098 dhlgth = fabs(hlgth);
9105 resk = wgk[25] * fc;
9106 *resabs = fabs(resk);
9107 for (j = 1; j <= 12; ++j) {
9109 absc = hlgth * xgk[jtw - 1];
9110 d__1 = centr - absc;
9112 d__1 = centr + absc;
9114 fv1[jtw - 1] = fval1;
9115 fv2[jtw - 1] = fval2;
9116 fsum = fval1 + fval2;
9117 resg += wg[j - 1] * fsum;
9118 resk += wgk[jtw - 1] * fsum;
9119 *resabs += wgk[jtw - 1] * (fabs(fval1) + fabs(fval2));
9122 for (j = 1; j <= 13; ++j) {
9123 jtwm1 = (j << 1) - 1;
9124 absc = hlgth * xgk[jtwm1 - 1];
9125 d__1 = centr - absc;
9127 d__1 = centr + absc;
9129 fv1[jtwm1 - 1] = fval1;
9130 fv2[jtwm1 - 1] = fval2;
9131 fsum = fval1 + fval2;
9132 resk += wgk[jtwm1 - 1] * fsum;
9133 *resabs += wgk[jtwm1 - 1] * (fabs(fval1) + fabs(fval2));
9137 *resasc = wgk[25] * (d__1 = fc - reskh, fabs(d__1));
9138 for (j = 1; j <= 25; ++j) {
9139 *resasc += wgk[j - 1] * ((d__1 = fv1[j - 1] - reskh, fabs(d__1)) + (
9140 d__2 = fv2[j - 1] - reskh, fabs(d__2)));
9143 *result = resk * hlgth;
9146 *abserr = (d__1 = (resk - resg) * hlgth, fabs(d__1));
9147 if (*resasc != 0. && *abserr != 0.) {
9149 d__3 = *abserr * 200. / *resasc;
9150 d__1 = 1., d__2 = pow(d__3,
c_b270);
9151 *abserr = *resasc *
min(d__1,d__2);
9153 if (*resabs > uflow / (epmach * 50.)) {
9155 d__1 = epmach * 50. * *resabs;
9156 *abserr =
max(d__1,*abserr);
9161 void dqk61_(
const D_fp& f,
const double *a,
const double *b,
double *
9162 result,
double *abserr,
double *resabs,
double *resasc)
9166 double wg[15] = { .007968192496166605615465883474674,
9167 .018466468311090959142302131912047,
9168 .028784707883323369349719179611292,
9169 .038799192569627049596801936446348,
9170 .048402672830594052902938140422808,
9171 .057493156217619066481721689402056,
9172 .065974229882180495128128515115962,
9173 .073755974737705206268243850022191,
9174 .08075589522942021535469493846053,
9175 .086899787201082979802387530715126,
9176 .092122522237786128717632707087619,
9177 .09636873717464425963946862635181,
9178 .099593420586795267062780282103569,
9179 .101762389748405504596428952168554,
9180 .102852652893558840341285636705415 };
9181 double xgk[31] = { .999484410050490637571325895705811,
9182 .996893484074649540271630050918695,
9183 .991630996870404594858628366109486,
9184 .983668123279747209970032581605663,
9185 .973116322501126268374693868423707,
9186 .960021864968307512216871025581798,
9187 .944374444748559979415831324037439,
9188 .926200047429274325879324277080474,
9189 .905573307699907798546522558925958,
9190 .882560535792052681543116462530226,
9191 .857205233546061098958658510658944,
9192 .829565762382768397442898119732502,
9193 .799727835821839083013668942322683,
9194 .767777432104826194917977340974503,
9195 .733790062453226804726171131369528,
9196 .69785049479331579693229238802664,
9197 .660061064126626961370053668149271,
9198 .620526182989242861140477556431189,
9199 .57934523582636169175602493217254,
9200 .536624148142019899264169793311073,
9201 .492480467861778574993693061207709,
9202 .447033769538089176780609900322854,
9203 .400401254830394392535476211542661,
9204 .352704725530878113471037207089374,
9205 .304073202273625077372677107199257,
9206 .254636926167889846439805129817805,
9207 .204525116682309891438957671002025,
9208 .153869913608583546963794672743256,
9209 .102806937966737030147096751318001,
9210 .051471842555317695833025213166723,0. };
9211 double wgk[31] = { .00138901369867700762455159122676,
9212 .003890461127099884051267201844516,
9213 .00663070391593129217331982636975,
9214 .009273279659517763428441146892024,
9215 .011823015253496341742232898853251,
9216 .01436972950704580481245143244358,
9217 .016920889189053272627572289420322,
9218 .019414141193942381173408951050128,
9219 .021828035821609192297167485738339,
9220 .024191162078080601365686370725232,
9221 .026509954882333101610601709335075,
9222 .028754048765041292843978785354334,
9223 .030907257562387762472884252943092,
9224 .032981447057483726031814191016854,
9225 .034979338028060024137499670731468,
9226 .036882364651821229223911065617136,
9227 .038678945624727592950348651532281,
9228 .040374538951535959111995279752468,
9229 .04196981021516424614714754128597,
9230 .043452539701356069316831728117073,
9231 .044814800133162663192355551616723,
9232 .046059238271006988116271735559374,
9233 .047185546569299153945261478181099,
9234 .048185861757087129140779492298305,
9235 .049055434555029778887528165367238,
9236 .049795683427074206357811569379942,
9237 .050405921402782346840893085653585,
9238 .050881795898749606492297473049805,
9239 .051221547849258772170656282604944,
9240 .051426128537459025933862879215781,
9241 .051494729429451567558340433647099 };
9244 double d__1, d__2, d__3;
9248 double fc, fv1[30], fv2[30];
9250 double resg, resk, fsum, fval1, fval2;
9252 double dabsc, hlgth, centr, reskh, uflow;
9253 double epmach, dhlgth;
9352 centr = (*b + *a) * .5;
9353 hlgth = (*b - *a) * .5;
9354 dhlgth = fabs(hlgth);
9362 resk = wgk[30] * fc;
9363 *resabs = fabs(resk);
9364 for (j = 1; j <= 15; ++j) {
9366 dabsc = hlgth * xgk[jtw - 1];
9367 d__1 = centr - dabsc;
9369 d__1 = centr + dabsc;
9371 fv1[jtw - 1] = fval1;
9372 fv2[jtw - 1] = fval2;
9373 fsum = fval1 + fval2;
9374 resg += wg[j - 1] * fsum;
9375 resk += wgk[jtw - 1] * fsum;
9376 *resabs += wgk[jtw - 1] * (fabs(fval1) + fabs(fval2));
9379 for (j = 1; j <= 15; ++j) {
9380 jtwm1 = (j << 1) - 1;
9381 dabsc = hlgth * xgk[jtwm1 - 1];
9382 d__1 = centr - dabsc;
9384 d__1 = centr + dabsc;
9386 fv1[jtwm1 - 1] = fval1;
9387 fv2[jtwm1 - 1] = fval2;
9388 fsum = fval1 + fval2;
9389 resk += wgk[jtwm1 - 1] * fsum;
9390 *resabs += wgk[jtwm1 - 1] * (fabs(fval1) + fabs(fval2));
9394 *resasc = wgk[30] * (d__1 = fc - reskh, fabs(d__1));
9395 for (j = 1; j <= 30; ++j) {
9396 *resasc += wgk[j - 1] * ((d__1 = fv1[j - 1] - reskh, fabs(d__1)) + (
9397 d__2 = fv2[j - 1] - reskh, fabs(d__2)));
9400 *result = resk * hlgth;
9403 *abserr = (d__1 = (resk - resg) * hlgth, fabs(d__1));
9404 if (*resasc != 0. && *abserr != 0.) {
9406 d__3 = *abserr * 200. / *resasc;
9407 d__1 = 1., d__2 = pow(d__3,
c_b270);
9408 *abserr = *resasc *
min(d__1,d__2);
9410 if (*resabs > uflow / (epmach * 50.)) {
9412 d__1 = epmach * 50. * *resabs;
9413 *abserr =
max(d__1,*abserr);
9418 void dqmomo_(
const double *alfa,
const double *beta,
double *
9419 ri,
double *rj,
double *rg,
double *rh,
const long *integr)
9425 double anm1, ralf, rbet, alfp1, alfp2, betp1, betp2;
9498 ralf = pow(
c_b320, alfp1);
9499 rbet = pow(
c_b320, betp1);
9503 ri[1] = ralf / alfp1;
9504 rj[1] = rbet / betp1;
9505 ri[2] = ri[1] * *alfa / alfp2;
9506 rj[2] = rj[1] * *beta / betp2;
9509 for (i__ = 3; i__ <= 25; ++i__) {
9510 ri[i__] = -(ralf + an * (an - alfp2) * ri[i__ - 1]) / (anm1 * (an +
9512 rj[i__] = -(rbet + an * (an - betp2) * rj[i__ - 1]) / (anm1 * (an +
9527 rg[1] = -ri[1] / alfp1;
9528 rg[2] = -(ralf + ralf) / (alfp2 * alfp2) - rg[1];
9532 for (i__ = 3; i__ <= 25; ++i__) {
9533 rg[i__] = -(an * (an - alfp2) * rg[im1] - an * ri[im1] + anm1 * ri[
9534 i__]) / (anm1 * (an + alfp1));
9547 rh[1] = -rj[1] / betp1;
9548 rh[2] = -(rbet + rbet) / (betp2 * betp2) - rh[1];
9552 for (i__ = 3; i__ <= 25; ++i__) {
9553 rh[i__] = -(an * (an - betp2) * rh[im1] - an * rj[im1] + anm1 * rj[
9554 i__]) / (anm1 * (an + betp1));
9560 for (i__ = 2; i__ <= 25; i__ += 2) {
9565 for (i__ = 2; i__ <= 25; i__ += 2) {
9573 void dqng_(
const D_fp& f,
const double *a,
const double *b,
const double *epsabs,
9574 const double *epsrel,
double *result,
double *abserr,
9575 long *neval,
long *ier)
9579 double x1[5] = { .973906528517171720077964012084452,
9580 .865063366688984510732096688423493,
9581 .679409568299024406234327365114874,
9582 .433395394129247190799265943165784,
9583 .14887433898163121088482600112972 };
9584 double w10[5] = { .066671344308688137593568809893332,
9585 .149451349150580593145776339657697,
9586 .219086362515982043995534934228163,
9587 .269266719309996355091226921569469,
9588 .295524224714752870173892994651338 };
9589 double x2[5] = { .995657163025808080735527280689003,
9590 .930157491355708226001207180059508,
9591 .780817726586416897063717578345042,
9592 .562757134668604683339000099272694,
9593 .294392862701460198131126603103866 };
9594 double w21a[5] = { .03255816230796472747881897245939,
9595 .07503967481091995276704314091619,
9596 .109387158802297641899210590325805,
9597 .134709217311473325928054001771707,
9598 .147739104901338491374841515972068 };
9599 double w21b[6] = { .011694638867371874278064396062192,
9600 .05475589657435199603138130024458,
9601 .093125454583697605535065465083366,
9602 .123491976262065851077958109831074,
9603 .142775938577060080797094273138717,
9604 .149445554002916905664936468389821 };
9605 double x3[11] = { .999333360901932081394099323919911,
9606 .987433402908088869795961478381209,
9607 .954807934814266299257919200290473,
9608 .900148695748328293625099494069092,
9609 .82519831498311415084706673258852,
9610 .732148388989304982612354848755461,
9611 .622847970537725238641159120344323,
9612 .499479574071056499952214885499755,
9613 .364901661346580768043989548502644,
9614 .222254919776601296498260928066212,
9615 .074650617461383322043914435796506 };
9616 double w43a[10] = { .016296734289666564924281974617663,
9617 .037522876120869501461613795898115,
9618 .054694902058255442147212685465005,
9619 .067355414609478086075553166302174,
9620 .073870199632393953432140695251367,
9621 .005768556059769796184184327908655,
9622 .027371890593248842081276069289151,
9623 .046560826910428830743339154433824,
9624 .061744995201442564496240336030883,
9625 .071387267268693397768559114425516 };
9626 double w43b[12] = { .001844477640212414100389106552965,
9627 .010798689585891651740465406741293,
9628 .021895363867795428102523123075149,
9629 .032597463975345689443882222526137,
9630 .042163137935191811847627924327955,
9631 .050741939600184577780189020092084,
9632 .058379395542619248375475369330206,
9633 .064746404951445885544689259517511,
9634 .069566197912356484528633315038405,
9635 .072824441471833208150939535192842,
9636 .074507751014175118273571813842889,
9637 .074722147517403005594425168280423 };
9638 double x4[22] = { .999902977262729234490529830591582,
9639 .99798989598667874542749632236596,
9640 .992175497860687222808523352251425,
9641 .981358163572712773571916941623894,
9642 .965057623858384619128284110607926,
9643 .943167613133670596816416634507426,
9644 .91580641468550720959182643072005,
9645 .883221657771316501372117548744163,
9646 .845710748462415666605902011504855,
9647 .803557658035230982788739474980964,
9648 .75700573068549555832894279343202,
9649 .70627320978732181982409427474084,
9650 .651589466501177922534422205016736,
9651 .593223374057961088875273770349144,
9652 .531493605970831932285268948562671,
9653 .46676362304202284487196678165927,
9654 .399424847859218804732101665817923,
9655 .329874877106188288265053371824597,
9656 .258503559202161551802280975429025,
9657 .185695396568346652015917141167606,
9658 .111842213179907468172398359241362,
9659 .037352123394619870814998165437704 };
9660 double w87a[21] = { .00814837738414917290000287844819,
9661 .018761438201562822243935059003794,
9662 .027347451050052286161582829741283,
9663 .033677707311637930046581056957588,
9664 .036935099820427907614589586742499,
9665 .002884872430211530501334156248695,
9666 .013685946022712701888950035273128,
9667 .023280413502888311123409291030404,
9668 .030872497611713358675466394126442,
9669 .035693633639418770719351355457044,
9670 9.15283345202241360843392549948e-4,
9671 .005399280219300471367738743391053,
9672 .010947679601118931134327826856808,
9673 .01629873169678733526266570322328,
9674 .02108156888920383511243306018819,
9675 .02537096976925382724346799983171,
9676 .02918969775647575250144615408492,
9677 .032373202467202789685788194889595,
9678 .034783098950365142750781997949596,
9679 .036412220731351787562801163687577,
9680 .037253875503047708539592001191226 };
9681 double w87b[23] = { 2.74145563762072350016527092881e-4,
9682 .001807124155057942948341311753254,
9683 .00409686928275916486445807068348,
9684 .006758290051847378699816577897424,
9685 .009549957672201646536053581325377,
9686 .01232944765224485369462663996378,
9687 .015010447346388952376697286041943,
9688 .0175489679862431910996653529259,
9689 .019938037786440888202278192730714,
9690 .022194935961012286796332102959499,
9691 .024339147126000805470360647041454,
9692 .026374505414839207241503786552615,
9693 .02828691078877120065996800298796,
9694 .030052581128092695322521110347341,
9695 .031646751371439929404586051078883,
9696 .033050413419978503290785944862689,
9697 .034255099704226061787082821046821,
9698 .035262412660156681033782717998428,
9699 .036076989622888701185500318003895,
9700 .036698604498456094498018047441094,
9701 .037120549269832576114119958413599,
9702 .037334228751935040321235449094698,
9703 .037361073762679023410321241766599 };
9706 double d__1, d__2, d__3, d__4;
9710 double fv1[5], fv2[5], fv3[5], fv4[5];
9712 double absc, fval, res10, res21=0., res43=0., res87, fval1, fval2,
9713 hlgth, centr, reskh, uflow;
9714 double epmach, dhlgth, resabs=0., resasc=0., fcentr, savfun[21];
9855 d__1 = epmach * 50.;
9856 if (*epsabs <= 0. && *epsrel <
max(d__1,5e-29)) {
9859 hlgth = (*b - *a) * .5;
9860 dhlgth = fabs(hlgth);
9861 centr = (*b + *a) * .5;
9868 for (l = 1; l <= 3; ++l) {
9876 res21 = w21b[5] * fcentr;
9877 resabs = w21b[5] * fabs(fcentr);
9878 for (k = 1; k <= 5; ++k) {
9879 absc = hlgth * x1[k - 1];
9880 d__1 = centr + absc;
9882 d__1 = centr - absc;
9884 fval = fval1 + fval2;
9885 res10 += w10[k - 1] * fval;
9886 res21 += w21a[k - 1] * fval;
9887 resabs += w21a[k - 1] * (fabs(fval1) + fabs(fval2));
9888 savfun[k - 1] = fval;
9894 for (k = 1; k <= 5; ++k) {
9896 absc = hlgth * x2[k - 1];
9897 d__1 = centr + absc;
9899 d__1 = centr - absc;
9901 fval = fval1 + fval2;
9902 res21 += w21b[k - 1] * fval;
9903 resabs += w21b[k - 1] * (fabs(fval1) + fabs(fval2));
9904 savfun[ipx - 1] = fval;
9912 *result = res21 * hlgth;
9915 resasc = w21b[5] * (d__1 = fcentr - reskh, fabs(d__1));
9916 for (k = 1; k <= 5; ++k) {
9917 resasc = resasc + w21a[k - 1] * ((d__1 = fv1[k - 1] - reskh, fabs(
9918 d__1)) + (d__2 = fv2[k - 1] - reskh, fabs(d__2))) + w21b[k
9919 - 1] * ((d__3 = fv3[k - 1] - reskh, fabs(d__3)) + (d__4 =
9920 fv4[k - 1] - reskh, fabs(d__4)));
9923 *abserr = (d__1 = (res21 - res10) * hlgth, fabs(d__1));
9930 res43 = w43b[11] * fcentr;
9932 for (k = 1; k <= 10; ++k) {
9933 res43 += savfun[k - 1] * w43a[k - 1];
9936 for (k = 1; k <= 11; ++k) {
9938 absc = hlgth * x3[k - 1];
9939 d__1 = absc + centr;
9940 d__2 = centr - absc;
9941 fval = f(d__1) + f(d__2);
9942 res43 += fval * w43b[k - 1];
9943 savfun[ipx - 1] = fval;
9949 *result = res43 * hlgth;
9950 *abserr = (d__1 = (res43 - res21) * hlgth, fabs(d__1));
9956 res87 = w87b[22] * fcentr;
9958 for (k = 1; k <= 21; ++k) {
9959 res87 += savfun[k - 1] * w87a[k - 1];
9962 for (k = 1; k <= 22; ++k) {
9963 absc = hlgth * x4[k - 1];
9964 d__1 = absc + centr;
9965 d__2 = centr - absc;
9966 res87 += w87b[k - 1] * (f(d__1) + f(d__2));
9969 *result = res87 * hlgth;
9970 *abserr = (d__1 = (res87 - res43) * hlgth, fabs(d__1));
9972 if (resasc != 0. && *abserr != 0.) {
9974 d__3 = *abserr * 200. / resasc;
9975 d__1 = 1., d__2 = pow(d__3,
c_b270);
9976 *abserr = resasc *
min(d__1,d__2);
9978 if (resabs > uflow / (epmach * 50.)) {
9980 d__1 = epmach * 50. * resabs;
9981 *abserr =
max(d__1,*abserr);
9984 d__1 = *epsabs, d__2 = *epsrel * fabs(*result);
9985 if (*abserr <=
max(d__1,d__2)) {
10000 void dqpsrt_(
const long *limit,
long *last,
long *maxerr,
10001 double *ermax,
double *elist,
long *iord,
long *nrmax)
10007 long i__, j, k, ido, ibeg, jbnd, isucc, jupbn;
10008 double errmin, errmax;
10086 errmax = elist[*maxerr];
10092 for (i__ = 1; i__ <= i__1; ++i__) {
10093 isucc = iord[*nrmax - 1];
10095 if (errmax <= elist[isucc]) {
10098 iord[*nrmax] = isucc;
10109 if (*last > *limit / 2 + 2) {
10110 jupbn = *limit + 3 - *last;
10112 errmin = elist[*last];
10123 for (i__ = ibeg; i__ <= i__1; ++i__) {
10126 if (errmax >= elist[isucc]) {
10129 iord[i__ - 1] = isucc;
10133 iord[jbnd] = *maxerr;
10134 iord[jupbn] = *last;
10140 iord[i__ - 1] = *maxerr;
10143 for (j = i__; j <= i__1; ++j) {
10146 if (errmin < elist[isucc]) {
10149 iord[k + 1] = isucc;
10156 iord[k + 1] = *last;
10161 *maxerr = iord[*nrmax];
10162 *ermax = elist[*maxerr];
10166 double dqwgtc_(
const double *x,
const double *c__,
const double *,
10167 const double *,
const double *,
const long *)
10184 ret_val = 1. / (*x - *c__);
10188 double dqwgtf_(
const double *x,
const double *omega,
const double *,
10189 const double *,
const double *,
const long *integr)
10213 ret_val = cos(omx);
10216 ret_val = sin(omx);
10221 double dqwgts_(
const double *x,
const double *a,
const double *b,
10222 const double *alfa,
const double *beta,
const long *integr)
10245 ret_val = pow(xma, *alfa) * pow(bmx, *beta);
10253 ret_val *= log(xma);
10256 ret_val *= log(bmx);
10259 ret_val = ret_val * log(xma) * log(bmx);
10268 sys_float *elist,
long *iord,
long *last)
10279 sys_float area1, area2, area12, erro12, defab1, defab2;
10282 long iroff1, iroff2;
10283 sys_float error1, error2, defabs, epmach, errbnd, resabs, errmax;
10488 r__1 = epmach * 50.f;
10489 if (*epsabs <= 0.f && *epsrel <
max(r__1,5e-15f)) {
10508 qk15_(f, a, b, result, abserr, &defabs, &resabs);
10511 qk21_(f, a, b, result, abserr, &defabs, &resabs);
10514 qk31_(f, a, b, result, abserr, &defabs, &resabs);
10517 qk41_(f, a, b, result, abserr, &defabs, &resabs);
10520 qk51_(f, a, b, result, abserr, &defabs, &resabs);
10523 qk61_(f, a, b, result, abserr, &defabs, &resabs);
10526 rlist[1] = *result;
10527 elist[1] = *abserr;
10533 r__1 = *epsabs, r__2 = *epsrel * fabs(*result);
10534 errbnd =
max(r__1,r__2);
10535 if (*abserr <= epmach * 50.f * defabs && *abserr > errbnd) {
10541 if (*ier != 0 || ( *abserr <= errbnd && *abserr != resabs ) || *abserr == 0.f)
10562 for (*last = 2; *last <= i__1; ++(*last)) {
10566 a1 = alist__[maxerr];
10567 b1 = (alist__[maxerr] + blist[maxerr]) * .5f;
10569 b2 = blist[maxerr];
10571 qk15_(f, &a1, &b1, &area1, &error1, &resabs, &defab1);
10574 qk21_(f, &a1, &b1, &area1, &error1, &resabs, &defab1);
10577 qk31_(f, &a1, &b1, &area1, &error1, &resabs, &defab1);
10580 qk41_(f, &a1, &b1, &area1, &error1, &resabs, &defab1);
10583 qk51_(f, &a1, &b1, &area1, &error1, &resabs, &defab1);
10586 qk61_(f, &a1, &b1, &area1, &error1, &resabs, &defab1);
10589 qk15_(f, &a2, &b2, &area2, &error2, &resabs, &defab2);
10592 qk21_(f, &a2, &b2, &area2, &error2, &resabs, &defab2);
10595 qk31_(f, &a2, &b2, &area2, &error2, &resabs, &defab2);
10598 qk41_(f, &a2, &b2, &area2, &error2, &resabs, &defab2);
10601 qk51_(f, &a2, &b2, &area2, &error2, &resabs, &defab2);
10604 qk61_(f, &a2, &b2, &area2, &error2, &resabs, &defab2);
10611 area12 = area1 + area2;
10612 erro12 = error1 + error2;
10613 errsum = errsum + erro12 - errmax;
10614 area = area + area12 - rlist[maxerr];
10615 if (defab1 == error1 || defab2 == error2) {
10618 if ((r__1 = rlist[maxerr] - area12, fabs(r__1)) <= fabs(area12) *
10619 1e-5f && erro12 >= errmax * .99f) {
10622 if (*last > 10 && erro12 > errmax) {
10626 rlist[maxerr] = area1;
10627 rlist[*last] = area2;
10629 r__1 = *epsabs, r__2 = *epsrel * fabs(area);
10630 errbnd =
max(r__1,r__2);
10631 if (errsum <= errbnd) {
10638 if (iroff1 >= 6 || iroff2 >= 20) {
10645 if (*last == *limit) {
10653 r__1 = fabs(a1), r__2 = fabs(b2);
10654 if (
max(r__1,r__2) <= (epmach * 100.f + 1.f) * (fabs(a2) + uflow *
10662 if (error2 > error1) {
10665 alist__[*last] =
a2;
10666 blist[maxerr] =
b1;
10668 elist[maxerr] = error1;
10669 elist[*last] = error2;
10672 alist__[maxerr] =
a2;
10673 alist__[*last] =
a1;
10675 rlist[maxerr] = area2;
10676 rlist[*last] = area1;
10677 elist[maxerr] = error2;
10678 elist[*last] = error1;
10686 qpsrt_(limit, last, &maxerr, &errmax, &elist[1], &iord[1], &nrmax);
10688 if (*ier != 0 || errsum <= errbnd) {
10700 for (k = 1; k <= i__1; ++k) {
10701 *result += rlist[k];
10707 *neval = (keyf * 10 + 1) * ((*neval << 1) + 1);
10710 *neval = *neval * 30 + 15;
10719 long *ier,
long *limit,
const long *lenw,
long *last,
long *iwork,
10722 long l1, l2, l3, lvl;
10888 if (*limit < 1 || *lenw < *limit << 2) {
10898 qage_(f, a, b, epsabs, epsrel, key, limit, result, abserr, neval,
10899 ier, &work[1], &work[l1], &work[l2], &work[l3], &iwork[1], last);
10909 xerror_(
"abnormal return from qag ", &
c__26, ier, &lvl, 26);
10919 sys_float *elist,
long *iord,
long *last)
10934 sys_float area1, area2, area12, small=0.f, erro12;
10940 long iroff1, iroff2, iroff3;
10941 sys_float res3la[3], error1, error2, rlist2[52];
10943 sys_float defabs, epmach, erlarg=0.f, abseps, correc=0.f, errbnd, resabs;
11180 r__1 = epmach * 50.f;
11181 if (*epsabs <= 0.f && *epsrel <
max(r__1,5e-15f)) {
11207 rlist[1] = *result;
11208 elist[1] = *abserr;
11210 dres = fabs(*result);
11212 r__1 = *epsabs, r__2 = *epsrel * dres;
11213 errbnd =
max(r__1,r__2);
11214 if (*abserr <= epmach * 100.f * defabs && *abserr > errbnd) {
11220 if (*ier != 0 || ( *abserr <= errbnd && *abserr != resabs ) || *abserr == 0.f)
11230 rlist2[0] = *result;
11247 if (dres >= (1.f - epmach * 50.f) * defabs) {
11255 for (*last = 2; *last <= i__1; ++(*last)) {
11260 a1 = alist__[maxerr];
11261 b1 = (alist__[maxerr] + blist[maxerr]) * .5f;
11263 b2 = blist[maxerr];
11265 qk15i_(f, &boun, inf, &a1, &b1, &area1, &error1, &resabs, &
11267 qk15i_(f, &boun, inf, &a2, &b2, &area2, &error2, &resabs, &
11273 area12 = area1 + area2;
11274 erro12 = error1 + error2;
11275 errsum = errsum + erro12 - errmax;
11276 area = area + area12 - rlist[maxerr];
11277 if (defab1 == error1 || defab2 == error2) {
11280 if ((r__1 = rlist[maxerr] - area12, fabs(r__1)) > fabs(area12) *
11281 1e-5f || erro12 < errmax * .99f) {
11291 if (*last > 10 && erro12 > errmax) {
11295 rlist[maxerr] = area1;
11296 rlist[*last] = area2;
11298 r__1 = *epsabs, r__2 = *epsrel * fabs(area);
11299 errbnd =
max(r__1,r__2);
11304 if (iroff1 + iroff2 >= 10 || iroff3 >= 20) {
11314 if (*last == *limit) {
11322 r__1 = fabs(a1), r__2 = fabs(b2);
11323 if (
max(r__1,r__2) <= (epmach * 100.f + 1.f) * (fabs(a2) + uflow *
11330 if (error2 > error1) {
11333 alist__[*last] =
a2;
11334 blist[maxerr] =
b1;
11336 elist[maxerr] = error1;
11337 elist[*last] = error2;
11340 alist__[maxerr] =
a2;
11341 alist__[*last] =
a1;
11343 rlist[maxerr] = area2;
11344 rlist[*last] = area1;
11345 elist[maxerr] = error2;
11346 elist[*last] = error1;
11354 qpsrt_(limit, last, &maxerr, &errmax, &elist[1], &iord[1], &nrmax);
11355 if (errsum <= errbnd) {
11368 if ((r__1 = b1 - a1, fabs(r__1)) > small) {
11378 if ((r__1 = blist[maxerr] - alist__[maxerr], fabs(r__1)) > small) {
11384 if (ierro == 3 || erlarg <= ertest) {
11395 if (*last > *limit / 2 + 2) {
11396 jupbnd = *limit + 3 - *last;
11399 for (k =
id; k <= i__2; ++k) {
11400 maxerr = iord[nrmax];
11401 errmax = elist[maxerr];
11402 if ((r__1 = blist[maxerr] - alist__[maxerr], fabs(r__1)) > small)
11414 rlist2[numrl2 - 1] = area;
11415 qelg_(&numrl2, rlist2, &reseps, &abseps, res3la, &nres);
11417 if (ktmin > 5 && *abserr < errsum * .001f) {
11420 if (abseps >= *abserr) {
11428 r__1 = *epsabs, r__2 = *epsrel * fabs(reseps);
11429 ertest =
max(r__1,r__2);
11430 if (*abserr <= ertest) {
11444 errmax = elist[maxerr];
11463 if (*abserr == oflow) {
11466 if (*ier + ierro == 0) {
11475 if (*result != 0.f && area != 0.f) {
11478 if (*abserr > errsum) {
11486 if (*abserr / fabs(*result) > errsum / fabs(area)) {
11494 r__1 = fabs(*result), r__2 = fabs(area);
11495 if (ksgn == -1 &&
max(r__1,r__2) <= defabs * .01f) {
11498 if (.01f > *result / area || *result / area > 100.f || errsum > fabs(area)
11509 for (k = 1; k <= i__1; ++k) {
11510 *result += rlist[k];
11515 *neval = *last * 30 - 15;
11529 long *ier,
const long *limit,
const long *lenw,
long *last,
11532 long l1, l2, l3, lvl;
11706 if (*limit < 1 || *lenw < *limit << 2) {
11716 qagie_(f, bound, inf, epsabs, epsrel, limit, result, abserr, neval,
11717 ier, &work[1], &work[l1], &work[l2], &work[l3], &iwork[1], last);
11727 xerror_(
"abnormal return from qagi", &
c__26, ier, &lvl, 26);
11737 sys_float *pts,
long *iord,
long *level,
long *ndin,
long *last)
11752 long nres,
nint, jlow, npts;
11753 sys_float area1, area2, area12, erro12;
11759 long iroff1, iroff2, iroff3;
11764 sys_float defabs, epmach, erlarg, abseps, correc=0.f, errbnd, resabs;
11769 long maxerr, levcur;
12047 r__1 = epmach * 50.f;
12048 if (*npts2 < 2 || *limit <= npts ||
12049 ( *epsabs <= 0.f && *epsrel <
max(r__1, 5e-15f)) ) {
12063 pts[1] =
min(*a,*b);
12068 for (i__ = 1; i__ <= i__1; ++i__) {
12069 pts[i__ + 1] = points[i__];
12073 pts[npts + 2] =
max(*a,*b);
12081 for (i__ = 1; i__ <= i__1; ++i__) {
12084 for (j = ip1; j <= i__2; ++j) {
12085 if (pts[i__] <= pts[j]) {
12095 if (pts[1] !=
min(*a,*b) || pts[nintp1] !=
max(*a,*b)) {
12108 for (i__ = 1; i__ <= i__2; ++i__) {
12110 qk21_(f, &a1, &b1, &area1, &error1, &defabs, &resa);
12114 if (error1 == resa && error1 != 0.f) {
12119 elist[i__] = error1;
12122 rlist[i__] = area1;
12129 for (i__ = 1; i__ <= i__2; ++i__) {
12130 if (ndin[i__] == 1) {
12131 elist[i__] = *abserr;
12133 errsum += elist[i__];
12140 *neval = nint * 21;
12141 dres = fabs(*result);
12143 r__1 = *epsabs, r__2 = *epsrel * dres;
12144 errbnd =
max(r__1,r__2);
12145 if (*abserr <= epmach * 100.f * resabs && *abserr > errbnd) {
12152 for (i__ = 1; i__ <= i__2; ++i__) {
12156 for (j = jlow; j <= i__1; ++j) {
12158 if (elist[ind1] > elist[ind2]) {
12166 if (ind1 == iord[i__]) {
12169 iord[k] = iord[i__];
12174 if (*limit < *npts2) {
12178 if (*ier != 0 || *abserr <= errbnd) {
12185 rlist2[0] = *result;
12187 errmax = elist[maxerr];
12206 if (dres >= (1.f - epmach * 50.f) * resabs) {
12214 for (*last = *npts2; *last <= i__2; ++(*last)) {
12219 levcur = level[maxerr] + 1;
12220 a1 = alist__[maxerr];
12221 b1 = (alist__[maxerr] + blist[maxerr]) * .5f;
12223 b2 = blist[maxerr];
12225 qk21_(f, &a1, &b1, &area1, &error1, &resa, &defab1);
12226 qk21_(f, &a2, &b2, &area2, &error2, &resa, &defab2);
12232 area12 = area1 + area2;
12233 erro12 = error1 + error2;
12234 errsum = errsum + erro12 - errmax;
12235 area = area + area12 - rlist[maxerr];
12236 if (defab1 == error1 || defab2 == error2) {
12239 if ((r__1 = rlist[maxerr] - area12, fabs(r__1)) > fabs(area12) *
12240 1e-5f || erro12 < errmax * .99f) {
12250 if (*last > 10 && erro12 > errmax) {
12254 level[maxerr] = levcur;
12255 level[*last] = levcur;
12256 rlist[maxerr] = area1;
12257 rlist[*last] = area2;
12259 r__1 = *epsabs, r__2 = *epsrel * fabs(area);
12260 errbnd =
max(r__1,r__2);
12265 if (iroff1 + iroff2 >= 10 || iroff3 >= 20) {
12275 if (*last == *limit) {
12283 r__1 = fabs(a1), r__2 = fabs(b2);
12284 if (
max(r__1,r__2) <= (epmach * 100.f + 1.f) * (fabs(a2) + uflow *
12291 if (error2 > error1) {
12294 alist__[*last] =
a2;
12295 blist[maxerr] =
b1;
12297 elist[maxerr] = error1;
12298 elist[*last] = error2;
12301 alist__[maxerr] =
a2;
12302 alist__[*last] =
a1;
12304 rlist[maxerr] = area2;
12305 rlist[*last] = area1;
12306 elist[maxerr] = error2;
12307 elist[*last] = error1;
12315 qpsrt_(limit, last, &maxerr, &errmax, &elist[1], &iord[1], &nrmax);
12317 if (errsum <= errbnd) {
12328 if (levcur + 1 <= levmax) {
12338 if (level[maxerr] + 1 <= levmax) {
12344 if (ierro == 3 || erlarg <= ertest) {
12355 if (*last > *limit / 2 + 2) {
12356 jupbnd = *limit + 3 - *last;
12359 for (k =
id; k <= i__1; ++k) {
12360 maxerr = iord[nrmax];
12361 errmax = elist[maxerr];
12363 if (level[maxerr] + 1 <= levmax) {
12374 rlist2[numrl2 - 1] = area;
12378 qelg_(&numrl2, rlist2, &reseps, &abseps, res3la, &nres);
12380 if (ktmin > 5 && *abserr < errsum * .001f) {
12383 if (abseps >= *abserr) {
12391 r__1 = *epsabs, r__2 = *epsrel * fabs(reseps);
12392 ertest =
max(r__1,r__2);
12394 if (*abserr < ertest) {
12409 errmax = elist[maxerr];
12423 if (*abserr == oflow) {
12426 if (*ier + ierro == 0) {
12435 if (*result != 0.f && area != 0.f) {
12438 if (*abserr > errsum) {
12446 if (*abserr / fabs(*result) > errsum / fabs(area)) {
12454 r__1 = fabs(*result), r__2 = fabs(area);
12455 if (ksgn == -1 &&
max(r__1,r__2) <= resabs * .01f) {
12458 if (.01f > *result / area || *result / area > 100.f || errsum > fabs(area)
12469 for (k = 1; k <= i__2; ++k) {
12470 *result += rlist[k];
12486 long *neval,
long *ier,
const long *leniw,
const long *lenw,
12487 long *last,
long *iwork,
sys_float *work)
12489 long l1, l2, l3, l4, lvl;
12694 if (*leniw < *npts2 * 3 - 2 || *lenw < (*leniw << 1) - *npts2 || *npts2 <
12701 limit = (*leniw - *npts2) / 2;
12707 qagpe_(f, a, b, npts2, &points[1], epsabs, epsrel, &limit, result,
12708 abserr, neval, ier, &work[1], &work[l1], &work[l2], &work[l3], &
12709 work[l4], &iwork[1], &iwork[l1], &iwork[l2], last);
12719 xerror_(
"abnormal return from qagp", &
c__26, ier, &lvl, 26);
12730 long *iord,
long *last)
12743 sys_float area1, area2, area12, small=0.f, erro12;
12749 long iroff1, iroff2, iroff3;
12750 sys_float res3la[3], error1, error2, rlist2[52];
12752 sys_float defabs, epmach, erlarg=0.f, abseps, correc=0.f, errbnd, resabs;
12985 r__1 = epmach * 50.f;
12986 if (*epsabs <= 0.f && *epsrel <
max(r__1,5e-15f)) {
12999 qk21_(f, a, b, result, abserr, &defabs, &resabs);
13003 dres = fabs(*result);
13005 r__1 = *epsabs, r__2 = *epsrel * dres;
13006 errbnd =
max(r__1,r__2);
13008 rlist[1] = *result;
13009 elist[1] = *abserr;
13011 if (*abserr <= epmach * 100.f * defabs && *abserr > errbnd) {
13017 if (*ier != 0 || ( *abserr <= errbnd && *abserr != resabs ) ||
13026 rlist2[0] = *result;
13042 if (dres >= (1.f - epmach * 50.f) * defabs) {
13050 for (*last = 2; *last <= i__1; ++(*last)) {
13055 a1 = alist__[maxerr];
13056 b1 = (alist__[maxerr] + blist[maxerr]) * .5f;
13058 b2 = blist[maxerr];
13060 qk21_(f, &a1, &b1, &area1, &error1, &resabs, &defab1);
13061 qk21_(f, &a2, &b2, &area2, &error2, &resabs, &defab2);
13066 area12 = area1 + area2;
13067 erro12 = error1 + error2;
13068 errsum = errsum + erro12 - errmax;
13069 area = area + area12 - rlist[maxerr];
13070 if (defab1 == error1 || defab2 == error2) {
13073 if ((r__1 = rlist[maxerr] - area12, fabs(r__1)) > fabs(area12) *
13074 1e-5f || erro12 < errmax * .99f) {
13084 if (*last > 10 && erro12 > errmax) {
13088 rlist[maxerr] = area1;
13089 rlist[*last] = area2;
13091 r__1 = *epsabs, r__2 = *epsrel * fabs(area);
13092 errbnd =
max(r__1,r__2);
13097 if (iroff1 + iroff2 >= 10 || iroff3 >= 20) {
13107 if (*last == *limit) {
13115 r__1 = fabs(a1), r__2 = fabs(b2);
13116 if (
max(r__1,r__2) <= (epmach * 100.f + 1.f) * (fabs(a2) + uflow *
13123 if (error2 > error1) {
13126 alist__[*last] =
a2;
13127 blist[maxerr] =
b1;
13129 elist[maxerr] = error1;
13130 elist[*last] = error2;
13133 alist__[maxerr] =
a2;
13134 alist__[*last] =
a1;
13136 rlist[maxerr] = area2;
13137 rlist[*last] = area1;
13138 elist[maxerr] = error2;
13139 elist[*last] = error1;
13147 qpsrt_(limit, last, &maxerr, &errmax, &elist[1], &iord[1], &nrmax);
13149 if (errsum <= errbnd) {
13163 if ((r__1 = b1 - a1, fabs(r__1)) > small) {
13173 if ((r__1 = blist[maxerr] - alist__[maxerr], fabs(r__1)) > small) {
13179 if (ierro == 3 || erlarg <= ertest) {
13190 if (*last > *limit / 2 + 2) {
13191 jupbnd = *limit + 3 - *last;
13194 for (k =
id; k <= i__2; ++k) {
13195 maxerr = iord[nrmax];
13196 errmax = elist[maxerr];
13198 if ((r__1 = blist[maxerr] - alist__[maxerr], fabs(r__1)) > small)
13210 rlist2[numrl2 - 1] = area;
13211 qelg_(&numrl2, rlist2, &reseps, &abseps, res3la, &nres);
13213 if (ktmin > 5 && *abserr < errsum * .001f) {
13216 if (abseps >= *abserr) {
13224 r__1 = *epsabs, r__2 = *epsrel * fabs(reseps);
13225 ertest =
max(r__1,r__2);
13227 if (*abserr <= ertest) {
13241 errmax = elist[maxerr];
13248 small = (r__1 = *b - *a, fabs(r__1)) * .375f;
13260 if (*abserr == oflow) {
13263 if (*ier + ierro == 0) {
13272 if (*result != 0.f && area != 0.f) {
13275 if (*abserr > errsum) {
13283 if (*abserr / fabs(*result) > errsum / fabs(area)) {
13291 r__1 = fabs(*result), r__2 = fabs(area);
13292 if (ksgn == -1 &&
max(r__1,r__2) <= defabs * .01f) {
13295 if (.01f > *result / area || *result / area > 100.f || errsum > fabs(area)
13306 for (k = 1; k <= i__1; ++k) {
13307 *result += rlist[k];
13316 *neval = *last * 42 - 21;
13324 const long *limit,
const long *lenw,
long *last,
long *iwork,
13327 long l1, l2, l3, lvl;
13500 if (*limit < 1 || *lenw < *limit << 2) {
13510 qagse_(f, a, b, epsabs, epsrel, limit, result, abserr, neval, ier, &
13511 work[1], &work[l1], &work[l2], &work[l3], &iwork[1], last);
13521 xerror_(
"abnormal return from qags", &
c__26, ier, &lvl, 26);
13528 const sys_float *epsrel,
const long *limit,
13542 sys_float area1, area2, area12, erro12;
13545 long iroff1, iroff2;
13546 sys_float error1, error2, epmach, errbnd, errmax;
13744 r__1 = epmach * 50.f;
13745 if (*c__ == *a || *c__ == *b ||
13746 ( *epsabs <= 0.f && *epsrel <
max(r__1, 5e-15f))) {
13763 qc25c_(f, &aa, &bb, c__, result, abserr, &krule, neval);
13765 rlist[1] = *result;
13766 elist[1] = *abserr;
13774 r__1 = *epsabs, r__2 = *epsrel * fabs(*result);
13775 errbnd =
max(r__1,r__2);
13780 r__1 = fabs(*result) * .01f;
13781 if (*abserr <
min(r__1,errbnd) || *ier == 1) {
13790 rlist[1] = *result;
13803 for (*last = 2; *last <= i__1; ++(*last)) {
13808 a1 = alist__[maxerr];
13809 b1 = (alist__[maxerr] + blist[maxerr]) * .5f;
13810 b2 = blist[maxerr];
13811 if (*c__ <= b1 && *c__ > a1) {
13812 b1 = (*c__ +
b2) * .5f;
13814 if (*c__ > b1 && *c__ < b2) {
13815 b1 = (a1 + *c__) * .5f;
13819 qc25c_(f, &a1, &b1, c__, &area1, &error1, &krule, &nev);
13821 qc25c_(f, &a2, &b2, c__, &area2, &error2, &krule, &nev);
13827 area12 = area1 + area2;
13828 erro12 = error1 + error2;
13829 errsum = errsum + erro12 - errmax;
13830 area = area + area12 - rlist[maxerr];
13831 if ((r__1 = rlist[maxerr] - area12, fabs(r__1)) < fabs(area12) *
13832 1e-5f && erro12 >= errmax * .99f && krule == 0) {
13835 if (*last > 10 && erro12 > errmax && krule == 0) {
13838 rlist[maxerr] = area1;
13839 rlist[*last] = area2;
13841 r__1 = *epsabs, r__2 = *epsrel * fabs(area);
13842 errbnd =
max(r__1,r__2);
13843 if (errsum <= errbnd) {
13850 if (iroff1 >= 6 && iroff2 > 20) {
13857 if (*last == *limit) {
13865 r__1 = fabs(a1), r__2 = fabs(b2);
13866 if (
max(r__1,r__2) <= (epmach * 100.f + 1.f) * (fabs(a2) + uflow *
13874 if (error2 > error1) {
13877 alist__[*last] =
a2;
13878 blist[maxerr] =
b1;
13880 elist[maxerr] = error1;
13881 elist[*last] = error2;
13884 alist__[maxerr] =
a2;
13885 alist__[*last] =
a1;
13887 rlist[maxerr] = area2;
13888 rlist[*last] = area1;
13889 elist[maxerr] = error2;
13890 elist[*last] = error1;
13898 qpsrt_(limit, last, &maxerr, &errmax, &elist[1], &iord[1], &nrmax);
13900 if (*ier != 0 || errsum <= errbnd) {
13912 for (k = 1; k <= i__1; ++k) {
13913 *result += rlist[k];
13919 *result = -(*result);
13928 sys_float *abserr,
long *neval,
long *ier,
13929 long *limit,
const long *lenw,
long *last,
long *iwork,
13932 long l1, l2, l3, lvl;
14094 if (*limit < 1 || *lenw < *limit << 2) {
14103 qawce_(f, a, b, c__, epsabs, epsrel, limit, result, abserr, neval,
14104 ier, &work[1], &work[l1], &work[l2], &work[l3], &iwork[1], last);
14114 xerror_(
"abnormal return from qawc", &
c__26, ier, &lvl, 26);
14120 const long *integr,
14121 const sys_float *epsabs,
const long *limlst,
const long *limit,
14134 int chebmo_dim1, chebmo_offset, i__1;
14389 chebmo_dim1 = *maxp1;
14390 chebmo_offset = 1 + chebmo_dim1;
14391 chebmo -= chebmo_offset;
14404 if ( ( *integr != 1 && *integr != 2 ) || *epsabs <= 0.f || *limlst < 3) {
14410 if (*omega != 0.f) {
14417 if (*integr == 1) {
14419 abserr, neval, ier, &alist__[1], &blist[1], &rlist[1], &elist[
14420 1], &iord[1], &last);
14422 rslst[1] = *result;
14423 erlst[1] = *abserr;
14434 cycle = dl * pi / fabs(*omega);
14445 if (*epsabs > uflow / p1) {
14446 eps = *epsabs * p1;
14458 for (*lst = 1; *lst <= i__1; ++(*lst)) {
14463 qawoe_(f, &c1, &c2, omega, integr, &epsa, &
c_b390, limit, lst,
14464 maxp1, &rslst[*lst], &erlst[*lst], &nev, &ierlst[*lst], &last,
14465 &alist__[1], &blist[1], &rlist[1], &elist[1], &iord[1], &
14466 nnlog[1], &momcom, &chebmo[chebmo_offset]);
14469 errsum += erlst[*lst];
14470 drl = (r__1 = rslst[*lst], fabs(r__1)) * 50.f;
14474 if (errsum + drl <= *epsabs && *lst >= 6) {
14478 r__1 = correc, r__2 = erlst[*lst];
14479 correc =
max(r__1,r__2);
14480 if (ierlst[*lst] != 0) {
14482 r__1 = ep, r__2 = correc * p1;
14483 eps =
max(r__1,r__2);
14485 if (ierlst[*lst] != 0) {
14488 if (*ier == 7 && errsum + drl <= correc * 10.f && *lst > 5) {
14495 psum[0] = rslst[1];
14498 psum[numrl2 - 1] = psum[ll - 1] + rslst[*lst];
14505 if (*lst == *limlst) {
14511 qelg_(&numrl2, psum, &reseps, &abseps, res3la, &nres);
14517 if (ktmin >= 15 && *abserr <= (errsum + drl) * .001f) {
14520 if (abseps > *abserr && *lst != 3) {
14531 if (*abserr + correc * 10.f <= *epsabs ||
14532 ( *abserr <= *epsabs && correc * 10.f >= *epsabs) ) {
14536 if (*ier != 0 && *ier != 7) {
14550 *abserr += correc * 10.f;
14554 if (*result != 0.f && psum[numrl2 - 1] != 0.f) {
14557 if (*abserr > errsum) {
14560 if (psum[numrl2 - 1] == 0.f) {
14564 if (*abserr / fabs(*result) > (errsum + drl) / (r__1 = psum[numrl2 - 1],
14568 if (*ier >= 1 && *ier != 7) {
14573 *result = psum[numrl2 - 1];
14574 *abserr = errsum + drl;
14580 const long *integr,
14582 long *neval,
long *ier,
const long *limlst,
long *lst,
14583 const long *leniw,
const long *maxp1,
14584 const long *lenw,
long *iwork,
sys_float *work)
14586 long l1, l2, l3, l4, l5, l6, ll2, lvl;
14794 if (*limlst < 3 || *leniw < *limlst + 2 || *maxp1 < 1 || *lenw < (*leniw
14795 << 1) + *maxp1 * 25) {
14801 limit = (*leniw - *limlst) / 2;
14809 qawfe_(f, a, omega, integr, epsabs, limlst, &limit, maxp1, result,
14810 abserr, neval, ier, &work[1], &work[l1], &iwork[1], lst, &work[l2]
14811 , &work[l3], &work[l4], &work[l5], &iwork[l1], &iwork[ll2], &work[
14822 xerror_(
"abnormal return from qawf", &
c__26, ier, &lvl, 26);
14828 const sys_float *omega,
const long *integr,
14830 const long *limit,
const long *icall,
const long *maxp1,
14833 elist,
long *iord,
long *nnlog,
long *momcom,
sys_float *chebmo)
14836 int chebmo_dim1, chebmo_offset, i__1, i__2;
14846 sys_float area1, area2, area12, small, erro12, defab1, defab2, width;
14849 long ktmin, nrmax, nrmom;
14852 long iroff1, iroff2, iroff3;
14853 sys_float res3la[3], error1, error2, rlist2[52];
14855 sys_float defabs, domega, epmach, erlarg=0.f, abseps, correc=0.f, errbnd,
15128 chebmo_dim1 = *maxp1;
15129 chebmo_offset = 1 + chebmo_dim1;
15130 chebmo -= chebmo_offset;
15150 r__1 = epmach * 50.f;
15151 if ( ( *integr != 1 && *integr != 2) ||
15152 (*epsabs <= 0.f && *epsrel <
max(r__1,5e-15f)) ||
15153 *icall < 1 || *maxp1 < 1) {
15163 domega = fabs(*omega);
15170 qc25f_(f, a, b, &domega, integr, &nrmom, maxp1, &
c__0, result,
15171 abserr, neval, &defabs, &resabs, momcom, &chebmo[chebmo_offset]);
15175 dres = fabs(*result);
15177 r__1 = *epsabs, r__2 = *epsrel * dres;
15178 errbnd =
max(r__1,r__2);
15179 rlist[1] = *result;
15180 elist[1] = *abserr;
15182 if (*abserr <= epmach * 100.f * defabs && *abserr > errbnd) {
15188 if (*ier != 0 || *abserr <= errbnd) {
15210 small = (r__1 = *b - *a, fabs(r__1)) * .75f;
15214 if ((r__1 = *b - *a, fabs(r__1)) * .5f * domega > 2.f) {
15219 rlist2[0] = *result;
15221 if ((r__1 = *b - *a, fabs(r__1)) * .25f * domega <= 2.f) {
15225 if (dres >= (1.f - epmach * 50.f) * defabs) {
15233 for (*last = 2; *last <= i__1; ++(*last)) {
15238 nrmom = nnlog[maxerr] + 1;
15239 a1 = alist__[maxerr];
15240 b1 = (alist__[maxerr] + blist[maxerr]) * .5f;
15242 b2 = blist[maxerr];
15244 qc25f_(f, &a1, &b1, &domega, integr, &nrmom, maxp1, &
c__0, &
15245 area1, &error1, &nev, &resabs, &defab1, momcom, &chebmo[
15248 qc25f_(f, &a2, &b2, &domega, integr, &nrmom, maxp1, &
c__1, &
15249 area2, &error2, &nev, &resabs, &defab2, momcom, &chebmo[
15256 area12 = area1 + area2;
15257 erro12 = error1 + error2;
15258 errsum = errsum + erro12 - errmax;
15259 area = area + area12 - rlist[maxerr];
15260 if (defab1 == error1 || defab2 == error2) {
15263 if ((r__1 = rlist[maxerr] - area12, fabs(r__1)) > fabs(area12) *
15264 1e-5f || erro12 < errmax * .99f) {
15274 if (*last > 10 && erro12 > errmax) {
15278 rlist[maxerr] = area1;
15279 rlist[*last] = area2;
15280 nnlog[maxerr] = nrmom;
15281 nnlog[*last] = nrmom;
15283 r__1 = *epsabs, r__2 = *epsrel * fabs(area);
15284 errbnd =
max(r__1,r__2);
15289 if (iroff1 + iroff2 >= 10 || iroff3 >= 20) {
15299 if (*last == *limit) {
15307 r__1 = fabs(a1), r__2 = fabs(b2);
15308 if (
max(r__1,r__2) <= (epmach * 100.f + 1.f) * (fabs(a2) + uflow *
15315 if (error2 > error1) {
15318 alist__[*last] =
a2;
15319 blist[maxerr] =
b1;
15321 elist[maxerr] = error1;
15322 elist[*last] = error2;
15325 alist__[maxerr] =
a2;
15326 alist__[*last] =
a1;
15328 rlist[maxerr] = area2;
15329 rlist[*last] = area1;
15330 elist[maxerr] = error2;
15331 elist[*last] = error1;
15339 qpsrt_(limit, last, &maxerr, &errmax, &elist[1], &iord[1], &nrmax);
15341 if (errsum <= errbnd) {
15347 if (*last == 2 && extall) {
15357 if ((r__1 = b1 - a1, fabs(r__1)) > small) {
15368 width = (r__1 = blist[maxerr] - alist__[maxerr], fabs(r__1));
15369 if (width > small) {
15382 if (width * .25f * domega > 2.f) {
15391 if (ierro == 3 || erlarg <= ertest) {
15401 if (*last > *limit / 2 + 2) {
15402 jupbnd = *limit + 3 - *last;
15406 for (k =
id; k <= i__2; ++k) {
15407 maxerr = iord[nrmax];
15408 errmax = elist[maxerr];
15409 if ((r__1 = blist[maxerr] - alist__[maxerr], fabs(r__1)) > small)
15421 rlist2[numrl2 - 1] = area;
15425 qelg_(&numrl2, rlist2, &reseps, &abseps, res3la, &nres);
15427 if (ktmin > 5 && *abserr < errsum * .001f) {
15430 if (abseps >= *abserr) {
15438 r__1 = *epsabs, r__2 = *epsrel * fabs(reseps);
15439 ertest =
max(r__1,r__2);
15441 if (*abserr <= ertest) {
15456 errmax = elist[maxerr];
15465 rlist2[numrl2 - 1] = area;
15477 if (*abserr == oflow || nres == 0) {
15480 if (*ier + ierro == 0) {
15489 if (*result != 0.f && area != 0.f) {
15492 if (*abserr > errsum) {
15500 if (*abserr / fabs(*result) > errsum / fabs(area)) {
15508 r__1 = fabs(*result), r__2 = fabs(area);
15509 if (ksgn == -1 &&
max(r__1,r__2) <= defabs * .01f) {
15512 if (.01f > *result / area || *result / area > 100.f || errsum >= fabs(
15523 for (k = 1; k <= i__1; ++k) {
15524 *result += rlist[k];
15533 if (*integr == 2 && *omega < 0.f) {
15534 *result = -(*result);
15541 const sys_float *omega,
const long *integr,
15544 long *neval,
long *ier,
const long *leniw,
const long *maxp1,
15545 const long *lenw,
long *last,
long *iwork,
sys_float *work)
15547 long l1, l2, l3, l4, lvl;
15548 long limit, momcom;
15752 if (*leniw < 2 || *maxp1 < 1 || *lenw < (*leniw << 1) + *maxp1 * 25) {
15758 limit = *leniw / 2;
15763 qawoe_(f, a, b, omega, integr, epsabs, epsrel, &limit, &
c__1, maxp1,
15764 result, abserr, neval, ier, last, &work[1], &work[l1], &work[l2],
15765 &work[l3], &iwork[1], &iwork[l1], &momcom, &work[l4]);
15775 xerror_(
"abnormal return from qawo", &
c__26, ier, &lvl, 26);
15782 const long *integr,
const sys_float *epsabs,
15796 sys_float area1, area2, area12, erro12;
15799 long iroff1, iroff2;
15800 sys_float resas1, resas2, error1, error2, epmach, errbnd, centre,
16018 r__1 = epmach * 50.f;
16019 if (*b <= *a || ( *epsabs == 0.f && *epsrel <
max(r__1,5e-15f)) || *alfa <=
16020 -1.f || *beta <= -1.f || *integr < 1 || *integr > 4 || *limit < 2)
16028 qmomo_(alfa, beta, ri, rj, rg, rh, integr);
16033 centre = (*b + *a) * .5f;
16034 qc25s_(f, a, b, a, ¢re, alfa, beta, ri, rj, rg, rh, &area1, &
16035 error1, &resas1, integr, &nev);
16037 qc25s_(f, a, b, ¢re, b, alfa, beta, ri, rj, rg, rh, &area2, &
16038 error2, &resas2, integr, &nev);
16041 *result = area1 + area2;
16042 *abserr = error1 + error2;
16047 r__1 = *epsabs, r__2 = *epsrel * fabs(*result);
16048 errbnd =
max(r__1,r__2);
16053 if (error2 > error1) {
16057 alist__[2] = centre;
16066 alist__[1] = centre;
16080 if (*abserr <= errbnd || *ier == 1) {
16095 for (*last = 3; *last <= i__1; ++(*last)) {
16099 a1 = alist__[maxerr];
16100 b1 = (alist__[maxerr] + blist[maxerr]) * .5f;
16102 b2 = blist[maxerr];
16104 qc25s_(f, a, b, &a1, &b1, alfa, beta, ri, rj, rg, rh, &area1, &
16105 error1, &resas1, integr, &nev);
16107 qc25s_(f, a, b, &a2, &b2, alfa, beta, ri, rj, rg, rh, &area2, &
16108 error2, &resas2, integr, &nev);
16114 area12 = area1 + area2;
16115 erro12 = error1 + error2;
16116 errsum = errsum + erro12 - errmax;
16117 area = area + area12 - rlist[maxerr];
16118 if (*a == a1 || *b == b2) {
16121 if (resas1 == error1 || resas2 == error2) {
16127 if ((r__1 = rlist[maxerr] - area12, fabs(r__1)) < fabs(area12) *
16128 1e-5f && erro12 >= errmax * .99f) {
16131 if (*last > 10 && erro12 > errmax) {
16135 rlist[maxerr] = area1;
16136 rlist[*last] = area2;
16141 r__1 = *epsabs, r__2 = *epsrel * fabs(area);
16142 errbnd =
max(r__1,r__2);
16143 if (errsum <= errbnd) {
16150 if (*last == *limit) {
16157 if (iroff1 >= 6 || iroff2 >= 20) {
16165 r__1 = fabs(a1), r__2 = fabs(b2);
16166 if (
max(r__1,r__2) <= (epmach * 100.f + 1.f) * (fabs(a2) + uflow *
16174 if (error2 > error1) {
16177 alist__[*last] =
a2;
16178 blist[maxerr] =
b1;
16180 elist[maxerr] = error1;
16181 elist[*last] = error2;
16184 alist__[maxerr] =
a2;
16185 alist__[*last] =
a1;
16187 rlist[maxerr] = area2;
16188 rlist[*last] = area1;
16189 elist[maxerr] = error2;
16190 elist[*last] = error1;
16198 qpsrt_(limit, last, &maxerr, &errmax, &elist[1], &iord[1], &nrmax);
16200 if (*ier != 0 || errsum <= errbnd) {
16212 for (k = 1; k <= i__1; ++k) {
16213 *result += rlist[k];
16223 const long *integr,
const sys_float *epsabs,
16225 long *neval,
long *ier,
const long *limit,
const long *lenw,
16226 long *last,
long *iwork,
sys_float *work)
16228 long l1, l2, l3, lvl;
16411 if (*limit < 2 || *lenw < *limit << 2) {
16421 qawse_(f, a, b, alfa, beta, integr, epsabs, epsrel, limit, result,
16422 abserr, neval, ier, &work[1], &work[l1], &work[l2], &work[l3], &
16433 xerror_(
"abnormal return from qaws", &
c__26, ier, &lvl, 26);
16440 sys_float *abserr,
long *krul,
long *neval)
16444 sys_float x[11] = { .9914448613738104f,.9659258262890683f,
16445 .9238795325112868f,.8660254037844386f,.7933533402912352f,
16446 .7071067811865475f,.6087614290087206f,.5f,.3826834323650898f,
16447 .2588190451025208f,.1305261922200516f };
16456 sys_float ak22, fval[25], res12, res24;
16458 sys_float amom0, amom1, amom2, cheb12[13], cheb24[25];
16545 cc = (2.f * *c__ - *b - *a) / (*b - *a);
16546 if (fabs(cc) < 1.1f) {
16553 qk15w_(f,
qwgtc_, c__, &p2, &p3, &p4, &kp, a, b, result,
16554 abserr, &resabs, &resasc);
16556 if (resasc == *abserr) {
16564 hlgth = (*b - *a) * .5f;
16565 centr = (*b + *a) * .5f;
16567 r__1 = hlgth + centr;
16568 fval[0] = f(r__1) * .5f;
16569 fval[12] = f(centr);
16570 r__1 = centr - hlgth;
16571 fval[24] = f(r__1) * .5f;
16572 for (i__ = 2; i__ <= 12; ++i__) {
16573 u = hlgth * x[i__ - 2];
16576 fval[i__ - 1] = f(r__1);
16578 fval[isym - 1] = f(r__1);
16584 qcheb_(x, fval, cheb12, cheb24);
16590 amom0 = log((r__1 = (1.f - cc) / (cc + 1.f), fabs(r__1)));
16591 amom1 = cc * amom0 + 2.f;
16592 res12 = cheb12[0] * amom0 + cheb12[1] * amom1;
16593 res24 = cheb24[0] * amom0 + cheb24[1] * amom1;
16594 for (k = 3; k <= 13; ++k) {
16595 amom2 = cc * 2.f * amom1 - amom0;
16596 ak22 = (
sys_float) ((k - 2) * (k - 2));
16597 if (k / 2 << 1 == k) {
16598 amom2 -= 4.f / (ak22 - 1.f);
16600 res12 += cheb12[k - 1] * amom2;
16601 res24 += cheb24[k - 1] * amom2;
16606 for (k = 14; k <= 25; ++k) {
16607 amom2 = cc * 2.f * amom1 - amom0;
16608 ak22 = (
sys_float) ((k - 2) * (k - 2));
16609 if (k / 2 << 1 == k) {
16610 amom2 -= 4.f / (ak22 - 1.f);
16612 res24 += cheb24[k - 1] * amom2;
16618 *abserr = (r__1 = res24 - res12, fabs(r__1));
16624 const sys_float *omega,
const long *integr,
16625 const long *nrmom,
const long *maxp1,
const long *ksave,
16632 sys_float x[11] = { .9914448613738104f,.9659258262890683f,
16633 .9238795325112868f,.8660254037844386f,.7933533402912352f,
16634 .7071067811865475f,.6087614290087206f,.5f,.3826834323650898f,
16635 .2588190451025208f,.1305261922200516f };
16638 int chebmo_dim1, chebmo_offset, i__1;
16643 long i__, j, k, m=0;
16644 sys_float v[28], d1[25], d2[25], p2, p3, p4, ac, an, as, an2, ass, par2,
16645 conc, asap, par22, fval[25], estc, cons;
16650 sys_float resc12, resc24, hlgth, centr, ress12, ress24, oflow;
16754 chebmo_dim1 = *maxp1;
16755 chebmo_offset = 1 + chebmo_dim1;
16756 chebmo -= chebmo_offset;
16792 centr = (*b + *a) * .5f;
16793 hlgth = (*b - *a) * .5f;
16794 parint = *omega * hlgth;
16800 if (fabs(parint) > 2.f) {
16803 qk15w_(f,
qwgtf_, omega, &p2, &p3, &p4, integr, a, b, result,
16804 abserr, resabs, resasc);
16812 conc = hlgth * cos(centr * *omega);
16813 cons = hlgth * sin(centr * *omega);
16820 if (*nrmom < *momcom || *ksave == 1) {
16827 par2 = parint * parint;
16828 par22 = par2 + 2.f;
16829 sinpar = sin(parint);
16830 cospar = cos(parint);
16834 v[0] = sinpar * 2.f / parint;
16835 v[1] = (cospar * 8.f + (par2 + par2 - 8.f) * sinpar / parint) / par2;
16836 v[2] = ((par2 - 12.f) * 32.f * cospar + ((par2 - 80.f) * par2 + 192.f) *
16837 2.f * sinpar / parint) / (par2 * par2);
16839 as = parint * 24.f *
sinpar;
16840 if (fabs(parint) > 24.f) {
16853 for (k = 1; k <= i__1; ++k) {
16855 d__[k - 1] = (an2 - 4.f) * -2.f * (par22 - an2 - an2);
16856 d2[k - 1] = (an - 1.f) * (an - 2.f) * par2;
16857 d1[k] = (an + 3.f) * (an + 4.f) * par2;
16858 v[k + 2] = as - (an2 - 4.f) * ac;
16863 d__[noequ - 1] = (an2 - 4.f) * -2.f * (par22 - an2 - an2);
16864 v[noequ + 2] = as - (an2 - 4.f) * ac;
16865 v[3] -= par2 * 56.f * v[2];
16867 asap = (((((par2 * 210.f - 1.f) * cospar - (par2 * 105.f - 63.f) * ass) /
16868 an2 - (1.f - par2 * 15.f) * cospar + ass * 15.f) / an2 - cospar +
16869 ass * 3.f) / an2 - cospar) / an2;
16870 v[noequ + 2] -= asap * 2.f * par2 * (an - 1.f) * (an - 2.f);
16875 sgtsl_(&noequ, d1, d__, d2, &v[3], &iers);
16883 for (i__ = 4; i__ <= 13; ++i__) {
16885 v[i__ - 1] = ((an2 - 4.f) * ((par22 - an2 - an2) * 2.f * v[i__ - 2] -
16886 ac) + as - par2 * (an + 1.f) * (an + 2.f) * v[i__ - 3]) / (
16887 par2 * (an - 1.f) * (an - 2.f));
16892 for (j = 1; j <= 13; ++j) {
16893 chebmo[m + ((j << 1) - 1) * chebmo_dim1] = v[j - 1];
16899 v[0] = (sinpar - parint * cospar) * 2.f / par2;
16900 v[1] = (18.f - 48.f / par2) * sinpar / par2 + (48.f / par2 - 2.f) *
16902 ac = parint * -24.f * cospar;
16903 as = sinpar * -8.f;
16904 if (fabs(parint) > 24.f) {
16915 for (k = 1; k <= i__1; ++k) {
16917 d__[k - 1] = (an2 - 4.f) * -2.f * (par22 - an2 - an2);
16918 d2[k - 1] = (an - 1.f) * (an - 2.f) * par2;
16919 d1[k] = (an + 3.f) * (an + 4.f) * par2;
16920 v[k + 1] = ac + (an2 - 4.f) * as;
16925 d__[noequ - 1] = (an2 - 4.f) * -2.f * (par22 - an2 - an2);
16926 v[noequ + 1] = ac + (an2 - 4.f) * as;
16927 v[2] -= par2 * 42.f * v[1];
16928 ass = parint * cospar;
16929 asap = (((((par2 * 105.f - 63.f) * ass + (par2 * 210.f - 1.f) *
sinpar) /
16930 an2 + (par2 * 15.f - 1.f) * sinpar - ass * 15.f) / an2 - ass *
16931 3.f - sinpar) / an2 -
sinpar) / an2;
16932 v[noequ + 1] -= asap * 2.f * par2 * (an - 1.f) * (an - 2.f);
16937 sgtsl_(&noequ, d1, d__, d2, &v[2], &iers);
16945 for (i__ = 3; i__ <= 12; ++i__) {
16947 v[i__ - 1] = ((an2 - 4.f) * ((par22 - an2 - an2) * 2.f * v[i__ - 2] +
16948 as) + ac - par2 * (an + 1.f) * (an + 2.f) * v[i__ - 3]) / (
16949 par2 * (an - 1.f) * (an - 2.f));
16954 for (j = 1; j <= 12; ++j) {
16955 chebmo[m + (j << 1) * chebmo_dim1] = v[j - 1];
16959 if (*nrmom < *momcom) {
16962 if (*momcom < *maxp1 - 1 && *nrmom >= *momcom) {
16969 r__1 = centr + hlgth;
16970 fval[0] = f(r__1) * .5f;
16971 fval[12] = f(centr);
16972 r__1 = centr - hlgth;
16973 fval[24] = f(r__1) * .5f;
16974 for (i__ = 2; i__ <= 12; ++i__) {
16976 r__1 = hlgth * x[i__ - 2] + centr;
16977 fval[i__ - 1] = f(r__1);
16978 r__1 = centr - hlgth * x[i__ - 2];
16979 fval[isym - 1] = f(r__1);
16982 qcheb_(x, fval, cheb12, cheb24);
16986 resc12 = cheb12[12] * chebmo[m + chebmo_dim1 * 13];
16989 for (j = 1; j <= 6; ++j) {
16990 resc12 += cheb12[k - 1] * chebmo[m + k * chebmo_dim1];
16991 ress12 += cheb12[k] * chebmo[m + (k + 1) * chebmo_dim1];
16995 resc24 = cheb24[24] * chebmo[m + chebmo_dim1 * 25];
16997 *resabs = fabs(cheb24[24]);
16999 for (j = 1; j <= 12; ++j) {
17000 resc24 += cheb24[k - 1] * chebmo[m + k * chebmo_dim1];
17001 ress24 += cheb24[k] * chebmo[m + (k + 1) * chebmo_dim1];
17002 *resabs = (r__1 = cheb24[k - 1], fabs(r__1)) + (r__2 = cheb24[k],
17007 estc = (r__1 = resc24 - resc12, fabs(r__1));
17008 ests = (r__1 = ress24 - ress12, fabs(r__1));
17009 *resabs *= fabs(hlgth);
17010 if (*integr == 2) {
17013 *result = conc * resc24 - cons * ress24;
17014 *abserr = (r__1 = conc * estc, fabs(r__1)) + (r__2 = cons * ests, fabs(
17018 *result = conc * ress24 + cons * resc24;
17019 *abserr = (r__1 = conc * ests, fabs(r__1)) + (r__2 = cons * estc, fabs(
17030 const long *integr,
long *nev)
17034 sys_float x[11] = { .9914448613738104f,.9659258262890683f,
17035 .9238795325112868f,.8660254037844386f,.7933533402912352f,
17036 .7071067811865475f,.6087614290087206f,.5f,.3826834323650898f,
17037 .2588190451025208f,.1305261922200516f };
17045 sys_float u, dc, fix, fval[25], res12, res24;
17164 if (*bl == *a && (*alfa != 0.f || *integr == 2 || *integr == 4)) {
17167 if (*br == *b && (*beta != 0.f || *integr == 3 || *integr == 4)) {
17175 qk15w_(f,
qwgts_, a, b, alfa, beta, integr, bl, br, result,
17176 abserr, &resabs, resasc);
17189 hlgth = (*br - *bl) * .5f;
17190 centr = (*br + *bl) * .5f;
17192 r__1 = hlgth + centr;
17193 d__1 = (double) (fix - hlgth);
17194 d__2 = (double) (*beta);
17195 fval[0] = f(r__1) * .5f * pow(d__1, d__2);
17196 d__1 = (double) fix;
17197 d__2 = (double) (*beta);
17198 fval[12] = f(centr) * pow(d__1, d__2);
17199 r__1 = centr - hlgth;
17200 d__1 = (double) (fix + hlgth);
17201 d__2 = (double) (*beta);
17202 fval[24] = f(r__1) * .5f * pow(d__1, d__2);
17203 for (i__ = 2; i__ <= 12; ++i__) {
17204 u = hlgth * x[i__ - 2];
17207 d__1 = (double) (fix - u);
17208 d__2 = (double) (*beta);
17209 fval[i__ - 1] = f(r__1) * pow(d__1, d__2);
17211 d__1 = (double) (fix + u);
17212 d__2 = (double) (*beta);
17213 fval[isym - 1] = f(r__1) * pow(d__1, d__2);
17216 d__1 = (double) hlgth;
17217 d__2 = (double) (*alfa + 1.f);
17218 factor = pow(d__1, d__2);
17226 qcheb_(x, fval, cheb12, cheb24);
17230 for (i__ = 1; i__ <= 13; ++i__) {
17231 res12 += cheb12[i__ - 1] * ri[i__];
17232 res24 += cheb24[i__ - 1] * ri[i__];
17235 for (i__ = 14; i__ <= 25; ++i__) {
17236 res24 += cheb24[i__ - 1] * ri[i__];
17239 if (*integr == 1) {
17245 dc = log(*br - *bl);
17246 *result = res24 * dc;
17247 *abserr = (r__1 = (res24 - res12) * dc, fabs(r__1));
17250 for (i__ = 1; i__ <= 13; ++i__) {
17251 res12 += cheb12[i__ - 1] * rg[i__];
17252 res24 = res12 + cheb24[i__ - 1] * rg[i__];
17255 for (i__ = 14; i__ <= 25; ++i__) {
17256 res24 += cheb24[i__ - 1] * rg[i__];
17266 fval[0] *= log(fix - hlgth);
17267 fval[12] *= log(fix);
17268 fval[24] *= log(fix + hlgth);
17269 for (i__ = 2; i__ <= 12; ++i__) {
17270 u = hlgth * x[i__ - 2];
17272 fval[i__ - 1] *= log(fix - u);
17273 fval[isym - 1] *= log(fix + u);
17276 qcheb_(x, fval, cheb12, cheb24);
17280 for (i__ = 1; i__ <= 13; ++i__) {
17281 res12 += cheb12[i__ - 1] * ri[i__];
17282 res24 += cheb24[i__ - 1] * ri[i__];
17285 for (i__ = 14; i__ <= 25; ++i__) {
17286 res24 += cheb24[i__ - 1] * ri[i__];
17289 if (*integr == 3) {
17295 dc = log(*br - *bl);
17296 *result = res24 * dc;
17297 *abserr = (r__1 = (res24 - res12) * dc, fabs(r__1));
17300 for (i__ = 1; i__ <= 13; ++i__) {
17301 res12 += cheb12[i__ - 1] * rg[i__];
17302 res24 += cheb24[i__ - 1] * rg[i__];
17305 for (i__ = 14; i__ <= 25; ++i__) {
17306 res24 += cheb24[i__ - 1] * rg[i__];
17310 *result = (*result + res24) * factor;
17311 *abserr = (*abserr + (r__1 = res24 - res12, fabs(r__1))) * factor;
17323 hlgth = (*br - *bl) * .5f;
17324 centr = (*br + *bl) * .5f;
17326 r__1 = hlgth + centr;
17327 d__1 = (double) (fix + hlgth);
17328 d__2 = (double) (*alfa);
17329 fval[0] = f(r__1) * .5f * pow(d__1, d__2);
17330 d__1 = (double) fix;
17331 d__2 = (double) (*alfa);
17332 fval[12] = f(centr) * pow(d__1, d__2);
17333 r__1 = centr - hlgth;
17334 d__1 = (double) (fix - hlgth);
17335 d__2 = (double) (*alfa);
17336 fval[24] = f(r__1) * .5f * pow(d__1, d__2);
17337 for (i__ = 2; i__ <= 12; ++i__) {
17338 u = hlgth * x[i__ - 2];
17341 d__1 = (double) (fix + u);
17342 d__2 = (double) (*alfa);
17343 fval[i__ - 1] = f(r__1) * pow(d__1, d__2);
17345 d__1 = (double) (fix - u);
17346 d__2 = (double) (*alfa);
17347 fval[isym - 1] = f(r__1) * pow(d__1, d__2);
17350 d__1 = (double) hlgth;
17351 d__2 = (double) (*beta + 1.f);
17352 factor = pow(d__1, d__2);
17357 if (*integr == 2 || *integr == 4) {
17363 qcheb_(x, fval, cheb12, cheb24);
17364 for (i__ = 1; i__ <= 13; ++i__) {
17365 res12 += cheb12[i__ - 1] * rj[i__];
17366 res24 += cheb24[i__ - 1] * rj[i__];
17369 for (i__ = 14; i__ <= 25; ++i__) {
17370 res24 += cheb24[i__ - 1] * rj[i__];
17373 if (*integr == 1) {
17379 dc = log(*br - *bl);
17380 *result = res24 * dc;
17381 *abserr = (r__1 = (res24 - res12) * dc, fabs(r__1));
17384 for (i__ = 1; i__ <= 13; ++i__) {
17385 res12 += cheb12[i__ - 1] * rh[i__];
17386 res24 += cheb24[i__ - 1] * rh[i__];
17389 for (i__ = 14; i__ <= 25; ++i__) {
17390 res24 += cheb24[i__ - 1] * rh[i__];
17400 fval[0] *= log(hlgth + fix);
17401 fval[12] *= log(fix);
17402 fval[24] *= log(fix - hlgth);
17403 for (i__ = 2; i__ <= 12; ++i__) {
17404 u = hlgth * x[i__ - 2];
17406 fval[i__ - 1] *= log(u + fix);
17407 fval[isym - 1] *= log(fix - u);
17410 qcheb_(x, fval, cheb12, cheb24);
17414 for (i__ = 1; i__ <= 13; ++i__) {
17415 res12 += cheb12[i__ - 1] * rj[i__];
17416 res24 += cheb24[i__ - 1] * rj[i__];
17419 for (i__ = 14; i__ <= 25; ++i__) {
17420 res24 += cheb24[i__ - 1] * rj[i__];
17423 if (*integr == 2) {
17426 dc = log(*br - *bl);
17427 *result = res24 * dc;
17428 *abserr = (r__1 = (res24 - res12) * dc, fabs(r__1));
17434 for (i__ = 1; i__ <= 13; ++i__) {
17435 res12 += cheb12[i__ - 1] * rh[i__];
17436 res24 += cheb24[i__ - 1] * rh[i__];
17439 for (i__ = 14; i__ <= 25; ++i__) {
17440 res24 += cheb24[i__ - 1] * rh[i__];
17444 *result = (*result + res24) * factor;
17445 *abserr = (*abserr + (r__1 = res24 - res12, fabs(r__1))) * factor;
17453 sys_float v[12], alam, alam1, alam2, part1, part2, part3;
17509 for (i__ = 1; i__ <= 12; ++i__) {
17511 v[i__ - 1] = fval[i__] - fval[j];
17512 fval[i__] += fval[j];
17515 alam1 = v[0] - v[8];
17516 alam2 = x[6] * (v[2] - v[6] - v[10]);
17517 cheb12[4] = alam1 + alam2;
17518 cheb12[10] = alam1 - alam2;
17519 alam1 = v[1] - v[7] - v[9];
17520 alam2 = v[3] - v[5] - v[11];
17521 alam = x[3] * alam1 + x[9] * alam2;
17522 cheb24[4] = cheb12[4] + alam;
17523 cheb24[22] = cheb12[4] - alam;
17524 alam = x[9] * alam1 - x[3] * alam2;
17525 cheb24[10] = cheb12[10] + alam;
17526 cheb24[16] = cheb12[10] - alam;
17527 part1 = x[4] * v[4];
17528 part2 = x[8] * v[8];
17529 part3 = x[6] * v[6];
17530 alam1 = v[0] + part1 + part2;
17531 alam2 = x[2] * v[2] + part3 + x[10] * v[10];
17532 cheb12[2] = alam1 + alam2;
17533 cheb12[12] = alam1 - alam2;
17534 alam = x[1] * v[1] + x[3] * v[3] + x[5] * v[5] + x[7] * v[7] + x[9] * v[9]
17536 cheb24[2] = cheb12[2] + alam;
17537 cheb24[24] = cheb12[2] - alam;
17538 alam = x[11] * v[1] - x[9] * v[3] + x[7] * v[5] - x[5] * v[7] + x[3] * v[
17540 cheb24[12] = cheb12[12] + alam;
17541 cheb24[14] = cheb12[12] - alam;
17542 alam1 = v[0] - part1 + part2;
17543 alam2 = x[10] * v[2] - part3 + x[2] * v[10];
17544 cheb12[6] = alam1 + alam2;
17545 cheb12[8] = alam1 - alam2;
17546 alam = x[5] * v[1] - x[9] * v[3] - x[1] * v[5] - x[11] * v[7] + x[3] * v[
17548 cheb24[6] = cheb12[6] + alam;
17549 cheb24[20] = cheb12[6] - alam;
17550 alam = x[7] * v[1] - x[3] * v[3] - x[11] * v[5] + x[1] * v[7] - x[9] * v[
17552 cheb24[8] = cheb12[8] + alam;
17553 cheb24[18] = cheb12[8] - alam;
17554 for (i__ = 1; i__ <= 6; ++i__) {
17556 v[i__ - 1] = fval[i__] - fval[j];
17557 fval[i__] += fval[j];
17560 alam1 = v[0] + x[8] * v[4];
17561 alam2 = x[4] * v[2];
17562 cheb12[3] = alam1 + alam2;
17563 cheb12[11] = alam1 - alam2;
17564 cheb12[7] = v[0] - v[4];
17565 alam = x[2] * v[1] + x[6] * v[3] + x[10] * v[5];
17566 cheb24[3] = cheb12[3] + alam;
17567 cheb24[23] = cheb12[3] - alam;
17568 alam = x[6] * (v[1] - v[3] - v[5]);
17569 cheb24[7] = cheb12[7] + alam;
17570 cheb24[19] = cheb12[7] - alam;
17571 alam = x[10] * v[1] - x[6] * v[3] + x[2] * v[5];
17572 cheb24[11] = cheb12[11] + alam;
17573 cheb24[15] = cheb12[11] - alam;
17574 for (i__ = 1; i__ <= 3; ++i__) {
17576 v[i__ - 1] = fval[i__] - fval[j];
17577 fval[i__] += fval[j];
17580 cheb12[5] = v[0] + x[8] * v[2];
17581 cheb12[9] = fval[1] - x[8] * fval[3];
17582 alam = x[4] * v[1];
17583 cheb24[5] = cheb12[5] + alam;
17584 cheb24[21] = cheb12[5] - alam;
17585 alam = x[8] * fval[2] - fval[4];
17586 cheb24[9] = cheb12[9] + alam;
17587 cheb24[17] = cheb12[9] - alam;
17588 cheb12[1] = fval[1] + fval[3];
17589 alam = fval[2] + fval[4];
17590 cheb24[1] = cheb12[1] + alam;
17591 cheb24[25] = cheb12[1] - alam;
17592 cheb12[13] = v[0] - v[2];
17593 cheb24[13] = cheb12[13];
17594 alam = .16666666666666666f;
17595 for (i__ = 2; i__ <= 12; ++i__) {
17596 cheb12[i__] *= alam;
17601 cheb12[13] *= alam;
17602 for (i__ = 2; i__ <= 24; ++i__) {
17603 cheb24[i__] *= alam;
17606 cheb24[1] = alam * .5f * cheb24[1];
17607 cheb24[25] = alam * .5f * cheb24[25];
17621 long k1, k2, k3, ib, ie;
17626 sys_float err1, err2, err3, tol1, tol2, tol3;
17628 sys_float e1abs, oflow, error, delta1, delta2, delta3;
17630 long newelm, limexp;
17716 *result = epstab[*n];
17721 epstab[*n + 2] = epstab[*n];
17722 newelm = (*n - 1) / 2;
17723 epstab[*n] = oflow;
17727 for (i__ = 1; i__ <= i__1; ++i__) {
17730 res = epstab[k1 + 2];
17736 err2 = fabs(delta2);
17739 tol2 =
max(r__1,e1abs) * epmach;
17741 err3 = fabs(delta3);
17743 r__1 = e1abs, r__2 = fabs(e0);
17744 tol3 =
max(r__1,r__2) * epmach;
17745 if (err2 > tol2 || err3 > tol3) {
17755 *abserr = err2 + err3;
17762 err1 = fabs(delta1);
17764 r__1 = e1abs, r__2 = fabs(e3);
17765 tol1 =
max(r__1,r__2) * epmach;
17770 if (err1 <= tol1 || err2 <= tol2 || err3 <= tol3) {
17773 ss = 1.f / delta1 + 1.f / delta2 - 1.f / delta3;
17774 epsinf = (r__1 = ss *
e1, fabs(r__1));
17780 if (epsinf > 1e-4f) {
17784 *n = i__ + i__ - 1;
17792 res = e1 + 1.f / ss;
17795 error = err2 + (r__1 = res -
e2, fabs(r__1)) + err3;
17796 if (error > *abserr) {
17808 if (*n == limexp) {
17809 *n = (limexp / 2 << 1) - 1;
17812 if (num / 2 << 1 == num) {
17817 for (i__ = 1; i__ <= i__1; ++i__) {
17819 epstab[ib] = epstab[ib2];
17826 indx = num - *n + 1;
17828 for (i__ = 1; i__ <= i__1; ++i__) {
17829 epstab[i__] = epstab[indx];
17837 res3la[*nres] = *result;
17844 *abserr = (r__1 = *result - res3la[3], fabs(r__1)) + (r__2 = *result -
17845 res3la[2], fabs(r__2)) + (r__3 = *result - res3la[1], fabs(r__3));
17846 res3la[1] = res3la[2];
17847 res3la[2] = res3la[3];
17848 res3la[3] = *result;
17851 r__1 = *abserr, r__2 = epmach * 5.f * fabs(*result);
17852 *abserr =
max(r__1,r__2);
17862 sys_float xgk[8] = { .9914553711208126f,.9491079123427585f,
17863 .8648644233597691f,.7415311855993944f,.5860872354676911f,
17864 .4058451513773972f,.2077849550078985f,0.f };
17865 sys_float wgk[8] = { .02293532201052922f,.06309209262997855f,
17866 .1047900103222502f,.1406532597155259f,.1690047266392679f,
17867 .1903505780647854f,.2044329400752989f,.2094821410847278f };
17868 sys_float wg[4] = { .1294849661688697f,.2797053914892767f,
17869 .3818300505051189f,.4179591836734694f };
17879 sys_float absc, resg, resk, fsum, fval1, fval2;
17975 centr = (*a + *b) * .5f;
17976 hlgth = (*b - *a) * .5f;
17977 dhlgth = fabs(hlgth);
17984 resk = fc * wgk[7];
17985 *resabs = fabs(resk);
17986 for (j = 1; j <= 3; ++j) {
17988 absc = hlgth * xgk[jtw - 1];
17989 r__1 = centr - absc;
17991 r__1 = centr + absc;
17993 fv1[jtw - 1] = fval1;
17994 fv2[jtw - 1] = fval2;
17995 fsum = fval1 + fval2;
17996 resg += wg[j - 1] * fsum;
17997 resk += wgk[jtw - 1] * fsum;
17998 *resabs += wgk[jtw - 1] * (fabs(fval1) + fabs(fval2));
18001 for (j = 1; j <= 4; ++j) {
18002 jtwm1 = (j << 1) - 1;
18003 absc = hlgth * xgk[jtwm1 - 1];
18004 r__1 = centr - absc;
18006 r__1 = centr + absc;
18008 fv1[jtwm1 - 1] = fval1;
18009 fv2[jtwm1 - 1] = fval2;
18010 fsum = fval1 + fval2;
18011 resk += wgk[jtwm1 - 1] * fsum;
18012 *resabs += wgk[jtwm1 - 1] * (fabs(fval1) + fabs(fval2));
18015 reskh = resk * .5f;
18016 *resasc = wgk[7] * (r__1 = fc - reskh, fabs(r__1));
18017 for (j = 1; j <= 7; ++j) {
18018 *resasc += wgk[j - 1] * ((r__1 = fv1[j - 1] - reskh, fabs(r__1)) + (
18019 r__2 = fv2[j - 1] - reskh, fabs(r__2)));
18022 *result = resk * hlgth;
18025 *abserr = (r__1 = (resk - resg) * hlgth, fabs(r__1));
18026 if (*resasc != 0.f && *abserr != 0.f) {
18028 d__1 = (double) (*abserr * 200.f / *resasc);
18029 r__1 = 1.f, r__2 = pow(d__1,
c_b270);
18030 *abserr = *resasc *
min(r__1,r__2);
18032 if (*resabs > uflow / (epmach * 50.f)) {
18034 r__1 = epmach * 50.f * *resabs;
18035 *abserr =
max(r__1,*abserr);
18046 sys_float xgk[8] = { .9914553711208126f,.9491079123427585f,
18047 .8648644233597691f,.7415311855993944f,.5860872354676911f,
18048 .4058451513773972f,.2077849550078985f,0.f };
18049 sys_float wgk[8] = { .02293532201052922f,.06309209262997855f,
18050 .1047900103222502f,.1406532597155259f,.1690047266392679f,
18051 .1903505780647854f,.2044329400752989f,.2094821410847278f };
18052 sys_float wg[8] = { 0.f,.1294849661688697f,0.f,.2797053914892767f,0.f,
18053 .3818300505051189f,0.f,.4179591836734694f };
18061 sys_float fc, fv1[7], fv2[7], absc, dinf, resg, resk, fsum, absc1,
18062 absc2, fval1, fval2, hlgth, centr, reskh, uflow;
18181 centr = (*a + *b) * .5f;
18182 hlgth = (*b - *a) * .5f;
18183 tabsc1 = *boun + dinf * (1.f - centr) / centr;
18189 fc = fval1 / centr / centr;
18195 resk = wgk[7] * fc;
18196 *resabs = fabs(resk);
18197 for (j = 1; j <= 7; ++j) {
18198 absc = hlgth * xgk[j - 1];
18199 absc1 = centr - absc;
18200 absc2 = centr + absc;
18201 tabsc1 = *boun + dinf * (1.f - absc1) / absc1;
18202 tabsc2 = *boun + dinf * (1.f - absc2) / absc2;
18213 fval1 = fval1 / absc1 / absc1;
18214 fval2 = fval2 / absc2 / absc2;
18215 fv1[j - 1] = fval1;
18216 fv2[j - 1] = fval2;
18217 fsum = fval1 + fval2;
18218 resg += wg[j - 1] * fsum;
18219 resk += wgk[j - 1] * fsum;
18220 *resabs += wgk[j - 1] * (fabs(fval1) + fabs(fval2));
18223 reskh = resk * .5f;
18224 *resasc = wgk[7] * (r__1 = fc - reskh, fabs(r__1));
18225 for (j = 1; j <= 7; ++j) {
18226 *resasc += wgk[j - 1] * ((r__1 = fv1[j - 1] - reskh, fabs(r__1)) + (
18227 r__2 = fv2[j - 1] - reskh, fabs(r__2)));
18230 *result = resk * hlgth;
18233 *abserr = (r__1 = (resk - resg) * hlgth, fabs(r__1));
18234 if (*resasc != 0.f && *abserr != 0.f) {
18236 d__1 = (double) (*abserr * 200.f / *resasc);
18237 r__1 = 1.f, r__2 = pow(d__1,
c_b270);
18238 *abserr = *resasc *
min(r__1,r__2);
18240 if (*resabs > uflow / (epmach * 50.f)) {
18242 r__1 = epmach * 50.f * *resabs;
18243 *abserr =
max(r__1,*abserr);
18255 sys_float xgk[8] = { .9914553711208126f,.9491079123427585f,
18256 .8648644233597691f,.7415311855993944f,.5860872354676911f,
18257 .4058451513773972f,.2077849550078985f,0.f };
18258 sys_float wgk[8] = { .02293532201052922f,.06309209262997855f,
18259 .1047900103222502f,.1406532597155259f,.1690047266392679f,
18260 .1903505780647854f,.2044329400752989f,.2094821410847278f };
18261 sys_float wg[4] = { .1294849661688697f,.2797053914892767f,
18262 .3818300505051889f,.4179591836734694f };
18272 sys_float absc, resg, resk, fsum, absc1, absc2, fval1, fval2;
18381 centr = (*a + *b) * .5f;
18382 hlgth = (*b - *a) * .5f;
18383 dhlgth = fabs(hlgth);
18388 fc = f(centr) * w(¢r, p1, p2, p3, p4, kp);
18390 resk = wgk[7] * fc;
18391 *resabs = fabs(resk);
18392 for (j = 1; j <= 3; ++j) {
18394 absc = hlgth * xgk[jtw - 1];
18395 absc1 = centr - absc;
18396 absc2 = centr + absc;
18397 fval1 = f(absc1) * w(&absc1, p1, p2, p3, p4, kp);
18398 fval2 = f(absc2) * w(&absc2, p1, p2, p3, p4, kp);
18399 fv1[jtw - 1] = fval1;
18400 fv2[jtw - 1] = fval2;
18401 fsum = fval1 + fval2;
18402 resg += wg[j - 1] * fsum;
18403 resk += wgk[jtw - 1] * fsum;
18404 *resabs += wgk[jtw - 1] * (fabs(fval1) + fabs(fval2));
18407 for (j = 1; j <= 4; ++j) {
18408 jtwm1 = (j << 1) - 1;
18409 absc = hlgth * xgk[jtwm1 - 1];
18410 absc1 = centr - absc;
18411 absc2 = centr + absc;
18412 fval1 = f(absc1) * w(&absc1, p1, p2, p3, p4, kp);
18413 fval2 = f(absc2) * w(&absc2, p1, p2, p3, p4, kp);
18414 fv1[jtwm1 - 1] = fval1;
18415 fv2[jtwm1 - 1] = fval2;
18416 fsum = fval1 + fval2;
18417 resk += wgk[jtwm1 - 1] * fsum;
18418 *resabs += wgk[jtwm1 - 1] * (fabs(fval1) + fabs(fval2));
18421 reskh = resk * .5f;
18422 *resasc = wgk[7] * (r__1 = fc - reskh, fabs(r__1));
18423 for (j = 1; j <= 7; ++j) {
18424 *resasc += wgk[j - 1] * ((r__1 = fv1[j - 1] - reskh, fabs(r__1)) + (
18425 r__2 = fv2[j - 1] - reskh, fabs(r__2)));
18428 *result = resk * hlgth;
18431 *abserr = (r__1 = (resk - resg) * hlgth, fabs(r__1));
18432 if (*resasc != 0.f && *abserr != 0.f) {
18434 d__1 = (double) (*abserr * 200.f / *resasc);
18435 r__1 = 1.f, r__2 = pow(d__1,
c_b270);
18436 *abserr = *resasc *
min(r__1,r__2);
18438 if (*resabs > uflow / (epmach * 50.f)) {
18440 r__1 = epmach * 50.f * *resabs;
18441 *abserr =
max(r__1,*abserr);
18451 sys_float xgk[11] = { .9956571630258081f,.9739065285171717f,
18452 .9301574913557082f,.8650633666889845f,.7808177265864169f,
18453 .6794095682990244f,.5627571346686047f,.4333953941292472f,
18454 .2943928627014602f,.1488743389816312f,0.f };
18455 sys_float wgk[11] = { .01169463886737187f,.03255816230796473f,
18456 .054755896574352f,.07503967481091995f,.09312545458369761f,
18457 .1093871588022976f,.1234919762620659f,.1347092173114733f,
18458 .1427759385770601f,.1477391049013385f,.1494455540029169f };
18459 sys_float wg[5] = { .06667134430868814f,.1494513491505806f,
18460 .219086362515982f,.2692667193099964f,.2955242247147529f };
18470 sys_float absc, resg, resk, fsum, fval1, fval2;
18569 centr = (*a + *b) * .5f;
18570 hlgth = (*b - *a) * .5f;
18571 dhlgth = fabs(hlgth);
18578 resk = wgk[10] * fc;
18579 *resabs = fabs(resk);
18580 for (j = 1; j <= 5; ++j) {
18582 absc = hlgth * xgk[jtw - 1];
18583 r__1 = centr - absc;
18585 r__1 = centr + absc;
18587 fv1[jtw - 1] = fval1;
18588 fv2[jtw - 1] = fval2;
18589 fsum = fval1 + fval2;
18590 resg += wg[j - 1] * fsum;
18591 resk += wgk[jtw - 1] * fsum;
18592 *resabs += wgk[jtw - 1] * (fabs(fval1) + fabs(fval2));
18595 for (j = 1; j <= 5; ++j) {
18596 jtwm1 = (j << 1) - 1;
18597 absc = hlgth * xgk[jtwm1 - 1];
18598 r__1 = centr - absc;
18600 r__1 = centr + absc;
18602 fv1[jtwm1 - 1] = fval1;
18603 fv2[jtwm1 - 1] = fval2;
18604 fsum = fval1 + fval2;
18605 resk += wgk[jtwm1 - 1] * fsum;
18606 *resabs += wgk[jtwm1 - 1] * (fabs(fval1) + fabs(fval2));
18609 reskh = resk * .5f;
18610 *resasc = wgk[10] * (r__1 = fc - reskh, fabs(r__1));
18611 for (j = 1; j <= 10; ++j) {
18612 *resasc += wgk[j - 1] * ((r__1 = fv1[j - 1] - reskh, fabs(r__1)) + (
18613 r__2 = fv2[j - 1] - reskh, fabs(r__2)));
18616 *result = resk * hlgth;
18619 *abserr = (r__1 = (resk - resg) * hlgth, fabs(r__1));
18620 if (*resasc != 0.f && *abserr != 0.f) {
18622 d__1 = (double) (*abserr * 200.f / *resasc);
18623 r__1 = 1.f, r__2 = pow(d__1,
c_b270);
18624 *abserr = *resasc *
min(r__1,r__2);
18626 if (*resabs > uflow / (epmach * 50.f)) {
18628 r__1 = epmach * 50.f * *resabs;
18629 *abserr =
max(r__1,*abserr);
18639 sys_float xgk[16] = { .9980022986933971f,.9879925180204854f,
18640 .9677390756791391f,.9372733924007059f,.8972645323440819f,
18641 .8482065834104272f,.7904185014424659f,.72441773136017f,
18642 .650996741297417f,.5709721726085388f,.4850818636402397f,
18643 .3941513470775634f,.2991800071531688f,.2011940939974345f,
18644 .1011420669187175f,0.f };
18645 sys_float wgk[16] = { .005377479872923349f,.01500794732931612f,
18646 .02546084732671532f,.03534636079137585f,.04458975132476488f,
18647 .05348152469092809f,.06200956780067064f,.06985412131872826f,
18648 .07684968075772038f,.08308050282313302f,.08856444305621177f,
18649 .09312659817082532f,.09664272698362368f,.09917359872179196f,
18650 .1007698455238756f,.1013300070147915f };
18651 sys_float wg[8] = { .03075324199611727f,.07036604748810812f,
18652 .1071592204671719f,.1395706779261543f,.1662692058169939f,
18653 .1861610000155622f,.1984314853271116f,.2025782419255613f };
18663 sys_float absc, resg, resk, fsum, fval1, fval2;
18757 centr = (*a + *b) * .5f;
18758 hlgth = (*b - *a) * .5f;
18759 dhlgth = fabs(hlgth);
18766 resk = wgk[15] * fc;
18767 *resabs = fabs(resk);
18768 for (j = 1; j <= 7; ++j) {
18770 absc = hlgth * xgk[jtw - 1];
18771 r__1 = centr - absc;
18773 r__1 = centr + absc;
18775 fv1[jtw - 1] = fval1;
18776 fv2[jtw - 1] = fval2;
18777 fsum = fval1 + fval2;
18778 resg += wg[j - 1] * fsum;
18779 resk += wgk[jtw - 1] * fsum;
18780 *resabs += wgk[jtw - 1] * (fabs(fval1) + fabs(fval2));
18783 for (j = 1; j <= 8; ++j) {
18784 jtwm1 = (j << 1) - 1;
18785 absc = hlgth * xgk[jtwm1 - 1];
18786 r__1 = centr - absc;
18788 r__1 = centr + absc;
18790 fv1[jtwm1 - 1] = fval1;
18791 fv2[jtwm1 - 1] = fval2;
18792 fsum = fval1 + fval2;
18793 resk += wgk[jtwm1 - 1] * fsum;
18794 *resabs += wgk[jtwm1 - 1] * (fabs(fval1) + fabs(fval2));
18797 reskh = resk * .5f;
18798 *resasc = wgk[15] * (r__1 = fc - reskh, fabs(r__1));
18799 for (j = 1; j <= 15; ++j) {
18800 *resasc += wgk[j - 1] * ((r__1 = fv1[j - 1] - reskh, fabs(r__1)) + (
18801 r__2 = fv2[j - 1] - reskh, fabs(r__2)));
18804 *result = resk * hlgth;
18807 *abserr = (r__1 = (resk - resg) * hlgth, fabs(r__1));
18808 if (*resasc != 0.f && *abserr != 0.f) {
18810 d__1 = (double) (*abserr * 200.f / *resasc);
18811 r__1 = 1.f, r__2 = pow(d__1,
c_b270);
18812 *abserr = *resasc *
min(r__1,r__2);
18814 if (*resabs > uflow / (epmach * 50.f)) {
18816 r__1 = epmach * 50.f * *resabs;
18817 *abserr =
max(r__1,*abserr);
18827 sys_float xgk[21] = { .9988590315882777f,.9931285991850949f,
18828 .9815078774502503f,.9639719272779138f,.9408226338317548f,
18829 .9122344282513259f,.878276811252282f,.8391169718222188f,
18830 .7950414288375512f,.7463319064601508f,.6932376563347514f,
18831 .636053680726515f,.5751404468197103f,.5108670019508271f,
18832 .4435931752387251f,.3737060887154196f,.301627868114913f,
18833 .2277858511416451f,.1526054652409227f,.07652652113349733f,0.f };
18834 sys_float wgk[21] = { .003073583718520532f,.008600269855642942f,
18835 .01462616925697125f,.02038837346126652f,.02588213360495116f,
18836 .0312873067770328f,.0366001697582008f,.04166887332797369f,
18837 .04643482186749767f,.05094457392372869f,.05519510534828599f,
18838 .05911140088063957f,.06265323755478117f,.06583459713361842f,
18839 .06864867292852162f,.07105442355344407f,.07303069033278667f,
18840 .07458287540049919f,.07570449768455667f,.07637786767208074f,
18841 .07660071191799966f };
18842 sys_float wg[10] = { .01761400713915212f,.04060142980038694f,
18843 .06267204833410906f,.08327674157670475f,.1019301198172404f,
18844 .1181945319615184f,.1316886384491766f,.1420961093183821f,
18845 .1491729864726037f,.1527533871307259f };
18855 sys_float absc, resg, resk, fsum, fval1, fval2;
18952 centr = (*a + *b) * .5f;
18953 hlgth = (*b - *a) * .5f;
18954 dhlgth = fabs(hlgth);
18961 resk = wgk[20] * fc;
18962 *resabs = fabs(resk);
18963 for (j = 1; j <= 10; ++j) {
18965 absc = hlgth * xgk[jtw - 1];
18966 r__1 = centr - absc;
18968 r__1 = centr + absc;
18970 fv1[jtw - 1] = fval1;
18971 fv2[jtw - 1] = fval2;
18972 fsum = fval1 + fval2;
18973 resg += wg[j - 1] * fsum;
18974 resk += wgk[jtw - 1] * fsum;
18975 *resabs += wgk[jtw - 1] * (fabs(fval1) + fabs(fval2));
18978 for (j = 1; j <= 10; ++j) {
18979 jtwm1 = (j << 1) - 1;
18980 absc = hlgth * xgk[jtwm1 - 1];
18981 r__1 = centr - absc;
18983 r__1 = centr + absc;
18985 fv1[jtwm1 - 1] = fval1;
18986 fv2[jtwm1 - 1] = fval2;
18987 fsum = fval1 + fval2;
18988 resk += wgk[jtwm1 - 1] * fsum;
18989 *resabs += wgk[jtwm1 - 1] * (fabs(fval1) + fabs(fval2));
18992 reskh = resk * .5f;
18993 *resasc = wgk[20] * (r__1 = fc - reskh, fabs(r__1));
18994 for (j = 1; j <= 20; ++j) {
18995 *resasc += wgk[j - 1] * ((r__1 = fv1[j - 1] - reskh, fabs(r__1)) + (
18996 r__2 = fv2[j - 1] - reskh, fabs(r__2)));
18999 *result = resk * hlgth;
19002 *abserr = (r__1 = (resk - resg) * hlgth, fabs(r__1));
19003 if (*resasc != 0.f && *abserr != 0.f) {
19005 d__1 = (double) (*abserr * 200.f / *resasc);
19006 r__1 = 1.f, r__2 = pow(d__1,
c_b270);
19007 *abserr = *resasc *
min(r__1,r__2);
19009 if (*resabs > uflow / (epmach * 50.f)) {
19011 r__1 = epmach * 50.f * *resabs;
19012 *abserr =
max(r__1,*abserr);
19022 sys_float xgk[26] = { .9992621049926098f,.9955569697904981f,
19023 .9880357945340772f,.9766639214595175f,.9616149864258425f,
19024 .9429745712289743f,.9207471152817016f,.8949919978782754f,
19025 .8658470652932756f,.833442628760834f,.7978737979985001f,
19026 .7592592630373576f,.7177664068130844f,.6735663684734684f,
19027 .6268100990103174f,.577662930241223f,.5263252843347192f,
19028 .473002731445715f,.4178853821930377f,.3611723058093878f,
19029 .3030895389311078f,.2438668837209884f,.1837189394210489f,
19030 .1228646926107104f,.06154448300568508f,0.f };
19031 sys_float wgk[26] = { .001987383892330316f,.005561932135356714f,
19032 .009473973386174152f,.01323622919557167f,.0168478177091283f,
19033 .02043537114588284f,.02400994560695322f,.02747531758785174f,
19034 .03079230016738749f,.03400213027432934f,.03711627148341554f,
19035 .04008382550403238f,.04287284502017005f,.04550291304992179f,
19036 .04798253713883671f,.05027767908071567f,.05236288580640748f,
19037 .05425112988854549f,.05595081122041232f,.05743711636156783f,
19038 .05868968002239421f,.05972034032417406f,.06053945537604586f,
19039 .06112850971705305f,.06147118987142532f,.06158081806783294f };
19040 sys_float wg[13] = { .01139379850102629f,.02635498661503214f,
19041 .04093915670130631f,.05490469597583519f,.06803833381235692f,
19042 .08014070033500102f,.09102826198296365f,.1005359490670506f,
19043 .1085196244742637f,.1148582591457116f,.1194557635357848f,
19044 .12224244299031f,.1231760537267155f };
19054 sys_float absc, resg, resk, fsum, fval1, fval2;
19150 centr = (*a + *b) * .5f;
19151 hlgth = (*b - *a) * .5f;
19152 dhlgth = fabs(hlgth);
19158 resg = wg[12] * fc;
19159 resk = wgk[25] * fc;
19160 *resabs = fabs(resk);
19161 for (j = 1; j <= 12; ++j) {
19163 absc = hlgth * xgk[jtw - 1];
19164 r__1 = centr - absc;
19166 r__1 = centr + absc;
19168 fv1[jtw - 1] = fval1;
19169 fv2[jtw - 1] = fval2;
19170 fsum = fval1 + fval2;
19171 resg += wg[j - 1] * fsum;
19172 resk += wgk[jtw - 1] * fsum;
19173 *resabs += wgk[jtw - 1] * (fabs(fval1) + fabs(fval2));
19176 for (j = 1; j <= 13; ++j) {
19177 jtwm1 = (j << 1) - 1;
19178 absc = hlgth * xgk[jtwm1 - 1];
19179 r__1 = centr - absc;
19181 r__1 = centr + absc;
19183 fv1[jtwm1 - 1] = fval1;
19184 fv2[jtwm1 - 1] = fval2;
19185 fsum = fval1 + fval2;
19186 resk += wgk[jtwm1 - 1] * fsum;
19187 *resabs += wgk[jtwm1 - 1] * (fabs(fval1) + fabs(fval2));
19190 reskh = resk * .5f;
19191 *resasc = wgk[25] * (r__1 = fc - reskh, fabs(r__1));
19192 for (j = 1; j <= 25; ++j) {
19193 *resasc += wgk[j - 1] * ((r__1 = fv1[j - 1] - reskh, fabs(r__1)) + (
19194 r__2 = fv2[j - 1] - reskh, fabs(r__2)));
19197 *result = resk * hlgth;
19200 *abserr = (r__1 = (resk - resg) * hlgth, fabs(r__1));
19201 if (*resasc != 0.f && *abserr != 0.f) {
19203 d__1 = (double) (*abserr * 200.f / *resasc);
19204 r__1 = 1.f, r__2 = pow(d__1,
c_b270);
19205 *abserr = *resasc *
min(r__1,r__2);
19207 if (*resabs > uflow / (epmach * 50.f)) {
19209 r__1 = epmach * 50.f * *resabs;
19210 *abserr =
max(r__1,*abserr);
19220 sys_float xgk[31] = { .9994844100504906f,.9968934840746495f,
19221 .9916309968704046f,.9836681232797472f,.9731163225011263f,
19222 .9600218649683075f,.94437444474856f,.9262000474292743f,
19223 .9055733076999078f,.8825605357920527f,.8572052335460611f,
19224 .8295657623827684f,.7997278358218391f,.7677774321048262f,
19225 .7337900624532268f,.6978504947933158f,.660061064126627f,
19226 .6205261829892429f,.5793452358263617f,.5366241481420199f,
19227 .4924804678617786f,.4470337695380892f,.4004012548303944f,
19228 .3527047255308781f,.3040732022736251f,.2546369261678898f,
19229 .2045251166823099f,.1538699136085835f,.102806937966737f,
19230 .0514718425553177f,0.f };
19231 sys_float wgk[31] = { .001389013698677008f,.003890461127099884f,
19232 .006630703915931292f,.009273279659517763f,.01182301525349634f,
19233 .0143697295070458f,.01692088918905327f,.01941414119394238f,
19234 .02182803582160919f,.0241911620780806f,.0265099548823331f,
19235 .02875404876504129f,.03090725756238776f,.03298144705748373f,
19236 .03497933802806002f,.03688236465182123f,.03867894562472759f,
19237 .04037453895153596f,.04196981021516425f,.04345253970135607f,
19238 .04481480013316266f,.04605923827100699f,.04718554656929915f,
19239 .04818586175708713f,.04905543455502978f,.04979568342707421f,
19240 .05040592140278235f,.05088179589874961f,.05122154784925877f,
19241 .05142612853745903f,.05149472942945157f };
19242 sys_float wg[15] = { .007968192496166606f,.01846646831109096f,
19243 .02878470788332337f,.03879919256962705f,.04840267283059405f,
19244 .05749315621761907f,.0659742298821805f,.07375597473770521f,
19245 .08075589522942022f,.08689978720108298f,.09212252223778613f,
19246 .09636873717464426f,.09959342058679527f,.1017623897484055f,
19247 .1028526528935588f };
19257 sys_float absc, resg, resk, fsum, fval1, fval2;
19353 centr = (*b + *a) * .5f;
19354 hlgth = (*b - *a) * .5f;
19355 dhlgth = fabs(hlgth);
19362 resk = wgk[30] * fc;
19363 *resabs = fabs(resk);
19364 for (j = 1; j <= 15; ++j) {
19366 absc = hlgth * xgk[jtw - 1];
19367 r__1 = centr - absc;
19369 r__1 = centr + absc;
19371 fv1[jtw - 1] = fval1;
19372 fv2[jtw - 1] = fval2;
19373 fsum = fval1 + fval2;
19374 resg += wg[j - 1] * fsum;
19375 resk += wgk[jtw - 1] * fsum;
19376 *resabs += wgk[jtw - 1] * (fabs(fval1) + fabs(fval2));
19379 for (j = 1; j <= 15; ++j) {
19380 jtwm1 = (j << 1) - 1;
19381 absc = hlgth * xgk[jtwm1 - 1];
19382 r__1 = centr - absc;
19384 r__1 = centr + absc;
19386 fv1[jtwm1 - 1] = fval1;
19387 fv2[jtwm1 - 1] = fval2;
19388 fsum = fval1 + fval2;
19389 resk += wgk[jtwm1 - 1] * fsum;
19390 *resabs += wgk[jtwm1 - 1] * (fabs(fval1) + fabs(fval2));
19393 reskh = resk * .5f;
19394 *resasc = wgk[30] * (r__1 = fc - reskh, fabs(r__1));
19395 for (j = 1; j <= 30; ++j) {
19396 *resasc += wgk[j - 1] * ((r__1 = fv1[j - 1] - reskh, fabs(r__1)) + (
19397 r__2 = fv2[j - 1] - reskh, fabs(r__2)));
19400 *result = resk * hlgth;
19403 *abserr = (r__1 = (resk - resg) * hlgth, fabs(r__1));
19404 if (*resasc != 0.f && *abserr != 0.f) {
19406 d__1 = (double) (*abserr * 200.f / *resasc);
19407 r__1 = 1.f, r__2 = pow(d__1,
c_b270);
19408 *abserr = *resasc *
min(r__1,r__2);
19410 if (*resabs > uflow / (epmach * 50.f)) {
19412 r__1 = epmach * 50.f * *resabs;
19413 *abserr =
max(r__1,*abserr);
19428 sys_float anm1, ralf, rbet, alfp1, alfp2, betp1, betp2;
19496 alfp1 = *alfa + 1.f;
19497 betp1 = *beta + 1.f;
19498 alfp2 = *alfa + 2.f;
19499 betp2 = *beta + 2.f;
19500 d__1 = (double) alfp1;
19501 ralf = pow(
c_b320, d__1);
19502 d__1 = (double) betp1;
19503 rbet = pow(
c_b320, d__1);
19507 ri[1] = ralf / alfp1;
19508 rj[1] = rbet / betp1;
19509 ri[2] = ri[1] * *alfa / alfp2;
19510 rj[2] = rj[1] * *beta / betp2;
19513 for (i__ = 3; i__ <= 25; ++i__) {
19514 ri[i__] = -(ralf + an * (an - alfp2) * ri[i__ - 1]) / (anm1 * (an +
19516 rj[i__] = -(rbet + an * (an - betp2) * rj[i__ - 1]) / (anm1 * (an +
19522 if (*integr == 1) {
19525 if (*integr == 3) {
19531 rg[1] = -ri[1] / alfp1;
19532 rg[2] = -(ralf + ralf) / (alfp2 * alfp2) - rg[1];
19536 for (i__ = 3; i__ <= 25; ++i__) {
19537 rg[i__] = -(an * (an - alfp2) * rg[im1] - an * ri[im1] + anm1 * ri[
19538 i__]) / (anm1 * (an + alfp1));
19544 if (*integr == 2) {
19551 rh[1] = -rj[1] / betp1;
19552 rh[2] = -(rbet + rbet) / (betp2 * betp2) - rh[1];
19556 for (i__ = 3; i__ <= 25; ++i__) {
19557 rh[i__] = -(an * (an - betp2) * rh[im1] - an * rj[im1] + anm1 * rj[
19558 i__]) / (anm1 * (an + betp1));
19564 for (i__ = 2; i__ <= 25; i__ += 2) {
19565 rh[i__] = -rh[i__];
19569 for (i__ = 2; i__ <= 25; i__ += 2) {
19570 rj[i__] = -rj[i__];
19582 sys_float x1[5] = { .9739065285171717f,.8650633666889845f,
19583 .6794095682990244f,.4333953941292472f,.1488743389816312f };
19584 sys_float x2[5] = { .9956571630258081f,.9301574913557082f,
19585 .7808177265864169f,.5627571346686047f,.2943928627014602f };
19586 sys_float x3[11] = { .9993333609019321f,.9874334029080889f,
19587 .9548079348142663f,.9001486957483283f,.8251983149831142f,
19588 .732148388989305f,.6228479705377252f,.4994795740710565f,
19589 .3649016613465808f,.2222549197766013f,.07465061746138332f };
19590 sys_float x4[22] = { .9999029772627292f,.9979898959866787f,
19591 .9921754978606872f,.9813581635727128f,.9650576238583846f,
19592 .9431676131336706f,.9158064146855072f,.8832216577713165f,
19593 .8457107484624157f,.803557658035231f,.7570057306854956f,
19594 .7062732097873218f,.6515894665011779f,.5932233740579611f,
19595 .5314936059708319f,.4667636230420228f,.3994248478592188f,
19596 .3298748771061883f,.2585035592021616f,.1856953965683467f,
19597 .1118422131799075f,.03735212339461987f };
19598 sys_float w10[5] = { .06667134430868814f,.1494513491505806f,
19599 .219086362515982f,.2692667193099964f,.2955242247147529f };
19600 sys_float w21a[5] = { .03255816230796473f,.07503967481091995f,
19601 .1093871588022976f,.1347092173114733f,.1477391049013385f };
19602 sys_float w21b[6] = { .01169463886737187f,.054755896574352f,
19603 .09312545458369761f,.1234919762620659f,.1427759385770601f,
19604 .1494455540029169f };
19605 sys_float w43a[10] = { .01629673428966656f,.0375228761208695f,
19606 .05469490205825544f,.06735541460947809f,.07387019963239395f,
19607 .005768556059769796f,.02737189059324884f,.04656082691042883f,
19608 .06174499520144256f,.0713872672686934f };
19609 sys_float w43b[12] = { .001844477640212414f,.01079868958589165f,
19610 .02189536386779543f,.03259746397534569f,.04216313793519181f,
19611 .05074193960018458f,.05837939554261925f,.06474640495144589f,
19612 .06956619791235648f,.07282444147183321f,.07450775101417512f,
19613 .07472214751740301f };
19614 sys_float w87a[21] = { .008148377384149173f,.01876143820156282f,
19615 .02734745105005229f,.03367770731163793f,.03693509982042791f,
19616 .002884872430211531f,.0136859460227127f,.02328041350288831f,
19617 .03087249761171336f,.03569363363941877f,9.152833452022414e-4f,
19618 .005399280219300471f,.01094767960111893f,.01629873169678734f,
19619 .02108156888920384f,.02537096976925383f,.02918969775647575f,
19620 .03237320246720279f,.03478309895036514f,.03641222073135179f,
19621 .03725387550304771f };
19622 sys_float w87b[23] = { 2.741455637620724e-4f,.001807124155057943f,
19623 .004096869282759165f,.006758290051847379f,.009549957672201647f,
19624 .01232944765224485f,.01501044734638895f,.01754896798624319f,
19625 .01993803778644089f,.02219493596101229f,.02433914712600081f,
19626 .02637450541483921f,.0282869107887712f,.0300525811280927f,
19627 .03164675137143993f,.0330504134199785f,.03425509970422606f,
19628 .03526241266015668f,.0360769896228887f,.03669860449845609f,
19629 .03712054926983258f,.03733422875193504f,.03736107376267902f };
19637 sys_float fv1[5], fv2[5], fv3[5], fv4[5];
19639 sys_float absc, fval, res10, res21=0., res43=0., res87, fval1, fval2, hlgth,
19640 centr, reskh, uflow;
19641 sys_float epmach, dhlgth, resabs=0., resasc=0., fcentr, savfun[21];
19781 r__1 = 5e-15f, r__2 = epmach * 50.f;
19782 if (*epsabs <= 0.f && *epsrel <
max(r__1,r__2)) {
19785 hlgth = (*b - *a) * .5f;
19786 dhlgth = fabs(hlgth);
19787 centr = (*b + *a) * .5f;
19794 for (l = 1; l <= 3; ++l) {
19802 res21 = w21b[5] * fcentr;
19803 resabs = w21b[5] * fabs(fcentr);
19804 for (k = 1; k <= 5; ++k) {
19805 absc = hlgth * x1[k - 1];
19806 r__1 = centr + absc;
19808 r__1 = centr - absc;
19810 fval = fval1 + fval2;
19811 res10 += w10[k - 1] * fval;
19812 res21 += w21a[k - 1] * fval;
19813 resabs += w21a[k - 1] * (fabs(fval1) + fabs(fval2));
19814 savfun[k - 1] = fval;
19815 fv1[k - 1] = fval1;
19816 fv2[k - 1] = fval2;
19820 for (k = 1; k <= 5; ++k) {
19822 absc = hlgth * x2[k - 1];
19823 r__1 = centr + absc;
19825 r__1 = centr - absc;
19827 fval = fval1 + fval2;
19828 res21 += w21b[k - 1] * fval;
19829 resabs += w21b[k - 1] * (fabs(fval1) + fabs(fval2));
19830 savfun[ipx - 1] = fval;
19831 fv3[k - 1] = fval1;
19832 fv4[k - 1] = fval2;
19838 *result = res21 * hlgth;
19840 reskh = res21 * .5f;
19841 resasc = w21b[5] * (r__1 = fcentr - reskh, fabs(r__1));
19842 for (k = 1; k <= 5; ++k) {
19843 resasc = resasc + w21a[k - 1] * ((r__1 = fv1[k - 1] - reskh, fabs(
19844 r__1)) + (r__2 = fv2[k - 1] - reskh, fabs(r__2))) + w21b[
19845 k - 1] * ((r__3 = fv3[k - 1] - reskh, fabs(r__3)) + (r__4
19846 = fv4[k - 1] - reskh, fabs(r__4)));
19849 *abserr = (r__1 = (res21 - res10) * hlgth, fabs(r__1));
19856 res43 = w43b[11] * fcentr;
19858 for (k = 1; k <= 10; ++k) {
19859 res43 += savfun[k - 1] * w43a[k - 1];
19862 for (k = 1; k <= 11; ++k) {
19864 absc = hlgth * x3[k - 1];
19865 r__1 = absc + centr;
19866 r__2 = centr - absc;
19867 fval = f(r__1) + f(r__2);
19868 res43 += fval * w43b[k - 1];
19869 savfun[ipx - 1] = fval;
19875 *result = res43 * hlgth;
19876 *abserr = (r__1 = (res43 - res21) * hlgth, fabs(r__1));
19882 res87 = w87b[22] * fcentr;
19884 for (k = 1; k <= 21; ++k) {
19885 res87 += savfun[k - 1] * w87a[k - 1];
19888 for (k = 1; k <= 22; ++k) {
19889 absc = hlgth * x4[k - 1];
19890 r__1 = absc + centr;
19891 r__2 = centr - absc;
19892 res87 += w87b[k - 1] * (f(r__1) + f(r__2));
19895 *result = res87 * hlgth;
19896 *abserr = (r__1 = (res87 - res43) * hlgth, fabs(r__1));
19898 if (resasc != 0.f && *abserr != 0.f) {
19900 d__1 = (double) (*abserr * 200.f / resasc);
19901 r__1 = 1.f, r__2 = pow(d__1,
c_b270);
19902 *abserr = resasc *
min(r__1,r__2);
19904 if (resabs > uflow / (epmach * 50.f)) {
19906 r__1 = epmach * 50.f * resabs;
19907 *abserr =
max(r__1,*abserr);
19910 r__1 = *epsabs, r__2 = *epsrel * fabs(*result);
19911 if (*abserr <=
max(r__1,r__2)) {
19926 void qpsrt_(
const long *limit,
long *last,
long *maxerr,
19933 long i__, j, k, ido, ibeg, jbnd, isucc, jupbn;
20018 errmax = elist[*maxerr];
20024 for (i__ = 1; i__ <= i__1; ++i__) {
20025 isucc = iord[*nrmax - 1];
20027 if (errmax <= elist[isucc]) {
20030 iord[*nrmax] = isucc;
20042 if (*last > *limit / 2 + 2) {
20043 jupbn = *limit + 3 - *last;
20045 errmin = elist[*last];
20056 for (i__ = ibeg; i__ <= i__1; ++i__) {
20059 if (errmax >= elist[isucc]) {
20062 iord[i__ - 1] = isucc;
20066 iord[jbnd] = *maxerr;
20067 iord[jupbn] = *last;
20073 iord[i__ - 1] = *maxerr;
20076 for (j = i__; j <= i__1; ++j) {
20079 if (errmin < elist[isucc]) {
20082 iord[k + 1] = isucc;
20089 iord[k + 1] = *last;
20094 *maxerr = iord[*nrmax];
20095 *ermax = elist[*maxerr];
20117 ret_val = 1.f / (*x - *c__);
20146 ret_val = cos(omx);
20149 ret_val = sin(omx);
20156 const long *integr)
20160 double d__1, d__2, d__3, d__4;
20180 d__1 = (double) xma;
20181 d__2 = (double) (*alfa);
20182 d__3 = (double) bmx;
20183 d__4 = (double) (*beta);
20184 ret_val = pow(d__1, d__2) * pow(d__3, d__4);
20192 ret_val *= log(xma);
20195 ret_val *= log(bmx);
20198 ret_val = ret_val * log(xma) * log(bmx);
20213 long kb, kp1, nm1, nm2;
20279 for (k = 1; k <= i__1; ++k) {
20284 if ((r__1 = c__[kp1], fabs(r__1)) < (r__2 = c__[k], fabs(r__2))) {
20306 if (c__[k] != 0.f) {
20313 t = -c__[kp1] / c__[k];
20314 c__[kp1] = d__[kp1] + t * d__[k];
20315 d__[kp1] = e[kp1] + t * e[k];
20317 b[kp1] += t * b[k];
20321 if (c__[*n] != 0.f) {
20335 b[nm1] = (b[nm1] - d__[nm1] * b[*n]) / c__[nm1];
20340 for (kb = 1; kb <= i__1; ++kb) {
20342 b[k] = (b[k] - d__[k] * b[k + 1] - e[k] * b[k + 2]) / c__[k];
20353 void dgtsl_(
const long *n,
double *c__,
double *d__,
20354 double *e,
double *b,
long *info)
20363 long kb, kp1, nm1, nm2;
20429 for (k = 1; k <= i__1; ++k) {
20434 if ((d__1 = c__[kp1], fabs(d__1)) < (d__2 = c__[k], fabs(d__2))) {
20456 if (c__[k] != 0.) {
20463 t = -c__[kp1] / c__[k];
20464 c__[kp1] = d__[kp1] + t * d__[k];
20465 d__[kp1] = e[kp1] + t * e[k];
20467 b[kp1] += t * b[k];
20471 if (c__[*n] != 0.) {
20485 b[nm1] = (b[nm1] - d__[nm1] * b[*n]) / c__[nm1];
20490 for (kb = 1; kb <= i__1; ++kb) {
20492 b[k] = (b[k] - d__[k] * b[k + 1] - e[k] * b[k + 2]) / c__[k];
void qk15i_(const E_fp &f, const sys_float *boun, const long *inf, const sys_float *a, const sys_float *b, sys_float *result, sys_float *abserr, sys_float *resabs, sys_float *resasc)
sys_float(* E_fp1)(const sys_float *, const sys_float *, const sys_float *, const sys_float *, const sys_float *, const long *)
void qc25f_(const E_fp &f, const sys_float *a, const sys_float *b, const sys_float *omega, const long *integr, const long *nrmom, const long *maxp1, const long *ksave, sys_float *result, sys_float *abserr, long *neval, sys_float *resabs, sys_float *resasc, long *momcom, sys_float *chebmo)
double dqwgtc_(const double *x, const double *c__, const double *, const double *, const double *, const long *)
void dqk61_(const D_fp &f, const double *a, const double *b, double *result, double *abserr, double *resabs, double *resasc)
static void sgtsl_(const long *, sys_float *, sys_float *, sys_float *, sys_float *, long *)
void dqk15w_(const D_fp &f, D_fp1 w, const double *p1, const double *p2, const double *p3, const double *p4, const long *kp, const double *a, const double *b, double *result, double *abserr, double *resabs, double *resasc)
void qk51_(const E_fp &f, const sys_float *a, const sys_float *b, sys_float *result, sys_float *abserr, sys_float *resabs, sys_float *resasc)
void qelg_(long *n, sys_float *epstab, sys_float *result, sys_float *abserr, sys_float *res3la, long *nres)
void qk15w_(const E_fp &f, E_fp1 w, const sys_float *p1, const sys_float *p2, const sys_float *p3, const sys_float *p4, const long *kp, const sys_float *a, const sys_float *b, sys_float *result, sys_float *abserr, sys_float *resabs, sys_float *resasc)
void qagse_(const E_fp &f, const sys_float *a, const sys_float *b, const sys_float *epsabs, const sys_float *epsrel, const long *limit, sys_float *result, sys_float *abserr, long *neval, long *ier, sys_float *alist__, sys_float *blist, sys_float *rlist, sys_float *elist, long *iord, long *last)
void qk21_(const E_fp &f, const sys_float *a, const sys_float *b, sys_float *result, sys_float *abserr, sys_float *resabs, sys_float *resasc)
void dqawoe_(const D_fp &f, const double *a, const double *b, const double *omega, const long *integr, const double *epsabs, const double *epsrel, const long *limit, const long *icall, const long *maxp1, double *result, double *abserr, long *neval, long *ier, long *last, double *alist__, double *blist, double *rlist, double *elist, long *iord, long *nnlog, long *momcom, double *chebmo)
static const double c_b20
double(* D_fp1)(const double *, const double *, const double *, const double *, const double *, const long *)
void dqagse_(const D_fp &f, const double *a, const double *b, const double *epsabs, const double *epsrel, const long *limit, double *result, double *abserr, long *neval, long *ier, double *alist__, double *blist, double *rlist, double *elist, long *iord, long *last)
void dqag_(const D_fp &f, const double *a, const double *b, const double *epsabs, const double *epsrel, const long *key, double *result, double *abserr, long *neval, long *ier, long *limit, const long *lenw, long *last, long *iwork, double *work)
void dqagie_(const D_fp &f, const double *bound, const long *inf, const double *epsabs, const double *epsrel, const long *limit, double *result, double *abserr, long *neval, long *ier, double *alist__, double *blist, double *rlist, double *elist, long *iord, long *last)
static const double c_b21
void qage_(const E_fp &f, const sys_float *a, const sys_float *b, const sys_float *epsabs, const sys_float *epsrel, const long *key, const long *limit, sys_float *result, sys_float *abserr, long *neval, long *ier, sys_float *alist__, sys_float *blist, sys_float *rlist, sys_float *elist, long *iord, long *last)
void dqk15_(const D_fp &f, const double *a, const double *b, double *result, double *abserr, double *resabs, double *resasc)
static const double c_b320
void dqmomo_(const double *alfa, const double *beta, double *ri, double *rj, double *rg, double *rh, const long *integr)
static double d1mach(long i)
void qk41_(const E_fp &f, const sys_float *a, const sys_float *b, sys_float *result, sys_float *abserr, sys_float *resabs, sys_float *resasc)
static sys_float r1mach(long i)
double dqwgtf_(const double *x, const double *omega, const double *, const double *, const double *, const long *integr)
void qawse_(const E_fp &f, const sys_float *a, const sys_float *b, const sys_float *alfa, const sys_float *beta, const long *integr, const sys_float *epsabs, const sys_float *epsrel, const long *limit, sys_float *result, sys_float *abserr, long *neval, long *ier, sys_float *alist__, sys_float *blist, sys_float *rlist, sys_float *elist, long *iord, long *last)
void qc25c_(const E_fp &f, const sys_float *a, const sys_float *b, const sys_float *c__, sys_float *result, sys_float *abserr, long *krul, long *neval)
void dqagp_(const D_fp &f, const double *a, const double *b, const long *npts2, const double *points, const double *epsabs, const double *epsrel, double *result, double *abserr, long *neval, long *ier, const long *leniw, const long *lenw, long *last, long *iwork, double *work)
sys_float qwgtf_(const sys_float *x, const sys_float *omega, const sys_float *, const sys_float *, const sys_float *, const long *integr)
void dqawc_(const D_fp &f, const double *a, const double *b, const double *c__, const double *epsabs, const double *epsrel, double *result, double *abserr, long *neval, long *ier, long *limit, const long *lenw, long *last, long *iwork, double *work)
void qagie_(const E_fp &f, const sys_float *bound, const long *inf, const sys_float *epsabs, const sys_float *epsrel, const long *limit, sys_float *result, sys_float *abserr, long *neval, long *ier, sys_float *alist__, sys_float *blist, sys_float *rlist, sys_float *elist, long *iord, long *last)
static void dgtsl_(const long *, double *, double *, double *, double *, long *)
void qk15_(const E_fp &f, const sys_float *a, const sys_float *b, sys_float *result, sys_float *abserr, sys_float *resabs, sys_float *resasc)
void dqage_(const D_fp &f, const double *a, const double *b, const double *epsabs, const double *epsrel, const long *key, const long *limit, double *result, double *abserr, long *neval, long *ier, double *alist__, double *blist, double *rlist, double *elist, long *iord, long *last)
sys_float qwgtc_(const sys_float *x, const sys_float *c__, const sys_float *, const sys_float *, const sys_float *, const long *)
void qaws_(const E_fp &f, const sys_float *a, const sys_float *b, const sys_float *alfa, const sys_float *beta, const long *integr, const sys_float *epsabs, const sys_float *epsrel, sys_float *result, sys_float *abserr, long *neval, long *ier, const long *limit, const long *lenw, long *last, long *iwork, sys_float *work)
STATIC void sinpar(double, double, double, double *, double *, double *, double *, double *, long *)
void dqawo_(const D_fp &f, const double *a, const double *b, const double *omega, const long *integr, const double *epsabs, const double *epsrel, double *result, double *abserr, long *neval, long *ier, const long *leniw, long *maxp1, const long *lenw, long *last, long *iwork, double *work)
void dqawfe_(const D_fp &f, const double *a, const double *omega, const long *integr, const double *epsabs, const long *limlst, const long *limit, const long *maxp1, double *result, double *abserr, long *neval, long *ier, double *rslst, double *erlst, long *ierlst, long *lst, double *alist__, double *blist, double *rlist, double *elist, long *iord, long *nnlog, double *chebmo)
static const sys_float c_b390
void qk61_(const E_fp &f, const sys_float *a, const sys_float *b, sys_float *result, sys_float *abserr, sys_float *resabs, sys_float *resasc)
void dqc25s_(const D_fp &f, const double *a, const double *b, const double *bl, const double *br, const double *alfa, const double *beta, const double *ri, const double *rj, const double *rg, const double *rh, double *result, double *abserr, double *resasc, const long *integr, long *nev)
void dqags_(const D_fp &f, const double *a, const double *b, const double *epsabs, const double *epsrel, double *result, double *abserr, long *neval, long *ier, const long *limit, const long *lenw, long *last, long *iwork, double *work)
void dqawf_(const D_fp &f, const double *a, const double *omega, const long *integr, const double *epsabs, double *result, double *abserr, long *neval, long *ier, long *limlst, long *lst, const long *leniw, const long *maxp1, const long *lenw, long *iwork, double *work)
void qc25s_(const E_fp &f, const sys_float *a, const sys_float *b, const sys_float *bl, const sys_float *br, const sys_float *alfa, const sys_float *beta, sys_float *ri, sys_float *rj, sys_float *rg, sys_float *rh, sys_float *result, sys_float *abserr, sys_float *resasc, const long *integr, long *nev)
void qawo_(const E_fp &f, const sys_float *a, const sys_float *b, const sys_float *omega, const long *integr, const sys_float *epsabs, const sys_float *epsrel, sys_float *result, sys_float *abserr, long *neval, long *ier, const long *leniw, const long *maxp1, const long *lenw, long *last, long *iwork, sys_float *work)
sys_float qwgts_(const sys_float *x, const sys_float *a, const sys_float *b, const sys_float *alfa, const sys_float *beta, const long *integr)
void dqaws_(const D_fp &f, const double *a, const double *b, const double *alfa, const double *beta, const long *integr, const double *epsabs, const double *epsrel, double *result, double *abserr, long *neval, long *ier, long *limit, const long *lenw, long *last, long *iwork, double *work)
void dqawce_(const D_fp &f, const double *a, const double *b, const double *c__, const double *epsabs, const double *epsrel, const long *limit, double *result, double *abserr, long *neval, long *ier, double *alist__, double *blist, double *rlist, double *elist, long *iord, long *last)
void qcheb_(sys_float *x, sys_float *fval, sys_float *cheb12, sys_float *cheb24)
void qawoe_(const E_fp &f, const sys_float *a, const sys_float *b, const sys_float *omega, const long *integr, const sys_float *epsabs, const sys_float *epsrel, const long *limit, const long *icall, const long *maxp1, sys_float *result, sys_float *abserr, long *neval, long *ier, long *last, sys_float *alist__, sys_float *blist, sys_float *rlist, sys_float *elist, long *iord, long *nnlog, long *momcom, sys_float *chebmo)
void dqng_(const D_fp &f, const double *a, const double *b, const double *epsabs, const double *epsrel, double *result, double *abserr, long *neval, long *ier)
void qawc_(const E_fp &f, const sys_float *a, const sys_float *b, const sys_float *c__, const sys_float *epsabs, const sys_float *epsrel, sys_float *result, sys_float *abserr, long *neval, long *ier, long *limit, const long *lenw, long *last, long *iwork, sys_float *work)
void qmomo_(const sys_float *alfa, const sys_float *beta, sys_float *ri, sys_float *rj, sys_float *rg, sys_float *rh, const long *integr)
void qagpe_(const E_fp &f, const sys_float *a, const sys_float *b, const long *npts2, const sys_float *points, const sys_float *epsabs, const sys_float *epsrel, const long *limit, sys_float *result, sys_float *abserr, long *neval, long *ier, sys_float *alist__, sys_float *blist, sys_float *rlist, sys_float *elist, sys_float *pts, long *iord, long *level, long *ndin, long *last)
void dqk51_(const D_fp &f, const double *a, const double *b, double *result, double *abserr, double *resabs, double *resasc)
void dqpsrt_(const long *limit, long *last, long *maxerr, double *ermax, double *elist, long *iord, long *nrmax)
void dqc25f_(const D_fp &f, const double *a, const double *b, const double *omega, const long *integr, const long *nrmom, const long *maxp1, const long *ksave, double *result, double *abserr, long *neval, double *resabs, double *resasc, long *momcom, double *chebmo)
void dqk15i_(const D_fp &f, const double *boun, const long *inf, const double *a, const double *b, double *result, double *abserr, double *resabs, double *resasc)
void dqelg_(long *n, double *epstab, double *result, double *abserr, double *res3la, long *nres)
double dqwgts_(const double *x, const double *a, const double *b, const double *alfa, const double *beta, const long *integr)
void qpsrt_(const long *limit, long *last, long *maxerr, sys_float *ermax, sys_float *elist, long *iord, long *nrmax)
void qags_(const E_fp &f, const sys_float *a, const sys_float *b, const sys_float *epsabs, const sys_float *epsrel, sys_float *result, sys_float *abserr, long *neval, long *ier, const long *limit, const long *lenw, long *last, long *iwork, sys_float *work)
void qagp_(const E_fp &f, const sys_float *a, const sys_float *b, const long *npts2, const sys_float *points, const sys_float *epsabs, const sys_float *epsrel, sys_float *result, sys_float *abserr, long *neval, long *ier, const long *leniw, const long *lenw, long *last, long *iwork, sys_float *work)
#define DEBUG_ENTRY(funcname)
static const double c_b270
static const sys_float c_b391
void dqk21_(const D_fp &f, const double *a, const double *b, double *result, double *abserr, double *resabs, double *resasc)
void qawf_(const E_fp &f, const sys_float *a, const sys_float *omega, const long *integr, const sys_float *epsabs, sys_float *result, sys_float *abserr, long *neval, long *ier, const long *limlst, long *lst, const long *leniw, const long *maxp1, const long *lenw, long *iwork, sys_float *work)
int fprintf(const Output &stream, const char *format,...)
void dqagi_(const D_fp &f, const double *bound, const long *inf, const double *epsabs, const double *epsrel, double *result, double *abserr, long *neval, long *ier, long *limit, const long *lenw, long *last, long *iwork, double *work)
void qag_(const E_fp &f, const sys_float *a, const sys_float *b, const sys_float *epsabs, const sys_float *epsrel, const long *key, sys_float *result, sys_float *abserr, long *neval, long *ier, long *limit, const long *lenw, long *last, long *iwork, sys_float *work)
void dqk41_(const D_fp &f, const double *a, const double *b, double *result, double *abserr, double *resabs, double *resasc)
void dqcheb_(const double *x, double *fval, double *cheb12, double *cheb24)
void dqk31_(const D_fp &f, const double *a, const double *b, double *result, double *abserr, double *resabs, double *resasc)
void qagi_(const E_fp &f, const sys_float *bound, const long *inf, const sys_float *epsabs, const sys_float *epsrel, sys_float *result, sys_float *abserr, long *neval, long *ier, const long *limit, const long *lenw, long *last, long *iwork, sys_float *work)
STATIC long xerror_(const char *mesg, const long *nmesg, const long *nerr, const long *level, long)
void qawfe_(const E_fp &f, const sys_float *a, const sys_float *omega, const long *integr, const sys_float *epsabs, const long *limlst, const long *limit, const long *maxp1, sys_float *result, sys_float *abserr, long *neval, long *ier, sys_float *rslst, sys_float *erlst, long *ierlst, long *lst, sys_float *alist__, sys_float *blist, sys_float *rlist, sys_float *elist, long *iord, long *nnlog, sys_float *chebmo)
void qawce_(const E_fp &f, const sys_float *a, const sys_float *b, const sys_float *c__, const sys_float *epsabs, const sys_float *epsrel, const long *limit, sys_float *result, sys_float *abserr, long *neval, long *ier, sys_float *alist__, sys_float *blist, sys_float *rlist, sys_float *elist, long *iord, long *last)
void dqc25c_(const D_fp &f, const double *a, const double *b, const double *c__, double *result, double *abserr, long *krul, long *neval)
void qng_(const E_fp &f, sys_float *a, sys_float *b, sys_float *epsabs, sys_float *epsrel, sys_float *result, sys_float *abserr, long *neval, long *ier)
void qk31_(const E_fp &f, const sys_float *a, const sys_float *b, sys_float *result, sys_float *abserr, sys_float *resabs, sys_float *resasc)
void dqagpe_(const D_fp &f, const double *a, const double *b, const long *npts2, const double *points, const double *epsabs, const double *epsrel, const long *limit, double *result, double *abserr, long *neval, long *ier, double *alist__, double *blist, double *rlist, double *elist, double *pts, long *iord, long *level, long *ndin, long *last)
void dqawse_(const D_fp &f, const double *a, const double *b, const double *alfa, const double *beta, const long *integr, const double *epsabs, const double *epsrel, const long *limit, double *result, double *abserr, long *neval, long *ier, double *alist__, double *blist, double *rlist, double *elist, long *iord, long *last)