26 #ifndef O2SCL_CERN_ADAPT_H
27 #define O2SCL_CERN_ADAPT_H
29 #include <o2scl/misc.h>
30 #include <o2scl/inte.h>
31 #include <o2scl/inte_gauss56_cern.h>
32 #include <o2scl/string_conv.h>
33 #include <o2scl/vector_derint.h>
35 #ifndef DOXYGEN_NO_O2NS
39 #ifdef O2SCL_NEVER_DEFINED
48 template<
class func_t=funct,
class fp_t=
double>
49 class inte_newton_cotes :
public inte<func_t,fp_t> {
56 virtual int integ_err(func_t &func, fp_t a, fp_t b,
57 fp_t &res, fp_t &err) {
59 std::vector<fp_t> x(8), y(8);
63 for(
size_t i=0;i<8;i++) {
67 res=o2scl::vector_integ_extended8<std::vector<fp_t>,
69 fp_t res2=o2scl::vector_integ_extended4<std::vector<fp_t>,
82 template<
class func_t=funct,
class fp_t=
double>
83 class inte_newton_cotes2 :
public inte<func_t,fp_t> {
90 virtual int integ_err(func_t &func, fp_t a, fp_t b,
91 fp_t &res, fp_t &err) {
93 std::vector<fp_t> x(11), y(11);
97 for(
size_t i=0;i<11;i++) {
109 res=dx*prenum/preden*(c1*(y[0]+y[10])+
118 for(
size_t i=0;i<10;i++) {
119 x[i]=((fp_t)i)*dx2+a;
129 fp_t res2=dx2*prenum/preden*(c1*(y[0]+y[9])+
145 template<
class func_t=funct,
class fp_t=
double>
146 class inte_newton_cotes_open :
public inte<func_t,fp_t> {
153 virtual int integ_err(func_t &func, fp_t a, fp_t b,
154 fp_t &res, fp_t &err) {
156 std::vector<fp_t> x(7), y(7);
161 for(
size_t i=0;i<7;i++) {
162 x[i]=((fp_t)i+1)*dx+a;
172 res=dx*prenum/preden*(c1*y[0]+c2*y[1]+c3*y[2]+c4*y[3]+
173 c3*y[4]+c2*y[5]+c1*y[6]);
178 for(
size_t i=0;i<6;i++) {
179 x[i]=((fp_t)i+1)*dx2+a;
188 fp_t res2=dx2*prenum/preden*(c1*(y[0]+y[5])+
239 template<
class func_t=
funct,
240 class def_inte_t=inte_gauss56_cern<
funct,double,
241 inte_gauss56_coeffs_double>,
242 size_t nsub=100,
class fp_t=
double>
245 #ifndef DOXYGEN_INTERNAL
284 fp_t &res, fp_t &err) {
286 fp_t tvals=0.0, terss, xlob, xhib, yhib=0.0, te,
root=0.0;
320 fp_t bin=(b-a)/((fp_t)nsubdivd);
321 for(i=0;i<nsubdivd;i++) {
325 if (i==nsubdivd-1)
xhi[i]=b;
327 it->integ_err(func,xlob,xhib,
tval[i],te);
332 for(
size_t iter=1;iter<=nsub;iter++) {
345 std::cout <<
"inte_adapt_cern Iter: " << iter;
346 std::cout.setf(std::ios::showpos);
347 std::cout <<
" Res: " << tvals;
348 std::cout.unsetf(std::ios::showpos);
349 std::cout <<
" Err: " << sqrt(two*terss);
351 std::cout <<
" Tol: " << this->
tol_abs << std::endl;
353 std::cout <<
" Tol: " << this->
tol_rel*abs(tvals)
358 std::cout <<
"Press a key and type enter to continue. " ;
364 root=sqrt(two*terss);
377 std::string s=
"Reached maximum number ("+
itos(nsub)+
378 ") of subdivisions in inte_adapt_cern::integ_err().";
394 fp_t xnew=(
xlo[ibig]+
xhi[ibig])/two;
446 fp_t &value, fp_t &errsq) {
457 template<
class vec_t>
473 #ifdef O2SCL_NEVER_DEFINED
480 template<
class func_t=
funct,
481 class def_inte_t=inte_gauss56_cern<
funct,double,
482 inte_gauss56_coeffs_double>,
485 class inte_qagil_cern :
public inte<func_t,fp_t> {
496 virtual fp_t transform(fp_t t) {
497 fp_t x=upper_limit-(1-t)/t, y;
513 fo=std::bind(std::mem_fn<fp_t(fp_t)>(&inte_qagil_cern::transform),
514 this,std::placeholders::_1);
524 typedef std::function<fp_t(fp_t)> internal_funct;
529 inte<internal_funct,fp_t> *it;
539 virtual int integ_err(func_t &func, fp_t a, fp_t b,
540 fp_t &res, fp_t &err) {
543 int ret=it->integ_err(fo,0.0,1.0,res,err);
550 int set_inte(inte<internal_funct,fp_t> &i) {
556 inte_adapt_cern<internal_funct,def_inte_t,nsub,fp_t> def_inte;
563 #ifndef DOXYGEN_NO_O2NS