inte_qawc_gsl.h
Go to the documentation of this file.
1  /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-2020, 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 #ifndef O2SCL_GSL_INTE_QAWC_H
24 #define O2SCL_GSL_INTE_QAWC_H
25 
26 /** \file inte_qawc_gsl.h
27  \brief File defining \ref o2scl::inte_qawc_gsl
28 */
29 
30 #include <o2scl/err_hnd.h>
31 #include <o2scl/inte.h>
32 #include <o2scl/inte_singular_gsl.h>
33 
34 #ifndef DOXYGEN_NO_O2NS
35 namespace o2scl {
36 #endif
37 
38  /** \brief Chebyshev integration base class (GSL)
39 
40  This class provides the basic Chebyshev integration functions
41  for use in the GSL-based integration classes which
42  require them. See \ref gslinte_subsect in the User's
43  guide for general information about the GSL integration classes.
44  */
45  template<class func_t> class inte_cheb_gsl :
46  public inte_transform_gsl<func_t> {
47 
48  protected:
49 
50  /// Compute the Chebyshev moments
51  void compute_moments(double cc, double *moment) {
52  size_t k;
53 
54  double a0 = log (fabs ((1.0 - cc) / (1.0 + cc)));
55  double a1 = 2 + a0 * cc;
56 
57  moment[0] = a0;
58  moment[1] = a1;
59 
60  for (k = 2; k < 25; k++) {
61  double a2;
62 
63  if ((k % 2) == 0) {
64  a2 = 2.0 * cc * a1 - a0;
65  } else {
66  const double km1 = k - 1.0;
67  a2 = 2.0 * cc * a1 - a0 - 4.0 / (km1 * km1 - 1.0);
68  }
69 
70  moment[k] = a2;
71 
72  a0 = a1;
73  a1 = a2;
74  }
75  }
76 
77  /** \brief Compute Chebyshev series expansion using a FFT method
78 
79  The Chebyshev coefficients for the truncated expansions,
80  \f[
81  f(x) =
82  \frac{a_0}{2}T_0(x) + \frac{a_d}{2}T_d(x) +
83  \sum_{k=}^{d-1} a_k^{(d)}T_k(x),
84  \f]
85  are computed for \f$ d=12 \f$ and \f$ d=24 \f$ using an FFT
86  algorithm from \ref Tolstov62 that is adapted so that the both
87  sets of coefficients are computed simultaneously.
88 
89  Given the function specified in \c f, this function computes
90  the 13 Chebyshev coefficients, \f$ C^{12}_{k} \f$ of degree 12
91  and 25 Chebyshev coefficients of degree 24, \f$ C^{24}_{k}
92  \f$, for the interval \f$ [a,b] \f$ using a FFT method.
93 
94  These coefficients are constructed to approximate
95  the original function with
96  \f[
97  f = \sum_{k=1}^{13} C^{12}_{k} T_{k-1}(x)
98  \f]
99  and
100  \f[
101  f = \sum_{k=1}^{25} C^{24}_{k} T_{k-1}(x)
102  \f]
103  where \f$ T_{k-1}(x) \f$ is the Chebyshev polynomial of
104  degree \f$ k-1 \f$ evaluated at the point \f$ x \f$.
105 
106  It is assumed that memory for \c cheb12 and \c cheb24 has
107  been allocated beforehand.
108 
109  Originally written in QUADPACK by R. Piessens and E. de
110  Doncker, translated into C for GSL by Brian Gough, and then
111  rewritten for \o2.
112  */
113  template<class func2_t>
114  void inte_cheb_series(func2_t &f, double a, double b,
115  double *cheb12, double *cheb24) {
116  size_t i;
117  double fval[25], v[12];
118 
119  /* These are the values of cos(pi*k/24) for k=1..11 needed for the
120  Chebyshev expansion of f(x) */
121 
122  const double x[11] = { 0.9914448613738104,
123  0.9659258262890683,
124  0.9238795325112868,
125  0.8660254037844386,
126  0.7933533402912352,
127  0.7071067811865475,
128  0.6087614290087206,
129  0.5000000000000000,
130  0.3826834323650898,
131  0.2588190451025208,
132  0.1305261922200516 };
133 
134  const double center = 0.5 * (b + a);
135  const double half_length = 0.5 * (b - a);
136 
137  double y1, y2, y3;
138  y1=f(b);
139  y2=f(center);
140  y3=f(a);
141  fval[0] = 0.5 * y1;
142  fval[12] = y2;
143  fval[24] = 0.5 * y3;
144 
145  for (i = 1; i < 12; i++) {
146  const size_t j = 24 - i;
147  const double u = half_length * x[i-1];
148  double yp, ym;
149  yp=f(center+u);
150  ym=f(center-u);
151  fval[i] = yp;
152  fval[j] = ym;
153  }
154 
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];
159  }
160 
161  {
162  const double alam1 = v[0] - v[8];
163  const double alam2 = x[5] * (v[2] - v[6] - v[10]);
164 
165  cheb12[3] = alam1 + alam2;
166  cheb12[9] = alam1 - alam2;
167  }
168 
169  {
170  const double alam1 = v[1] - v[7] - v[9];
171  const double alam2 = v[3] - v[5] - v[11];
172  {
173  const double alam = x[2] * alam1 + x[8] * alam2;
174 
175  cheb24[3] = cheb12[3] + alam;
176  cheb24[21] = cheb12[3] - alam;
177  }
178 
179  {
180  const double alam = x[8] * alam1 - x[2] * alam2;
181  cheb24[9] = cheb12[9] + alam;
182  cheb24[15] = cheb12[9] - alam;
183  }
184  }
185 
186  {
187  const double part1 = x[3] * v[4];
188  const double part2 = x[7] * v[8];
189  const double part3 = x[5] * v[6];
190 
191  {
192  const double alam1 = v[0] + part1 + part2;
193  const double alam2 = x[1] * v[2] + part3 + x[9] * v[10];
194 
195  cheb12[1] = alam1 + alam2;
196  cheb12[11] = alam1 - alam2;
197  }
198 
199  {
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;
204  }
205  }
206 
207  {
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;
212  }
213 
214  {
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;
219  }
220 
221  {
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;
226  }
227 
228  {
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;
233  }
234 
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];
239  }
240 
241  {
242  const double alam1 = v[0] + x[7] * v[4];
243  const double alam2 = x[3] * v[2];
244 
245  cheb12[2] = alam1 + alam2;
246  cheb12[10] = alam1 - alam2;
247  }
248 
249  cheb12[6] = v[0] - v[4];
250 
251  {
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;
255  }
256 
257  {
258  const double alam = x[5] * (v[1] - v[3] - v[5]);
259  cheb24[6] = cheb12[6] + alam;
260  cheb24[18] = cheb12[6] - alam;
261  }
262 
263  {
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;
267  }
268 
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];
273  }
274 
275  cheb12[4] = v[0] + x[7] * v[2];
276  cheb12[8] = fval[0] - x[7] * fval[2];
277 
278  {
279  const double alam = x[3] * v[1];
280  cheb24[4] = cheb12[4] + alam;
281  cheb24[20] = cheb12[4] - alam;
282  }
283 
284  {
285  const double alam = x[7] * fval[1] - fval[3];
286  cheb24[8] = cheb12[8] + alam;
287  cheb24[16] = cheb12[8] - alam;
288  }
289 
290  cheb12[0] = fval[0] + fval[2];
291 
292  {
293  const double alam = fval[1] + fval[3];
294  cheb24[0] = cheb12[0] + alam;
295  cheb24[24] = cheb12[0] - alam;
296  }
297 
298  cheb12[12] = v[0] - v[2];
299  cheb24[12] = cheb12[12];
300 
301  for (i = 1; i < 12; i++) {
302  cheb12[i] *= 1.0 / 6.0;
303  }
304 
305  cheb12[0] *= 1.0 / 12.0;
306  cheb12[12] *= 1.0 / 12.0;
307 
308  for (i = 1; i < 24; i++) {
309  cheb24[i] *= 1.0 / 12.0;
310  }
311 
312  cheb24[0] *= 1.0 / 24.0;
313  cheb24[24] *= 1.0 / 24.0;
314  }
315 
316  };
317 
318  /** \brief Adaptive Cauchy principal value integration (GSL)
319 
320  The Cauchy principal value of the integral of
321  \f[
322  \int_a^b \frac{f(x)}{x-c}~dx =
323  \lim_{\epsilon\to 0^+}
324  \left\{ \int_a^{c-\epsilon} \frac{f(x)}{x-c}~dx +
325  \int_{c+\epsilon}^b \frac{f(x)}{x-c}~dx \right\}.
326  \f]
327  over \f$ (a,b), \f$ with a singularity at \f$ c, \f$ is
328  computed. The adaptive refinement algorithm described for
329  inte_qag_gsl is used with modifications to ensure that
330  subdivisions do not occur at the singular point \f$ x = c\f$ .
331  When a subinterval contains the point \f$ x = c \f$ or is close
332  to it, a special 25-point modified Clenshaw-Curtis rule is used
333  to control the singularity. Further away from the singularity
334  the algorithm uses a Gauss-Kronrod integration rule.
335 
336  The location of the singularity must be specified before-hand in
337  inte_qawc_gsl::s, and the singularity must not be at one of the
338  endpoints. Note that when integrating a function of the form \f$
339  \frac{f(x)}{(x-s)} \f$, the denominator \f$ (x-s) \f$ must not
340  be specified in the argument \c func to integ(). Note that this
341  is different from how the \ref inte_cauchy_cern operates.
342 
343  See \ref gslinte_subsect in the User's guide for general
344  information about the GSL integration classes.
345 
346  \future Make inte_cauchy_cern and this class consistent in the
347  way which they require the user to provide the denominator
348  in the integrand
349  */
350  template<class func_t> class inte_qawc_gsl :
351  public inte_cheb_gsl<func_t> {
352 
353  public:
354 
355  inte_qawc_gsl() {
356  }
357 
358  virtual ~inte_qawc_gsl() {}
359 
360  /// The singularity
361  double s;
362 
363  /** \brief Integrate function \c func from \c a to \c b and place
364  the result in \c res and the error in \c err
365  */
366  virtual int integ_err(func_t &func, double a, double b,
367  double &res, double &err) {
368 
369  return this->qawc(func,a,b,s,this->tol_abs,this->tol_rel,&res,&err);
370  }
371 
372 #ifndef DOXYGEN_INTERNAL
373 
374  protected:
375 
376  /** \brief The full GSL integration routine called by integ_err()
377  */
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) {
381 
382  double area, errsum;
383  double result0, abserr0;
384  double tolerance;
385  size_t iteration = 0;
386  int roundoff_type1 = 0, roundoff_type2 = 0, error_type = 0;
387  int err_reliable;
388  int sign = 1;
389  double lower, higher;
390 
391  /* Initialize results */
392 
393  *result = 0;
394  *abserr = 0;
395 
396  size_t limit=this->w->limit;
397 
398  if (b < a) {
399  lower = b;
400  higher = a;
401  sign = -1;
402  } else {
403  lower = a;
404  higher = b;
405  }
406 
407  this->w->initialise(lower,higher);
408 
409  double dbl_eps=std::numeric_limits<double>::epsilon();
410 
411  if (epsabs <= 0 && (epsrel < 50 * dbl_eps || epsrel < 0.5e-28)) {
412  this->last_iter=0;
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().";
416  O2SCL_ERR(estr.c_str(),exc_ebadtol);
417  }
418 
419  if (c == a || c == b) {
420  this->last_iter=0;
421  std::string estr="Cannot integrate with singularity on endpoint ";
422  estr+="in inte_qawc_gsl::qawc().";
423  O2SCL_ERR(estr.c_str(),exc_einval);
424  }
425 
426  /* perform the first integration */
427 
428  this->qc25c(func,lower,higher,c,&result0,&abserr0, &err_reliable);
429 
430  this->w->set_initial_result (result0, abserr0);
431 
432  /* Test on accuracy, use 0.01 relative error as an extra safety
433  margin on the first iteration (ignored for subsequent iterations)
434  */
435  tolerance = GSL_MAX_DBL (epsabs, epsrel * fabs (result0));
436 
437  if (abserr0 < tolerance && abserr0 < 0.01 * fabs(result0)) {
438 
439  this->last_iter=1;
440  *result = sign * result0;
441  *abserr = abserr0;
442  return success;
443 
444  } else if (limit == 1) {
445 
446  *result = sign * result0;
447  *abserr = abserr0;
448 
449  this->last_iter=1;
450  std::string estr="A maximum of 1 iteration was insufficient ";
451  estr+="in inte_qawc_gsl::qawc().";
452  O2SCL_CONV_RET(estr.c_str(),exc_emaxiter,this->err_nonconv);
453  }
454 
455  area = result0;
456  errsum = abserr0;
457 
458  iteration = 1;
459 
460  do {
461 
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;
467 
468  /* Bisect the subinterval with the largest error estimate */
469 
470  this->w->retrieve (&a_i, &b_i, &r_i, &e_i);
471 
472  a1 = a_i;
473  b1 = 0.5 * (a_i + b_i);
474  a2 = b1;
475  b2 = b_i;
476 
477  if (c > a1 && c <= b1) {
478  b1 = 0.5 * (c + b2) ;
479  a2 = b1;
480  } else if (c > b1 && c < b2) {
481  b1 = 0.5 * (a1 + c) ;
482  a2 = b1;
483  }
484 
485  qc25c (func, a1, b1, c, &area1, &error1, &err_reliable1);
486  qc25c (func, a2, b2, c, &area2, &error2, &err_reliable2);
487 
488  area12 = area1 + area2;
489  error12 = error1 + error2;
490 
491  errsum += (error12 - e_i);
492  area += area12 - r_i;
493 
494  if (err_reliable1 && err_reliable2) {
495  double delta = r_i - area12;
496 
497  if (fabs (delta) <= 1.0e-5 * fabs (area12) &&
498  error12 >= 0.99 * e_i) {
499  roundoff_type1++;
500  }
501  if (iteration >= 10 && error12 > e_i) {
502  roundoff_type2++;
503  }
504  }
505 
506  tolerance = GSL_MAX_DBL (epsabs, epsrel * fabs (area));
507 
508  if (errsum > tolerance) {
509  if (roundoff_type1 >= 6 || roundoff_type2 >= 20) {
510  error_type = 2; /* round off error */
511  }
512 
513  /* set error flag in the case of bad integrand behaviour at
514  a point of the integration range */
515 
516  if (this->w->subinterval_too_small (a1, a2, b2)) {
517  error_type = 3;
518  }
519  }
520 
521  this->w->update (a1, b1, area1, error1, a2, b2, area2, error2);
522 
523  this->w->retrieve (&a_i, &b_i, &r_i, &e_i);
524 
525  if (this->verbose>0) {
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;
532  if (this->verbose>1) {
533  char ch;
534  std::cout << "Press a key and type enter to continue. " ;
535  std::cin >> ch;
536  }
537  }
538 
539  iteration++;
540 
541  } while (iteration < limit && !error_type && errsum > tolerance);
542 
543  *result = sign * this->w->sum_results();
544  *abserr = errsum;
545 
546  this->last_iter=iteration;
547  if (errsum <= tolerance) {
548  return success;
549  } else if (error_type == 2) {
550  std::string estr="Roundoff error prevents tolerance ";
551  estr+="from being achieved in inte_qawc_gsl::qawc().";
552  O2SCL_CONV_RET(estr.c_str(),exc_eround,this->err_nonconv);
553  } else if (error_type == 3) {
554  std::string estr="Bad integrand behavior ";
555  estr+=" in inte_qawc_gsl::qawc().";
556  O2SCL_CONV_RET(estr.c_str(),exc_esing,this->err_nonconv);
557  } else if (iteration == limit) {
558  std::string estr="Maximum number of subdivisions ("+itos(iteration);
559  estr+=") reached in inte_qawc_gsl::qawc().";
560  O2SCL_CONV_RET(estr.c_str(),exc_emaxiter,this->err_nonconv);
561  } else {
562  std::string estr="Could not integrate function in inte_qawc_gsl::";
563  estr+="qawc() (it may have returned a non-finite result).";
564  O2SCL_ERR(estr.c_str(),exc_efailed);
565  }
566 
567  // No return statement needed since the above if statement
568  // always forces a return, but some compilers like having one
569  // anyway.
570  return o2scl::success;
571  }
572 
573  /// 25-point quadrature for Cauchy principal values
574  void qc25c(func_t &func, double a, double b, double c,
575  double *result, double *abserr, int *err_reliable) {
576 
577  double cc = (2 * c - b - a) / (b - a);
578 
579  if (fabs (cc) > 1.1) {
580  double resabs, resasc;
581 
582  this->gauss_kronrod(func,a,b,result,abserr,&resabs,&resasc);
583 
584  if (*abserr == resasc) {
585  *err_reliable = 0;
586  } else {
587  *err_reliable = 1;
588  }
589 
590  return;
591 
592  } else {
593 
594  double cheb12[13], cheb24[25], moment[25];
595  double res12 = 0, res24 = 0;
596  size_t i;
597  this->inte_cheb_series(func, a, b, cheb12, cheb24);
598  this->compute_moments (cc, moment);
599 
600  for (i = 0; i < 13; i++) {
601  res12 += cheb12[i] * moment[i];
602  }
603 
604  for (i = 0; i < 25; i++) {
605  res24 += cheb24[i] * moment[i];
606  }
607 
608  *result = res24;
609  *abserr = fabs(res24 - res12) ;
610  *err_reliable = 0;
611 
612  return;
613  }
614  }
615 
616  /// Add the singularity to the function
617  virtual double transform(double t, func_t &func) {
618  double y;
619  y=func(t);
620  return y/(t-s);
621  }
622 
623 #endif
624 
625  /// Return string denoting type ("inte_qawc_gsl")
626  const char *type() { return "inte_qawc_gsl"; }
627 
628  };
629 
630 #ifndef DOXYGEN_NO_O2NS
631 }
632 #endif
633 
634 #endif
o2scl::exc_esing
@ exc_esing
apparent singularity detected
Definition: err_hnd.h:93
o2scl::inte_cheb_gsl::inte_cheb_series
void inte_cheb_series(func2_t &f, double a, double b, double *cheb12, double *cheb24)
Compute Chebyshev series expansion using a FFT method.
Definition: inte_qawc_gsl.h:114
o2scl::dtos
std::string dtos(const fp_t &x, int prec=6, bool auto_prec=false)
Convert a double to a string.
Definition: string_conv.h:77
o2scl::inte_kronrod_gsl< funct >::w
inte_workspace_gsl * w
The integration workspace.
Definition: inte_kronrod_gsl.h:620
o2scl::inte_workspace_gsl::retrieve
int retrieve(double *a, double *b, double *r, double *e) const
Retrieve the ith result from the workspace stack.
o2scl::inte_qawc_gsl
Adaptive Cauchy principal value integration (GSL)
Definition: inte_qawc_gsl.h:350
o2scl::inte_workspace_gsl::set_initial_result
int set_initial_result(double result, double error)
Update the workspace with the result and error from the first integration.
o2scl::exc_efailed
@ exc_efailed
generic failure
Definition: err_hnd.h:61
o2scl
The main O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl names...
Definition: anneal.h:42
o2scl::inte_qawc_gsl::qawc
int qawc(func_t &func, const double a, const double b, const double c, const double epsabs, const double epsrel, double *result, double *abserr)
The full GSL integration routine called by integ_err()
Definition: inte_qawc_gsl.h:378
o2scl::exc_emaxiter
@ exc_emaxiter
exceeded max number of iterations
Definition: err_hnd.h:73
o2scl::inte_workspace_gsl::initialise
int initialise(double a, double b)
Initialize the workspace for an integration with limits a and b.
o2scl::inte_qawc_gsl::type
const char * type()
Return string denoting type ("inte_qawc_gsl")
Definition: inte_qawc_gsl.h:626
o2scl::inte_transform_gsl::gauss_kronrod
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.
Definition: inte_singular_gsl.h:761
o2scl::inte_cheb_gsl::compute_moments
void compute_moments(double cc, double *moment)
Compute the Chebyshev moments.
Definition: inte_qawc_gsl.h:51
o2scl::success
@ success
Success.
Definition: err_hnd.h:47
o2scl::inte_workspace_gsl::sum_results
double sum_results()
Add up all of the contributions to construct the final result.
o2scl::exc_einval
@ exc_einval
invalid argument supplied by user
Definition: err_hnd.h:59
o2scl::inte_cheb_gsl
Chebyshev integration base class (GSL)
Definition: inte_qawc_gsl.h:45
o2scl::inte< funct >::last_iter
size_t last_iter
The most recent number of iterations taken.
Definition: inte.h:70
O2SCL_CONV_RET
#define O2SCL_CONV_RET(d, n, b)
Set a "convergence" error and return the error value.
Definition: err_hnd.h:297
o2scl::inte_qawc_gsl::qc25c
void qc25c(func_t &func, double a, double b, double c, double *result, double *abserr, int *err_reliable)
25-point quadrature for Cauchy principal values
Definition: inte_qawc_gsl.h:574
o2scl::inte_qawc_gsl::s
double s
The singularity.
Definition: inte_qawc_gsl.h:361
o2scl::inte_workspace_gsl::update
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.
o2scl::inte_transform_gsl
Integrate a function with a singularity (GSL) [abstract base].
Definition: inte_singular_gsl.h:749
o2scl::inte< funct >::tol_abs
double tol_abs
The maximum absolute uncertainty in the value of the integral (default )
Definition: inte.h:80
o2scl::inte< funct >::verbose
int verbose
Verbosity.
Definition: inte.h:67
o2scl::itos
std::string itos(int x)
Convert an integer to a string.
o2scl::inte< funct >::tol_rel
double tol_rel
The maximum relative uncertainty in the value of the integral (default )
Definition: inte.h:75
o2scl::exc_eround
@ exc_eround
failed because of roundoff error
Definition: err_hnd.h:87
O2SCL_ERR
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
o2scl::inte_workspace_gsl::limit
size_t limit
Maximum number of subintervals allocated.
Definition: inte_kronrod_gsl.h:507
o2scl::inte_qawc_gsl::integ_err
virtual int integ_err(func_t &func, double a, double b, double &res, double &err)
Integrate function func from a to b and place the result in res and the error in err.
Definition: inte_qawc_gsl.h:366
o2scl::exc_ebadtol
@ exc_ebadtol
user specified an invalid tolerance
Definition: err_hnd.h:77
o2scl::inte_workspace_gsl::subinterval_too_small
int subinterval_too_small(double a1, double a2, double b2)
Test whether the proposed subdivision falls before floating-point precision.
o2scl::inte_qawc_gsl::transform
virtual double transform(double t, func_t &func)
Add the singularity to the function.
Definition: inte_qawc_gsl.h:617

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