inte_gauss_cern.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-2020, 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 /** \file inte_gauss_cern.h
24  \brief File defining \ref o2scl::inte_gauss_cern
25 */
26 #ifndef O2SCL_CERN_GAUSS_H
27 #define O2SCL_CERN_GAUSS_H
28 
29 #include <o2scl/misc.h>
30 #include <o2scl/inte.h>
31 
32 #ifdef O2SCL_LD_TYPES
33 #include <boost/multiprecision/cpp_dec_float.hpp>
34 #endif
35 
36 #ifndef DOXYGEN_NO_O2NS
37 namespace o2scl {
38 #endif
39 
40  /** \brief Integration weights and abcissas for
41  o2scl::inte_gauss_cern and \ref o2scl::inte_cauchy_cern in
42  double precision
43  */
45 
46  public:
47 
48  /** \brief Integration abscissas for \ref o2scl::inte_gauss_cern and
49  \ref o2scl::inte_cauchy_cern in double precision
50  */
51  double x[12];
52 
53  /** \brief Integration weights for \ref o2scl::inte_gauss_cern and
54  \ref o2scl::inte_cauchy_cern in double precision
55  */
56  double w[12];
57 
59  x[0]=0.96028985649753623;
60  x[1]=0.79666647741362674;
61  x[2]=0.52553240991632899;
62  x[3]=0.18343464249564980;
63  x[4]=0.98940093499164993;
64  x[5]=0.94457502307323258;
65  x[6]=0.86563120238783175;
66  x[7]=0.75540440835500303;
67  x[8]=0.61787624440264375;
68  x[9]=0.45801677765722739;
69  x[10]=0.28160355077925891;
70  x[11]=0.95012509837637440e-1;
71 
72  w[0]=0.10122853629037626;
73  w[1]=0.22238103445337447;
74  w[2]=0.31370664587788729;
75  w[3]=0.36268378337836198;
76  w[4]=0.27152459411754095e-1;
77  w[5]=0.62253523938647893e-1;
78  w[6]=0.95158511682492785e-1;
79  w[7]=0.12462897125553387;
80  w[8]=0.14959598881657673;
81  w[9]=0.16915651939500254;
82  w[10]=0.18260341504492359;
83  w[11]=0.18945061045506850;
84  }
85  };
86 
87  /** \brief Integration weights and abcissas for
88  o2scl::inte_gauss_cern and \ref o2scl::inte_cauchy_cern in
89  long double precision
90 
91  \warning The long double type doesn't work uniformly across
92  systems and so the accuracy when using this class varies.
93  */
95 
96  public:
97 
98  /** \brief Integration abscissas for \ref o2scl::inte_gauss_cern and
99  \ref o2scl::inte_cauchy_cern in long double precision
100  */
101  long double x[12];
102 
103  /** \brief Integration weights for \ref o2scl::inte_gauss_cern and
104  \ref o2scl::inte_cauchy_cern in long double precision
105  */
106  long double w[12];
107 
109 
110  x[0]=0.96028985649753623168356086856947299L;
111  x[1]=0.79666647741362673959155393647583044L;
112  x[2]=0.52553240991632898581773904918924635L;
113  x[3]=0.18343464249564980493947614236018398L;
114  x[4]=0.98940093499164993259615417345033263L;
115  x[5]=0.94457502307323257607798841553460835L;
116  x[6]=0.86563120238783174388046789771239313L;
117  x[7]=0.75540440835500303389510119484744227L;
118  x[8]=0.61787624440264374844667176404879102L;
119  x[9]=0.45801677765722738634241944298357757L;
120  x[10]=0.28160355077925891323046050146049611L;
121  x[11]=0.095012509837637440185319335424958063L;
122 
123  w[0]=0.10122853629037625915253135430996219L;
124  w[1]=0.22238103445337447054435599442624088L;
125  w[2]=0.31370664587788728733796220198660131L;
126  w[3]=0.36268378337836198296515044927719561L;
127  w[4]=0.027152459411754094851780572456018104L;
128  w[5]=0.062253523938647892862843836994377694L;
129  w[6]=0.095158511682492784809925107602246226L;
130  w[7]=0.12462897125553387205247628219201642L;
131  w[8]=0.14959598881657673208150173054747855L;
132  w[9]=0.16915651939500253818931207903035996L;
133  w[10]=0.18260341504492358886676366796921994L;
134  w[11]=0.18945061045506849628539672320828311L;
135  }
136 
137  };
138 
139 #if defined(O2SCL_LD_TYPES) || defined(DOXYGEN)
140 
141  /** \brief Integration weights and abcissas for
142  o2scl::inte_gauss_cern and \ref o2scl::inte_cauchy_cern for
143  the cpp_dec_float_50 type
144 
145  \note Experimental, and only included if
146  O2SCL_LD_TYPES is defined during the library configuration.
147 
148  \comment
149  Weights and abcissas originally generated using cpp_dec_float_100
150  numbers by AWS using code in ~/wcs/int5/sbox on 10/7/19.
151  \endcomment
152  */
154 
155  public:
156 
157  typedef boost::multiprecision::cpp_dec_float_50 cpp_dec_float_50;
158 
159  /** \brief Integration abscissas for \ref o2scl::inte_gauss_cern and
160  \ref o2scl::inte_cauchy_cern in cpp_dec_float_50 precision
161  */
162  cpp_dec_float_50 x[12];
163 
164  /** \brief Integration weights for \ref o2scl::inte_gauss_cern and
165  \ref o2scl::inte_cauchy_cern in cpp_dec_float_50 precision
166  */
167  cpp_dec_float_50 w[12];
168 
170 
171  // We convert from strings to ensure the full accuracy
172  x[0]=cpp_dec_float_50
173  ("9.60289856497536231683560868569472990428235234301452e-01");
174  x[1]=cpp_dec_float_50
175  ("7.96666477413626739591553936475830436837171731615965e-01");
176  x[2]=cpp_dec_float_50
177  ("5.25532409916328985817739049189246349041964243120393e-01");
178  x[3]=cpp_dec_float_50
179  ("1.83434642495649804939476142360183980666757812912974e-01");
180  x[4]=cpp_dec_float_50
181  ("9.89400934991649932596154173450332627426274071657645e-01");
182  x[5]=cpp_dec_float_50
183  ("9.44575023073232576077988415534608345091139272591073e-01");
184  x[6]=cpp_dec_float_50
185  ("8.65631202387831743880467897712393132387335384847527e-01");
186  x[7]=cpp_dec_float_50
187  ("7.55404408355003033895101194847442268353813656457503e-01");
188  x[8]=cpp_dec_float_50
189  ("6.17876244402643748446671764048791018991882217765658e-01");
190  x[9]=cpp_dec_float_50
191  ("4.58016777657227386342419442983577573540031613035523e-01");
192  x[10]=cpp_dec_float_50
193  ("2.81603550779258913230460501460496106486069490770600e-01");
194  x[11]=cpp_dec_float_50
195  ("9.50125098376374401853193354249580631303530556890655e-02");
196 
197  w[0]=cpp_dec_float_50
198  ("1.01228536290376259152531354309962190115394091051685e-01");
199  w[1]=cpp_dec_float_50
200  ("2.22381034453374470544355994426240884430130870051250e-01");
201  w[2]=cpp_dec_float_50
202  ("3.13706645877887287337962201986601313260328999002735e-01");
203  w[3]=cpp_dec_float_50
204  ("3.62683783378361982965150449277195612194146039894331e-01");
205  w[4]=cpp_dec_float_50
206  ("2.71524594117540948517805724560181035122673755667608e-02");
207  w[5]=cpp_dec_float_50
208  ("6.22535239386478928628438369943776942749865083529069e-02");
209  w[6]=cpp_dec_float_50
210  ("9.51585116824927848099251076022462263552635031837127e-02");
211  w[7]=cpp_dec_float_50
212  ("1.24628971255533872052476282192016420144886859222203e-01");
213  w[8]=cpp_dec_float_50
214  ("1.49595988816576732081501730547478548970491068207836e-01");
215  w[9]=cpp_dec_float_50
216  ("1.69156519395002538189312079030359962211639473416028e-01");
217  w[10]=cpp_dec_float_50
218  ("1.82603415044923588866763667969219939383556223654649e-01");
219  w[11]=cpp_dec_float_50
220  ("1.89450610455068496285396723208283105146908988395903e-01");
221  }
222 
223  };
224 
225 #endif
226 
227  /** \brief Gaussian quadrature (CERNLIB)
228 
229  For any interval \f$ (a,b) \f$, we define \f$ g_8(a,b) \f$ and
230  \f$ g_{16}(a,b) \f$ to be the 8- and 16-point Gaussian
231  quadrature approximations to
232  \f[
233  I=\int_a^b f(x)~dx
234  \f]
235  and define
236  \f[
237  r(a,b)=\frac{|g_{16}(a,b)-g_{8}(a,b)|}{1+g_{16}(a,b)}
238  \f]
239  The function integ() returns \f$ G \f$ given by
240  \f[
241  G=\sum_{i=1}^{k} g_{16}(x_{i-1},x_i)
242  \f]
243  where \f$ x_0=a \f$ and \f$ x_k=b \f$ and the subdivision
244  points \f$ x_i \f$ are given by
245  \f[
246  x_i=x_{i-1}+\lambda(B-x_{i-1})
247  \f]
248  where \f$ \lambda \f$ is the first number in the sequence \f$
249  1,\frac{1}{2},\frac{1}{4},... \f$ for which
250  \f[
251  r(x_{i-1},x_i)<\mathrm{eps}.
252  \f]
253  If, at any stage, the ratio
254  \f[
255  q=\left| \frac{x_i-x_{i-1}}{b-a} \right|
256  \f]
257  is so small so that \f$ 1+0.005 q \f$ is indistinguishable from
258  unity, then the accuracy is required is not reachable and
259  the error handler is called.
260 
261  Unless there is severe cancellation, inte::tol_rel may be
262  considered as specifying a bound on the relative error of the
263  integral in the case that \f$ |I|>1 \f$ and an absolute error if
264  \f$ |I|<1 \f$. More precisely, if \f$ k \f$ is the number of
265  subintervals from above, and if
266  \f[
267  I_{abs} = \int_a^b |f(x)|~dx
268  \f]
269  then
270  \f[
271  \frac{|G-I|}{I_{abs}+k}<\mathrm{tol_rel}
272  \f]
273  will nearly always be true when no error is returned. For
274  functions with no singualarities in the interval, the accuracy
275  will usually be higher than this.
276 
277  If the desired accuracy is not achieved, the integration
278  functions will call the error handler and return the best guess,
279  unless \ref inte::err_nonconv is false, in which case the error
280  handler is not called.
281 
282  This function is based on the CERNLIB routines GAUSS and
283  DGAUSS which are documented at
284  http://wwwasdoc.web.cern.ch/wwwasdoc/shortwrupsdir/d103/top.html
285 
286  \note Currently \o2 supports only types \c double and
287  \c long \c double for the floating point type \c fp_t .
288 
289  \future Allow user to change \c cst?
290  */
291  template<class func_t=funct, class fp_t=double,
292  class weights_t=inte_gauss_coeffs_double>
293  class inte_gauss_cern : public inte<func_t,fp_t> {
294 
295  public:
296 
297  const fp_t *w, *x;
298  weights_t wgts;
299 
300  inte_gauss_cern() {
301  w=&(wgts.w[0]);
302  x=&(wgts.x[0]);
303  }
304 
305  virtual ~inte_gauss_cern() {
306  }
307 
308  /** \brief Integrate function \c func from \c a to \c b.
309  */
310  virtual int integ_err(func_t &func, fp_t a, fp_t b,
311  fp_t &res, fp_t &err) {
312 
313  fp_t y1, y2;
314  err=0.0;
315 
316  size_t itx=0;
317 
318  int i;
319  bool loop=true, loop2=false;
320  static const fp_t cst=0.005;
321  fp_t h=0.0;
322  if (b==a) {
323  res=0.0;
324  return o2scl::success;
325  }
326  fp_t cnst=cst/(b-a);
327  fp_t aa=0.0, bb=a;
328  while (loop==true || loop2==true) {
329  itx++;
330  if (loop==true) {
331  aa=bb;
332  bb=b;
333  }
334  fp_t c1=(bb+aa)/2.0;
335  fp_t c2=(bb-aa)/2.0;
336  fp_t s8=0.0;
337  for(i=0;i<4;i++) {
338  fp_t u=c2*x[i];
339  y1=func(c1+u);
340  y2=func(c1-u);
341  s8+=w[i]*(y1+y2);
342  }
343  fp_t s16=0.0;
344  for(i=4;i<12;i++) {
345  fp_t u=c2*x[i];
346  y1=func(c1+u);
347  y2=func(c1-u);
348  s16+=w[i]*(y1+y2);
349  }
350  s16*=c2;
351 
352  loop=false;
353  loop2=false;
354 
355  if (o2scl::o2abs(s16-c2*s8)<this->tol_rel*
356  (1.0+o2scl::o2abs(s16))) {
357  h+=s16;
358  if (bb!=b) loop=true;
359  } else {
360  bb=c1;
361  fp_t one=1;
362  if (one+cnst*o2scl::o2abs(c2)!=one) {
363  loop2=true;
364  } else {
365  this->last_iter=itx;
366  O2SCL_CONV2_RET("Failed to reach required accuracy in cern_",
367  "gauss::integ().",exc_efailed,this->err_nonconv);
368  }
369  }
370  }
371  this->last_iter=itx;
372  res=h;
373  return o2scl::success;
374  }
375 
376  };
377 
378 #ifndef DOXYGEN_NO_O2NS
379 }
380 #endif
381 
382 #endif
o2scl::inte_gauss_coeffs_double::w
double w[12]
Integration weights for o2scl::inte_gauss_cern and o2scl::inte_cauchy_cern in double precision.
Definition: inte_gauss_cern.h:56
O2SCL_CONV2_RET
#define O2SCL_CONV2_RET(d, d2, n, b)
Set an error and return the error value, two-string version.
Definition: err_hnd.h:303
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_gauss_cern::integ_err
virtual int integ_err(func_t &func, fp_t a, fp_t b, fp_t &res, fp_t &err)
Integrate function func from a to b.
Definition: inte_gauss_cern.h:310
o2scl::inte_gauss_coeffs_double::x
double x[12]
Integration abscissas for o2scl::inte_gauss_cern and o2scl::inte_cauchy_cern in double precision.
Definition: inte_gauss_cern.h:51
o2scl::inte_gauss_coeffs_long_double::w
long double w[12]
Integration weights for o2scl::inte_gauss_cern and o2scl::inte_cauchy_cern in long double precision.
Definition: inte_gauss_cern.h:106
o2scl::success
@ success
Success.
Definition: err_hnd.h:47
o2scl::inte_gauss_coeffs_long_double::x
long double x[12]
Integration abscissas for o2scl::inte_gauss_cern and o2scl::inte_cauchy_cern in long double precision...
Definition: inte_gauss_cern.h:101
o2scl::inte_gauss_coeffs_cpp_dec_float_50::x
cpp_dec_float_50 x[12]
Integration abscissas for o2scl::inte_gauss_cern and o2scl::inte_cauchy_cern in cpp_dec_float_50 prec...
Definition: inte_gauss_cern.h:162
o2scl::inte< funct, double >::last_iter
size_t last_iter
The most recent number of iterations taken.
Definition: inte.h:70
o2scl::inte_gauss_coeffs_long_double
Integration weights and abcissas for o2scl::inte_gauss_cern and o2scl::inte_cauchy_cern in long doubl...
Definition: inte_gauss_cern.h:94
o2scl::inte< funct, double >::tol_rel
double tol_rel
The maximum relative uncertainty in the value of the integral (default )
Definition: inte.h:75
o2scl::inte_gauss_coeffs_cpp_dec_float_50
Integration weights and abcissas for o2scl::inte_gauss_cern and o2scl::inte_cauchy_cern for the cpp_d...
Definition: inte_gauss_cern.h:153
o2scl::inte
Base integration class [abstract base].
Definition: inte.h:51
o2scl::inte_gauss_cern
Gaussian quadrature (CERNLIB)
Definition: inte_gauss_cern.h:293
o2scl::funct
std::function< double(double)> funct
One-dimensional function typedef in src/base/funct.h.
Definition: funct.h:48
o2scl::o2abs
float o2abs(const float x)
Absolute value for single precision numbers.
o2scl::inte< funct, double >::err_nonconv
bool err_nonconv
If true, call the error handler if the routine does not converge or reach the desired tolerance (defa...
Definition: inte.h:88
o2scl::inte_gauss_coeffs_cpp_dec_float_50::w
cpp_dec_float_50 w[12]
Integration weights for o2scl::inte_gauss_cern and o2scl::inte_cauchy_cern in cpp_dec_float_50 precis...
Definition: inte_gauss_cern.h:167
o2scl::inte_gauss_coeffs_double
Integration weights and abcissas for o2scl::inte_gauss_cern and o2scl::inte_cauchy_cern in double pre...
Definition: inte_gauss_cern.h:44

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