23 #ifndef O2SCL_GSL_INTE_SINGULAR_H
24 #define O2SCL_GSL_INTE_SINGULAR_H
32 #include <gsl/gsl_integration.h>
34 #include <o2scl/string_conv.h>
35 #include <o2scl/inte.h>
36 #include <o2scl/inte_kronrod_gsl.h>
38 #ifndef DOXYGEN_NO_O2NS
94 double dbl_eps=std::numeric_limits<double>::epsilon();
95 int status=(fabs(result) >= (1-50*dbl_eps)*resabs);
140 void qelg(
struct extrapolation_table *
table,
double *result,
143 double *epstab =
table->rlist2;
144 double *res3la =
table->res3la;
145 const size_t n =
table->n - 1;
147 const double current = epstab[n];
149 double absolute = GSL_DBL_MAX;
150 double relative = 5 * GSL_DBL_EPSILON * fabs (current);
152 const size_t newelm = n / 2;
153 const size_t n_orig = n;
157 const size_t nres_orig =
table->nres;
160 *abserr = GSL_DBL_MAX;
164 *abserr = GSL_MAX_DBL (absolute, relative);
168 epstab[n + 2] = epstab[n];
169 epstab[n] = GSL_DBL_MAX;
171 for (i = 0; i < newelm; i++) {
172 double res = epstab[n - 2 * i + 2];
173 double e0 = epstab[n - 2 * i - 2];
174 double e1 = epstab[n - 2 * i - 1];
177 double e1abs = fabs (e1);
178 double delta2 = e2 - e1;
179 double err2 = fabs (delta2);
180 double tol2 = GSL_MAX_DBL (fabs (e2), e1abs) * GSL_DBL_EPSILON;
181 double delta3 = e1 - e0;
182 double err3 = fabs (delta3);
183 double tol3 = GSL_MAX_DBL (e1abs, fabs (e0)) * GSL_DBL_EPSILON;
185 double e3, delta1, err1, tol1, ss;
187 if (err2 <= tol2 && err3 <= tol3) {
192 absolute = err2 + err3;
193 relative = 5 * GSL_DBL_EPSILON * fabs (res);
194 *abserr = GSL_MAX_DBL (absolute, relative);
198 e3 = epstab[n - 2 * i];
199 epstab[n - 2 * i] = e1;
201 err1 = fabs (delta1);
202 tol1 = GSL_MAX_DBL (e1abs, fabs (e3)) * GSL_DBL_EPSILON;
207 if (err1 <= tol1 || err2 <= tol2 || err3 <= tol3) {
212 ss = (1 / delta1 + 1 / delta2) - 1 / delta3;
217 if (fabs (ss * e1) <= 0.0001) {
225 epstab[n - 2 * i] = res;
228 const double error = err2 + fabs (res - e2) + err3;
230 if (error <= *abserr) {
240 const size_t limexp = 50 - 1;
242 if (n_final == limexp) {
243 n_final = 2 * (limexp / 2);
247 if (n_orig % 2 == 1) {
248 for (i = 0; i <= newelm; i++) {
249 epstab[1 + i * 2] = epstab[i * 2 + 3];
252 for (i = 0; i <= newelm; i++) {
253 epstab[i * 2] = epstab[i * 2 + 2];
256 if (n_orig != n_final) {
257 for (i = 0; i <= n_final; i++) {
258 epstab[i] = epstab[n_orig - n_final + i];
262 table->n = n_final + 1;
266 res3la[nres_orig] = *result;
267 *abserr = GSL_DBL_MAX;
271 *abserr = (fabs (*result - res3la[2]) + fabs (*result - res3la[1])
272 + fabs (*result - res3la[0]));
274 res3la[0] = res3la[1];
275 res3la[1] = res3la[2];
285 table->nres = nres_orig + 1;
287 *abserr = GSL_MAX_DBL (*abserr, 5 * GSL_DBL_EPSILON * fabs (*result));
294 size_t i = workspace->
i ;
295 const size_t * level = workspace->
level;
306 workspace->
nrmax = 0;
307 workspace->
i = workspace->
order[0] ;
313 int id = workspace->
nrmax;
316 const size_t * level = workspace->
level;
317 const size_t * order = workspace->
order;
319 size_t limit = workspace->
limit ;
320 size_t last = workspace->
size - 1 ;
322 if (last > (1 + limit / 2)) {
323 jupbnd = limit + 1 - last;
328 for (k =
id; k <= jupbnd; k++) {
329 size_t i_max = order[workspace->
nrmax];
331 workspace->
i = i_max ;
348 int qags(func_t &func,
const double a,
const double b,
349 const double l_epsabs,
const double l_epsrel,
350 double *result,
double *abserr) {
353 double res_ext, err_ext;
354 double result0, abserr0, resabs0, resasc0;
358 double error_over_large_intervals = 0;
359 double reseps = 0, abseps = 0, correc = 0;
361 int roundoff_type1 = 0, roundoff_type2 = 0, roundoff_type3 = 0;
362 int error_type = 0, error_type2 = 0;
364 size_t iteration = 0;
366 int positive_integrand = 0;
368 int disallow_extrapolation = 0;
370 struct extrapolation_table
table;
379 size_t limit=this->
w->
limit;
383 if (this->
tol_abs <= 0 && (this->
tol_rel < 50 * GSL_DBL_EPSILON ||
387 std::string estr=
"Tolerance cannot be achieved with given ";
388 estr+=
"value of tol_abs, "+
dtos(l_epsabs)+
", and tol_rel, "+
389 dtos(l_epsrel)+
", in inte_singular_gsl::qags().";
395 this->
gauss_kronrod(func,a,b,&result0,&abserr0,&resabs0,&resasc0);
399 tolerance = GSL_MAX_DBL (this->
tol_abs, this->
tol_rel * fabs (result0));
401 if (abserr0 <= 100 * GSL_DBL_EPSILON * resabs0 &&
402 abserr0 > tolerance) {
409 std::string estr=
"Cannot reach tolerance because of roundoff error ";
410 estr+=
"on first attempt in inte_singular_gsl::qags().";
413 }
else if ((abserr0 <= tolerance &&
414 abserr0 != resasc0) || abserr0 == 0.0) {
421 }
else if (limit == 1) {
428 "in inte_singular_gsl::qags().",
441 err_ext = GSL_DBL_MAX;
451 std::cout << this->
type();
452 std::cout <<
" Iter: " << iteration;
453 std::cout.setf(std::ios::showpos);
454 std::cout <<
" Res: " << area;
455 std::cout.unsetf(std::ios::showpos);
456 std::cout <<
" Err: " << errsum
457 <<
" Tol: " << tolerance << std::endl;
460 std::cout <<
"Press a key and type enter to continue. " ;
465 size_t current_level;
466 double a1, b1, a2, b2;
467 double a_i, b_i, r_i, e_i;
468 double area1 = 0, area2 = 0, area12 = 0;
469 double error1 = 0, error2 = 0, error12 = 0;
470 double resasc1, resasc2;
471 double resabs1, resabs2;
476 this->
w->
retrieve (&a_i, &b_i, &r_i, &e_i);
478 current_level = this->
w->
level[this->
w->
i] + 1;
481 b1 = 0.5 * (a_i + b_i);
487 this->
gauss_kronrod(func,a1,b1,&area1,&error1,&resabs1,&resasc1);
488 this->
gauss_kronrod(func,a2,b2,&area2,&error2,&resabs2,&resasc2);
490 area12 = area1 + area2;
491 error12 = error1 + error2;
502 errsum = errsum + error12 - e_i;
503 area = area + area12 - r_i;
505 tolerance = GSL_MAX_DBL (this->
tol_abs, this->
tol_rel * fabs (area));
506 if (resasc1 != error1 && resasc2 != error2) {
507 double delta = r_i - area12;
509 if (fabs (delta) <= 1.0e-5 * fabs (area12) &&
510 error12 >= 0.99 * e_i) {
517 if (iteration > 10 && error12 > e_i) {
524 if (roundoff_type1 + roundoff_type2 >= 10 ||
525 roundoff_type3 >= 20) {
530 if (roundoff_type2 >= 5) {
543 this->
w->
update(a1,b1,area1,error1,a2,b2,area2,error2);
545 if (errsum <= tolerance) {
549 std::cout << this->
type();
550 std::cout <<
" Iter: " << iteration;
551 std::cout.setf(std::ios::showpos);
552 std::cout <<
" Res: " << area;
553 std::cout.unsetf(std::ios::showpos);
554 std::cout <<
" Err: " << errsum
555 <<
" Tol: " << tolerance << std::endl;
558 std::cout <<
"Press a key and type enter to continue. " ;
570 if (iteration >= limit - 1) {
575 if (iteration == 2) {
576 error_over_large_intervals = errsum;
582 if (disallow_extrapolation) {
586 error_over_large_intervals += -last_e_i;
589 error_over_large_intervals += error12;
605 if (!error_type2 && error_over_large_intervals > ertest) {
619 if (ktmin > 5 && err_ext < 0.001 * errsum) {
623 if (abseps < err_ext) {
627 correc = error_over_large_intervals;
628 ertest = GSL_MAX_DBL (this->
tol_abs,
629 this->
tol_rel * fabs (reseps));
630 if (err_ext <= ertest) {
638 disallow_extrapolation = 1;
641 if (error_type == 5) {
649 error_over_large_intervals = errsum;
651 }
while (iteration < limit);
656 if (err_ext == GSL_DBL_MAX)
659 if (error_type || error_type2) {
667 if (res_ext != 0.0 && area != 0.0) {
668 if (err_ext / fabs (res_ext) > errsum / fabs (area))
670 }
else if (err_ext > errsum) {
672 }
else if (area == 0.0) {
680 double max_area = GSL_MAX_DBL (fabs (res_ext),fabs (area));
682 if (!positive_integrand && max_area < 0.01 * resabs0)
687 double ratio = res_ext / area;
689 if (ratio < 0.01 || ratio > 100.0 || errsum > fabs (area))
702 if (error_type > 2) {
708 if (error_type == 0) {
710 }
else if (error_type == 1) {
711 std::string estr=
"Maximum number of subdivisions ("+
itos(iteration);
712 estr+=
") reached in inte_singular_gsl::qags().";
714 }
else if (error_type == 2) {
715 std::string estr=
"Roundoff error prevents tolerance ";
716 estr+=
"from being achieved in inte_singular_gsl::qags().";
718 }
else if (error_type == 3) {
719 std::string estr=
"Bad integrand behavior ";
720 estr+=
"in inte_singular_gsl::qags().";
722 }
else if (error_type == 4) {
723 std::string estr=
"Roundoff error detected in extrapolation table ";
724 estr+=
"in inte_singular_gsl::qags().";
726 }
else if (error_type == 5) {
727 std::string estr=
"Integral is divergent or slowly convergent ";
728 estr+=
"in inte_singular_gsl::qags().";
732 std::string estr=
"Could not integrate function in inte_kronrod_gsl";
733 estr+=
"::qags() (it may have returned a non-finite result).";
755 virtual double transform(
double t, func_t &func)=0;
761 (func_t &func,
double a,
double b,
762 double *result,
double *abserr,
double *resabs,
double *resasc) {
764 funct fmp=std::bind(std::mem_fn<
double(
double,func_t &)>
766 this,std::placeholders::_1,std::ref(func));
769 (fmp,a,b,result,abserr,resabs,resasc);
774 #ifndef DOXYGEN_NO_O2NS