21 typedef mpfr::mpreal
mpREAL;
22 typedef const mpREAL & mpREAL_t;
24 typedef std::function<int(
unsigned ydim,
mpREAL *y,
mpREAL *e,
const mpREAL & x,
void *fdata)> f1Type;
25 typedef const f1Type & f1Type_t;
26 typedef std::function<int(
unsigned ydim,
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 *);
37 int TanhSinh1(
mpREAL *oval,
mpREAL *oerr, f1Type_t f,
unsigned ydim, mpREAL_t epsrel, PrintHookerType PrintHooker,
void *fdata) {
41 const mpREAL tol = 10*epsrel;
42 mpREAL s[ydim], e_s[ydim], val[ydim], err[ydim], v[ydim], h = 2;
43 if(f(ydim, s, e_s, c, fdata))
return FAILURE;
46 mpREAL p[ydim], fp[ydim], fm[ydim], q, t, eh;
47 for(
int j=0; j<ydim; j++) p[j] = fp[j] = fm[j] = 0;
53 if(omp_in_parallel() && ( !omp_get_nested() || omp_get_active_level()>1 )) parallel =
false;
56 int total = CPUCORES/2;
57 mpREAL xs[total], ws[total], fps[total*ydim], e_fps[total*ydim], fms[total*ydim], e_fms[total*ydim];
58 int RC1[total], RC2[total];
59 for(
int i=0; i<total; i++) {
62 ws[i] = (t+1/t)*r/(1+u);
67 auto prec = mpfr::mpreal::get_default_prec();
68 auto rnd = mpfr::mpreal::get_default_rnd();
69 #pragma omp parallel for
70 for(
int ii=0; ii<2*total; ii++) {
71 int i = ii/2, i2 = ii%2;
72 mpfr::mpreal::set_default_prec(prec);
73 mpfr::mpreal::set_default_rnd(rnd);
75 if(i2==0) {
if (
a+x >
a) RC1[i] = f(ydim,fps+ydim*i,e_fps+ydim*i,
a+x,fdata); }
76 else if(i2==1) {
if (b-x < b) RC2[i] = f(ydim,fms+ydim*i,e_fms+ydim*i,b-x,fdata); }
80 for(
int i=0; i<total; i++) {
81 if(RC1[i]!=0 || RC2[i]!=0)
return FAILURE;
84 for(
int j=0; j<ydim; j++) {
85 if (
a+x >
a)
if (isfinite(fps[i])) {
86 fp[j] = fps[i*ydim+j];
87 if(e_s[j] < e_fps[i*ydim+j]) e_s[j] = e_fps[i*ydim+j];
89 if (b-x < b)
if (isfinite(fms[i])) {
90 fm[j] = fms[i*ydim+j];
91 if(e_s[j] < e_fms[i*ydim+j]) e_s[j] = e_fms[i*ydim+j];
93 q = ws[i]*(fp[j]+fm[j]);
95 ok = ok && (fabs(q)<=epsrel*fabs(p[j]));
108 mpREAL y[ydim], e_y[ydim];
110 if(f(ydim,y,e_y,
a+x,fdata))
return FAILURE;
111 for(
int j=0; j<ydim; j++)
if (isfinite(y[j])) {
113 if(e_s[j] < e_y[j]) e_s[j] = e_y[j];
117 if(f(ydim,y,e_y,b-x,fdata))
return FAILURE;
118 for(
int j=0; j<ydim; j++)
if (isfinite(y[j])) {
120 if(e_s[j] < e_y[j]) e_s[j] = e_y[j];
124 for(
int j=0; j<ydim; j++) {
127 ok = ok && (fabs(q)<=epsrel*fabs(p[j]));
132 for(
int j=0; j<ydim; j++) {
136 err[j] = fabs(
d*h*v[j]) + (b-
a)*e_s[j];
139 cout <<
RED <<
"Large Level: " <<
RESET << k << endl;
140 for(
int j=0; j<ydim; j++) {
141 cout << val[j] <<
" +- " << err[j].toString(3) << endl;
146 if(PrintHooker) PrintHooker(val, err, &kk, fdata);
148 for(
int j=0; j<ydim; j++) {
155 for(
int j=0; j<ydim; j++) {
156 ok = ok && (fabs(v[j])<=tol*fabs(s[j]));
162 for(
int j=0; j<ydim; j++) {
169 int TanhSinhN(
mpREAL *val,
mpREAL *err, fnType_t f,
unsigned xdim,
unsigned ydim, mpREAL_t
eps, PrintHookerType PrintHooker,
void *fdata) {
171 auto f1 = [f](
unsigned ydim,
mpREAL *y,
mpREAL *e, mpREAL_t x,
void *fdata)->
int {
174 if(f(ydim,y,e,1,xs,fdata))
return FAILURE;
177 return TanhSinh1(val,err,f1,ydim,
eps,PrintHooker,fdata);
179 auto f1 = [f,xdim,
eps](
unsigned ydim,
mpREAL *y,
mpREAL *e, mpREAL_t x,
void *fdata)->
int {
180 auto f2 =[f,x](
unsigned ydim,
mpREAL *y1,
mpREAL *e1,
unsigned xdim1,
const mpREAL *xs1,
void *fdata)->
int {
184 for(
int i=1; i<=xdim1; i++) xs[i] = xs1[i-1];
186 for(
int i=0; i<xdim1; i++) xs[i] = xs1[i];
189 if(f(ydim,y1,e1,xdim1+1,xs,fdata))
return FAILURE;
192 return TanhSinhN(y,e,f2,xdim-1,ydim,
eps/xdim,NULL,fdata);
194 return TanhSinh1(val,err,f1,ydim,
eps,PrintHooker,fdata);
205 TanhSinhMP::TanhSinhMP(
size_t k) { K = k; }
207 ex TanhSinhMP::mp2ex(mpREAL_t num) {
209 oss.precision(MPDigits);
211 string sn = oss.str();
212 numeric ret(sn.c_str());
216 int TanhSinhMP::Wrapper(
unsigned ydim,
mpREAL *y,
mpREAL *e,
unsigned xdim,
const mpREAL *x,
void *fdata) {
217 for(
int j=0; j<ydim; j++) e[j] = 0;
219 self->IntegrandMP(xdim, x, ydim, y, self->mpParameter, self->mpLambda);
221 for(
int j=0; j<ydim; j++) {
222 if(!isfinite(y[j])) { ok =
false;
break; }
225 mpfr::mpreal::set_default_prec(mpfr::digits2bits(self->MPDigits*10));
226 self->IntegrandMP(xdim, x, ydim, y, self->mpParameter, self->mpLambda);
227 mpfr::mpreal::set_default_prec(mpfr::digits2bits(self->MPDigits));
232 void TanhSinhMP::DefaultPrintHooker(
mpREAL* result,
mpREAL* epsabs,
size_t * nrun,
void *fdata) {
234 if(*nrun == self->K + 1979)
return;
235 if(self->RunTime>0) {
236 auto cur_timer = time(NULL);
237 auto used_time = difftime(cur_timer,self->StartTimer);
238 if(used_time>self->RunTime) {
240 *nrun =
self->K + 1979;
248 auto e0 = epsabs[0].toString(3);
249 auto e1 = epsabs[1].toString(3);
250 cout <<
" L: " << (*nrun) <<
", ";
251 if(self->ReIm==3 || self->ReIm==1) cout <<
"[" << r0 <<
", " << e0 <<
"]";
252 if(self->ReIm==3 || self->ReIm==2) cout <<
"+I*[" << r1 <<
", " << e1 <<
"]";
257 bool rExit = (epsabs[0] <
self->EpsAbs+1E-50Q) || (epsabs[0] < fabs(result[0])*
self->EpsRel+1E-50Q);
258 bool iExit = (epsabs[1] <
self->EpsAbs+1E-50Q) || (epsabs[1] < fabs(result[1])*
self->EpsRel+1E-50Q);
260 *nrun =
self->K + 1979;
266 fn << pid <<
".int.done";
269 *nrun =
self->K + 1979;
271 cmd <<
"rm " << fn.str();
272 system(cmd.str().c_str());
273 if(
Verbose>10) cout <<
" Exit: " << fn.str() << endl;
277 ex TanhSinhMP::Integrate(
size_t n) {
279 CPUCORES = omp_get_num_procs();
280 if(mpfr_buildopt_tls_p()<=0)
throw Error(
"Integrate: mpfr_buildopt_tls_p()<=0.");
282 mpfr::mpreal::set_default_prec(mpfr::digits2bits(MPDigits));
283 mpPi = mpfr::const_pi();
286 unsigned int xdim = XDim;
287 unsigned int ydim = 2;
288 mpREAL result[ydim], estabs[ydim];
291 int nok = TanhSinhN(result, estabs, Wrapper, xdim, ydim, epsrel, NULL,
this);
292 epsrel = min(EpsAbs/(fabs(result[0])+EpsAbs), EpsAbs/(fabs(result[1])+EpsAbs));
294 StartTimer = time(NULL);
295 nok = TanhSinhN(result, estabs, Wrapper, xdim, ydim, epsrel, n==0 ? PrintHooker : NULL,
this);
299 mpREAL abs_res = sqrt(result[0]*result[0]+result[1]*result[1]);
300 mpREAL abs_est = sqrt(estabs[0]*estabs[0]+estabs[1]*estabs[1]);
301 mpREAL mpfr_eps = 10*mpfr::machine_epsilon();
302 if( (abs_res < mpfr_eps) && (abs_est < mpfr_eps) ) {
303 cout <<
ErrColor <<
"TanhSinhMP Failed with 0 result returned!" <<
RESET << endl;
309 if(!isfinite(result[0]) || !isfinite(result[1])) FResult +=
NaN;
312 FResult += VE(mp2ex(result[0]), mp2ex(estabs[0]));
313 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