inte_qaws_gsl.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-2017, Jerry Gagelman and Andrew W. Steiner
5 
6  This file is part of O2scl.
7 
8  O2scl is free software; you can redistribute it and/or modify
9  it under the terms of the GNU General Public License as published by
10  the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version.
12 
13  O2scl is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with O2scl. If not, see <http://www.gnu.org/licenses/>.
20 
21  -------------------------------------------------------------------
22 */
23 /*
24  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Brian Gough
25  *
26  * This program is free software; you can redistribute it and/or modify
27  * it under the terms of the GNU General Public License as published by
28  * the Free Software Foundation; either version 3 of the License, or (at
29  * your option) any later version.
30  *
31  * This program is distributed in the hope that it will be useful, but
32  * WITHOUT ANY WARRANTY; without even the implied warranty of
33  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
34  * General Public License for more details.
35  *
36  * You should have received a copy of the GNU General Public License
37  * along with this program; if not, write to the Free Software
38  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
39  * 02110-1301, USA.
40  */
41 #ifndef GSL_INTE_QAWS_H
42 #define GSL_INTE_QAWS_H
43 
44 /** \file inte_qaws_gsl.h
45  \brief File defining \ref o2scl::inte_qaws_gsl
46 */
47 
48 #include <o2scl/inte_qawc_gsl.h>
49 
50 #ifndef DOXYGEN_NO_O2NS
51 namespace o2scl {
52 #endif
53 
54  /** \brief Adaptive integration with with algebraic-logarithmic
55  singularities at the end-points (GSL)
56 
57  This class computes the weighted integral
58  \f[
59  \int_a^b f(x)(x - a)^\alpha (b - x)^\beta \log^\mu(x - a)
60  \log^\nu(b - x)~dx
61  \f]
62  where the parameters of the weight function must satisfy
63  \f[
64  \alpha > -1, \quad \beta > -1, \quad
65  \mu \in \{0, 1\}, \quad \nu \in \{0, 1\},
66  \f]
67  and which are set by \ref set_weight(). Note that setting \f$
68  \mu=0 \f$ or \f$ \nu=0 \f$ removes the respective factor \f$
69  \log^mu(\ldots) \f$ or \f$ \log^\nu(\ldots) \f$ from the weight.
70 
71  The adaptive refinement algorithm described for \ref
72  inte_qag_gsl is used. When a subinterval contains one of the
73  endpoints, a special 25-point modified Clenshaw-Curtis rule is
74  used to control the singularities. For subintervals which do not
75  include the endpoints, a Gauss-Kronrod integration rule is used.
76 
77  See \ref gslinte_subsect in the User's guide for general
78  information about the GSL integration classes.
79  */
80  template<class func_t=funct11> class inte_qaws_gsl :
81  public inte_cheb_gsl<func_t> {
82 
83 #ifndef DOXYGEN_INTERNAL
84 
85  protected:
86 
87  /** \name Data from \c gsl_integration_qaws_table
88  */
89  //@{
90  double alpha;
91  double beta;
92  int mu;
93  int nu;
94 
95  double ri[25];
96  double rj[25];
97  double rg[25];
98  double rh[25];
99  //@}
100 
101  /** \brief Set the array values \c ri, \c rj, \c rg, \c rh from the
102  current values \c alpha and \c beta.
103 
104  This is the function from the GSL source code \c integration/qmomo.c
105  that initializes \c gsl_integration_qaws_table.
106  */
108 
109  const double alpha_p1 = this->alpha + 1.0;
110  const double beta_p1 = this->beta + 1.0;
111 
112  const double alpha_p2 = this->alpha + 2.0;
113  const double beta_p2 = this->beta + 2.0;
114 
115  const double r_alpha = pow (2.0, alpha_p1);
116  const double r_beta = pow (2.0,beta_p1);
117 
118  size_t i;
119 
120  double an, anm1;
121 
122  this->ri[0] = r_alpha / alpha_p1;
123  this->ri[1] = this->ri[0] * this->alpha / alpha_p2;
124 
125  an = 2.0;
126  anm1 = 1.0;
127 
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));
131  anm1 = an;
132  an = an + 1.0;
133  }
134 
135  rj[0] = r_beta / beta_p1;
136  rj[1] = rj[0] * this->beta / beta_p2;
137 
138  an = 2.0;
139  anm1 = 1.0;
140 
141  for (i = 2; i < 25; i++) {
142  rj[i] = (-(r_beta + an * (an - beta_p2) * rj[i - 1])
143  / (anm1 * (an + beta_p1)));
144  anm1 = an;
145  an = an + 1.0;
146  }
147 
148  this->rg[0] = -this->ri[0] / alpha_p1;
149  this->rg[1] = -this->rg[0] - 2.0 * r_alpha / (alpha_p2 * alpha_p2);
150 
151  an = 2.0;
152  anm1 = 1.0;
153 
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)));
159  anm1 = an;
160  an = an + 1.0;
161  }
162 
163  this->rh[0] = -this->rj[0] / beta_p1;
164  this->rh[1] = -this->rh[0] - 2.0 * r_beta / (beta_p2 * beta_p2);
165 
166  an = 2.0;
167  anm1 = 1.0;
168 
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)));
174  anm1 = an;
175  an = an + 1.0;
176  }
177 
178  for (i = 1; i < 25; i += 2) {
179  this->rj[i] *= -1;
180  this->rh[i] *= -1;
181  }
182 
183  return;
184  }
185 
186  /** \brief True if algebraic-logarithmic singularity is present at the
187  right endpoint in the definition \c f_trans.
188  */
189  bool fn_qaws_R;
190 
191  /** \brief True if algebraic-logarithmic singularity is present at the
192  left endpoint in the definition \c f_trans.
193  */
194  bool fn_qaws_L;
195 
196  /// Left endpoint in definition of \c f_trans
198 
199  /// Right endpoint in definition of \c f_trans.
201 
202  /** \brief Weighted integrand.
203  */
204  virtual double transform(double t, func_t &func) {
205 
206  double factor = 1.0,y;
207 
208  if (fn_qaws_L) {
209  if (alpha != 0.0) {
210  factor *= pow(t - left_endpoint,alpha);
211  }
212  if (mu == 1) {
213  factor *= log(t - left_endpoint);
214  }
215  }
216 
217  if (fn_qaws_R) {
218  if (beta != 0.0) {
219  factor *= pow(right_endpoint - t,beta);
220  }
221  if (nu == 1) {
222  factor *= log(right_endpoint - t);
223  }
224  }
225 
226  return func(t)*factor;
227  }
228 
229  /** \brief Clenshaw-Curtis 25-point integration and error estimator
230  for functions with an algebraic-logarithmic singularity at the
231  endpoint(s).
232  */
233  void qc25s(func_t &func, double a, double b, double a1, double b1,
234  double &result, double &abserr, int &err_reliable) {
235 
236  // Transformed function object for inte_cheb_series()
237  funct11 fmp=
238  std::bind(std::mem_fn<double(double,func_t &)>
240  this,std::placeholders::_1,func);
241 
242  this->left_endpoint = a;
243  this->right_endpoint = b;
244 
245  if (a1 == a && (this->alpha != 0.0 || this->mu != 0)) {
246 
247  double cheb12[13], cheb24[25];
248 
249  double factor = pow(0.5 * (b1 - a1),this->alpha + 1.0);
250 
251  // weighted_function.function = &fn_qaws_R;
252  this->fn_qaws_R = true;
253  this->fn_qaws_L = false;
254 
255  this->inte_cheb_series(fmp,a1,b1,cheb12,cheb24);
256 
257  if (this->mu == 0) {
258  double res12 = 0,res24 = 0;
259  double u = factor;
260 
261  this->compute_result(this->ri,cheb12,cheb24,res12,res24);
262 
263  result = u * res24;
264  abserr = fabs(u * (res24 - res12));
265 
266  } else {
267 
268  double res12a = 0,res24a = 0;
269  double res12b = 0,res24b = 0;
270 
271  double u = factor * log(b1 - a1);
272  double v = factor;
273 
274  this->compute_result(this->ri,cheb12,cheb24,res12a,res24a);
275  this->compute_result(this->rg,cheb12,cheb24,res12b,res24b);
276 
277  result = u * res24a + v * res24b;
278  abserr = fabs(u*(res24a - res12a)) + fabs(v*(res24b - res12b));
279  }
280 
281  err_reliable = 0;
282  return;
283 
284  } else if (b1 == b && (this->beta != 0.0 || this->nu != 0)) {
285 
286  double cheb12[13], cheb24[25];
287  double factor = pow(0.5 * (b1 - a1), this->beta + 1.0);
288 
289  // weighted_function.function = &fn_qaws_L;
290  this->fn_qaws_L = true;
291  this->fn_qaws_R = false;
292 
293  this->inte_cheb_series(fmp,a1,b1,cheb12,cheb24);
294 
295  if (this->nu == 0) {
296 
297  double res12 = 0, res24 = 0;
298  double u = factor;
299 
300  this->compute_result(this->rj, cheb12,cheb24,res12,res24);
301 
302  result = u * res24;
303  abserr = fabs(u * (res24 - res12));
304 
305  } else {
306 
307  double res12a = 0, res24a = 0;
308  double res12b = 0, res24b = 0;
309 
310  double u = factor * log(b1 - a1);
311  double v = factor;
312 
313  this->compute_result(this->rj,cheb12,cheb24,res12a,res24a);
314  this->compute_result(this->rh,cheb12,cheb24,res12b,res24b);
315 
316  result = u * res24a + v * res24b;
317  abserr = fabs(u*(res24a - res12a)) + fabs(v*(res24b - res12b));
318  }
319 
320  err_reliable = 0;
321  return;
322 
323  } else {
324 
325  double resabs, resasc;
326 
327  // weighted_function.function = &fn_qaws;
328  this->fn_qaws_R = true;
329  this->fn_qaws_L = true;
330 
331  this->gauss_kronrod(func,a1,b1,&result,&abserr,&resabs,&resasc);
332 
333  if (abserr == resasc) {
334  err_reliable = 0;
335  } else {
336  err_reliable = 1;
337  }
338 
339  return;
340  }
341  }
342 
343  /** \brief Compute the 13-point and 25-point approximations from
344  the Chebyshev moments and coefficients.
345  */
346  void compute_result(double *r, double *cheb12, double *cheb24,
347  double &result12, double &result24) {
348 
349  result12=0.0;
350  result24=0.0;
351 
352  size_t i;
353  for (i = 0; i < 13; i++) {
354  result12 += r[i] * cheb12[i];
355  }
356 
357  for (i = 0; i < 25; i++) {
358  result24 += r[i] * cheb24[i];
359  }
360  }
361 
362 #endif
363 
364  public:
365 
366  /** \brief Initialize the adptive workspace as with the constructor
367  \ref inte_qag_gsl::inte_qag_gsl.
368 
369  The default paramters \f$ \alpha, \beta, \mu, \nu \f$ of the weight
370  function are all zero.
371  */
372  inte_qaws_gsl() : inte_cheb_gsl<func_t>() {
373  set_weight(0.0,0.0,0,0);
374  }
375 
376  ~inte_qaws_gsl() {}
377 
378  /** \brief Sets the exponents of singularites of the weight function.
379 
380  The parameters determine the exponents of the weight function
381  \f[
382  W(x) = (x-a)^\alpha (b-x)^\beta \log^\mu(x-a) \log^\nu(b-x),
383  \f]
384  and must satsify
385  \f[
386  \alpha > -1, \quad \beta > -1, \quad
387  \mu \in \{0, 1\}, \quad \nu \in \{0, 1\}.
388  \f]
389  In order for the adaptive algorithm to run quickly, a table of
390  Chebyshev weights for the particular parameters are computed in
391  advance.
392  */
393  int set_weight(double u_alpha, double u_beta, int u_mu, int u_nu) {
394 
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().";
398  O2SCL_ERR(estr.c_str(),exc_einval);
399  }
400  if (u_beta < -1.0) {
401  std::string estr=((std::string)"Variable beta must be ")+
402  "greater than -1.0 in inte_qaws_gsl().";
403  O2SCL_ERR(estr.c_str(),exc_einval);
404  }
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().";
408  O2SCL_ERR(estr.c_str(),exc_einval);
409  }
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().";
413  O2SCL_ERR(estr.c_str(),exc_einval);
414  }
415 
416  this->alpha = u_alpha;
417  this->beta = u_beta;
418  this->mu = u_mu;
419  this->nu = u_nu;
420 
422  return success;
423  }
424 
425  /** \brief Returns the current values (via reference) of the
426  weight-function's parameters.
427  */
428  void get_weight(double &u_alpha, double &u_beta, int &u_mu, int &u_nu) {
429  u_alpha = this->alpha;
430  u_beta = this->beta;
431  u_mu = this->mu;
432  u_nu = this->nu;
433  }
434 
435  /** \brief Integrate the function \c func on the interval (\c a, \c b)
436  returning the \c result and error estimate \c abserr.
437  */
438  virtual int integ_err(func_t &func, double a, double b,
439  double &result, double &abserr) {
440 
441  double area, errsum;
442  double result0, abserr0;
443  double tolerance;
444  this->last_iter = 0;
445  int roundoff_type1 = 0, roundoff_type2 = 0, error_type = 0;
446 
447  /* Initialize results */
448 
449  this->w->initialise(a,b);
450 
451  result = 0;
452  abserr = 0;
453 
454  size_t limit=this->w->limit;
455 
456  if (b <= a) {
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().";
460  O2SCL_ERR(estr.c_str(),exc_einval);
461  }
462 
463  double dbl_eps=std::numeric_limits<double>::epsilon();
464 
465  if (this->tol_abs <= 0 && (this->tol_rel < 50 * dbl_eps ||
466  this->tol_rel < 0.5e-28)) {
467  this->last_iter=0;
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().";
471  O2SCL_ERR(estr.c_str(),exc_ebadtol);
472  }
473 
474  /* perform the first integration */
475 
476  {
477  double area1, area2;
478  double error1, error2;
479  int err_reliable1, err_reliable2;
480  double a1 = a;
481  double b1 = 0.5 * (a + b);
482  double a2 = b1;
483  double b2 = b;
484 
485  this->qc25s(func, a, b, a1, b1, area1, error1, err_reliable1);
486  this->qc25s(func, a, b, a2, b2, area2, error2, err_reliable2);
487 
488  this->last_iter = 2;
489 
490  if (error1 > error2) {
491  this->w->append_interval(a1, b1, area1, error1);
492  this->w->append_interval(a2, b2, area2, error2);
493  } else {
494  this->w->append_interval(a2, b2, area2, error2);
495  this->w->append_interval(a1, b1, area1, error1);
496  }
497 
498  result0 = area1 + area2;
499  abserr0 = error1 + error2;
500  }
501 
502  /* Test on accuracy */
503 
504  tolerance = GSL_MAX_DBL (this->tol_abs, this->tol_rel * fabs (result0));
505 
506  /* Test on accuracy, use 0.01 relative error as an extra safety
507  margin on the first iteration (ignored for subsequent iterations) */
508 
509  if (abserr0 < tolerance && abserr0 < 0.01 * fabs(result0)) {
510  result = result0;
511  abserr = abserr0;
512  return success;
513  } else if (limit == 1) {
514  result = result0;
515  abserr = abserr0;
516 
517  std::string estr = "A maximum of 1 iteration was insufficient ";
518  estr += "in inte_qaws_gsl::gsl_qaws().";
519  O2SCL_CONV_RET(estr.c_str(), exc_emaxiter, this->err_nonconv);
520  }
521 
522  area = result0;
523  errsum = abserr0;
524 
525  do {
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;
531 
532  /* Bisect the subinterval with the largest error estimate */
533  this->w->retrieve(&a_i,&b_i,&r_i,&e_i);
534 
535  a1 = a_i;
536  b1 = 0.5 * (a_i + b_i);
537  a2 = b1;
538  b2 = b_i;
539 
540  qc25s(func, a, b, a1, b1, area1, error1, err_reliable1);
541  qc25s(func, a, b, a2, b2, area2, error2, err_reliable2);
542 
543  area12 = area1 + area2;
544  error12 = error1 + error2;
545 
546  errsum += (error12 - e_i);
547  area += area12 - r_i;
548 
549  if (err_reliable1 && err_reliable2) {
550 
551  double delta = r_i - area12;
552 
553  if (fabs(delta) <= 1.0e-5 * fabs (area12)
554  && error12 >= 0.99 * e_i) {
555  roundoff_type1++;
556  }
557  if (this->last_iter >= 10 && error12 > e_i) {
558  roundoff_type2++;
559  }
560  }
561 
562  tolerance = GSL_MAX_DBL (this->tol_abs, this->tol_rel * fabs (area));
563 
564  if (errsum > tolerance) {
565  if (roundoff_type1 >= 6 || roundoff_type2 >= 20) {
566  // round off error
567  error_type = 2;
568  }
569 
570  /* set error flag in the case of bad integrand behaviour at
571  a point of the integration range */
572  if (this->w->subinterval_too_small (a1, a2, b2)) {
573  error_type = 3;
574  }
575  }
576 
577  this->w->update(a1, b1, area1, error1, a2, b2, area2, error2);
578  this->w->retrieve(&a_i,&b_i,&r_i,&e_i);
579  this->last_iter++;
580 
581  } while (this->last_iter < this->w->limit
582  && !error_type && errsum > tolerance);
583 
584  result = this->w->sum_results();
585  abserr = errsum;
586  this->interror = abserr;
587 
588  if (errsum <= tolerance) {
589  return success;
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().";
593  O2SCL_CONV_RET(estr.c_str(),exc_eround,this->err_nonconv);
594  } else if (error_type == 3) {
595  std::string estr="Bad integrand behavior ";
596  estr+=" in inte_qaws_gsl::gsl_qaws().";
597  O2SCL_CONV_RET(estr.c_str(),exc_esing,this->err_nonconv);
598  } else if (this->last_iter == limit) {
599  std::string estr="Maximum number of subdivisions ("+itos(limit);
600  estr+=") reached in inte_qaws_gsl::gsl_qaws().";
601  O2SCL_CONV_RET(estr.c_str(),exc_emaxiter,this->err_nonconv);
602  } else {
603  std::string estr="Could not integrate function in ";
604  estr+="inte_qaws_gsl::gsl_qaws().";
605  O2SCL_ERR(estr.c_str(),exc_efailed);
606  }
607 
608  // No return statement needed since the above if statement
609  // always forces a return, but some compilers like having one
610  // anyway.
611  return o2scl::success;
612  }
613 
614  /// Return string denoting type ("inte_qaws_gsl")
615  const char *type() { return "inte_qaws_gsl"; }
616 
617  };
618 
619 #ifndef DOXYGEN_NO_O2NS
620 }
621 #endif
622 
623 #endif
virtual int integ_err(func_t &func, double a, double b, double &result, double &abserr)
Integrate the function func on the interval (a, b) returning the result and error estimate abserr...
bool fn_qaws_L
True if algebraic-logarithmic singularity is present at the left endpoint in the definition f_trans...
The main O<span style=&#39;position: relative; top: 0.3em; font-size: 0.8em&#39;>2</span>scl O$_2$scl names...
Definition: anneal.h:42
bool fn_qaws_R
True if algebraic-logarithmic singularity is present at the right endpoint in the definition f_trans...
void get_weight(double &u_alpha, double &u_beta, int &u_mu, int &u_nu)
Returns the current values (via reference) of the weight-function&#39;s parameters.
#define O2SCL_CONV_RET(d, n, b)
Set a "convergence" error and return the error value.
Definition: err_hnd.h:292
std::function< double(double)> funct11
One-dimensional function typedef.
Definition: funct.h:44
double left_endpoint
Left endpoint in definition of f_trans.
void inte_cheb_series(func2_t &f, double a, double b, double *cheb12, double *cheb24)
Compute Chebyshev series expansion using a FFT method.
invalid argument supplied by user
Definition: err_hnd.h:59
apparent singularity detected
Definition: err_hnd.h:93
exceeded max number of iterations
Definition: err_hnd.h:73
Adaptive integration with with algebraic-logarithmic singularities at the end-points (GSL) ...
Definition: inte_qaws_gsl.h:80
size_t last_iter
The most recent number of iterations taken.
Definition: inte.h:63
bool err_nonconv
If true, call the error handler if the routine does not converge or reach the desired tolerance (defa...
Definition: inte.h:81
double right_endpoint
Right endpoint in definition of f_trans.
int update(double a1, double b1, double area1, double error1, double a2, double b2, double area2, double error2)
Determine which new subinterval to add to the workspace stack and perform update. ...
generic failure
Definition: err_hnd.h:61
Chebyshev integration base class (GSL)
Definition: inte_qawc_gsl.h:46
int subinterval_too_small(double a1, double a2, double b2)
Test whether the proposed subdivision falls before floating-point precision.
double tol_abs
The maximum absolute uncertainty in the value of the integral (default )
Definition: inte.h:73
virtual void gauss_kronrod(func_t &func, double a, double b, double *result, double *abserr, double *resabs, double *resasc)
Integration wrapper for internal transformed function type.
int set_weight(double u_alpha, double u_beta, int u_mu, int u_nu)
Sets the exponents of singularites of the weight function.
int initialise(double a, double b)
Initialize the workspace for an integration with limits a and b.
std::string dtos(double x, int prec=6, bool auto_prec=false)
Convert a double to a string.
inte_qaws_gsl()
Initialize the adptive workspace as with the constructor inte_qag_gsl::inte_qag_gsl.
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
virtual double transform(double t, func_t &func)
Weighted integrand.
void compute_result(double *r, double *cheb12, double *cheb24, double &result12, double &result24)
Compute the 13-point and 25-point approximations from the Chebyshev moments and coefficients.
double interror
The uncertainty for the last integration computation.
Definition: inte.h:117
Integrate a function with a singularity (GSL) [abstract base].
void append_interval(double a1, double b1, double area1, double error1)
Push a new interval to the workspace stack.
inte_workspace_gsl * w
The integration workspace.
double sum_results()
Add up all of the contributions to construct the final result.
void qc25s(func_t &func, double a, double b, double a1, double b1, double &result, double &abserr, int &err_reliable)
Clenshaw-Curtis 25-point integration and error estimator for functions with an algebraic-logarithmic ...
double tol_rel
The maximum relative uncertainty in the value of the integral (default )
Definition: inte.h:68
size_t limit
Maximum number of subintervals allocated.
int retrieve(double *a, double *b, double *r, double *e) const
Retrieve the ith result from the workspace stack.
const char * type()
Return string denoting type ("inte_qaws_gsl")
std::string itos(int x)
Convert an integer to a string.
user specified an invalid tolerance
Definition: err_hnd.h:77
Success.
Definition: err_hnd.h:47
failed because of roundoff error
Definition: err_hnd.h:87
void initialise_qaws_table()
Set the array values ri, rj, rg, rh from the current values alpha and beta.

Documentation generated with Doxygen. Provided under the GNU Free Documentation License (see License Information).