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

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