HepLib
Loading...
Searching...
No Matches
Integrates.cpp
Go to the documentation of this file.
1
6#include "SD.h"
7#include <math.h>
8#include <cmath>
9
10namespace 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::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::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::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::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::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:2323
HepLib::q2ex
ex q2ex(__float128 num)
__float128 to ex
Definition BASIC.cpp:867
HepLib::NNDigits
int NNDigits
Definition Init.cpp:158
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:144
HepLib::w
ex w
Definition Init.cpp:93
HepLib::vs
const Symbol vs
HepLib::iEpsilonN
const ex iEpsilonN
Definition Init.cpp:141
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:142
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:2314
HepLib::w2
ex w2
Definition BASIC.h:499
HepLib::PRE
string PRE
Definition Init.cpp:143
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.8