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[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[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[kid-1] = NaN;
645 else lstRE.append(NaN);
646 break;
647 } else {
648 if(kid>0) {
649 lstRE[kid-1] = co * res;
650 break;
651 } else lstRE.append(co * res);
652 }
653 }
654 if(Verbose>0 && Verbose<=5) cout << endl;
655 //----------------------------------------------------------------
656
657 if(use_dlclose) {
658 for(auto module : ex_modules) dlclose(module);
659 dlclose(main_module);
660 }
661
662 if(!ResultError.is_equal(NaN)) {
663 ResultError = 0;
664 for(auto item : lstRE) ResultError += item;
665 ResultError = VESimplify(ResultError);
666 }
667
668 if(key != "") {
669 ostringstream garfn;
670 garfn << key;
671 if(pkey != "") garfn << "-" << pkey;
672 garfn << ".res.gar";
673 archive ar;
674 ar.archive_ex(ResultError, "res");
675 ar.archive_ex(lstRE, "relst");
676 ar.archive_ex(19790923, "c");
677 ofstream out(garfn.str());
678 out << ar;
679 out.close();
680 }
681 }
682
689 void SecDec::ReIntegrates(const string & key, const string & pkey, qREAL err) {
690 if(IsZero) return;
691 if(Verbose>0) cout << Color_HighLight << " Integrates @ " << now() << RESET << endl;
692
693 lst lstRE;
694 auto pid = getpid();
695 ostringstream fsofn, sofn, cmd;
696 if(true) {
697 set_precision(NNDigits);
698 sofn << key << ".so";
699 ostringstream garfn;
700 garfn << key << ".ci.gar";
701 archive ar;
702 ifstream in(garfn.str());
703 in >> ar;
704 in.close();
705 auto c = ar.unarchive_ex("c");
706 if(c!=19790923) throw Error("Integrates: *.ci.gar error!");
707 auto gl = ar.unarchive_ex("soLimit");
708 eps_lst = ex_to<lst>(ar.unarchive_ex("eps_lst"));
709 soLimit = ex2int(gl);
710 auto res = ar.unarchive_ex("res");
711 ciResult.clear();
712 for(auto item : ex_to<lst>(res)) ciResult.push_back(ex_to<lst>(item));
713 garfn.clear();
714 garfn.str("");
715 garfn << key;
716 if(pkey != "") garfn << "-" << pkey;
717 garfn << ".las.gar";
718 if(file_exists(garfn.str().c_str())) {
719 archive la_ar;
720 ifstream la_in(garfn.str());
721 la_in >> la_ar;
722 la_in.close();
723 auto la_c = la_ar.unarchive_ex("c");
724 auto la_res = la_ar.unarchive_ex("res");
725 if(la_c!=19790923) throw Error("Integrates: *.ci.gar error!");
726 for(auto item : ex_to<lst>(la_res)) {
727 LambdaMap[item.op(0)] = item.op(1);
728 }
729 }
730
731 if(true) {
732 garfn.clear();
733 garfn.str("");
734 garfn << key;
735 if(pkey != "") garfn << "-" << pkey;
736 garfn << ".res.gar";
737 if(!file_exists(garfn.str().c_str())) {
738 throw Error("Integrates: File Not Found: " + garfn.str());
739 }
740
741 archive res_ar;
742 ifstream res_in(garfn.str());
743 res_in >> res_ar;
744 res_in.close();
745 auto res_c = res_ar.unarchive_ex("c");
746 auto relst = res_ar.unarchive_ex("relst");
747 if(res_c!=19790923) throw Error("*.res.gar error with ReIntegrates!");
748 lstRE = ex_to<lst>(relst);
749 }
750 reset_precision();
751 }
752
753 void* main_module = dlopen(sofn.str().c_str(), RTLD_NOW);
754 if(main_module == nullptr) {
755 main_module = dlopen(("./"+sofn.str()).c_str(), RTLD_NOW);
756 if(main_module == nullptr) {
757 cout << "dlerror(): " << dlerror() << endl;
758 throw Error("Integrates: could not open main module!");
759 }
760 }
761
762 vector<void*> ex_modules;
763 for(int n=1; true; n++) {
764 ostringstream ex_sofn;
765 ex_sofn << key << "X" << n << ".so";
766 if(file_exists(ex_sofn.str().c_str())) {
767 void* module = dlopen(ex_sofn.str().c_str(), RTLD_NOW);
768 if(module == nullptr) {
769 module = dlopen(("./"+ex_sofn.str()).c_str(), RTLD_NOW);
770 if(module == nullptr) {
771 cout << "dlerror(): " << dlerror() << endl;
772 throw Error("Integrates: could not open ex-module!");
773 }
774 }
775 ex_modules.push_back(module);
776 if(!Debug && key == "") remove(ex_sofn.str().c_str());
777 } else break;
778 }
779
780 int npara = 0;
781 lst plRepl;
782 for(auto kv : Parameter) {
783 plRepl.append(PL(kv.first)==kv.second);
784 if(kv.first>npara) npara = kv.first;
785 }
786 plRepl.sort();
787 plRepl.unique();
788
789 int total = ciResult.size(), current = 0;
790 qREAL stot = sqrtq(total*1.Q);
791
792 ResultError = 0;
793 //----------------------------------------------------------------
794 for(auto &item : ciResult) {
795 current++;
796 auto cmerr = ex2q(VEMaxErr(lstRE.op(current-1)));
797 if(cmerr < err) continue;
798 if(Verbose>10) {
799 char es[64];
800 quadmath_snprintf(es, sizeof es, "%.10QG", cmerr);
801 cout << PRE << "\\--Current Err: " << es << endl;
802 }
803 if(Verbose>0) {
804 cout << "\r \r";
805 cout << PRE << "\\--Integrating [" <<current<<"/"<<total<< "] " << flush;
806 }
807
808 unsigned int xsize = 0;
809 ex co;
810 exvector xs, fxs;
811 xsize = ex_to<numeric>(item.op(1)).to_int();
812 co = item.op(2).subs(plRepl).subs(iEpsilon==iEpsilonN);
813 if(co.has(WRA(w))) co = ContinuousWRA(co);
814
815 if(xsize<1) {
816 // { expr, xs.size(), kvf.op(0), -1}
817 ex exint = item.op(0).subs(plRepl).subs(iEpsilon==iEpsilonN);
818 if(exint.has(WRA(w))) exint = ContinuousWRA(exint);
819 ex res = VE(NN(exint),0);
820 if(Verbose>5) {
821 cout << "XDim=" << xsize << endl;
822 cout << Color_HighLight << " IRes = "<< HepLib::SD::VEResult(VESimplify(res)) << RESET << endl;
823 }
824 lstRE[current-1] = res*co;
825 continue;
826 }
827
828 //{ idx, xs.size(), kvf.op(0), ft_n }
829 ex intid = item.op(0);
830 ex ftid = item.op(3);
831
832 if(co.is_zero()) continue;
833 co = collect_ex(co, eps);
834 if(co.is_zero()) continue;
835 if(co.has(PL(w))) throw Error("Integrates: PL found @ " + ex2str(co));
836 qREAL cmax = -1;
837 int reim = 0;
838 if(ReIm==3) reim = 3;
839
840 for(int si=co.ldegree(eps); si<=co.degree(eps); si++) {
841 auto tmp = co.coeff(eps, si);
842 if(tmp.has(eps)) throw Error("Integrates: eps found @ " + ex2str(tmp));
843 tmp = collect_ex(tmp, ep);
844 for(int i=tmp.ldegree(ep); i<=tmp.degree(ep); i++) {
845 auto ccRes = NN(tmp.coeff(ep, i)).expand();
846 lst css;
847 css.append(ccRes);
848 if(is_a<add>(ccRes)) {
849 for(auto item : ccRes) css.append(item);
850 }
851
852 exmap eps_map;
853 ex epn = ex(1)/111;
854 for(auto epi : eps_lst) {
855 eps_map[epi.op(0)] = epn;
856 epn = epn / 13;
857 }
858
859 for(int ci=0; ci<css.nops(); ci++) {
860 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);
861 if(nt.has(ep)) throw Error("Integrates: ep found @ nt = "+ex2str(nt));
862
863 lst nt_lst;
864 if(!is_a<numeric>(nt)) {
865 auto cv_lst = collect_lst(nt,[](const ex &e)->bool{return Symbol::has(e);});
866 for(auto nti : cv_lst) {
867 auto nnt = nti.op(0);
868 if(!is_a<numeric>(nnt)) throw Error("Integrates: Not a number with nt = "+ex2str(nnt));
869 nt_lst.append(nnt);
870 }
871 } else nt_lst.append(nt);
872
873 for(auto nnt : nt_lst) {
874 if(ReIm!=3 && reim!=3) {
875 if(ex_to<numeric>(nnt).imag()==0) {
876 if(reim==2) reim = 3;
877 else reim = 1;
878 } else if(ex_to<numeric>(nnt).real()==0) {
879 if(reim==1) reim = 3;
880 else reim = 2;
881 } else {
882 reim = 3;
883 }
884 }
885 nnt = NN(abs(nnt)); // no PL here, then nReplacement
886
887 qREAL qnt = ex2q(nnt);
888 if(qnt > cmax) cmax = qnt;
889 }
890 }
891 }
892 }
893 if(cmax<=0) {
894 while(true) {
895 auto co2 = subs(co,lst{exp(w1*log(w2)+w3)==pow(w2,w1)*exp(w3),exp(w1*log(w2))==pow(w2,w1)});
896 if(is_zero(co2-co)) break;
897 co = co2;
898 }
899 if(normal(co).is_zero()) continue;
900 throw Error("Integrates: cmax<=0 with co = "+ex2str(co));
901 }
902 if(reim!=3 && ReIm!=3) {
903 if(ReIm==2) {
904 if(reim==1) reim=2;
905 else if(reim==2) reim=1;
906 }
907 }
908
909 char d1[20], d2[20];
910 sprintf(d1, "%.6G", (double)(EpsAbs/cmax/stot));
911 sprintf(d2, "%.6G", (double)cmax);
912 if(Verbose>5) cout << "XDim=" << xsize << ", EpsAbs=" << d1 << "/" << d2 << endl;
913 auto las = LambdaMap[ftid];
914 bool hasF = (ftid>0);
915 if(hasF && las.is_zero()) throw Error("Integrates: lambda with the key(ft_n=" + ex2str(ftid) + ") is NOT found!");
916 if(hasF && !is_a<lst>(las)) {
917 if(!is_zero(las-ex(1979))) { // the convention for xPositive or explict real mode
918 throw Error("Integrates: something is wrong with the F-term @ ft_n = "+ex2str(ftid) + ", las=" + ex2str(las));
919 } else {
920 hasF = false;
921 }
922 }
923
924 IntegratorBase::SDD_Type fpD = nullptr;
925 IntegratorBase::SDQ_Type fpQ = nullptr;
926 IntegratorBase::SDMP_Type fpMP = nullptr;
927 IntegratorBase::FT_Type ftp = nullptr;
928 int idx = ex_to<numeric>(intid).to_int();
929 auto module = main_module;
930 if(idx>=soLimit) module = ex_modules[idx/soLimit-1];
931 ostringstream fname;
932 if(hasF) fname << "C";
933 fname << "SDD_" << idx;
934 fpD = (IntegratorBase::SDD_Type)dlsym(module, fname.str().c_str());
935 if(fpD==NULL) {
936 cout << "dlerror(): " << dlerror() << endl;
937 throw Error("Integrates: fpD==NULL");
938 }
939
940 fname.clear();
941 fname.str("");
942 if(hasF) fname << "C";
943 fname << "SDQ_" << idx;
944 fpQ = (IntegratorBase::SDQ_Type)dlsym(module, fname.str().c_str());
945 if(fpQ==NULL) {
946 cout << "dlerror(): " << dlerror() << endl;
947 throw Error("Integrates: fpQ==NULL");
948 }
949
950 fname.clear();
951 fname.str("");
952 if(hasF) fname << "C";
953 fname << "SDMP_" << idx;
954 fpMP = (IntegratorBase::SDMP_Type)dlsym(module, fname.str().c_str());
955
956 if(is_a<lst>(las)) {
957 fname.clear();
958 fname.str("");
959 fname << "FT_" << idx;
960 ftp = (IntegratorBase::FT_Type)dlsym(module, fname.str().c_str());
961 if(ftp==NULL) {
962 cout << "dlerror(): " << dlerror() << endl;
963 throw Error("Integrates: ftp==NULL.");
964 }
965 }
966
967 dREAL dlas[las.nops()], dpl[npara];
968 qREAL qlas[las.nops()], qpl[npara];
969 mpREAL mplas[las.nops()], mppl[npara];
970 for(auto kv : Parameter) {
971 qpl[kv.first] = ex2q(kv.second);
972 dpl[kv.first] = (dREAL)qpl[kv.first];
973 mppl[kv.first] = qpl[kv.first];
974 }
975
976 if(Integrator==NULL) Integrator = new HCubature();
977 Integrator->ReIm = reim;
978 if(MPDigits>0) Integrator->MPDigits = MPDigits;
979 Integrator->IntegrandD = fpD;
980 Integrator->IntegrandQ = fpQ;
981 Integrator->IntegrandMP = fpMP;
982 Integrator->FT = ftp;
983 Integrator->dParameter = dpl;
984 Integrator->dLambda = dlas;
985 Integrator->qParameter = qpl;
986 Integrator->qLambda = qlas;
987 Integrator->mpParameter = mppl;
988 Integrator->mpLambda = mplas;
989 Integrator->XDim = xsize;
990
991 if(hasF) {
992 qREAL lamax = ex2q(las.op(las.nops()-1));
993 if(lamax > IntLaMax) lamax = IntLaMax;
994
995 int ctryR = 0, ctry = 0, ctryL = 0;
996 int smin = -1;
997 ex min_err, min_res;
998 size_t min_eval;
999 qREAL log_lamax = log10q(lamax);
1000 qREAL log_lamin = log_lamax-1.Q;
1001
1002 ostringstream las_fn;
1003 las_fn << key;
1004 if(pkey != "") las_fn << "-" << pkey;
1005 las_fn << "-" << current << ".las";
1006 if(use_las && file_exists(las_fn.str().c_str())) {
1007 std::ifstream las_ifs;
1008 las_ifs.open(las_fn.str(), ios::in);
1009 if (!las_ifs) throw Error("Integrates: failed to open *.las file!");
1010 for(int i=0; i<las.nops()-1; i++) {
1011 dREAL la_tmp;
1012 las_ifs >> la_tmp;
1013 qlas[i] = la_tmp;
1014 dlas[i] = (dREAL)qlas[i];
1015 mplas[i] = qlas[i];
1016 }
1017 las_ifs.close();
1018 auto res = Integrator->Integrate();
1019 auto res_tmp = res.subs(VE(w1, w2)==w2);
1020 auto err = real_part(res_tmp);
1021 if(err < imag_part(res_tmp)) err = imag_part(res_tmp);
1022 min_err = err;
1023 min_res = res;
1024 } else {
1025 // ---------------------------------------
1026 while(true) {
1027 smin = -1;
1028 min_err = 0;
1029 ex lastResErr = 0;
1030 bool err_break = false;
1031 for(int s=0; s<=LambdaSplit; s++) {
1032 if(Verbose>10 && s==0) {
1033 if(ctryR>0 || ctry>0 || ctryL>0)
1034 cout << " ------------------------------" << endl;
1035 }
1036 auto log_cla = (log_lamin + s * (log_lamax-log_lamin) / LambdaSplit);
1037 auto cla = powq(10.Q, log_cla);
1038 if(cla < 1E-10) throw Error("NIntegrate: too small lambda.");
1039 for(int i=0; i<las.nops()-1; i++) {
1040 qlas[i] = ex2q(las.op(i)) * cla;
1041 dlas[i] = (dREAL)qlas[i];
1042 mplas[i] = qlas[i];
1043 }
1044
1045 auto res = Integrator->Integrate(CTryI);
1046 if(Verbose>10) {
1047 cout << "\r \r";
1048 if(res.has(NaN)) cout << " λ=" << (double)cla << "/" << Integrator->NEval << ": " << NaN << endl;
1049 else cout << " λ=" << (double)cla << "/" << Integrator->NEval << ": " << VEResult2(VESimplify(res)) << endl;
1050 }
1051
1052 if(res.has(NaN) && s==0) continue;
1053 else if(res.has(NaN)) break;
1054 ex res_abs = NN(abs(res.subs(VE(w1,w2)==w1)));
1055 if(lastResErr.is_zero()) lastResErr = res;
1056 auto diff = VESimplify(lastResErr - res);
1057 diff = diff.subs(VE(0,0)==0);
1058 exset ves;
1059 diff.find(VE(w0, w1), ves);
1060 for(auto ve : ves) {
1061 auto ve0 = abs(ve.op(0));
1062 if(ve0>ve.op(1)) {
1063 if(numeric("1.E10")*ve0<res_abs) continue; // avoid fluctuation aroud 0
1064 if(numeric("1.E10")*ve0<q2ex(EpsAbs)) continue; // avoid fluctuation aroud 0
1065 err_break = true;
1066 break;
1067 }
1068 }
1069 if(err_break) {
1070 if(Verbose>10) cout << Color_HighLight << " Error Break ..." << RESET << endl;
1071 break;
1072 }
1073 lastResErr = res;
1074
1075 auto res_tmp = res.subs(VE(w1, w2)==w2);
1076 auto err = real_part(res_tmp);
1077 if(err < imag_part(res_tmp)) err = imag_part(res_tmp);
1078 if(smin<0 || err < min_err) {
1079 min_err = err;
1080 min_res = res;
1081 min_eval = Integrator->NEval;
1082 smin = s;
1083 }
1084 if(s>0 && min_err < q2ex(EpsAbs/cmax/stot)) {
1085 // s>0 make sure at least 2 λs compatiable
1086 if(Verbose>5) {
1087 cout << Color_HighLight << " λ=" << (double)cla << "/" << min_eval << ": " << HepLib::SD::VEResult(VESimplify(min_res)) << RESET << endl;
1088 }
1089
1090 smin = -2;
1091 lstRE[current-1] = co * min_res;
1092 break;
1093 }
1094 if(err > 100 * min_err) break; // s>0 make sure at least 2 λs compatiable
1095 }
1096 if(smin == -2) break;
1097 if(smin == -1) throw Error("Integrates: smin = -1, too small lambda (<1E-10)!");
1098
1099 if(smin <= 0) {
1100 if((!err_break) && (ctryL >= CTryL || ctryR>0)) break;
1101 log_lamax = log_lamin;
1102 log_lamin -= 1.Q;
1103 if(!err_break) ctryL++;
1104 } else if(smin >= LambdaSplit) {
1105 if(ctryR >= CTryR || ctryL>0) break;
1106 log_lamin = log_lamax;
1107 log_lamax += log10q(CTryRRatio);
1108 ctryR++;
1109 } else {
1110 if(ctry >= CTryM) break;
1111 auto la1 = log_lamin + (smin-1) * (log_lamax-log_lamin) / LambdaSplit;
1112 auto la2 = log_lamin + (smin+1) * (log_lamax-log_lamin) / LambdaSplit;
1113 log_lamin = la1;
1114 log_lamax = la2;
1115 ctry++;
1116 }
1117 }
1118
1119 if(smin == -2) continue;
1120
1121 auto log_cla = (log_lamin + smin * (log_lamax-log_lamin) / LambdaSplit);
1122 auto cla = powq(10.Q, log_cla);
1123 if(Verbose>5) cout << Color_HighLight << " Final λ = " << (double)cla << " / " << las.op(las.nops()-1) << RESET << endl;
1124 for(int i=0; i<las.nops()-1; i++) {
1125 qlas[i] = ex2q(las.op(i)) * cla;
1126 dlas[i] = (dREAL)qlas[i];
1127 mplas[i] = qlas[i];
1128 }
1129 // ---------------------------------------
1130 }
1131
1132 // ---------------------------------------
1133 // try HookeJeeves
1134 // ---------------------------------------
1135 if( use_ErrMin && (min_err > q2ex(1E5 * EpsAbs/cmax/stot)) ) {
1136 ErrMin::lastResErr = min_res;
1137 auto miner = new HookeJeeves();
1138 ErrMin::miner = miner;
1139 ErrMin::Integrator = Integrator;
1140 dREAL oo[las.nops()-1], ip[las.nops()-1];
1141 for(int i=0; i<las.nops()-1; i++) ip[i] = oo[i] = (dREAL)qlas[i];
1142 ErrMin::lambda = oo;
1143 ErrMin::err_max = 1E100;
1144 auto oerrmin = ErrMin::err_min;
1145 ErrMin::err_min = (dREAL)(oerrmin < 0 ? -oerrmin * ex2q(min_err) : oerrmin/cmax);
1146 ErrMin::RunRND = 0;
1147 miner->Minimize(las.nops()-1, ErrMin::IntError, ip);
1148 delete miner;
1149 ErrMin::err_min = oerrmin;
1150 for(int i=0; i<las.nops()-1; i++) {
1151 qlas[i] = ErrMin::lambda[i];
1152 dlas[i] = (dREAL)qlas[i];
1153 mplas[i] = qlas[i];
1154 }
1155
1156 Integrator->dLambda = dlas;
1157 Integrator->qLambda = qlas; // Integrator->Lambda changed in ErrMin
1158 Integrator->mpLambda = mplas;
1159 if(Verbose>5) {
1160 cout << Color_HighLight << " Final λs: " << RESET;
1161 for(int i=0; i<xsize; i++) {
1162 char buffer[128];
1163 quadmath_snprintf(buffer, sizeof buffer, "%.6QG", qlas[i]);
1164 cout << buffer << " ";
1165 }
1166 cout << endl << " ------------------------------" << endl;
1167 }
1168 }
1169 // ---------------------------------------
1170
1171 if(save_las) {
1172 std::ofstream las_ofs;
1173 las_ofs.open(las_fn.str(), ios::out);
1174 if (las_ofs) {
1175 for(int i=0; i<las.nops()-1; i++) {
1176 dREAL la_tmp = (dREAL)qlas[i];
1177 las_ofs << la_tmp << " ";
1178 }
1179 las_ofs << endl;
1180 las_ofs.close();
1181 }
1182 }
1183 }
1184
1185 auto res = Integrator->Integrate();
1186 if(Verbose>5) {
1187 cout << Color_HighLight << " IRes = "<< HepLib::SD::VEResult(VESimplify(res)) << RESET << endl;
1188 }
1189 if(res.has(NaN)) {
1190 ResultError = NaN;
1191 lstRE[current-1] = NaN;
1192 } else {
1193 lstRE[current-1] = co * res;
1194 }
1195 }
1196 if(Verbose>0 && Verbose<=5) cout << endl;
1197
1198 //----------------------------------------------------------------
1199
1200 if(use_dlclose) {
1201 for(auto module : ex_modules) dlclose(module);
1202 dlclose(main_module);
1203 }
1204
1205 if(!ResultError.is_equal(NaN)) {
1206 ResultError = 0;
1207 for(auto item : lstRE) ResultError += item;
1208 ResultError = VESimplify(ResultError);
1209 }
1210
1211 if(key != "") {
1212 ostringstream garfn;
1213 garfn << key;
1214 if(pkey != "") garfn << "-" << pkey;
1215 garfn << ".res.gar";
1216 archive ar;
1217 ar.archive_ex(ResultError, "res");
1218 ar.archive_ex(lstRE, "relst");
1219 ar.archive_ex(19790923, "c");
1220 ofstream out(garfn.str());
1221 out << ar;
1222 out.close();
1223 }
1224 }
1225
1226}
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:689
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