HepLib
Integrates.cpp
Go to the documentation of this file.
1 
6 #include "SD.h"
7 #include <math.h>
8 #include <cmath>
9 
10 namespace HepLib::SD {
11 
18  ex SecDec::ContinuousWRA(ex expr_in, int nc) {
19  auto expr = expr_in;
20  static symbol xwra("xwra");
21  expr = xyz_pow_simplify(expr);
22  exset pows_set;
23  expr.find(sqrt(w), pows_set);
24  expr.find(pow(w1,w2), pows_set);
25  exmap pow_map;
26  for(auto item : pows_set) {
27  if(item.has(WRA(w))) {
28  if(item.match(pow(w1,w2)) && !item.op(1).info(info_flags::integer)) {
29  pow_map[item] = exp(log(item.op(0)) * item.op(1));
30  } else if(item.match(sqrt(w))) {
31  pow_map[item] = exp(log(item.op(0))/2);
32  }
33  }
34  }
35 
36  expr = expr.subs(pow_map);
37  expr = exp_simplify(expr);;
38 
39  int log_id = 0;
40  exset logs_set;
41  expr.find(log(w), logs_set);
42  exmap log_map;
43  lst id_logz_lst;
44  for(auto item : logs_set){
45  if(item.has(WRA(w))) {
46  log_id++;
47  id_logz_lst.append(iWF(log_id,item.op(0)));
48  log_map[item] = iWF(log_id);
49  }
50  }
51  if(log_id<1) return expr.subs(WRA(w)==w);
52 
53  expr = expr.subs(log_map).subs(WRA(w)==w);
54 
55  exmap log_map2;
56  for(auto id_logz : id_logz_lst) {
57  auto id = id_logz.op(0);
58  auto zz = id_logz.op(1);
59  exset wra_set;
60  zz.find(WRA(w), wra_set);
61  if(wra_set.size()<1) return expr;
62  if(wra_set.size()!=1) {
63  cout << zz << endl;
64  throw Error("ContinuousWRA: too many WRA.");
65  }
66  zz = zz.subs(WRA(w)==xwra);
67  ex wra = (*(wra_set.begin())).op(0);
68 
69  ex ret = log(zz.subs(xwra==wra));
70  if(nc<1) {
71  log_map2[iWF(id)] = ret;
72  break;
73  }
74  int total=0;
75  int ReIm[nc][2];
76  for(int k=0; k<=nc; k++) {
77  ex zzk;
78  if(k==0) zzk = NN(zz.subs(xwra==wra/(25*nc)));
79  else zzk = NN(zz.subs(xwra==k*wra/nc));
80  if(!is_a<numeric>(zzk)) throw Error("ContinuousWRA: zzk is not numeric: "+ex2str(zzk));
81  auto nzzk = ex_to<numeric>(zzk);
82  auto curR = real(nzzk);
83  auto curI = imag(nzzk);
84  if(is_zero(curI) && k==0 && curR<0) ReIm[total][1] = -1;
85  else if(is_zero(curR) || is_zero(curR)) continue;
86  else ReIm[total][1] = curI>0 ? 1 : -1;
87  ReIm[total][0] = curR>0 ? 1 : -1;
88  total++;
89  }
90 
91  int cutN = 0;
92  for(int k=0; k<total-1; k++) {
93  if(ReIm[k][0]*ReIm[k+1][0]<0 && ReIm[k][1]*ReIm[k+1][1]<0) throw Error("ContinuousWRA: 1<->3 or 2<->4 happened.");
94  if(ReIm[k][0]<0 && ReIm[k+1][0]<0 && ReIm[k][1]*ReIm[k+1][1]<0) {
95  if(ReIm[k][1]>0) cutN++;
96  else cutN--;
97  }
98  }
99  if(cutN!=0) ret += I * cutN * 2 * Pi;
100  log_map2[iWF(id)] = ret;
101  }
102  expr = expr.subs(log_map2);
103  if(expr.has(iWF(w))) throw Error("ContinuousWRA: iWF(w) still exists in final result.");
104  return expr;
105  }
106 
113  void SecDec::Integrates(const string & key, const string & pkey, int kid) {
114  if(MPDigits>0) mpfr::mpreal::set_default_prec(mpfr::digits2bits(MPDigits));
115  if(IsZero) return;
116 
117  if(Verbose>0) cout << Color_HighLight << " Integrates @ " << now() << RESET << endl;
118 
119  lst lstRE;
120  auto pid = getpid();
121  ostringstream fsofn, sofn, cmd;
122  if(key == "") {
123  sofn << pid << ".so";
124  fsofn << pid << "F.so";
125  } else {
126  set_precision(NNDigits);
127  sofn << key << ".so";
128  ostringstream garfn;
129  garfn << key << ".ci.gar";
130  archive ar;
131  ifstream in(garfn.str());
132  in >> ar;
133  in.close();
134  auto c = ar.unarchive_ex("c");
135  if(c!=19790923) throw Error("Integrates: *.ci.gar error!");
136  auto gl = ar.unarchive_ex("soLimit");
137  eps_lst = ex_to<lst>(ar.unarchive_ex("eps_lst"));
138  soLimit = ex2int(gl);
139  auto res = ar.unarchive_ex("res");
140  ciResult.clear();
141  for(auto item : ex_to<lst>(res)) ciResult.push_back(ex_to<lst>(item));
142  garfn.clear();
143  garfn.str("");
144  garfn << key;
145  if(pkey != "") garfn << "-" << pkey;
146  garfn << ".las.gar";
147  if(file_exists(garfn.str().c_str())) {
148  archive la_ar;
149  ifstream la_in(garfn.str());
150  la_in >> la_ar;
151  la_in.close();
152  auto la_c = la_ar.unarchive_ex("c");
153  auto la_res = la_ar.unarchive_ex("res");
154  if(la_c!=19790923) throw Error("Integrates: *.ci.gar error!");
155  for(auto item : ex_to<lst>(la_res)) {
156  LambdaMap[item.op(0)] = item.op(1);
157  }
158  }
159 
160  if(kid>0) {
161  garfn.clear();
162  garfn.str("");
163  garfn << key;
164  if(pkey != "") garfn << "-" << pkey;
165  garfn << ".res.gar";
166  if(!file_exists(garfn.str().c_str())) {
167  throw Error("Integrates: File Not Found: " + garfn.str());
168  }
169 
170  archive res_ar;
171  ifstream res_in(garfn.str());
172  res_in >> res_ar;
173  res_in.close();
174  auto res_c = res_ar.unarchive_ex("c");
175  auto relst = res_ar.unarchive_ex("relst");
176  if(res_c!=19790923) throw Error("*.res.gar error with kid!");
177  lstRE = ex_to<lst>(relst);
178  }
179  reset_precision();
180  }
181 
182  void* main_module = dlopen(sofn.str().c_str(), RTLD_NOW);
183  if(main_module == nullptr) {
184  main_module = dlopen(("./"+sofn.str()).c_str(), RTLD_NOW);
185  if(main_module == nullptr) {
186  cout << "dlerror(): " << dlerror() << endl;
187  throw Error("Integrates: could not open main module!");
188  }
189  }
190 
191  if(!Debug && key == "") {
192  if(file_exists(fsofn.str().c_str())) remove(fsofn.str().c_str());
193  remove(sofn.str().c_str());
194  }
195  vector<void*> ex_modules;
196  for(int n=1; true; n++) {
197  ostringstream ex_sofn;
198  if(key == "") {
199  ex_sofn << pid << "X" << n << ".so";
200  } else {
201  ex_sofn << key << "X" << n << ".so";
202  }
203  if(file_exists(ex_sofn.str().c_str())) {
204  void* module = dlopen(ex_sofn.str().c_str(), RTLD_NOW);
205  if(module == nullptr) {
206  module = dlopen(("./"+ex_sofn.str()).c_str(), RTLD_NOW);
207  if(module == nullptr) {
208  cout << "dlerror(): " << dlerror() << endl;
209  throw Error("Integrates: could not open ex-module!");
210  }
211  }
212  ex_modules.push_back(module);
213  if(!Debug && key == "") remove(ex_sofn.str().c_str());
214  } else break;
215  }
216 
217  int npara = 0;
218  lst plRepl;
219  for(auto kv : Parameter) {
220  plRepl.append(PL(kv.first)==kv.second);
221  if(kv.first>npara) npara = kv.first;
222  }
223  plRepl.sort();
224  plRepl.unique();
225 
226  int total = ciResult.size(), current = 0;
227  qREAL stot = sqrtq(total*1.Q);
228 
229  ResultError = 0;
230  //----------------------------------------------------------------
231  for(auto &item : ciResult) {
232  current++;
233  if(kid>0 && current != kid) continue;
234  if(Verbose>0) {
235  cout << "\r \r";
236  cout << PRE << "\\--Integrating [" <<current<<"/"<<total<< "] " << flush;
237  }
238 
239  unsigned int xsize = 0;
240  ex co;
241  exvector xs, fxs;
242  xsize = ex_to<numeric>(item.op(1)).to_int();
243  co = item.op(2).subs(plRepl).subs(iEpsilon==iEpsilonN);
244  if(co.has(WRA(w))) co = ContinuousWRA(co);
245 
246  if(xsize<1) {
247  // { expr, xs.size(), kvf.op(0), -1}
248  ex exint = item.op(0).subs(plRepl).subs(iEpsilon==iEpsilonN);
249  if(exint.has(WRA(w))) exint = ContinuousWRA(exint);
250  ex res = VE(NN(exint),0);
251  if(Verbose>5) {
252  cout << "XDim=" << xsize << endl;
253  cout << Color_HighLight << " IRes = "<< HepLib::SD::VEResult(VESimplify(res)) << RESET << endl;
254  }
255  if(kid>0) {
256  lstRE.let_op(kid-1) = res*co;
257  break;
258  }
259  lstRE.append(res*co);
260  continue;
261  }
262 
263  //{ idx, xs.size(), kvf.op(0), ft_n }
264  ex intid = item.op(0);
265  ex ftid = item.op(3);
266 
267  if(co.is_zero()) continue;
268  co = collect_ex(co, eps);
269  if(co.is_zero()) continue;
270  if(co.has(PL(w))) throw Error("Integrates: PL found @ " + ex2str(co));
271  qREAL cmax = -1;
272  int reim = 0;
273  if(ReIm==3) reim = 3;
274 
275  for(int si=co.ldegree(eps); si<=co.degree(eps); si++) {
276  auto tmp = co.coeff(eps, si);
277  if(tmp.has(eps)) throw Error("Integrates: eps found @ " + ex2str(tmp));
278  tmp = collect_ex(tmp, ep);
279  for(int i=tmp.ldegree(ep); i<=tmp.degree(ep); i++) {
280  auto ccRes = NN(tmp.coeff(ep, i)).expand();
281  lst css;
282  css.append(ccRes);
283  if(is_a<add>(ccRes)) {
284  for(auto item : ccRes) css.append(item);
285  }
286 
287  exmap eps_map;
288  ex epn = ex(1)/111;
289  for(auto epi : eps_lst) {
290  eps_map[epi.op(0)] = epn;
291  epn = epn / 13;
292  }
293 
294  for(int ci=0; ci<css.nops(); ci++) {
295  auto nt = NN(css.op(ci).subs(epz==1).subs(log(vs)==1).subs(vs==1).subs(nReplacement).subs(lst{CV(w1,w2)==w2}).subs(eps_map),30);
296  if(nt.has(ep)) throw Error("Integrates: ep found @ nt = "+ex2str(nt));
297 
298  lst nt_lst;
299  if(!is_a<numeric>(nt)) {
300  auto cv_lst = collect_lst(nt,[](const ex &e)->bool{return Symbol::has(e);});
301  for(auto nti : cv_lst) {
302  auto nnt = nti.op(0);
303  if(!is_a<numeric>(nnt)) throw Error("Integrates: Not a number with nt = "+ex2str(nnt));
304  nt_lst.append(nnt);
305  }
306  } else nt_lst.append(nt);
307 
308  for(auto nnt : nt_lst) {
309  if(ReIm!=3 && reim!=3) {
310  if(ex_to<numeric>(nnt).imag()==0) {
311  if(reim==2) reim = 3;
312  else reim = 1;
313  } else if(ex_to<numeric>(nnt).real()==0) {
314  if(reim==1) reim = 3;
315  else reim = 2;
316  } else {
317  reim = 3;
318  }
319  }
320  nnt = NN(abs(nnt)); // no PL here, then nReplacement
321 
322  qREAL qnt = ex2q(nnt);
323  if(qnt > cmax) cmax = qnt;
324  }
325  }
326  }
327  }
328  if(cmax<=0) {
329  while(true) {
330  auto co2 = subs(co,lst{exp(w1*log(w2)+w3)==pow(w2,w1)*exp(w3),exp(w1*log(w2))==pow(w2,w1)});
331  if(is_zero(co2-co)) break;
332  co = co2;
333  }
334  if(normal(co).is_zero()) continue;
335  throw Error("Integrates: cmax<=0 with co = "+ex2str(co));
336  }
337  if(reim!=3 && ReIm!=3) {
338  if(ReIm==2) {
339  if(reim==1) reim=2;
340  else if(reim==2) reim=1;
341  }
342  }
343 
344  char d1[20], d2[20];
345  sprintf(d1, "%.6G", (double)(EpsAbs/cmax/stot));
346  sprintf(d2, "%.6G", (double)cmax);
347  if(Verbose>5) cout << "XDim=" << xsize << ", EpsAbs=" << d1 << "/" << d2 << endl;
348 
349  auto las = LambdaMap[ftid];
350  bool hasF = (ftid>0);
351  if(hasF && las.is_zero()) throw Error("Integrates: lambda with the key(ft_n=" + ex2str(ftid) + ") is NOT found!");
352 
353  if(hasF && !is_a<lst>(las)) {
354  if(!is_zero(las-ex(1979))) { // the convention for xPositive or explict real mode
355  throw Error("Integrates: something is wrong with the F-term @ ft_n = "+ex2str(ftid) + ", las=" + ex2str(las));
356  } else {
357  hasF = false;
358  }
359  }
360 
361  IntegratorBase::SDD_Type fpD = nullptr;
362  IntegratorBase::SDQ_Type fpQ = nullptr;
363  IntegratorBase::SDMP_Type fpMP = nullptr;
364  IntegratorBase::FT_Type ftp = nullptr;
365  int idx = ex_to<numeric>(intid).to_int();
366  auto module = main_module;
367  if(idx>=soLimit) module = ex_modules[idx/soLimit-1];
368  ostringstream fname;
369  if(hasF) fname << "C";
370  fname << "SDD_" << idx;
371  fpD = (IntegratorBase::SDD_Type)dlsym(module, fname.str().c_str());
372 
373  if(fpD==NULL) {
374  cout << "dlerror(): " << dlerror() << endl;
375  throw Error("Integrates: fpD==NULL");
376  }
377 
378  fname.clear();
379  fname.str("");
380  if(hasF) fname << "C";
381  fname << "SDQ_" << idx;
382  fpQ = (IntegratorBase::SDQ_Type)dlsym(module, fname.str().c_str());
383  if(fpQ==NULL) {
384  cout << "dlerror(): " << dlerror() << endl;
385  throw Error("Integrates: fpQ==NULL");
386  }
387 
388  fname.clear();
389  fname.str("");
390  if(hasF) fname << "C";
391  fname << "SDMP_" << idx;
392  fpMP = (IntegratorBase::SDMP_Type)dlsym(module, fname.str().c_str());
393  if(fpMP==NULL) {
394  cout << "dlerror(): " << dlerror() << endl;
395  throw Error("Integrates: fpMP==NULL");
396  }
397 
398  if(is_a<lst>(las)) {
399  fname.clear();
400  fname.str("");
401  fname << "FT_" << idx;
402  ftp = (IntegratorBase::FT_Type)dlsym(module, fname.str().c_str());
403  if(ftp==NULL) {
404  cout << "dlerror(): " << dlerror() << endl;
405  throw Error("Integrates: ftp==NULL.");
406  }
407  }
408 
409  dREAL dlas[las.nops()], dpl[npara];
410  qREAL qlas[las.nops()], qpl[npara];
411  mpREAL mplas[las.nops()], mppl[npara];
412  for(auto kv : Parameter) {
413  qpl[kv.first] = ex2q(kv.second);
414  dpl[kv.first] = (dREAL)qpl[kv.first];
415  mppl[kv.first] = qpl[kv.first];
416  }
417 
418  if(Integrator==NULL) Integrator = new HCubature();
419  Integrator->ReIm = reim;
420  if(MPDigits>0) Integrator->MPDigits = MPDigits;
421  Integrator->IntegrandD = fpD;
422  Integrator->IntegrandQ = fpQ;
423  Integrator->IntegrandMP = fpMP;
424  Integrator->FT = ftp;
425  Integrator->dParameter = dpl;
426  Integrator->dLambda = dlas;
427  Integrator->qParameter = qpl;
428  Integrator->qLambda = qlas;
429  Integrator->mpParameter = mppl;
430  Integrator->mpLambda = mplas;
431  Integrator->XDim = xsize;
432 
433  if(hasF) {
434  Integrator->EpsAbs = EpsAbs/cmax/stot/2;
435  Integrator->EpsRel = 0;
436 
437  qREAL lamax = ex2q(las.op(las.nops()-1));
438  if(lamax > IntLaMax) lamax = IntLaMax;
439  int ctryR = 0, ctry = 0, ctryL = 0;
440  int smin = -1;
441  ex min_err, min_res;
442  size_t min_eval;
443  qREAL log_lamax = log10q(lamax);
444  qREAL log_lamin = log_lamax-1.Q;
445 
446  ostringstream las_fn;
447  las_fn << key;
448  if(pkey != "") las_fn << "-" << pkey;
449  las_fn << "-" << current << ".las";
450  if(use_las && file_exists(las_fn.str().c_str())) {
451  std::ifstream las_ifs;
452  las_ifs.open(las_fn.str(), ios::in);
453  if (!las_ifs) throw Error("Integrates: failed to open *.las file!");
454  for(int i=0; i<las.nops()-1; i++) {
455  dREAL la_tmp;
456  las_ifs >> la_tmp;
457  dlas[i] = la_tmp;
458  qlas[i] = la_tmp;
459  mplas[i] = la_tmp;
460  }
461  las_ifs.close();
462  auto res = Integrator->Integrate();
463  auto res_tmp = res.subs(VE(w1, w2)==w2);
464  auto err = real_part(res_tmp);
465  if(err < imag_part(res_tmp)) err = imag_part(res_tmp);
466  min_err = err;
467  min_res = res;
468  } else {
469  // ---------------------------------------
470  auto prec = cout.precision();
471  cout.precision(5);
472  while(true) {
473  smin = -1;
474  min_err = 0;
475  ex lastResErr = 0;
476  bool err_break = false;
477  for(int s=0; s<=LambdaSplit; s++) {
478  if(Verbose>10 && s==0) {
479  if(ctryR>0 || ctry>0 || ctryL>0)
480  cout << " ------------------------------" << endl;
481  }
482  auto log_cla = (log_lamin + s * (log_lamax-log_lamin) / LambdaSplit);
483  auto cla = powq(10.Q, log_cla);
484  if(cla < 1E-10) throw Error("NIntegrate: too small lambda.");
485  for(int i=0; i<las.nops()-1; i++) {
486  qlas[i] = ex2q(las.op(i)) * cla;
487  mplas[i] = qlas[i];
488  dlas[i] = (dREAL)qlas[i];
489  }
490 
491  auto res = Integrator->Integrate(CTryI);
492  if(Verbose>10) {
493  cout << "\r \r";
494  if(res.has(NaN)) cout << " λ=" << (double)cla << "/" << Integrator->NEval << ": " << NaN << endl;
495  else cout << " λ=" << (double)cla << "/" << Integrator->NEval << ": " << VEResult2(VESimplify(res)) << endl;
496  }
497 
498  if(res.has(NaN) && s==0) continue;
499  else if(res.has(NaN)) break;
500  ex res_abs = NN(abs(res.subs(VE(w1,w2)==w1)));
501  if(lastResErr.is_zero()) lastResErr = res;
502  auto diff = VESimplify(lastResErr - res);
503  diff = diff.subs(VE(0,0)==0);
504  exset ves;
505  diff.find(VE(w0, w1), ves);
506  for(auto ve : ves) {
507  auto ve0 = abs(ve.op(0));
508  if(ve0>ve.op(1)) {
509  if(numeric("1.E10")*ve0<res_abs) continue; // avoid fluctuation aroud 0
510  if(numeric("1.E10")*ve0<q2ex(EpsAbs)) continue; // avoid fluctuation aroud 0
511  err_break = true;
512  break;
513  }
514  }
515  if(err_break) {
516  if(Verbose>10) cout << Color_HighLight << " Error Break ..." << RESET << endl;
517  break;
518  }
519  lastResErr = res;
520 
521  auto res_tmp = res.subs(VE(w1, w2)==w2);
522  auto err = real_part(res_tmp);
523  if(err < imag_part(res_tmp)) err = imag_part(res_tmp);
524  if(smin<0 || err < min_err) {
525  min_err = err;
526  min_res = res;
527  min_eval = Integrator->NEval;
528  smin = s;
529  }
530  if(s>0 && min_err < q2ex(EpsAbs/cmax/stot)) {
531  // s>0 make sure at least 2 λs compatiable
532  if(Verbose>5) {
533  cout << Color_HighLight << " λ=" << (double)cla << "/" << min_eval << ": " << HepLib::SD::VEResult(VESimplify(min_res)) << RESET << endl;
534  }
535 
536  smin = -2;
537  if(kid>0) lstRE.let_op(kid-1) = co * min_res;
538  else lstRE.append(co * min_res);
539  break;
540  }
541  if(err > 100 * min_err) break; // s>0 make sure at least 2 λs compatiable
542  }
543  if(smin == -2) break;
544  if(smin == -1) throw Error("Integrates: smin = -1, too small lambda (<1E-10)!");
545 
546  if(smin <= 0) {
547  if((!err_break) && (ctryL >= CTryL || ctryR>0)) break;
548  log_lamax = log_lamin;
549  log_lamin -= 1.Q;
550  if(!err_break) ctryL++;
551  } else if(smin >= LambdaSplit) {
552  if(ctryR >= CTryR || ctryL>0) break;
553  log_lamin = log_lamax;
554  log_lamax += log10q(CTryRRatio);
555  ctryR++;
556  } else {
557  if(ctry >= CTryM) break;
558  auto la1 = log_lamin + (smin-1) * (log_lamax-log_lamin) / LambdaSplit;
559  auto la2 = log_lamin + (smin+1) * (log_lamax-log_lamin) / LambdaSplit;
560  log_lamin = la1;
561  log_lamax = la2;
562  ctry++;
563  }
564  }
565  cout.precision(prec);
566 
567  if(smin == -2) {
568  if(kid>0) break;
569  continue;
570  }
571 
572  auto log_cla = (log_lamin + smin * (log_lamax-log_lamin) / LambdaSplit);
573  auto cla = powq(10.Q, log_cla);
574  if(Verbose>5) cout << Color_HighLight << " Final λ = " << (double)cla << " / " << las.op(las.nops()-1) << RESET << endl;
575  for(int i=0; i<las.nops()-1; i++) {
576  qlas[i] = ex2q(las.op(i)) * cla;
577  dlas[i] = (dREAL)qlas[i];
578  mplas[i] = qlas[i];
579  }
580  // ---------------------------------------
581  }
582 
583  // ---------------------------------------
584  // try HookeJeeves
585  // ---------------------------------------
586  if( use_ErrMin && (min_err > q2ex(1E5 * EpsAbs/cmax/stot)) ) {
587  ErrMin::lastResErr = min_res;
588  auto miner = new HookeJeeves();
589  ErrMin::miner = miner;
590  ErrMin::Integrator = Integrator;
591  dREAL oo[las.nops()-1], ip[las.nops()-1];
592  for(int i=0; i<las.nops()-1; i++) ip[i] = oo[i] = (dREAL)qlas[i];
593  ErrMin::lambda = oo;
594  ErrMin::err_max = 1E100;
595  auto oerrmin = ErrMin::err_min;
596  ErrMin::err_min = (dREAL)(oerrmin < 0 ? -oerrmin * ex2q(min_err) : oerrmin/cmax);
597  ErrMin::RunRND = 0;
598  miner->Minimize(las.nops()-1, ErrMin::IntError, ip);
599  delete miner;
600  ErrMin::err_min = oerrmin;
601  for(int i=0; i<las.nops()-1; i++) {
602  qlas[i] = ErrMin::lambda[i];
603  dlas[i] = (dREAL)qlas[i];
604  mplas[i] = qlas[i];
605  }
606 
607  Integrator->dLambda = dlas;
608  Integrator->qLambda = qlas; // Integrator->Lambda changed in ErrMin
609  Integrator->mpLambda = mplas;
610  if(Verbose>5) {
611  cout << Color_HighLight << " Final λs: " << RESET;
612  for(int i=0; i<xsize; i++) {
613  char buffer[128];
614  quadmath_snprintf(buffer, sizeof buffer, "%.6QG", qlas[i]);
615  cout << buffer << " ";
616  }
617  cout << endl << " ------------------------------" << endl;
618  }
619  }
620  // ---------------------------------------
621 
622  if(save_las) {
623  std::ofstream las_ofs;
624  las_ofs.open(las_fn.str(), ios::out);
625  if (las_ofs) {
626  for(int i=0; i<las.nops()-1; i++) {
627  dREAL la_tmp = (dREAL)qlas[i];
628  las_ofs << la_tmp << " ";
629  }
630  las_ofs << endl;
631  las_ofs.close();
632  }
633  }
634  }
635 
636  Integrator->EpsAbs = EpsAbs/cmax/stot;
637  Integrator->EpsRel = 0;
638  auto res = Integrator->Integrate();
639  if(Verbose>5) {
640  cout << Color_HighLight << " IRes = "<< HepLib::SD::VEResult(VESimplify(res)) << RESET << endl;
641  }
642  if(res.has(NaN)) {
643  ResultError = NaN;
644  if(kid>0) lstRE.let_op(kid-1) = NaN;
645  else lstRE.append(NaN);
646  break;
647  } else {
648  if(kid>0) {
649  lstRE.let_op(kid-1) = co * res;
650  break;
651  } else lstRE.append(co * res);
652  }
653  }
654  //----------------------------------------------------------------
655 
656  if(use_dlclose) {
657  for(auto module : ex_modules) dlclose(module);
658  dlclose(main_module);
659  }
660 
661  if(!ResultError.is_equal(NaN)) {
662  ResultError = 0;
663  for(auto item : lstRE) ResultError += item;
664  ResultError = VESimplify(ResultError);
665  }
666 
667  if(key != "") {
668  ostringstream garfn;
669  garfn << key;
670  if(pkey != "") garfn << "-" << pkey;
671  garfn << ".res.gar";
672  archive ar;
673  ar.archive_ex(ResultError, "res");
674  ar.archive_ex(lstRE, "relst");
675  ar.archive_ex(19790923, "c");
676  ofstream out(garfn.str());
677  out << ar;
678  out.close();
679  }
680  }
681 
688  void SecDec::ReIntegrates(const string & key, const string & pkey, qREAL err) {
689  if(IsZero) return;
690  if(Verbose>0) cout << Color_HighLight << " Integrates @ " << now() << RESET << endl;
691 
692  lst lstRE;
693  auto pid = getpid();
694  ostringstream fsofn, sofn, cmd;
695  if(true) {
696  set_precision(NNDigits);
697  sofn << key << ".so";
698  ostringstream garfn;
699  garfn << key << ".ci.gar";
700  archive ar;
701  ifstream in(garfn.str());
702  in >> ar;
703  in.close();
704  auto c = ar.unarchive_ex("c");
705  if(c!=19790923) throw Error("Integrates: *.ci.gar error!");
706  auto gl = ar.unarchive_ex("soLimit");
707  eps_lst = ex_to<lst>(ar.unarchive_ex("eps_lst"));
708  soLimit = ex2int(gl);
709  auto res = ar.unarchive_ex("res");
710  ciResult.clear();
711  for(auto item : ex_to<lst>(res)) ciResult.push_back(ex_to<lst>(item));
712  garfn.clear();
713  garfn.str("");
714  garfn << key;
715  if(pkey != "") garfn << "-" << pkey;
716  garfn << ".las.gar";
717  if(file_exists(garfn.str().c_str())) {
718  archive la_ar;
719  ifstream la_in(garfn.str());
720  la_in >> la_ar;
721  la_in.close();
722  auto la_c = la_ar.unarchive_ex("c");
723  auto la_res = la_ar.unarchive_ex("res");
724  if(la_c!=19790923) throw Error("Integrates: *.ci.gar error!");
725  for(auto item : ex_to<lst>(la_res)) {
726  LambdaMap[item.op(0)] = item.op(1);
727  }
728  }
729 
730  if(true) {
731  garfn.clear();
732  garfn.str("");
733  garfn << key;
734  if(pkey != "") garfn << "-" << pkey;
735  garfn << ".res.gar";
736  if(!file_exists(garfn.str().c_str())) {
737  throw Error("Integrates: File Not Found: " + garfn.str());
738  }
739 
740  archive res_ar;
741  ifstream res_in(garfn.str());
742  res_in >> res_ar;
743  res_in.close();
744  auto res_c = res_ar.unarchive_ex("c");
745  auto relst = res_ar.unarchive_ex("relst");
746  if(res_c!=19790923) throw Error("*.res.gar error with ReIntegrates!");
747  lstRE = ex_to<lst>(relst);
748  }
749  reset_precision();
750  }
751 
752  void* main_module = dlopen(sofn.str().c_str(), RTLD_NOW);
753  if(main_module == nullptr) {
754  main_module = dlopen(("./"+sofn.str()).c_str(), RTLD_NOW);
755  if(main_module == nullptr) {
756  cout << "dlerror(): " << dlerror() << endl;
757  throw Error("Integrates: could not open main module!");
758  }
759  }
760 
761  vector<void*> ex_modules;
762  for(int n=1; true; n++) {
763  ostringstream ex_sofn;
764  ex_sofn << key << "X" << n << ".so";
765  if(file_exists(ex_sofn.str().c_str())) {
766  void* module = dlopen(ex_sofn.str().c_str(), RTLD_NOW);
767  if(module == nullptr) {
768  module = dlopen(("./"+ex_sofn.str()).c_str(), RTLD_NOW);
769  if(module == nullptr) {
770  cout << "dlerror(): " << dlerror() << endl;
771  throw Error("Integrates: could not open ex-module!");
772  }
773  }
774  ex_modules.push_back(module);
775  if(!Debug && key == "") remove(ex_sofn.str().c_str());
776  } else break;
777  }
778 
779  int npara = 0;
780  lst plRepl;
781  for(auto kv : Parameter) {
782  plRepl.append(PL(kv.first)==kv.second);
783  if(kv.first>npara) npara = kv.first;
784  }
785  plRepl.sort();
786  plRepl.unique();
787 
788  int total = ciResult.size(), current = 0;
789  qREAL stot = sqrtq(total*1.Q);
790 
791  ResultError = 0;
792  //----------------------------------------------------------------
793  for(auto &item : ciResult) {
794  current++;
795  auto cmerr = ex2q(VEMaxErr(lstRE.op(current-1)));
796  if(cmerr < err) continue;
797  if(Verbose>10) {
798  char es[64];
799  quadmath_snprintf(es, sizeof es, "%.10QG", cmerr);
800  cout << PRE << "\\--Current Err: " << es << endl;
801  }
802  if(Verbose>0) {
803  cout << "\r \r";
804  cout << PRE << "\\--Integrating [" <<current<<"/"<<total<< "] " << flush;
805  }
806 
807  unsigned int xsize = 0;
808  ex co;
809  exvector xs, fxs;
810  xsize = ex_to<numeric>(item.op(1)).to_int();
811  co = item.op(2).subs(plRepl).subs(iEpsilon==iEpsilonN);
812  if(co.has(WRA(w))) co = ContinuousWRA(co);
813 
814  if(xsize<1) {
815  // { expr, xs.size(), kvf.op(0), -1}
816  ex exint = item.op(0).subs(plRepl).subs(iEpsilon==iEpsilonN);
817  if(exint.has(WRA(w))) exint = ContinuousWRA(exint);
818  ex res = VE(NN(exint),0);
819  if(Verbose>5) {
820  cout << "XDim=" << xsize << endl;
821  cout << Color_HighLight << " IRes = "<< HepLib::SD::VEResult(VESimplify(res)) << RESET << endl;
822  }
823  lstRE.let_op(current-1) = res*co;
824  continue;
825  }
826 
827  //{ idx, xs.size(), kvf.op(0), ft_n }
828  ex intid = item.op(0);
829  ex ftid = item.op(3);
830 
831  if(co.is_zero()) continue;
832  co = collect_ex(co, eps);
833  if(co.is_zero()) continue;
834  if(co.has(PL(w))) throw Error("Integrates: PL found @ " + ex2str(co));
835  qREAL cmax = -1;
836  int reim = 0;
837  if(ReIm==3) reim = 3;
838 
839  for(int si=co.ldegree(eps); si<=co.degree(eps); si++) {
840  auto tmp = co.coeff(eps, si);
841  if(tmp.has(eps)) throw Error("Integrates: eps found @ " + ex2str(tmp));
842  tmp = collect_ex(tmp, ep);
843  for(int i=tmp.ldegree(ep); i<=tmp.degree(ep); i++) {
844  auto ccRes = NN(tmp.coeff(ep, i)).expand();
845  lst css;
846  css.append(ccRes);
847  if(is_a<add>(ccRes)) {
848  for(auto item : ccRes) css.append(item);
849  }
850 
851  exmap eps_map;
852  ex epn = ex(1)/111;
853  for(auto epi : eps_lst) {
854  eps_map[epi.op(0)] = epn;
855  epn = epn / 13;
856  }
857 
858  for(int ci=0; ci<css.nops(); ci++) {
859  auto nt = NN(css.op(ci).subs(epz==1).subs(log(vs)==1).subs(vs==1).subs(nReplacement).subs(lst{CV(w1,w2)==w2}).subs(eps_map),30);
860  if(nt.has(ep)) throw Error("Integrates: ep found @ nt = "+ex2str(nt));
861 
862  lst nt_lst;
863  if(!is_a<numeric>(nt)) {
864  auto cv_lst = collect_lst(nt,[](const ex &e)->bool{return Symbol::has(e);});
865  for(auto nti : cv_lst) {
866  auto nnt = nti.op(0);
867  if(!is_a<numeric>(nnt)) throw Error("Integrates: Not a number with nt = "+ex2str(nnt));
868  nt_lst.append(nnt);
869  }
870  } else nt_lst.append(nt);
871 
872  for(auto nnt : nt_lst) {
873  if(ReIm!=3 && reim!=3) {
874  if(ex_to<numeric>(nnt).imag()==0) {
875  if(reim==2) reim = 3;
876  else reim = 1;
877  } else if(ex_to<numeric>(nnt).real()==0) {
878  if(reim==1) reim = 3;
879  else reim = 2;
880  } else {
881  reim = 3;
882  }
883  }
884  nnt = NN(abs(nnt)); // no PL here, then nReplacement
885 
886  qREAL qnt = ex2q(nnt);
887  if(qnt > cmax) cmax = qnt;
888  }
889  }
890  }
891  }
892  if(cmax<=0) {
893  while(true) {
894  auto co2 = subs(co,lst{exp(w1*log(w2)+w3)==pow(w2,w1)*exp(w3),exp(w1*log(w2))==pow(w2,w1)});
895  if(is_zero(co2-co)) break;
896  co = co2;
897  }
898  if(normal(co).is_zero()) continue;
899  throw Error("Integrates: cmax<=0 with co = "+ex2str(co));
900  }
901  if(reim!=3 && ReIm!=3) {
902  if(ReIm==2) {
903  if(reim==1) reim=2;
904  else if(reim==2) reim=1;
905  }
906  }
907 
908  char d1[20], d2[20];
909  sprintf(d1, "%.6G", (double)(EpsAbs/cmax/stot));
910  sprintf(d2, "%.6G", (double)cmax);
911  if(Verbose>5) cout << "XDim=" << xsize << ", EpsAbs=" << d1 << "/" << d2 << endl;
912  auto las = LambdaMap[ftid];
913  bool hasF = (ftid>0);
914  if(hasF && las.is_zero()) throw Error("Integrates: lambda with the key(ft_n=" + ex2str(ftid) + ") is NOT found!");
915  if(hasF && !is_a<lst>(las)) {
916  if(!is_zero(las-ex(1979))) { // the convention for xPositive or explict real mode
917  throw Error("Integrates: something is wrong with the F-term @ ft_n = "+ex2str(ftid) + ", las=" + ex2str(las));
918  } else {
919  hasF = false;
920  }
921  }
922 
923  IntegratorBase::SDD_Type fpD = nullptr;
924  IntegratorBase::SDQ_Type fpQ = nullptr;
925  IntegratorBase::SDMP_Type fpMP = nullptr;
926  IntegratorBase::FT_Type ftp = nullptr;
927  int idx = ex_to<numeric>(intid).to_int();
928  auto module = main_module;
929  if(idx>=soLimit) module = ex_modules[idx/soLimit-1];
930  ostringstream fname;
931  if(hasF) fname << "C";
932  fname << "SDD_" << idx;
933  fpD = (IntegratorBase::SDD_Type)dlsym(module, fname.str().c_str());
934  if(fpD==NULL) {
935  cout << "dlerror(): " << dlerror() << endl;
936  throw Error("Integrates: fpD==NULL");
937  }
938 
939  fname.clear();
940  fname.str("");
941  if(hasF) fname << "C";
942  fname << "SDQ_" << idx;
943  fpQ = (IntegratorBase::SDQ_Type)dlsym(module, fname.str().c_str());
944  if(fpQ==NULL) {
945  cout << "dlerror(): " << dlerror() << endl;
946  throw Error("Integrates: fpQ==NULL");
947  }
948 
949  fname.clear();
950  fname.str("");
951  if(hasF) fname << "C";
952  fname << "SDMP_" << idx;
953  fpMP = (IntegratorBase::SDMP_Type)dlsym(module, fname.str().c_str());
954 
955  if(is_a<lst>(las)) {
956  fname.clear();
957  fname.str("");
958  fname << "FT_" << idx;
959  ftp = (IntegratorBase::FT_Type)dlsym(module, fname.str().c_str());
960  if(ftp==NULL) {
961  cout << "dlerror(): " << dlerror() << endl;
962  throw Error("Integrates: ftp==NULL.");
963  }
964  }
965 
966  dREAL dlas[las.nops()], dpl[npara];
967  qREAL qlas[las.nops()], qpl[npara];
968  mpREAL mplas[las.nops()], mppl[npara];
969  for(auto kv : Parameter) {
970  qpl[kv.first] = ex2q(kv.second);
971  dpl[kv.first] = (dREAL)qpl[kv.first];
972  mppl[kv.first] = qpl[kv.first];
973  }
974 
975  if(Integrator==NULL) Integrator = new HCubature();
976  Integrator->ReIm = reim;
977  if(MPDigits>0) Integrator->MPDigits = MPDigits;
978  Integrator->IntegrandD = fpD;
979  Integrator->IntegrandQ = fpQ;
980  Integrator->IntegrandMP = fpMP;
981  Integrator->FT = ftp;
982  Integrator->dParameter = dpl;
983  Integrator->dLambda = dlas;
984  Integrator->qParameter = qpl;
985  Integrator->qLambda = qlas;
986  Integrator->mpParameter = mppl;
987  Integrator->mpLambda = mplas;
988  Integrator->XDim = xsize;
989 
990  if(hasF) {
991  qREAL lamax = ex2q(las.op(las.nops()-1));
992  if(lamax > IntLaMax) lamax = IntLaMax;
993 
994  int ctryR = 0, ctry = 0, ctryL = 0;
995  int smin = -1;
996  ex min_err, min_res;
997  size_t min_eval;
998  qREAL log_lamax = log10q(lamax);
999  qREAL log_lamin = log_lamax-1.Q;
1000 
1001  ostringstream las_fn;
1002  las_fn << key;
1003  if(pkey != "") las_fn << "-" << pkey;
1004  las_fn << "-" << current << ".las";
1005  if(use_las && file_exists(las_fn.str().c_str())) {
1006  std::ifstream las_ifs;
1007  las_ifs.open(las_fn.str(), ios::in);
1008  if (!las_ifs) throw Error("Integrates: failed to open *.las file!");
1009  for(int i=0; i<las.nops()-1; i++) {
1010  dREAL la_tmp;
1011  las_ifs >> la_tmp;
1012  qlas[i] = la_tmp;
1013  dlas[i] = (dREAL)qlas[i];
1014  mplas[i] = qlas[i];
1015  }
1016  las_ifs.close();
1017  auto res = Integrator->Integrate();
1018  auto res_tmp = res.subs(VE(w1, w2)==w2);
1019  auto err = real_part(res_tmp);
1020  if(err < imag_part(res_tmp)) err = imag_part(res_tmp);
1021  min_err = err;
1022  min_res = res;
1023  } else {
1024  // ---------------------------------------
1025  while(true) {
1026  smin = -1;
1027  min_err = 0;
1028  ex lastResErr = 0;
1029  bool err_break = false;
1030  for(int s=0; s<=LambdaSplit; s++) {
1031  if(Verbose>10 && s==0) {
1032  if(ctryR>0 || ctry>0 || ctryL>0)
1033  cout << " ------------------------------" << endl;
1034  }
1035  auto log_cla = (log_lamin + s * (log_lamax-log_lamin) / LambdaSplit);
1036  auto cla = powq(10.Q, log_cla);
1037  if(cla < 1E-10) throw Error("NIntegrate: too small lambda.");
1038  for(int i=0; i<las.nops()-1; i++) {
1039  qlas[i] = ex2q(las.op(i)) * cla;
1040  dlas[i] = (dREAL)qlas[i];
1041  mplas[i] = qlas[i];
1042  }
1043 
1044  auto res = Integrator->Integrate(CTryI);
1045  if(Verbose>10) {
1046  cout << "\r \r";
1047  if(res.has(NaN)) cout << " λ=" << (double)cla << "/" << Integrator->NEval << ": " << NaN << endl;
1048  else cout << " λ=" << (double)cla << "/" << Integrator->NEval << ": " << VEResult2(VESimplify(res)) << endl;
1049  }
1050 
1051  if(res.has(NaN) && s==0) continue;
1052  else if(res.has(NaN)) break;
1053  ex res_abs = NN(abs(res.subs(VE(w1,w2)==w1)));
1054  if(lastResErr.is_zero()) lastResErr = res;
1055  auto diff = VESimplify(lastResErr - res);
1056  diff = diff.subs(VE(0,0)==0);
1057  exset ves;
1058  diff.find(VE(w0, w1), ves);
1059  for(auto ve : ves) {
1060  auto ve0 = abs(ve.op(0));
1061  if(ve0>ve.op(1)) {
1062  if(numeric("1.E10")*ve0<res_abs) continue; // avoid fluctuation aroud 0
1063  if(numeric("1.E10")*ve0<q2ex(EpsAbs)) continue; // avoid fluctuation aroud 0
1064  err_break = true;
1065  break;
1066  }
1067  }
1068  if(err_break) {
1069  if(Verbose>10) cout << Color_HighLight << " Error Break ..." << RESET << endl;
1070  break;
1071  }
1072  lastResErr = res;
1073 
1074  auto res_tmp = res.subs(VE(w1, w2)==w2);
1075  auto err = real_part(res_tmp);
1076  if(err < imag_part(res_tmp)) err = imag_part(res_tmp);
1077  if(smin<0 || err < min_err) {
1078  min_err = err;
1079  min_res = res;
1080  min_eval = Integrator->NEval;
1081  smin = s;
1082  }
1083  if(s>0 && min_err < q2ex(EpsAbs/cmax/stot)) {
1084  // s>0 make sure at least 2 λs compatiable
1085  if(Verbose>5) {
1086  cout << Color_HighLight << " λ=" << (double)cla << "/" << min_eval << ": " << HepLib::SD::VEResult(VESimplify(min_res)) << RESET << endl;
1087  }
1088 
1089  smin = -2;
1090  lstRE.let_op(current-1) = co * min_res;
1091  break;
1092  }
1093  if(err > 100 * min_err) break; // s>0 make sure at least 2 λs compatiable
1094  }
1095  if(smin == -2) break;
1096  if(smin == -1) throw Error("Integrates: smin = -1, too small lambda (<1E-10)!");
1097 
1098  if(smin <= 0) {
1099  if((!err_break) && (ctryL >= CTryL || ctryR>0)) break;
1100  log_lamax = log_lamin;
1101  log_lamin -= 1.Q;
1102  if(!err_break) ctryL++;
1103  } else if(smin >= LambdaSplit) {
1104  if(ctryR >= CTryR || ctryL>0) break;
1105  log_lamin = log_lamax;
1106  log_lamax += log10q(CTryRRatio);
1107  ctryR++;
1108  } else {
1109  if(ctry >= CTryM) break;
1110  auto la1 = log_lamin + (smin-1) * (log_lamax-log_lamin) / LambdaSplit;
1111  auto la2 = log_lamin + (smin+1) * (log_lamax-log_lamin) / LambdaSplit;
1112  log_lamin = la1;
1113  log_lamax = la2;
1114  ctry++;
1115  }
1116  }
1117 
1118  if(smin == -2) continue;
1119 
1120  auto log_cla = (log_lamin + smin * (log_lamax-log_lamin) / LambdaSplit);
1121  auto cla = powq(10.Q, log_cla);
1122  if(Verbose>5) cout << Color_HighLight << " Final λ = " << (double)cla << " / " << las.op(las.nops()-1) << RESET << endl;
1123  for(int i=0; i<las.nops()-1; i++) {
1124  qlas[i] = ex2q(las.op(i)) * cla;
1125  dlas[i] = (dREAL)qlas[i];
1126  mplas[i] = qlas[i];
1127  }
1128  // ---------------------------------------
1129  }
1130 
1131  // ---------------------------------------
1132  // try HookeJeeves
1133  // ---------------------------------------
1134  if( use_ErrMin && (min_err > q2ex(1E5 * EpsAbs/cmax/stot)) ) {
1135  ErrMin::lastResErr = min_res;
1136  auto miner = new HookeJeeves();
1137  ErrMin::miner = miner;
1138  ErrMin::Integrator = Integrator;
1139  dREAL oo[las.nops()-1], ip[las.nops()-1];
1140  for(int i=0; i<las.nops()-1; i++) ip[i] = oo[i] = (dREAL)qlas[i];
1141  ErrMin::lambda = oo;
1142  ErrMin::err_max = 1E100;
1143  auto oerrmin = ErrMin::err_min;
1144  ErrMin::err_min = (dREAL)(oerrmin < 0 ? -oerrmin * ex2q(min_err) : oerrmin/cmax);
1145  ErrMin::RunRND = 0;
1146  miner->Minimize(las.nops()-1, ErrMin::IntError, ip);
1147  delete miner;
1148  ErrMin::err_min = oerrmin;
1149  for(int i=0; i<las.nops()-1; i++) {
1150  qlas[i] = ErrMin::lambda[i];
1151  dlas[i] = (dREAL)qlas[i];
1152  mplas[i] = qlas[i];
1153  }
1154 
1155  Integrator->dLambda = dlas;
1156  Integrator->qLambda = qlas; // Integrator->Lambda changed in ErrMin
1157  Integrator->mpLambda = mplas;
1158  if(Verbose>5) {
1159  cout << Color_HighLight << " Final λs: " << RESET;
1160  for(int i=0; i<xsize; i++) {
1161  char buffer[128];
1162  quadmath_snprintf(buffer, sizeof buffer, "%.6QG", qlas[i]);
1163  cout << buffer << " ";
1164  }
1165  cout << endl << " ------------------------------" << endl;
1166  }
1167  }
1168  // ---------------------------------------
1169 
1170  if(save_las) {
1171  std::ofstream las_ofs;
1172  las_ofs.open(las_fn.str(), ios::out);
1173  if (las_ofs) {
1174  for(int i=0; i<las.nops()-1; i++) {
1175  dREAL la_tmp = (dREAL)qlas[i];
1176  las_ofs << la_tmp << " ";
1177  }
1178  las_ofs << endl;
1179  las_ofs.close();
1180  }
1181  }
1182  }
1183 
1184  auto res = Integrator->Integrate();
1185  if(Verbose>5) {
1186  cout << Color_HighLight << " IRes = "<< HepLib::SD::VEResult(VESimplify(res)) << RESET << endl;
1187  }
1188  if(res.has(NaN)) {
1189  ResultError = NaN;
1190  lstRE.let_op(current-1) = NaN;
1191  } else {
1192  lstRE.let_op(current-1) = co * res;
1193  }
1194  }
1195  //----------------------------------------------------------------
1196 
1197  if(use_dlclose) {
1198  for(auto module : ex_modules) dlclose(module);
1199  dlclose(main_module);
1200  }
1201 
1202  if(!ResultError.is_equal(NaN)) {
1203  ResultError = 0;
1204  for(auto item : lstRE) ResultError += item;
1205  ResultError = VESimplify(ResultError);
1206  }
1207 
1208  if(key != "") {
1209  ostringstream garfn;
1210  garfn << key;
1211  if(pkey != "") garfn << "-" << pkey;
1212  garfn << ".res.gar";
1213  archive ar;
1214  ar.archive_ex(ResultError, "res");
1215  ar.archive_ex(lstRE, "relst");
1216  ar.archive_ex(19790923, "c");
1217  ofstream out(garfn.str());
1218  out << ar;
1219  out.close();
1220  }
1221  }
1222 
1223 }
RESET
#define RESET
Definition: BASIC.h:79
SD.h
SecDec header file.
HepLib::Error
class used to wrap error message
Definition: BASIC.h:242
HepLib::SD::ErrMin::lastResErr
static ex lastResErr
Definition: SD.h:380
HepLib::SD::ErrMin::Integrator
static IntegratorBase * Integrator
Definition: SD.h:371
HepLib::SD::ErrMin::RunRND
static size_t RunRND
Definition: SD.h:376
HepLib::SD::ErrMin::err_max
static dREAL err_max
Definition: SD.h:373
HepLib::SD::ErrMin::miner
static MinimizeBase * miner
Definition: SD.h:377
HepLib::SD::ErrMin::IntError
static dREAL IntError(int nvars, dREAL *las, dREAL *n1, dREAL *n2)
Definition: ErrMin.cpp:22
HepLib::SD::ErrMin::lambda
static dREAL * lambda
Definition: SD.h:378
HepLib::SD::ErrMin::err_min
static dREAL err_min
Definition: SD.h:374
HepLib::SD::HCubature
numerical integrator using HCubature
Definition: SD.h:166
HepLib::SD::HookeJeeves
class to minimize a function using HookeJeeves
Definition: SD.h:280
HepLib::SD::IntegratorBase::IntegrandMP
SDMP_Type IntegrandMP
Definition: SD.h:144
HepLib::SD::IntegratorBase::dParameter
const dREAL * dParameter
Definition: SD.h:146
HepLib::SD::IntegratorBase::IntegrandD
SDD_Type IntegrandD
Definition: SD.h:142
HepLib::SD::IntegratorBase::FT_Type
qREAL(* FT_Type)(const qREAL xx[], const qREAL pl[])
Definition: SD.h:138
HepLib::SD::IntegratorBase::SDD_Type
int(* SDD_Type)(const unsigned int xn, const dREAL x[], const unsigned int yn, dREAL y[], const dREAL pl[], const dREAL las[])
Definition: SD.h:135
HepLib::SD::IntegratorBase::dLambda
const dREAL * dLambda
Definition: SD.h:145
HepLib::SD::IntegratorBase::IntegrandQ
SDQ_Type IntegrandQ
Definition: SD.h:143
HepLib::SD::IntegratorBase::mpLambda
const mpREAL * mpLambda
Definition: SD.h:149
HepLib::SD::IntegratorBase::SDQ_Type
int(* SDQ_Type)(const unsigned int xn, const qREAL x[], const unsigned int yn, qREAL y[], const qREAL pl[], const qREAL las[])
Definition: SD.h:136
HepLib::SD::IntegratorBase::qParameter
const qREAL * qParameter
Definition: SD.h:148
HepLib::SD::IntegratorBase::Integrate
virtual ex Integrate(size_t n=0)=0
HepLib::SD::IntegratorBase::SDMP_Type
int(* SDMP_Type)(const unsigned int xn, const mpREAL x[], const unsigned int yn, mpREAL y[], const mpREAL pl[], const mpREAL las[])
Definition: SD.h:137
HepLib::SD::IntegratorBase::mpParameter
const mpREAL * mpParameter
Definition: SD.h:150
HepLib::SD::IntegratorBase::qLambda
const qREAL * qLambda
Definition: SD.h:147
HepLib::SD::SecDec::LambdaSplit
size_t LambdaSplit
Definition: SD.h:452
HepLib::SD::SecDec::nReplacement
exmap nReplacement
Definition: SD.h:425
HepLib::SD::SecDec::MPDigits
int MPDigits
Definition: SD.h:441
HepLib::SD::SecDec::use_ErrMin
bool use_ErrMin
Definition: SD.h:435
HepLib::SD::SecDec::CTryRRatio
dREAL CTryRRatio
Definition: SD.h:458
HepLib::SD::SecDec::IntLaMax
qREAL IntLaMax
Definition: SD.h:453
HepLib::SD::SecDec::use_las
bool use_las
Definition: SD.h:436
HepLib::SD::SecDec::Parameter
map< int, numeric > Parameter
Definition: SD.h:444
HepLib::SD::SecDec::EpsAbs
qREAL EpsAbs
Definition: SD.h:461
HepLib::SD::SecDec::save_las
bool save_las
Definition: SD.h:437
HepLib::SD::SecDec::Integrator
IntegratorBase * Integrator
Definition: SD.h:430
HepLib::SD::SecDec::use_dlclose
static bool use_dlclose
Definition: SD.h:416
HepLib::SD::SecDec::ResultError
ex ResultError
Definition: SD.h:432
HepLib::SD::SecDec::IsZero
bool IsZero
Definition: SD.h:433
HepLib::SD::SecDec::Integrates
void Integrates(const string &key="", const string &pkey="", int kid=0)
Contours, note that here we need to provide the specific Parameter.
Definition: Integrates.cpp:113
HepLib::SD::SecDec::ContinuousWRA
static ex ContinuousWRA(ex expr_in, int nc=15)
ContinuousWRA, note that here we need to provide the specific Parameter.
Definition: Integrates.cpp:18
HepLib::SD::SecDec::ReIntegrates
void ReIntegrates(const string &key, const string &pkey, qREAL err)
Contours, note that here we need to provide the specific Parameter.
Definition: Integrates.cpp:688
HepLib::SD::SecDec::CTryI
size_t CTryI
Definition: SD.h:457
HepLib::Symbol::has
static bool has(const ex &e)
HepLib::SD
namespace for Numerical integration with Sector Decomposition method
Definition: AsyMB.cpp:10
HepLib::SD::xyz_pow_simplify
ex xyz_pow_simplify(const ex expr_in)
Definition: Helpers.cpp:336
HepLib::SD::exp_simplify
ex exp_simplify(const ex expr_in)
Definition: Helpers.cpp:302
HepLib::SD::qREAL
__float128 qREAL
Definition: SD.h:123
HepLib::SD::VEResult
ex VEResult(ex expr)
Definition: Helpers.cpp:123
HepLib::SD::dREAL
long double dREAL
Definition: SD.h:121
HepLib::SD::VESimplify
ex VESimplify(ex expr)
Definition: Helpers.cpp:92
HepLib::SD::VEMaxErr
ex VEMaxErr(ex expr)
Definition: Helpers.cpp:131
HepLib::SD::mpREAL
mpfr::mpreal mpREAL
Definition: SD.h:125
HepLib::SD::VEResult2
ex VEResult2(ex expr)
Definition: Helpers.cpp:127
HepLib::Color_HighLight
const char * Color_HighLight
Definition: BASIC.cpp:248
HepLib::ex2q
__float128 ex2q(ex num)
ex of numeric to __float128
Definition: BASIC.cpp:879
HepLib::iEpsilon
const iSymbol iEpsilon
HepLib::reset_precision
void reset_precision()
Definition: BASIC.cpp:2313
HepLib::q2ex
ex q2ex(__float128 num)
__float128 to ex
Definition: BASIC.cpp:867
HepLib::NNDigits
int NNDigits
Definition: Init.cpp:155
HepLib::NN
ex NN(ex expr, int digits)
the nuerical evaluation
Definition: BASIC.cpp:1278
HepLib::w0
ex w0
Definition: BASIC.h:499
HepLib::ep
const Symbol ep
HepLib::collect_ex
ex collect_ex(ex const &expr_in, std::function< bool(const ex &)> has_func, int opt)
the collect function like Mathematica
Definition: BASIC.cpp:1202
HepLib::eps
const Symbol eps
HepLib::file_exists
bool file_exists(string fn)
Definition: BASIC.h:289
HepLib::now
string now(bool use_date)
date/time string
Definition: BASIC.cpp:525
HepLib::Debug
bool Debug
Definition: Init.cpp:141
HepLib::w
ex w
Definition: Init.cpp:90
HepLib::vs
const Symbol vs
HepLib::iEpsilonN
const ex iEpsilonN
Definition: Init.cpp:138
HepLib::epz
const Symbol epz
HepLib::collect_lst
lst collect_lst(ex const &expr_in, std::function< bool(const ex &)> has_func, int opt)
the collect function like Mathematica, reture the lst { {c1,v1}, {c2,v2}, ... }
Definition: BASIC.cpp:1222
HepLib::Verbose
int Verbose
Definition: Init.cpp:139
HepLib::ex2str
string ex2str(const ex &expr)
convert ex to output string, the defalut printer format will be used
Definition: BASIC.cpp:715
HepLib::w1
ex w1
Definition: BASIC.h:499
HepLib::w3
ex w3
Definition: BASIC.h:499
HepLib::set_precision
void set_precision(long prec, bool push)
Definition: BASIC.cpp:2304
HepLib::w2
ex w2
Definition: BASIC.h:499
HepLib::PRE
string PRE
Definition: Init.cpp:140
HepLib::ex2int
int ex2int(ex num)
ex to integer
Definition: BASIC.cpp:893
HepLib::NaN
const Symbol NaN
HepLib::subs
ex subs(const ex &e, init_list sl)
Definition: BASIC.h:51
Generated by doxygen 1.9.1