21 typedef mpfr::mpreal
mpREAL;
22 typedef const mpREAL & mpREAL_t;
24 typedef std::function<int(
unsigned yn,
mpREAL *y,
mpREAL *e,
const mpREAL & x,
void *fdata)> f1Type;
25 typedef const f1Type & f1Type_t;
26 typedef std::function<int(
unsigned yn,
mpREAL *y,
mpREAL *e,
unsigned xdim,
const mpREAL *x,
void *fdata)> fnType;
27 typedef const fnType & fnType_t;
28 typedef void (*PrintHookerType) (
mpREAL *,
mpREAL *,
size_t *,
void *);
39 int QuadPack1(
unsigned yn,
mpREAL *oval,
mpREAL *oerr, f1Type_t f, mpREAL_t epsabs, PrintHookerType PrintHooker,
void *fdata) {
41 Function F(yn,f,fdata);
43 return gk.QAG(F, 0, 1, epsabs, 0, oval, oerr, PrintHooker);
44 }
catch (
const char* reason) {
45 cout << reason << endl;
51 int QuadPackN(
unsigned yn,
mpREAL *oval,
mpREAL *oerr,
unsigned xdim, fnType_t f, mpREAL_t
eps, PrintHookerType PrintHooker,
void *fdata) {
53 auto f1 = [f,
eps](
unsigned yn,
mpREAL *y,
mpREAL *e, mpREAL_t x,
void *fdata)->
int {
56 return f(yn,y,e,1,xs,fdata);
58 return QuadPack1(yn, oval, oerr, f1,
eps, PrintHooker, fdata);
61 auto f1 = [f,
eps,xdim](
unsigned yn,
mpREAL *y,
mpREAL *e, mpREAL_t x,
void *fdata)->
int {
62 auto f2 =[f,x](
unsigned yn,
mpREAL *y1,
mpREAL *e1,
unsigned xdim1,
const mpREAL *xs1,
void *fdata1)->
int {
66 for(
int i=0; i<xdim1; i++) xs[i+1] = xs1[i];
68 for(
int i=0; i<xdim1; i++) xs[i] = xs1[i];
71 return f(yn,y1,e1,xdim1+1,xs,fdata1);
73 return QuadPackN(yn,y,e,xdim-1,f2,
eps/xdim,NULL,fdata);
75 return QuadPack1(yn, oval, oerr, f1,
eps, PrintHooker, fdata);
86 ex QuadMP::mp2ex(
const mpREAL & num) {
88 oss.precision(MPDigits);
90 string sn = oss.str();
91 numeric ret(sn.c_str());
95 int QuadMP::Wrapper(
unsigned ydim,
mpREAL *y,
mpREAL *e,
unsigned xdim,
const mpREAL *x,
void *fdata) {
96 auto self = (
QuadMP*)fdata;
97 self->IntegrandMP(xdim, x, ydim, y, self->mpParameter, self->mpLambda);
100 for(
int j=0; j<ydim; j++) {
101 if(!isfinite(y[j])) { ok =
false;
break; }
104 mpfr::mpreal::set_default_prec(mpfr::digits2bits(self->MPDigits*10));
105 self->IntegrandMP(xdim, x, ydim, y, self->mpParameter, self->mpLambda);
106 mpfr::mpreal::set_default_prec(mpfr::digits2bits(self->MPDigits));
108 for(
int j=0; j<ydim; j++) e[j] = 0;
112 void QuadMP::DefaultPrintHooker(
mpREAL* result,
mpREAL* epsabs,
size_t * nrun,
void *fdata) {
113 auto self = (
QuadMP*)fdata;
114 if(*nrun == self->nGK + 1979)
return;
115 if(self->RunTime>0) {
116 auto cur_timer = time(NULL);
117 auto used_time = difftime(cur_timer,self->StartTimer);
118 if(used_time>self->RunTime) {
120 *nrun =
self->nGK + 1979;
128 auto e0 = epsabs[0].toString(3);
129 auto e1 = epsabs[1].toString(3);
130 cout <<
" L: " << (*nrun) <<
", ";
131 if(self->ReIm==3 || self->ReIm==1) cout <<
"[" << r0 <<
", " << e0 <<
"]";
132 if(self->ReIm==3 || self->ReIm==2) cout <<
"+I*[" << r1 <<
", " << e1 <<
"]";
137 bool rExit = (epsabs[0] <
self->EpsAbs+1E-50Q) || (epsabs[0] < fabs(result[0])*
self->EpsRel+1E-50Q);
138 bool iExit = (epsabs[1] <
self->EpsAbs+1E-50Q) || (epsabs[1] < fabs(result[1])*
self->EpsRel+1E-50Q);
140 *nrun =
self->nGK + 1979;
146 fn << pid <<
".int.done";
149 *nrun =
self->nGK + 1979;
151 cmd <<
"rm " << fn.str();
152 auto rc = system(cmd.str().c_str());
153 if(
Verbose>10) cout <<
" Exit: " << fn.str() << endl;
157 ex QuadMP::Integrate(
size_t n) {
158 if(mpfr_buildopt_tls_p()<=0)
throw Error(
"Integrate: mpfr_buildopt_tls_p()<=0.");
160 mpfr::mpreal::set_default_prec(mpfr::digits2bits(MPDigits));
161 mpPi = mpfr::const_pi();
163 mpiEpsilon = complex<mpREAL>(0,mpfr::machine_epsilon()*100);
165 unsigned int xdim = XDim;
166 unsigned int ydim = 2;
167 mpREAL result[ydim], estabs[ydim];
171 _nGK_ = n==0 ? nGK : n;
172 StartTimer = time(NULL);
173 auto nok = QuadPackN(ydim, result, estabs, xdim, Wrapper, EpsAbs, n==0 ? PrintHooker : NULL,
this);
176 mpREAL abs_res = sqrt(result[0]*result[0]+result[1]*result[1]);
177 mpREAL abs_est = sqrt(estabs[0]*estabs[0]+estabs[1]*estabs[1]);
178 mpREAL mpfr_eps = 10*mpfr::machine_epsilon();
179 if( (abs_res < mpfr_eps) && (abs_est < mpfr_eps) ) {
180 cout <<
ErrColor <<
"QuadMP Failed with 0 result returned!" <<
RESET << endl;
186 if(isnan(result[0]) || isnan(result[1])) FResult +=
NaN;
189 FResult += VE(mp2ex(result[0]), mp2ex(estabs[0]));
190 FResult += VE(mp2ex(result[1]), mp2ex(estabs[1])) * I;
complex< mpREAL > mpCOMPLEX
bool file_exists(const char *fn)
class used to wrap error message
numerical integrator using TanhSinhMP
namespace for Numerical integration with Sector Decomposition method