23 #ifndef O2SCL_GSL_INTE_QAWO_H
24 #define O2SCL_GSL_INTE_QAWO_H
31 #include <o2scl/inte.h>
32 #include <o2scl/inte_qawc_gsl.h>
34 #ifndef DOXYGEN_NO_O2NS
79 virtual int integ_err(func_t &func,
double a,
double b,
80 double &res,
double &err) {
82 otable=gsl_integration_qawo_table_alloc
86 this->
w,this->
otable,&res,&err);
88 gsl_integration_qawo_table_free(
otable);
93 #ifndef DOXYGEN_INTERNAL
104 int qawo(func_t &func,
const double a,
const double epsabs,
106 gsl_integration_qawo_table *wf,
double *result,
double *abserr) {
109 double res_ext, err_ext;
110 double result0=0.0, abserr0=0.0, resabs0=0.0, resasc0=0.0;
114 double error_over_large_intervals = 0;
115 double reseps = 0, abseps = 0, correc = 0;
117 int roundoff_type1 = 0, roundoff_type2 = 0, roundoff_type3 = 0;
118 int error_type = 0, error_type2 = 0;
120 size_t iteration = 0;
122 int positive_integrand = 0;
125 int disallow_extrapolation = 0;
129 double b = a + wf->L ;
130 double abs_omega = fabs (wf->omega) ;
139 size_t limit=this->
w->
limit;
143 double dbl_eps=std::numeric_limits<double>::epsilon();
145 if (epsabs <= 0 && (epsrel < 50 * dbl_eps || epsrel < 0.5e-28)) {
147 std::string estr=
"Tolerance cannot be achieved with given ";
148 estr+=
"value of tol_abs, "+
dtos(epsabs)+
", and tol_rel, "+
149 dtos(epsrel)+
", in inte_qawo_gsl_sin::qawo().";
155 this->
qc25f(func, a, b, wf, 0,
156 &result0, &abserr0, &resabs0, &resasc0);
160 tolerance = GSL_MAX_DBL (epsabs, epsrel * fabs (result0));
163 std::cout <<
"inte_qawo_gsl Iter: " << 1;
164 std::cout.setf(std::ios::showpos);
165 std::cout <<
" Res: " << result0;
166 std::cout.unsetf(std::ios::showpos);
167 std::cout <<
" Err: " << abserr0
168 <<
" Tol: " << tolerance << std::endl;
171 std::cout <<
"Press a key and type enter to continue. " ;
176 if (abserr0 <= 100 * GSL_DBL_EPSILON * resabs0 &&
177 abserr0 > tolerance) {
182 std::string estr=
"Cannot reach tolerance because of roundoff error ";
183 estr+=
"on first attempt in inte_qawo_gsl_sin::qawo().";
185 }
else if ((abserr0 <= tolerance && abserr0 != resasc0) ||
192 }
else if (limit == 1) {
197 std::string estr=
"A maximum of 1 iteration was insufficient ";
198 estr+=
"in inte_qawo_gsl_sin::qawo().";
206 if (0.5 * abs_omega * fabs(b - a) <= 2) {
215 err_ext = GSL_DBL_MAX;
223 size_t current_level;
224 double a1, b1, a2, b2;
225 double a_i, b_i, r_i, e_i;
226 double area1 = 0, area2 = 0, area12 = 0;
227 double error1 = 0, error2 = 0, error12 = 0;
228 double resasc1=0.0, resasc2=0.0;
229 double resabs1=0.0, resabs2=0.0;
234 loc_w->
retrieve (&a_i, &b_i, &r_i, &e_i);
236 current_level = loc_w->
level[loc_w->
i] + 1;
238 if (current_level >= wf->n) {
244 b1 = 0.5 * (a_i + b_i);
250 qc25f(func, a1, b1, wf, current_level,
251 &area1, &error1, &resabs1, &resasc1);
252 qc25f(func, a2, b2, wf, current_level,
253 &area2, &error2, &resabs2, &resasc2);
255 area12 = area1 + area2;
256 error12 = error1 + error2;
266 errsum = errsum + error12 - e_i;
267 area = area + area12 - r_i;
269 tolerance = GSL_MAX_DBL (epsabs, epsrel * fabs (area));
272 std::cout <<
"inte_qawo_gsl Iter: " << iteration;
273 std::cout.setf(std::ios::showpos);
274 std::cout <<
" Res: " << area;
275 std::cout.unsetf(std::ios::showpos);
276 std::cout <<
" Err: " << errsum
277 <<
" Tol: " << tolerance << std::endl;
280 std::cout <<
"Press a key and type enter to continue. " ;
285 if (resasc1 != error1 && resasc2 != error2) {
287 double delta = r_i - area12;
289 if (fabs (delta) <= 1.0e-5 * fabs (area12) &&
290 error12 >= 0.99 * e_i) {
299 if (iteration > 10 && error12 > e_i) {
306 if (roundoff_type1 + roundoff_type2 >= 10 || roundoff_type3 >= 20) {
310 if (roundoff_type2 >= 5) {
323 loc_w->
update (a1, b1, area1, error1, a2, b2, area2, error2);
325 if (errsum <= tolerance) {
333 if (iteration >= limit - 1) {
340 if (iteration == 2 && extall) {
341 error_over_large_intervals = errsum;
347 if (disallow_extrapolation) {
352 error_over_large_intervals += -last_e_i;
354 if (current_level < loc_w->maximum_level)
356 error_over_large_intervals += error12;
375 double width = loc_w->
blist[i] - loc_w->
alist[i];
377 if (0.25 * fabs(width) * abs_omega > 2) {
382 error_over_large_intervals = errsum;
389 if (!error_type2 && error_over_large_intervals > ertest) {
402 error_over_large_intervals = errsum;
406 this->
qelg (&table, &reseps, &abseps);
410 if (ktmin > 5 && err_ext < 0.001 * errsum) {
414 if (abseps < err_ext) {
418 correc = error_over_large_intervals;
419 ertest = GSL_MAX_DBL (epsabs, epsrel * fabs (reseps));
420 if (err_ext <= ertest) {
428 disallow_extrapolation = 1;
431 if (error_type == 5) {
439 error_over_large_intervals = errsum;
441 }
while (iteration < limit);
446 if (err_ext == GSL_DBL_MAX)
449 if (error_type || error_type2) {
458 if (result != 0 && area != 0) {
459 if (err_ext / fabs (res_ext) > errsum / fabs (area)) {
462 }
else if (err_ext > errsum) {
464 }
else if (area == 0.0) {
472 double max_area = GSL_MAX_DBL (fabs (res_ext), fabs (area));
474 if (!positive_integrand && max_area < 0.01 * resabs0) {
480 double ratio = res_ext / area;
482 if (ratio < 0.01 || ratio > 100 || errsum > fabs (area)) {
496 if (error_type > 2) {
502 if (error_type == 0) {
504 }
else if (error_type == 1) {
505 std::string estr=
"Number of iterations was insufficient ";
506 estr+=
" in inte_qawo_gsl_sin::qawo().";
508 }
else if (error_type == 2) {
509 std::string estr=
"Roundoff error prevents tolerance ";
510 estr+=
"from being achieved in inte_qawo_gsl_sin::qawo().";
512 }
else if (error_type == 3) {
513 std::string estr=
"Bad integrand behavior ";
514 estr+=
" in inte_qawo_gsl_sin::qawo().";
516 }
else if (error_type == 4) {
517 std::string estr=
"Roundoff error detected in extrapolation table ";
518 estr+=
"in inte_qawo_gsl_sin::qawo().";
520 }
else if (error_type == 5) {
521 std::string estr=
"Integral is divergent or slowly convergent ";
522 estr+=
"in inte_qawo_gsl_sin::qawo().";
524 }
else if (error_type == -1) {
525 std::string estr=
"Exceeded limit of trigonometric table ";
526 estr+=
"inte_qawo_gsl_sin::qawo()";
529 std::string estr=
"Could not integrate function in inte_qawo_gsl";
530 estr+=
"::qawo() (it may have returned a non-finite result).";
541 void qc25f(func_t &func,
double a,
double b,
542 gsl_integration_qawo_table *wf,
size_t level,
543 double *result,
double *abserr,
double *resabs,
546 const double center = 0.5 * (a + b);
547 const double half_length = 0.5 * (b - a);
549 const double par =
omega * half_length;
551 if (fabs (par) < 2) {
560 double cheb12[13], cheb24[25];
561 double result_abs, res12_cos, res12_sin, res24_cos, res24_sin;
562 double est_cos, est_sin;
568 if (level >= wf->n) {
570 O2SCL_ERR(
"Table overflow in inte_qawo_gsl::qc25f().",
577 moment = wf->chebmo + 25 * level;
579 res12_cos = cheb12[12] * moment[12];
582 for (i = 0; i < 6; i++) {
583 size_t k = 10 - 2 * i;
584 res12_cos += cheb12[k] * moment[k];
585 res12_sin += cheb12[k + 1] * moment[k + 1];
588 res24_cos = cheb24[24] * moment[24];
591 result_abs = fabs(cheb24[24]) ;
593 for (i = 0; i < 12; i++) {
594 size_t k = 22 - 2 * i;
595 res24_cos += cheb24[k] * moment[k];
596 res24_sin += cheb24[k + 1] * moment[k + 1];
597 result_abs += fabs(cheb24[k]) + fabs(cheb24[k+1]);
600 est_cos = fabs(res24_cos - res12_cos);
601 est_sin = fabs(res24_sin - res12_sin);
603 c = half_length * cos(center *
omega);
604 s = half_length * sin(center *
omega);
606 if (wf->sine == GSL_INTEG_SINE) {
607 *result = c * res24_sin + s * res24_cos;
608 *abserr = fabs(c * est_sin) + fabs(s * est_cos);
610 *result = c * res24_cos - s * res24_sin;
611 *abserr = fabs(c * est_cos) + fabs(s * est_sin);
614 *resabs = result_abs * half_length;
615 *resasc = GSL_DBL_MAX;
623 return func(t)*sin(this->omega*t);
629 const char *
type() {
return "inte_qawo_gsl_sin"; }
659 double &res,
double &err) {
661 this->
otable=gsl_integration_qawo_table_alloc
665 this->
w,this->
otable,&res,&err);
667 gsl_integration_qawo_table_free(this->
otable);
672 #ifndef DOXYGEN_INTERNAL
678 return func(t)*cos(this->
omega*t);
684 const char *
type() {
return "inte_qawo_gsl_cos"; }
688 #ifndef DOXYGEN_NO_O2NS