27 #include "Lib3_HCubature.h"
34 int HCubature::Wrapper(
unsigned int xdim,
size_t npts,
const qREAL *x,
void *fdata,
unsigned int ydim,
qREAL *y) {
38 unsigned int nthreads =
self->Threads>0 ?
self->Threads : omp_get_num_procs();
39 #pragma omp parallel for num_threads(nthreads) schedule(dynamic, 1)
40 for(
int i=0; i<npts; i++) {
42 mpfr::mpreal::set_default_prec(mpfr::digits2bits(self->MPDigits));
43 int iDQMP =
self->inDQMP(x+i*xdim);
44 if( (self->IntegrandMP!=NULL) && (self->DQMP>2 || iDQMP>2) ) {
45 mpREAL mpx[xdim], mpy[ydim];
46 for(
int j=0; j<xdim; j++) mpx[j] = x[i*xdim+j];
47 self->IntegrandMP(xdim, mpx, ydim, mpy, self->mpParameter, self->mpLambda);
48 for(
int j=0; j<ydim; j++) y[i*ydim+j] = mpy[j].toFloat128();
49 }
else if(self->DQMP>1 || iDQMP>1) {
50 self->IntegrandQ(xdim, x+i*xdim, ydim, y+i*ydim, self->qParameter, self->qLambda);
52 for(
int j=0; j<ydim; j++) {
53 qREAL ytmp = y[i*ydim+j];
54 if(isnanq(ytmp) || isinfq(ytmp)) { ok =
false;
break; }
57 if(!ok && (self->IntegrandMP!=NULL)) {
58 mpREAL mpx[xdim], mpy[ydim];
59 for(
int j=0; j<xdim; j++) mpx[j] = x[i*xdim+j];
60 self->IntegrandMP(xdim, mpx, ydim, mpy, self->mpParameter, self->mpLambda);
61 for(
int j=0; j<ydim; j++) y[i*ydim+j] = mpy[j].toFloat128();
64 dREAL dx[xdim], dy[ydim];
65 for(
int j=0; j<xdim; j++) dx[j] = (
dREAL)x[i*xdim+j];
66 self->IntegrandD(xdim, dx, ydim, dy, self->dParameter, self->dLambda);
67 for(
int j=0; j<ydim; j++) y[i*ydim+j] = dy[j];
69 for(
int j=0; j<ydim; j++) {
70 qREAL ytmp = y[i*ydim+j];
71 if(isnanq(ytmp) || isinfq(ytmp)) { ok =
false;
break; }
74 if(!ok)
self->IntegrandQ(xdim, x+i*xdim, ydim, y+i*ydim, self->qParameter, self->qLambda);
77 for(
int j=0; j<ydim; j++) {
78 qREAL ytmp = y[i*ydim+j];
79 if(isnanq(ytmp) || isinfq(ytmp)) { ok =
false;
break; }
81 if(!ok && (self->IntegrandMP!=NULL)) {
82 mpREAL mpx[xdim], mpy[ydim];
83 for(
int j=0; j<xdim; j++) mpx[j] = x[i*xdim+j];
84 self->IntegrandMP(xdim, mpx, ydim, mpy, self->mpParameter, self->mpLambda);
85 for(
int j=0; j<ydim; j++) y[i*ydim+j] = mpy[j].toFloat128();
91 for(
int j=0; j<ydim; j++) {
92 qREAL ytmp = y[i*ydim+j];
93 if(isnanq(ytmp) || isinfq(ytmp)) { ok =
false;
break; }
95 if(!ok && (self->IntegrandMP!=NULL)) {
97 mpfr::mpreal::set_default_prec(mpfr::digits2bits(self->MPDigits*10));
98 mpREAL mpx[xdim], mpy[ydim];
99 for(
int j=0; j<xdim; j++) mpx[j] = x[i*xdim+j];
100 self->IntegrandMP(xdim, mpx, ydim, mpy, self->mpParameter, self->mpLambda);
101 for(
int j=0; j<ydim; j++) y[i*ydim+j] = mpy[j].toFloat128();
102 mpfr::mpreal::set_default_prec(mpfr::digits2bits(self->MPDigits));
106 for(
int j=0; j<ydim; j++) {
107 qREAL ytmp = y[i*ydim+j];
108 if(isnanq(ytmp) || isinfq(ytmp)) {
111 if(self->nNAN > self->NANMax) { NaNQ =
true;
break; }
112 else y[i*ydim+j] = 0;
115 if(self->ReIm == 1) y[i*ydim+1] = 0;
116 else if(self->ReIm == 2) y[i*ydim+0] = 0;
123 void HCubature::DefaultPrintHooker(
qREAL *result,
qREAL *epsabs,
size_t *nrun,
void *fdata) {
125 if(*nrun == self->MaxPTS + 1979)
return;
126 if(self->RunTime>0) {
127 auto cur_timer = time(NULL);
128 auto used_time = difftime(cur_timer,self->StartTimer);
129 if(used_time>self->RunTime) {
131 *nrun =
self->MaxPTS + 1979;
137 if((*nrun-self->NEval) >= self->RunPTS) {
139 char r0[64], r1[64], e0[32], e1[32];
140 quadmath_snprintf(r0,
sizeof r0,
"%.10QG", result[0]);
141 quadmath_snprintf(r1,
sizeof r1,
"%.10QG", result[1]);
142 quadmath_snprintf(e0,
sizeof e0,
"%.5QG", epsabs[0]);
143 quadmath_snprintf(e1,
sizeof e1,
"%.5QG", epsabs[1]);
144 cout <<
" N: " << (*nrun) <<
", ";
145 if(self->ReIm==3 || self->ReIm==1) cout <<
"[" << r0 <<
", " << e0 <<
"]";
146 if(self->ReIm==3 || self->ReIm==2) cout <<
"+I*[" << r1 <<
", " << e1 <<
"]";
152 if((isnanq(result[0]) || isnanq(result[1]) || isnanq(epsabs[0]) || isnanq(epsabs[1])) || (isinfq(result[0]) || isinfq(result[1]) || isinfq(epsabs[0]) || isinfq(epsabs[1]))) {
154 *nrun =
self->MaxPTS + 1979;
155 if(self->LastState>0)
self->LastState = -1;
160 if(epsabs[0] > 1E30*self->EpsAbs || epsabs[1] > 1E30*self->EpsAbs) {
162 *nrun =
self->MaxPTS + 1979;
163 if(self->LastState>0)
self->LastState = -1;
168 if((self->LastState == 0) || (epsabs[0]<=2*self->LastAbsErr[0] && epsabs[1]<=2*self->LastAbsErr[1])) {
169 self->LastResult[0] = result[0];
170 self->LastResult[1] = result[1];
171 self->LastAbsErr[0] = epsabs[0];
172 self->LastAbsErr[1] = epsabs[1];
174 self->lastNRUN = *nrun;
175 self->lastnNAN =
self->nNAN;
178 bool rExit = (epsabs[0] <
self->EpsAbs+1E-50Q) || (epsabs[0] < fabsq(result[0])*
self->EpsRel+1E-50Q);
179 bool iExit = (epsabs[1] <
self->EpsAbs+1E-50Q) || (epsabs[1] < fabsq(result[1])*
self->EpsRel+1E-50Q);
180 if(rExit && iExit && (*nrun)>
self->MinPTS) {
182 *nrun =
self->MaxPTS + 1979;
188 fn << pid <<
".int.done";
191 *nrun =
self->MaxPTS + 1979;
193 cmd <<
"rm " << fn.str();
194 auto rc = system(cmd.str().c_str());
195 if(
Verbose>10) cout <<
" Exit: " << fn.str() << endl;
199 ex HCubature::Integrate(
size_t tn) {
200 if(mpfr_buildopt_tls_p()<=0)
throw Error(
"Integrate: mpfr_buildopt_tls_p()<=0.");
202 mpfr::mpreal::set_default_prec(mpfr::digits2bits(MPDigits));
203 mpPi = mpfr::const_pi();
205 mpiEpsilon = complex<mpREAL>(0,mpfr::machine_epsilon()*100);
207 unsigned int xdim = XDim;
208 unsigned int ydim = 2;
209 qREAL result[ydim], estabs[ydim];
211 qREAL xmin[xdim], xmax[xdim];
212 for(
int i=0; i<xdim; i++) {
220 size_t _MinPTS_, _RunPTS_;
223 MaxPTS = RunPTS * RunMAX;
224 _MinPTS_ = MinPTS>0 ? MinPTS : _RunPTS_/10;
227 if(MaxPTS<10000) MaxPTS = 10000;
229 _MinPTS_ = MinPTS>0 ? MinPTS : _RunPTS_/10;
231 StartTimer = time(NULL);
233 int nok = hcubature_v(ydim, Wrapper,
this, xdim, xmin, xmax, _MinPTS_, _RunPTS_, MaxPTS, EpsAbs, EpsRel, result, estabs, tn==0 ? PrintHooker : NULL);
236 if( (cabsq(result[0]+result[1]*1.Qi) < FLT128_EPSILON) && (cabsq(estabs[0]+estabs[1]*1.Qi) < FLT128_EPSILON) ) {
237 cout <<
ErrColor <<
"HCubature Failed with 0 result returned!" <<
RESET << endl;
242 if(LastState==-1 && use_last) {
243 result[0] = LastResult[0];
244 result[1] = LastResult[1];
245 estabs[0] = LastAbsErr[0];
246 estabs[1] = LastAbsErr[1];
252 if(isnanq(result[0]) || isnanq(result[1])) FResult +=
NaN;
255 FResult += VE(
q2ex(result[0]),
q2ex(estabs[0]));
256 FResult += VE(
q2ex(result[1]),
q2ex(estabs[1])) * I;
263 if(nNAN * 1000 > NEval && NEval>0) {
264 cout <<
ErrColor <<
"NAN=" << nNAN <<
" v.s. RUN=" << NEval <<
RESET << endl;
270 int HCubature::inDQMP(
qREAL const *x) {
271 unsigned int xdim = XDim;
273 if(xdim<=MPXDim)
return 3;
275 for(
int i=0; i<xdim; i++) {
276 if(x[i] < xmin) xmin = x[i];
278 if(xmin < MPXLimit)
return 3;
282 static FT_Type last_ft = NULL;
286 for(
int i=0; i<xdim; i++) x0[i]=0.521Q;
287 qREAL ft0 = fabsq(FT(x0, qParameter));
288 if(ft0<1E-50) ft0 = 1;
291 ft = fabsq(FT(x, qParameter));
293 if(ft<MPFLimit)
return 3;
294 else if(ft<QFLimit)
return 2;
297 if(xdim <= QXDim || xmin < QXLimit)
return 2;
complex< mpREAL > mpCOMPLEX
complex< dREAL > dCOMPLEX
bool file_exists(const char *fn)
class used to wrap error message
numerical integrator using HCubature
namespace for Numerical integration with Sector Decomposition method
ex q2ex(__float128 num)
__float128 to ex