41 #ifndef GSL_INTE_QAWS_H
42 #define GSL_INTE_QAWS_H
48 #include <o2scl/inte_qawc_gsl.h>
50 #ifndef DOXYGEN_NO_O2NS
83 #ifndef DOXYGEN_INTERNAL
109 const double alpha_p1 = this->alpha + 1.0;
110 const double beta_p1 = this->beta + 1.0;
112 const double alpha_p2 = this->alpha + 2.0;
113 const double beta_p2 = this->beta + 2.0;
115 const double r_alpha = pow (2.0, alpha_p1);
116 const double r_beta = pow (2.0,beta_p1);
122 this->ri[0] = r_alpha / alpha_p1;
123 this->ri[1] = this->ri[0] * this->alpha / alpha_p2;
128 for (i = 2; i < 25; i++) {
129 this->ri[i] = -(r_alpha + an * (an - alpha_p2) * this->ri[i - 1])
130 / (anm1 * (an + alpha_p1));
135 rj[0] = r_beta / beta_p1;
136 rj[1] = rj[0] * this->beta / beta_p2;
141 for (i = 2; i < 25; i++) {
142 rj[i] = (-(r_beta + an * (an - beta_p2) * rj[i - 1])
143 / (anm1 * (an + beta_p1)));
148 this->rg[0] = -this->ri[0] / alpha_p1;
149 this->rg[1] = -this->rg[0] - 2.0 * r_alpha / (alpha_p2 * alpha_p2);
154 for (i = 2; i < 25; i++) {
155 this->rg[i] = (-(an * (an - alpha_p2) *
156 this->rg[i - 1] - an * this->ri[i - 1]
157 + anm1 * this->ri[i])
158 / (anm1 * (an + alpha_p1)));
163 this->rh[0] = -this->rj[0] / beta_p1;
164 this->rh[1] = -this->rh[0] - 2.0 * r_beta / (beta_p2 * beta_p2);
169 for (i = 2; i < 25; i++) {
170 this->rh[i] = (-(an * (an - beta_p2) *
171 this->rh[i - 1] - an * this->rj[i - 1]
172 + anm1 * this->rj[i])
173 / (anm1 * (an + beta_p1)));
178 for (i = 1; i < 25; i += 2) {
206 double factor = 1.0,y;
226 return func(t)*factor;
233 void qc25s(func_t &func,
double a,
double b,
double a1,
double b1,
234 double &result,
double &abserr,
int &err_reliable) {
238 std::bind(std::mem_fn<
double(
double,func_t &)>
240 this,std::placeholders::_1,std::ref(func));
242 this->left_endpoint = a;
243 this->right_endpoint = b;
245 if (a1 == a && (this->alpha != 0.0 || this->mu != 0)) {
247 double cheb12[13], cheb24[25];
249 double factor = pow(0.5 * (b1 - a1),this->alpha + 1.0);
252 this->fn_qaws_R =
true;
253 this->fn_qaws_L =
false;
258 double res12 = 0,res24 = 0;
264 abserr = fabs(u * (res24 - res12));
268 double res12a = 0,res24a = 0;
269 double res12b = 0,res24b = 0;
271 double u = factor * log(b1 - a1);
277 result = u * res24a + v * res24b;
278 abserr = fabs(u*(res24a - res12a)) + fabs(v*(res24b - res12b));
284 }
else if (b1 == b && (this->beta != 0.0 || this->nu != 0)) {
286 double cheb12[13], cheb24[25];
287 double factor = pow(0.5 * (b1 - a1), this->beta + 1.0);
290 this->fn_qaws_L =
true;
291 this->fn_qaws_R =
false;
297 double res12 = 0, res24 = 0;
303 abserr = fabs(u * (res24 - res12));
307 double res12a = 0, res24a = 0;
308 double res12b = 0, res24b = 0;
310 double u = factor * log(b1 - a1);
316 result = u * res24a + v * res24b;
317 abserr = fabs(u*(res24a - res12a)) + fabs(v*(res24b - res12b));
325 double resabs, resasc;
328 this->fn_qaws_R =
true;
329 this->fn_qaws_L =
true;
331 this->
gauss_kronrod(func,a1,b1,&result,&abserr,&resabs,&resasc);
333 if (abserr == resasc) {
347 double &result12,
double &result24) {
353 for (i = 0; i < 13; i++) {
354 result12 += r[i] * cheb12[i];
357 for (i = 0; i < 25; i++) {
358 result24 += r[i] * cheb24[i];
393 int set_weight(
double u_alpha,
double u_beta,
int u_mu,
int u_nu) {
395 if (u_alpha < -1.0) {
396 std::string estr=((std::string)
"Variable alpha must be ")+
397 "greater than -1.0 in inte_qaws_gsl().";
401 std::string estr=((std::string)
"Variable beta must be ")+
402 "greater than -1.0 in inte_qaws_gsl().";
405 if (u_mu != 0 && u_mu != 1) {
406 std::string estr=((std::string)
"Variable mu must be 0 or 1 ")+
407 "in inte_qaws_gsl().";
410 if (u_nu != 0 && u_nu != 1) {
411 std::string estr=((std::string)
"Variable nu must be 0 or 1 ")+
412 "in inte_qaws_gsl().";
416 this->alpha = u_alpha;
428 void get_weight(
double &u_alpha,
double &u_beta,
int &u_mu,
int &u_nu) {
429 u_alpha = this->alpha;
439 double &result,
double &abserr) {
442 double result0, abserr0;
445 int roundoff_type1 = 0, roundoff_type2 = 0, error_type = 0;
454 size_t limit=this->
w->
limit;
457 std::string estr=
"Integration limits, a="+
dtos(a);
458 estr+=
" and b="+
dtos(b)+
", must satisfy a < b";
459 estr+=
" in inte_qaws_gsl::gsl_qaws().";
463 double dbl_eps=std::numeric_limits<double>::epsilon();
468 std::string estr=
"Tolerance cannot be achieved with given ";
469 estr+=
"value of tol_abs, "+
dtos(this->
tol_abs)+
", and tol_rel, "+
470 dtos(this->
tol_rel)+
", in inte_qaws_gsl::integ_err().";
478 double error1, error2;
479 int err_reliable1, err_reliable2;
481 double b1 = 0.5 * (a + b);
485 this->
qc25s(func, a, b, a1, b1, area1, error1, err_reliable1);
486 this->
qc25s(func, a, b, a2, b2, area2, error2, err_reliable2);
490 if (error1 > error2) {
498 result0 = area1 + area2;
499 abserr0 = error1 + error2;
504 tolerance = GSL_MAX_DBL (this->
tol_abs, this->
tol_rel * fabs (result0));
509 if (abserr0 < tolerance && abserr0 < 0.01 * fabs(result0)) {
513 }
else if (limit == 1) {
517 std::string estr =
"A maximum of 1 iteration was insufficient ";
518 estr +=
"in inte_qaws_gsl::gsl_qaws().";
526 double a1, b1, a2, b2;
527 double a_i, b_i, r_i, e_i;
528 double area1 = 0, area2 = 0, area12 = 0;
529 double error1 = 0, error2 = 0, error12 = 0;
530 int err_reliable1, err_reliable2;
536 b1 = 0.5 * (a_i + b_i);
540 qc25s(func, a, b, a1, b1, area1, error1, err_reliable1);
541 qc25s(func, a, b, a2, b2, area2, error2, err_reliable2);
543 area12 = area1 + area2;
544 error12 = error1 + error2;
546 errsum += (error12 - e_i);
547 area += area12 - r_i;
549 if (err_reliable1 && err_reliable2) {
551 double delta = r_i - area12;
553 if (fabs(delta) <= 1.0e-5 * fabs (area12)
554 && error12 >= 0.99 * e_i) {
557 if (this->
last_iter >= 10 && error12 > e_i) {
562 tolerance = GSL_MAX_DBL (this->
tol_abs, this->
tol_rel * fabs (area));
564 if (errsum > tolerance) {
565 if (roundoff_type1 >= 6 || roundoff_type2 >= 20) {
577 this->
w->
update(a1, b1, area1, error1, a2, b2, area2, error2);
581 }
while (this->last_iter < this->
w->
limit
582 && !error_type && errsum > tolerance);
588 if (errsum <= tolerance) {
590 }
else if (error_type == 2) {
591 std::string estr=
"Round-off error prevents tolerance ";
592 estr+=
"from being achieved in inte_qaws_gsl::gsl_qaws().";
594 }
else if (error_type == 3) {
595 std::string estr=
"Bad integrand behavior ";
596 estr+=
" in inte_qaws_gsl::gsl_qaws().";
599 std::string estr=
"Maximum number of subdivisions ("+
itos(limit);
600 estr+=
") reached in inte_qaws_gsl::gsl_qaws().";
603 std::string estr=
"Could not integrate function in ";
604 estr+=
"inte_qaws_gsl::gsl_qaws().";
615 const char *
type() {
return "inte_qaws_gsl"; }
619 #ifndef DOXYGEN_NO_O2NS