23 #ifndef O2SCL_GSL_INTE_QAWC_H
24 #define O2SCL_GSL_INTE_QAWC_H
30 #include <o2scl/err_hnd.h>
31 #include <o2scl/inte.h>
32 #include <o2scl/inte_singular_gsl.h>
34 #ifndef DOXYGEN_NO_O2NS
54 double a0 = log (fabs ((1.0 - cc) / (1.0 + cc)));
55 double a1 = 2 + a0 * cc;
60 for (k = 2; k < 25; k++) {
64 a2 = 2.0 * cc * a1 - a0;
66 const double km1 = k - 1.0;
67 a2 = 2.0 * cc * a1 - a0 - 4.0 / (km1 * km1 - 1.0);
113 template<
class func2_t>
115 double *cheb12,
double *cheb24) {
117 double fval[25], v[12];
122 const double x[11] = { 0.9914448613738104,
132 0.1305261922200516 };
134 const double center = 0.5 * (b + a);
135 const double half_length = 0.5 * (b - a);
145 for (i = 1; i < 12; i++) {
146 const size_t j = 24 - i;
147 const double u = half_length * x[i-1];
155 for (i = 0; i < 12; i++) {
156 const size_t j = 24 - i;
157 v[i] = fval[i] - fval[j];
158 fval[i] = fval[i] + fval[j];
162 const double alam1 = v[0] - v[8];
163 const double alam2 = x[5] * (v[2] - v[6] - v[10]);
165 cheb12[3] = alam1 + alam2;
166 cheb12[9] = alam1 - alam2;
170 const double alam1 = v[1] - v[7] - v[9];
171 const double alam2 = v[3] - v[5] - v[11];
173 const double alam = x[2] * alam1 + x[8] * alam2;
175 cheb24[3] = cheb12[3] + alam;
176 cheb24[21] = cheb12[3] - alam;
180 const double alam = x[8] * alam1 - x[2] * alam2;
181 cheb24[9] = cheb12[9] + alam;
182 cheb24[15] = cheb12[9] - alam;
187 const double part1 = x[3] * v[4];
188 const double part2 = x[7] * v[8];
189 const double part3 = x[5] * v[6];
192 const double alam1 = v[0] + part1 + part2;
193 const double alam2 = x[1] * v[2] + part3 + x[9] * v[10];
195 cheb12[1] = alam1 + alam2;
196 cheb12[11] = alam1 - alam2;
200 const double alam1 = v[0] - part1 + part2;
201 const double alam2 = x[9] * v[2] - part3 + x[1] * v[10];
202 cheb12[5] = alam1 + alam2;
203 cheb12[7] = alam1 - alam2;
208 const double alam = (x[0] * v[1] + x[2] * v[3] + x[4] * v[5]
209 + x[6] * v[7] + x[8] * v[9] + x[10] * v[11]);
210 cheb24[1] = cheb12[1] + alam;
211 cheb24[23] = cheb12[1] - alam;
215 const double alam = (x[10] * v[1] - x[8] * v[3] + x[6] * v[5]
216 - x[4] * v[7] + x[2] * v[9] - x[0] * v[11]);
217 cheb24[11] = cheb12[11] + alam;
218 cheb24[13] = cheb12[11] - alam;
222 const double alam = (x[4] * v[1] - x[8] * v[3] - x[0] * v[5]
223 - x[10] * v[7] + x[2] * v[9] + x[6] * v[11]);
224 cheb24[5] = cheb12[5] + alam;
225 cheb24[19] = cheb12[5] - alam;
229 const double alam = (x[6] * v[1] - x[2] * v[3] - x[10] * v[5]
230 + x[0] * v[7] - x[8] * v[9] - x[4] * v[11]);
231 cheb24[7] = cheb12[7] + alam;
232 cheb24[17] = cheb12[7] - alam;
235 for (i = 0; i < 6; i++) {
236 const size_t j = 12 - i;
237 v[i] = fval[i] - fval[j];
238 fval[i] = fval[i] + fval[j];
242 const double alam1 = v[0] + x[7] * v[4];
243 const double alam2 = x[3] * v[2];
245 cheb12[2] = alam1 + alam2;
246 cheb12[10] = alam1 - alam2;
249 cheb12[6] = v[0] - v[4];
252 const double alam = x[1] * v[1] + x[5] * v[3] + x[9] * v[5];
253 cheb24[2] = cheb12[2] + alam;
254 cheb24[22] = cheb12[2] - alam;
258 const double alam = x[5] * (v[1] - v[3] - v[5]);
259 cheb24[6] = cheb12[6] + alam;
260 cheb24[18] = cheb12[6] - alam;
264 const double alam = x[9] * v[1] - x[5] * v[3] + x[1] * v[5];
265 cheb24[10] = cheb12[10] + alam;
266 cheb24[14] = cheb12[10] - alam;
269 for (i = 0; i < 3; i++) {
270 const size_t j = 6 - i;
271 v[i] = fval[i] - fval[j];
272 fval[i] = fval[i] + fval[j];
275 cheb12[4] = v[0] + x[7] * v[2];
276 cheb12[8] = fval[0] - x[7] * fval[2];
279 const double alam = x[3] * v[1];
280 cheb24[4] = cheb12[4] + alam;
281 cheb24[20] = cheb12[4] - alam;
285 const double alam = x[7] * fval[1] - fval[3];
286 cheb24[8] = cheb12[8] + alam;
287 cheb24[16] = cheb12[8] - alam;
290 cheb12[0] = fval[0] + fval[2];
293 const double alam = fval[1] + fval[3];
294 cheb24[0] = cheb12[0] + alam;
295 cheb24[24] = cheb12[0] - alam;
298 cheb12[12] = v[0] - v[2];
299 cheb24[12] = cheb12[12];
301 for (i = 1; i < 12; i++) {
302 cheb12[i] *= 1.0 / 6.0;
305 cheb12[0] *= 1.0 / 12.0;
306 cheb12[12] *= 1.0 / 12.0;
308 for (i = 1; i < 24; i++) {
309 cheb24[i] *= 1.0 / 12.0;
312 cheb24[0] *= 1.0 / 24.0;
313 cheb24[24] *= 1.0 / 24.0;
367 double &res,
double &err) {
372 #ifndef DOXYGEN_INTERNAL
378 int qawc(func_t &func,
const double a,
const double b,
const double c,
379 const double epsabs,
const double epsrel,
380 double *result,
double *abserr) {
383 double result0, abserr0;
385 size_t iteration = 0;
386 int roundoff_type1 = 0, roundoff_type2 = 0, error_type = 0;
389 double lower, higher;
396 size_t limit=this->
w->
limit;
409 double dbl_eps=std::numeric_limits<double>::epsilon();
411 if (epsabs <= 0 && (epsrel < 50 * dbl_eps || epsrel < 0.5e-28)) {
413 std::string estr=
"Tolerance cannot be achieved with given ";
414 estr+=
"value of tol_abs, "+
dtos(epsabs)+
", and tol_rel, "+
415 dtos(epsrel)+
", in inte_qawc_gsl::qawc().";
419 if (c == a || c == b) {
421 std::string estr=
"Cannot integrate with singularity on endpoint ";
422 estr+=
"in inte_qawc_gsl::qawc().";
428 this->
qc25c(func,lower,higher,c,&result0,&abserr0, &err_reliable);
435 tolerance = GSL_MAX_DBL (epsabs, epsrel * fabs (result0));
437 if (abserr0 < tolerance && abserr0 < 0.01 * fabs(result0)) {
440 *result = sign * result0;
444 }
else if (limit == 1) {
446 *result = sign * result0;
450 std::string estr=
"A maximum of 1 iteration was insufficient ";
451 estr+=
"in inte_qawc_gsl::qawc().";
462 double a1, b1, a2, b2;
463 double a_i, b_i, r_i, e_i;
464 double area1 = 0, area2 = 0, area12 = 0;
465 double error1 = 0, error2 = 0, error12 = 0;
466 int err_reliable1, err_reliable2;
470 this->
w->
retrieve (&a_i, &b_i, &r_i, &e_i);
473 b1 = 0.5 * (a_i + b_i);
477 if (c > a1 && c <= b1) {
478 b1 = 0.5 * (c + b2) ;
480 }
else if (c > b1 && c < b2) {
481 b1 = 0.5 * (a1 + c) ;
485 qc25c (func, a1, b1, c, &area1, &error1, &err_reliable1);
486 qc25c (func, a2, b2, c, &area2, &error2, &err_reliable2);
488 area12 = area1 + area2;
489 error12 = error1 + error2;
491 errsum += (error12 - e_i);
492 area += area12 - r_i;
494 if (err_reliable1 && err_reliable2) {
495 double delta = r_i - area12;
497 if (fabs (delta) <= 1.0e-5 * fabs (area12) &&
498 error12 >= 0.99 * e_i) {
501 if (iteration >= 10 && error12 > e_i) {
506 tolerance = GSL_MAX_DBL (epsabs, epsrel * fabs (area));
508 if (errsum > tolerance) {
509 if (roundoff_type1 >= 6 || roundoff_type2 >= 20) {
521 this->
w->
update (a1, b1, area1, error1, a2, b2, area2, error2);
523 this->
w->
retrieve (&a_i, &b_i, &r_i, &e_i);
526 std::cout <<
"inte_qawc_gsl Iter: " << iteration;
527 std::cout.setf(std::ios::showpos);
528 std::cout <<
" Res: " << area;
529 std::cout.unsetf(std::ios::showpos);
530 std::cout <<
" Err: " << errsum
531 <<
" Tol: " << tolerance << std::endl;
534 std::cout <<
"Press a key and type enter to continue. " ;
541 }
while (iteration < limit && !error_type && errsum > tolerance);
547 if (errsum <= tolerance) {
549 }
else if (error_type == 2) {
550 std::string estr=
"Roundoff error prevents tolerance ";
551 estr+=
"from being achieved in inte_qawc_gsl::qawc().";
553 }
else if (error_type == 3) {
554 std::string estr=
"Bad integrand behavior ";
555 estr+=
" in inte_qawc_gsl::qawc().";
557 }
else if (iteration == limit) {
558 std::string estr=
"Maximum number of subdivisions ("+
itos(iteration);
559 estr+=
") reached in inte_qawc_gsl::qawc().";
562 std::string estr=
"Could not integrate function in inte_qawc_gsl::";
563 estr+=
"qawc() (it may have returned a non-finite result).";
574 void qc25c(func_t &func,
double a,
double b,
double c,
575 double *result,
double *abserr,
int *err_reliable) {
577 double cc = (2 * c - b - a) / (b - a);
579 if (fabs (cc) > 1.1) {
580 double resabs, resasc;
584 if (*abserr == resasc) {
594 double cheb12[13], cheb24[25], moment[25];
595 double res12 = 0, res24 = 0;
600 for (i = 0; i < 13; i++) {
601 res12 += cheb12[i] * moment[i];
604 for (i = 0; i < 25; i++) {
605 res24 += cheb24[i] * moment[i];
609 *abserr = fabs(res24 - res12) ;
626 const char *
type() {
return "inte_qawc_gsl"; }
630 #ifndef DOXYGEN_NO_O2NS