HepLib
SecDec.cpp
Go to the documentation of this file.
1 
6 #include "SD.h"
7 #include <math.h>
8 #include <cmath>
9 #include<algorithm>
10 
11 namespace HepLib::SD {
12 
18  ex SecDecBase::XMonomials(const ex & expr_in) {
19  auto expr = expr_in;
20  auto xs = get_x_from(expr);
21  auto nx = xs.size();
22  if(nx<1) return expr;
23  if(!is_polynomial(expr, vec2lst(xs))) return expr;
24  auto cvs = collect_lst(expr,x(w));
25  vector<exvector> degs_vec;
26  for(auto cv : cvs) {
27  if(is_zero(cv.op(0))) continue;
28  exvector degs;
29  for(auto xi : xs) degs.push_back(cv.op(1).degree(xi));
30  degs_vec.push_back(degs);
31  }
32 
33  map<int,bool> flag;
34  auto nn = degs_vec.size();
35 
36  for(int i=0; i<nn; i++) flag[i] = true;
37  try {
38  matrix mat(nn,nx);
39  for(int ri=0; ri<nn; ri++) for(int ci=0; ci<nx; ci++) mat(ri,ci) = degs_vec[ri][ci];
40  auto vvi = SecDecG::RunQHull(mat);
41  for(auto vi : vvi) for(auto ii : vi) flag[ii] = false;
42  } catch(...) {
43  for(int i=0; i<nn; i++) flag[i] = false;
44  }
45 
46  for(int i=0; i<nn; i++) {
47  if(flag[i]) continue;
48  exvector & ti = degs_vec[i];
49  for(int j=i+1; j<nn; j++) {
50  if(flag[j]) continue;
51  exvector & tj = degs_vec[j];
52  int comp = 0;
53  for(int c=0; c<nx; c++) {
54  if(ti[c]<tj[c]) {
55  if(comp>0) { comp=0; break; }
56  comp = -1;
57  } else if(ti[c]>tj[c]) {
58  if(comp<0) { comp=0; break; }
59  comp = 1;
60  }
61  }
62  if(comp>0) {
63  flag[i]=true;
64  break;
65  } else if(comp<0) {
66  flag[j]=true;
67  }
68  }
69  }
70 
71  ex ret = 0;
72  for(int i=0; i<nn; i++) {
73  if(flag[i]) continue;
74  ex res = 1;
75  for(int c=0; c<nx; c++) res *= pow(xs[c],degs_vec[i][c]);
76  ret += res;
77  }
78 
79  return ret;
80  }
81 
88  bool SecDecBase::VerifySD(vector<exmap> map_vec, bool quick) {
89  lst xs;
90  exmap nxs;
91  for(auto kv : map_vec[0]) {
92  auto xi = kv.first;
93  if(is_zero(xi-x(-1))) continue;
94  xs.append(xi);
95  if(quick) nxs[xi] = xi.subs(x(w)==ex(w+1)/ex(w+13));
96  else nxs[xi];
97  }
98  ex chk_total = 1;
99  for(auto xi : xs) chk_total *= 1/(nxs[xi]+1);
100 
101  ex int_total = 0;
102  for(auto imap : map_vec) {
103  ex intg = 1;
104  for(auto kv : imap) {
105  auto xi = kv.first;
106  auto nxi = kv.first.subs(x(w)==(13+w));
107  if(is_zero(xi-x(-1))) intg *= kv.second;
108  else intg *= pow(kv.second, nxs[kv.first]);
109  }
110  while(true) {
111  auto tmp = intg;
112  tmp = tmp.subs(w0*pow(w1*w2,w3)==w0*pow(w1,w3)*pow(w2,w3));
113  tmp = xyz_pow_simplify(tmp);
114  if(is_zero(tmp-intg)) break;
115  intg = tmp;
116  }
117  intg = intg.subs(pow(y(w1),w2)==1/(w2+1)).subs(y(w)==1/ex(2));
118  int_total += intg;
119  }
120  int_total = int_total.normal();
121  return normal(chk_total-int_total).is_zero();
122  }
123  bool SecDec::VerifySD(vector<exmap> map_vec, bool quick) {
124  return SecDecBase::VerifySD(map_vec, quick);
125  }
126 
132  vector<exmap> SecDecBase::x2y(const lst & in_xpols) {
133  if(in_xpols.has(y(w))) throw Error("SecDecBase::x2y: y(w) found @ " + ex2str(in_xpols));
134  lst xpols_lst = in_xpols;
135  while(true) {
136  ex xpols = 1;
137  for(auto item : xpols_lst) {
138  if(item.match(x(w))) continue;
139  else if(item.match(pow(x(w1),w2))) continue;
140  else if(item.match(pow(w1,w2))) item = item.op(0);
141  //xpols *= factor_form(item, false);
142  xpols *= exfactor(item);
143  }
144  if(!is_a<mul>(xpols)) xpols = lst{ xpols };
145 
146  xpols_lst.remove_all();
147  for(auto item : xpols) {
148  if(item.match(x(w))) continue;
149  else if(item.match(pow(x(w1),w2))) continue;
150  else if(item.match(pow(w1,w2))) xpols_lst.append(item.op(0));
151  else xpols_lst.append(item);
152  }
153  xpols_lst.sort();
154  xpols_lst.unique();
155 
156  ex xpols2 = 1;
157  for(auto item : xpols_lst) xpols2 *= item;
158  if(is_a<lst>(xpols)) xpols = xpols.op(0);
159  if(is_zero(xpols2-xpols)) {
160  if(use_XMonomials) xpols = XMonomials(xpols);
161  return x2y(xpols);
162  }
163  }
164  throw Error("SecDecBase::x2y: unexpected region reached.");
165  return vector<exmap>();
166  }
167 
174  bool SecDec::IsBad(ex f, vector<exmap> vmap) {
175  for(auto &vi : vmap) {
176  auto ft = f.subs(vi);
177  auto xs_tmp = get_x_from(ft);
178  auto ys_tmp = get_y_from(ft);
179  int ysn = ys_tmp.size();
180  vector<int> y_free_indexs;
181  for(int j=0; j<ys_tmp.size()+10; j++) {
182  if(!ft.has(y(j))) y_free_indexs.push_back(j);
183  }
184  if(xs_tmp.size()>y_free_indexs.size()) throw Error("IsBad: xsize>ysize");
185  for(int i=0; i<xs_tmp.size(); i++) {
186  vi[xs_tmp[i]] = y(y_free_indexs[i]);
187  }
188  if(xs_tmp.size()>0) ft = ft.subs(vi);
189 
190  // need collect_common_factors
191  ft = collect_common_factors(expand_ex(ft,y(w)));
192  if(is_exactly_a<mul>(ft)) {
193  ex ret = 1;
194  for (auto item : ft) {
195  if( !(item.match(y(w)) || item.match(pow(y(w), w1))) ) {
196  ret *= item;
197  }
198  }
199  ft = ret;
200  } else if ( ft.match(y(w)) || ft.match(pow(y(w), w1)) ) {
201  ft = 1;
202  } // else ft not changed
203 
204  ys_tmp = get_y_from(ft);
205  for(int pi=1; pi<std::pow(2, ys_tmp.size()); pi++) {
206  lst yrepl;
207  int ci = pi;
208  for(int i=0; i<ys_tmp.size(); i++) {
209  yrepl.append(ys_tmp[i] == ((ci%2)==1 ? 1 : 0));
210  ci /= 2;
211  }
212  if(normal(ft.subs(yrepl)).is_zero()) return true;
213  }
214  }
215  return false;
216  }
217 
223  exvector SecDec::AutoEnd(ex po_ex) {
224  if(po_ex.nops()>2) throw Error("AutoEnd: Deltas found @ " + ex2str(po_ex));
225  lst const exlist = ex_to<lst>(po_ex.op(1));
226  if(!(exlist.op(0)-1).is_zero()) throw Error("AutoEnd: (!(exlist.op(0)-1).is_zero())");
227  auto xs = get_x_from(po_ex.op(0));
228  if(xs.size()<1) xs = get_y_from(po_ex.op(0));
229  int nx = xs.size();
230  for(int nn=0; nn<=nx; nn++) {
231  for(int pi=0; pi<std::pow(2, nx); pi++) {
232  int cpi = pi, cn1 = 0;
233  for(int i=0; i<nx; i++) {
234  if((cpi % 2) == 1) cn1++;
235  cpi /= 2;
236  }
237  if(cn1 != nn) continue;
238 
239  lst polists = lst{ po_ex.op(0) };
240  int bi = 0, bs = BisectionPoints.nops();
241  cpi = pi;
242  for(int i=0; i<nx; i++) {
243  if((cpi % 2) == 1) {
244  lst polists2;
245  auto x0 = BisectionPoints.op(bi % bs);
246  for(auto item : polists) {
247  symbol xy;
248  auto tmp = ex_to<lst>(subs(subs(item, lst{xs[i]==x0*xy}), lst{xy==xs[i]}));
249  tmp.let_op(0) = tmp.op(0) * x0;
250  polists2.append(tmp);
251 
252  tmp = ex_to<lst>(subs(subs(item, lst{xs[i]==(x0-1)*xy+1}), lst{xy==xs[i]}));
253  tmp.let_op(0) = tmp.op(0) * (1-x0);
254  polists2.append(tmp);
255  }
256  polists = polists2;
257  bi++;
258  }
259  cpi /= 2;
260  }
261  bool OK = true;
262  for(auto polist : polists) {
263  lst sdList;
264  for(int i=0; i<polist.nops(); i++) {
265  auto tmp = polist.op(i);
266  auto ntmp = NN(exlist.op(i));
267  if(!tmp.subs(lst{x(w)==0, y(w)==0}).normal().is_zero()) continue;
268  if( (!tmp.has(x(w)) && !tmp.has(y(w))) || (is_a<numeric>(ntmp) && ntmp>0) ) continue;
269  sdList.append(tmp);
270  }
271 
272  vector<exmap> vmap = SecDec->x2y(sdList);
273 
274  for(int ni=0; ni<polist.nops(); ni++) {
275  auto po = polist.op(ni);
276  auto expo = exlist.op(ni); // note that expo may decrease when using diff
277  if(expo.info(info_flags::posint)) continue;
278  if(IsBad(po, vmap)) {
279  OK = false;
280  break;
281  }
282  }
283  if(!OK) break;
284  }
285 
286  if(OK) {
287  exvector res;
288  for(auto item : polists) res.push_back(lst{ex_to<lst>(item), exlist});
289  return exvector(std::move(res));
290  }
291  }}
292 
293  cerr << ErrColor << "polynormial list: " << po_ex.op(0) << RESET << endl;
294  throw Error("AutoEnd Failed @ ALL possible bisections!");
295  return exvector();
296  }
297 
298  /*-----------------------------------------------------*/
299  /* SD */
300  /*-----------------------------------------------------*/
301  //return a lst, element pattern: { {{x1,n1}, {x2,n2}, ...}, {{e1, n1},{e2,n2}, ...} }
302  //e1 is a const term, e2 still the F-term
303  exvector SecDec::DS(const ex po_ex) {
304  // 1st element in input polist is the constant term, guess NOT necessary
305  // 2nd element in input polist is the F-term, required!
306  lst const polist = ex_to<lst>(po_ex.op(0));
307  lst const exlist = ex_to<lst>(po_ex.op(1));
308  lst sdList;
309  for(int i=0; i<polist.nops(); i++) {
310  auto tmp = polist.op(i);
311  auto ntmp = exlist.op(i);
312  if(!tmp.subs(lst{x(w)==0, y(w)==0}).normal().is_zero()) continue;
313  ntmp = NN(ntmp);
314  if( (!tmp.has(x(w)) && !tmp.has(y(w))) || (is_a<numeric>(ntmp) && ntmp>0) ) continue;
315  sdList.append(tmp);
316  }
317 
318  vector<exmap> vmap = SecDec->x2y(sdList);
319  if(!VerifySD(vmap)) {
320  for(auto vi : vmap) cout << vi << endl;
321  throw Error("VerifySD Failed!");
322  }
323 
324  exvector sd;
325  for(auto vi : vmap) {
326  auto ypolist = polist.subs(vi);
327  auto xs_tmp = get_x_from(ypolist);
328  auto ys_tmp = get_y_from(ypolist);
329  int ysn = ys_tmp.size();
330  vector<int> y_free_indexs;
331  for(int j=0; j<xs_tmp.size()+ys_tmp.size()+10; j++) {
332  if(!ypolist.has(y(j))) y_free_indexs.push_back(j);
333  }
334  if(xs_tmp.size()>y_free_indexs.size()) throw Error("DS: xsize>ysize");
335  for(int i=0; i<xs_tmp.size(); i++) {
336  vi[xs_tmp[i]] = y(y_free_indexs[i]);
337  }
338  if(xs_tmp.size()>0) ypolist = ypolist.subs(vi);
339  if(ypolist.has(x(w))) throw Error("DS: x(w) found @ " + ex2str(ypolist));
340 
341  // need collect_common_factors
342  auto ft = collect_common_factors(expand_ex(ypolist.op(1),y(w)));
343 
344  ex ct = 1, fsgin = 1;
345  if(is_a<mul>(ft)) {
346  ex ret = 1;
347  for (auto item : ft) {
348  if( !(item.match(y(w)) || item.match(pow(y(w), w1))) ) {
349  ret *= item;
350  }
351  }
352  ft = ret;
353  } else if( ft.match(y(w)) || ft.match(pow(y(w), w1)) ) {
354  ft = 1;
355  } // else ft not changed
356  bool hasWRA = false;
357  if(ft.has(WRA(w))) {
358  ft = 1;
359  hasWRA = true;
360  }
361 
362  exmap eps_map;
363  ex epn = ex(1)/111;
364  for(auto epi : eps_lst) {
365  eps_map[epi.op(0)] = epn;
366  epn = epn / 13;
367  }
368 
369  bool need_contour_deformation = ft.has(PL(w));
370  if(ft.has(y(w)) && !need_contour_deformation) {
371  auto tmp = ft.subs(nReplacement).subs(lst{ CV(w1,w2)==w2 }).subs(eps_map).expand();
372  if(is_a<add>(tmp)) {
373  need_contour_deformation = false;
374  auto first = tmp.op(0).subs(y(w)==1);
375  for(auto item : tmp) {
376  auto chk = item.subs(y(w)==1)*first;
377  auto nchk = NN(chk,30);
378  if(!is_a<numeric>(nchk)) {
379  throw Error("DS: Not a number: " + ex2str(chk));
380  }
381  if(nchk<0) {
382  need_contour_deformation = true;
383  break;
384  }
385  }
386  }
387  }
388 
389  if(disable_Contour || !need_contour_deformation) {
390  auto tmp = NN(ft.subs(y(w)==1).subs(nReplacement).subs(lst{ CV(w1,w2)==w2 }).subs(eps_map),30);
391 
392  if(!is_a<numeric>(tmp)) throw Error("DS: NOT a numeric with " + ex2str(tmp));
393  if(tmp<0) {
394  ct = exp(-I * Pi * exlist.op(1));
395  fsgin = -1;
396  }
397  ft = 1;
398  }
399 
400  exmap ymol;
401 
402  auto det = exfactor(vi[x(-1)]); // need factor
403  if(is_a<add>(det)) throw Error("DS: det is add " + ex2str(det));
404  auto ys = get_xy_from(det);
405  ex det1 = 1;
406  for(int i=0; i<ys.size(); i++) {
407  auto ndeg = det.degree(ys[i]);
408  if(ndeg!=det.ldegree(ys[i])) throw Error("DS: det is not a monomial with det = "+ex2str(det));
409  ymol[ys[i]] = ymol[ys[i]] + ndeg;
410  det1 *= pow(ys[i], ndeg);
411  }
412  if(!(is_a<numeric>(det/det1) && ex_to<numeric>(det/det1).is_integer())) {
413  throw Error("DS: det=" + ex2str(det) + ", det1=" + ex2str(det1));
414  }
415 
416  lst ftxlst;
417  for(auto xi : get_xy_from(ft)) ftxlst.append(xi);
418 
419  lst pol_exp_lst;
420  pol_exp_lst.append(lst{ subs(FTX(ft,ftxlst)*CT(ct*det/det1), y(w)==x(w)), 1 });
421 
422  for(int i=0; i<ypolist.nops(); i++) {
423  auto tmp = ypolist.op(i);
424  auto nexp = exlist.op(i).normal();
425  bool nchk = (!is_a<numeric>(nexp) || (is_a<numeric>(nexp) && nexp<0));
426  if(nexp.has(x(w)) || nexp.has(y(w))) throw Error("DS: x or y found in exp: "+ex2str(nexp));
427 
428  if(tmp.has(y(w))) {
429  lst ys;
430  for(auto yi : get_y_from(tmp)) ys.append(yi);
431  if(!tmp.is_polynomial(ys)) throw Error("DS: NOT a polynomial with "+ex2str(tmp));
432  }
433 
434  if(i==1 && fsgin<0) tmp = ex(0)-tmp; // check complete negtive F-term
435  auto tmp1 = tmp;
436  while(true) {
437  tmp=tmp1.subs(pow(y(w1)*w2, w3)==pow(y(w1),w3) * pow(w2, w3));
438  tmp = tmp.subs(pow(pow(y(w1),w2),w3)==pow(y(w1),w2*w3));
439  tmp = tmp.subs(pow(sqrt(y(w1)),w2)==pow(y(w1),w2/2));
440  tmp = tmp.subs(sqrt(pow(y(w1),w2))==pow(y(w1),w2/2));
441  if((tmp1-tmp).is_zero()) break;
442  tmp1 = tmp;
443  }
444 
445  if(tmp.has(y(w))) tmp = exfactor(tmp); // need factor
446 
447  ex rem = 1;
448  ex ct = 1;
449  ex tmps = tmp;
450  if(!is_a<mul>(tmps)) tmps = lst{tmps};
451  for (auto item : tmps) {
452  if(item.match(y(w))) {
453  auto yi = item;
454  ymol[yi] = ymol[yi] + nexp;
455  } else if(item.match(pow(y(w1), w2))) {
456  auto yi = item.op(0);
457  ymol[yi] = ymol[yi] + item.op(1) * nexp;
458  } else if(!item.has(y(w)) && !item.has(x(w))) {
459  if(is_a<numeric>(nexp) && ex_to<numeric>(nexp).is_integer()) {
460  ct *= item;
461  } else if(!item.has(PL(w)) && !item.has(WRA(w))) {
462  auto tr = NN(item.subs(nReplacement).subs(lst{ CV(w1,w2)==w2 }).subs(eps_map),30);
463  if(!is_a<numeric>(tr)) {
464  throw Error("DS: not numeric - item: " + ex2str(tr) + " ; " + ex2str(item));
465  }
466  auto nitem = ex_to<numeric>(tr);
467  if( nitem.is_real() && nitem<0 ) {
468  ct *= ex(-1)*item;
469  rem *= ex(-1);
470  } else {
471  ct *= item;
472  }
473  } else {
474  rem *= item;
475  }
476  } else {
477  if(nchk && item.subs(lst{y(w)==0,iEpsilon==0}).normal().is_zero()) {
478  throw Error("DS: zero item: " + ex2str(item) + " and exlist.op(i) = " + ex2str(nexp));
479  }
480  rem *= item;
481  }
482  }
483 
484  ex pnp = rem-(i==1 && (ft!=1 || hasWRA) ? iEpsilon : ex(0));
485  pnp = pnp.subs(y(w)==x(w));
486 
487  pol_exp_lst.append(lst{pnp, nexp});
488  pol_exp_lst.let_op(0) = pol_exp_lst.op(0).subs(CT(w)==CT(w*pow(ct,nexp))).subs(CT(0)==0);
489  }
490 
491  lst x_n_lst;
492  for(auto & kv : ymol) {
493  auto k = kv.first.subs(y(w)==x(w));
494  auto v = kv.second;
495  if(is_a<numeric>(v)) {
496  auto nv = ex_to<numeric>(v);
497  if(nv<=-1) {
498  throw Error("DS: " + ex2str(kv.first) + "^(" + ex2str(nv) + ") found, " + ex2str(vi));
499  }
500  }
501  x_n_lst.append(lst{k, v});
502  }
503  x_n_lst.sort();
504  sd.push_back(lst{x_n_lst, pol_exp_lst});
505  }
506 
507  return exvector(std::move(sd));
508  }
509 
510  // 1st element in output is the constant term
511  // 2nd element in both input and output is the F-term
512  lst SecDec::Normalize(const ex &input) {
513  ex const_term = 1;
514  lst plst, nlst;
515  lst in_plst = ex_to<lst>(input.op(0));
516  lst in_nlst = ex_to<lst>(input.op(1));
517 
518  auto vars = vec2lst(get_xy_from(input.op(0)));
519  int pnN = in_plst.nops();
520  for(int i=0; i<pnN; i++) {
521  auto pi = in_plst.op(i);
522  auto ni = in_nlst.op(i);
523 
524  if(!pi.is_polynomial(vars)) {
525  if(ni.info(info_flags::integer)) {
526  auto nd = numer_denom(exfactor(pi,o_flintf));
527  if(nd.op(0).is_polynomial(vars) && nd.op(1).is_polynomial(vars)) {
528  in_plst.append(nd.op(1));
529  in_nlst.append(ex(0)-ni);
530  in_plst.let_op(i) = nd.op(0);
531  goto ok;
532  }
533  } else {
534  auto nd = numer_denom(exfactor(pi,o_flintf));
535  if(nd.op(0).is_polynomial(vars) && nd.op(1).is_polynomial(vars)) {
536  if(xPositive(nd.op(1))) {
537  in_plst.append(nd.op(1));
538  in_nlst.append(ex(0)-ni);
539  in_plst.let_op(i) = nd.op(0);
540  goto ok;
541  } else if(xPositive(nd.op(0))) {
542  in_plst.append(nd.op(0));
543  in_nlst.append(ni);
544  in_plst.let_op(i) = nd.op(1);
545  in_nlst.let_op(i) = ex(0)-ni;
546  goto ok;
547  }
548  }
549  }
550 
551  if(true) {
552  ex pis;
553  if(!is_a<mul>(pi)) pis = lst{ pi };
554  else pis = pi;
555  ex rem = 1;
556  for(auto item : pis) {
557  if(item.is_polynomial(vars)) rem *= item;
558  else if(item.match(pow(w0,w1))) {
559  ex pp = item.op(0);
560  ex nn = item.op(1);
561  if(pp.is_polynomial(vars) && xPositive(pp)) {
562  in_plst.append(pp);
563  in_nlst.append(nn*ni);
564  } else {
565  auto nd = numer_denom(exfactor(pp,o_flintf));
566  if(nd.op(0).is_polynomial(vars) && nd.op(1).is_polynomial(vars)) {
567  if(xPositive(nd.op(0)) && xPositive(nd.op(1))) {
568  in_plst.append(nd.op(0));
569  in_plst.append(nd.op(1));
570  in_nlst.append(nn*ni);
571  in_nlst.append(ex(0)-nn*ni);
572  } else { cout << item << endl; goto nok; }
573  } else { cout << item << endl; goto nok; }
574  }
575  } else goto nok;
576  }
577  in_plst.let_op(i) = rem;
578  goto ok;
579  }
580 
581  nok: throw Error("SecDec::Normalize, NOT polynomial found: "+ex2str(pi));
582  ok: continue;
583  }
584  }
585 
586  for(int i=0; i<in_plst.nops(); i++) {
587  if(i!=1 && (in_nlst.op(i).is_zero() || in_plst.op(i)==ex(1))) continue;
588  if(i!=1 && !in_plst.op(i).has(x(w)) && !in_plst.op(i).has(y(w)) && !in_plst.op(i).has(z(w))) {
589  if(in_nlst.op(i).info(info_flags::integer) || xPositive(in_plst.op(i))) {
590  const_term *= pow(in_plst.op(i), in_nlst.op(i));
591  } else const_term *= exp(log(in_plst.op(i)) * in_nlst.op(i));
592  } else {
593  //auto ptmp = collect_common_factors(expand_ex(in_plst.op(i),{x(w),y(w),z(w)}));
594  auto ptmp = in_plst.op(i);
595  auto ntmp = in_nlst.op(i);
596  if(!is_a<mul>(ptmp)) ptmp = lst{ ptmp };
597 
598  exmap eps_map;
599  ex epn = ex(1)/111;
600  for(auto epi : eps_lst) {
601  eps_map[epi.op(0)] = epn;
602  epn = epn / 13;
603  }
604  ex tmul = 1;
605  for(int j=0; j<ptmp.nops(); j++) {
606  auto tmp = ptmp.op(j);
607  if(!tmp.has(x(w)) && !tmp.has(y(w)) && !tmp.has(z(w))) { // constant terms
608  if(ntmp.info(info_flags::integer)) {
609  const_term *= pow(tmp,ntmp);
610  } else if((tmp-vs).is_zero() || tmp.match(pow(vs,w))) {
611  const_term *= pow(tmp,ntmp);
612  } else if(!tmp.has(PL(w)) && !tmp.has(vs) && !tmp.has(WRA(w))) {
613  auto tr = NN(tmp.subs(nReplacement).subs(lst{ CV(w1,w2)==w2 }).subs(eps_map),30);
614  if(!is_a<numeric>(tr) || !tr.info(info_flags::real)) {
615  cerr << "tmp: " << tmp << endl;
616  cerr << "tr: " << tr << endl;
617  cerr << "nReplacement: " << nReplacement << endl;
618  throw Error("Normalize: tr is NOT numeric with nReplacement.");
619  }
620 
621  if(tr>0) {
622  const_term *= pow(tmp,ntmp);
623  } else {
624  const_term *= pow(-tmp,ntmp);
625  tmul = ex(0) - tmul;
626  }
627  } else {
628  tmul *= tmp;
629  }
630  } else if(tmp.match(pow(w1,w2)) && ntmp.info(info_flags::integer) && tmp.op(1).info(info_flags::integer)) {
631  plst.append(tmp.op(0));
632  nlst.append(ntmp * tmp.op(1));
633  } else if(tmp.match(pow(w1,w2)) && xPositive(tmp.op(0))) {
634  plst.append(tmp.op(0));
635  nlst.append(ntmp * tmp.op(1));
636  } else if(xSign(tmp)!=0 || xSign(tmp.subs(lst{y(w)==x(w),z(w)==x(w)}))!=0) {
637  if(tmp.subs(lst{x(w)==1, y(w)==1, z(w)==1})>0) {
638  plst.append(tmp);
639  nlst.append(ntmp);
640  } else {
641  plst.append(ex(0)-tmp);
642  nlst.append(ntmp);
643  tmul = ex(0) - tmul;
644  }
645  } else {
646  tmul *= tmp;
647  }
648  }
649 
650  if(i==1) {
651  plst.prepend(tmul);
652  nlst.prepend(ntmp);
653  } else if(tmul != 1) {
654  plst.append(tmul);
655  nlst.append(ntmp);
656  }
657 
658  }
659  }
660  plst.prepend(const_term);
661  nlst.prepend(1);
662 
663  lst plst_comb, nlst_comb;
664  exmap np;
665  for(int i=0; i<nlst.nops(); i++) {
666  if(i!=1) np[plst[i]] += nlst[i];
667  }
668  map<ex,int,ex_is_less> inp;
669  for(int i=0; i<nlst.nops(); i++) {
670  if(i==1) {
671  plst_comb.append(plst[1]);
672  nlst_comb.append(nlst[1]);
673  } else {
674  if(inp[plst[i]]==1) continue;
675  plst_comb.append(plst[i]);
676  nlst_comb.append(np[plst[i]]);
677  inp[plst[i]] = 1;
678  }
679  }
680 
681  if(nlst_comb.op(0)!=1) {
682  if(nlst_comb.op(0).info(info_flags::integer) || xPositive(plst_comb.op(0))) {
683  plst_comb.let_op(0) = pow(plst_comb.op(0),nlst_comb.op(0));
684  } else plst_comb.let_op(0) = exp(log(plst_comb.op(0))*nlst_comb.op(0));
685  nlst_comb.let_op(0) = 1;
686  }
687 
688  return (input.nops()>2) ? lst{plst_comb, nlst_comb, input.op(2)} : lst{plst_comb, nlst_comb};
689  }
690 
691  /*-----------------------------------------------------*/
692  /* 's Funtions in SD */
693  /*-----------------------------------------------------*/
694  // working with or without Deltas
696  if(IsZero || !use_XReOrders) return;
697  if(Integrands.size()<1) {
698  for(auto &fe : FunExp) {
699  auto xs = get_x_from(fe);
700  lst x2y;
701  for(int i=0; i<xs.size(); i++) x2y.append(xs[i]==y(i));
702  fe = ex_to<lst>(subs(fe, x2y).subs(y(w)==x(w)));
703  }
704  } else {
705  for(auto &vint : Integrands) {
706  auto xs = get_x_from(vint);
707  lst x2y;
708  for(int i=0; i<xs.size(); i++) x2y.append(xs[i]==y(i));
709  vint = vint.subs(x2y).subs(y(w)==x(w));
710  }
711  }
712  }
713 
714  // working with or without Deltas
716  for(int ri=0; ri<2; ri++) { // run twice, needs to check in more details
717  if(IsZero) return;
718  exvector funexp;
719  if(FunExp.size()<10) {
720  for(auto fe : FunExp) funexp.push_back(Normalize(fe));
721  } else {
722  funexp = GiNaC_Parallel(FunExp.size(), [&](int idx)->ex { return Normalize(FunExp[idx]); }, "Normalize");
723  }
724  FunExp.clear();
725  FunExp.shrink_to_fit();
726 
727  exmap fn;
728  for(auto fe : funexp) {
729  ex key = 1;
730  if(fe.nops()>2) key = iWF(fe.op(2)); // deltas
731  if(fe.op(1).op(0)!=1) {
732  cout << fe << endl;
733  throw Error("Normalizes: fe.op(0).op(0) is NOT 1.");
734  }
735  for(int i=1; i<fe.op(0).nops(); i++) key *= pow(fe.op(0).op(i), fe.op(1).op(i));
736  fn[key] += fe.op(0).op(0);
737  }
738 
739  map<ex,int,ex_is_less> ifn;
740  for(auto fe : funexp) {
741  ex key = 1;
742  if(fe.nops()>2) key = iWF(fe.op(2));
743  for(int i=1; i<fe.op(0).nops(); i++) key *= pow(fe.op(0).op(i), fe.op(1).op(i));
744  if(ifn[key]==1) continue;
745  lst funs, exps;
746  funs.append(fn[key]);
747  exps.append(1);
748  for(int i=1; i<fe.op(0).nops(); i++) {
749  funs.append(fe.op(0).op(i));
750  exps.append(fe.op(1).op(i));
751  }
752  if(fe.nops()>2) FunExp.push_back(lst{funs, exps, fe.op(2)});
753  else FunExp.push_back(lst{funs, exps});
754  ifn[key] = 1;
755  }
756  }
757  }
758 
759  // working with or without Deltas
761  exvector funexp;
762  for(auto fe : FunExp) {
763  funexp.push_back(fe);
764  }
765  FunExp.clear();
766  FunExp.shrink_to_fit();
767 
768  map<ex, ex, ex_is_less> fe_cc;
769  for(auto fe : funexp) {
770  lst fun, exp;
771  fun.append(1);
772  exp.append(1);
773  fun.append(fe.op(0).op(1));
774  exp.append(fe.op(1).op(1));
775  ex rem = 1;
776  for(int i=0; i<fe.op(1).nops(); i++) {
777  if(i==1) continue;
778  auto expi = fe.op(1).op(i);
779  if(expi.info(info_flags::nonnegint)) {
780  rem *= pow(fe.op(0).op(i), expi);
781  } else {
782  fun.append(fe.op(0).op(i));
783  exp.append(fe.op(1).op(i));
784  }
785  }
786  auto key = fe;
787  key.let_op(0) = fun;
788  key.let_op(1) = exp;
789  fe_cc[key] += rem;
790  }
791 
792  for(auto kv : fe_cc) {
793  lst fe = ex_to<lst>(kv.first);
794  let_op_append(fe, 0, kv.second);
795  let_op_append(fe, 1, 1);
796  FunExp.push_back(fe);
797  }
798  }
799 
800  // working with or without Deltas
802  exvector funexp;
803  for(auto fe : FunExp) {
804  funexp.push_back(fe);
805  }
806  FunExp.clear();
807  FunExp.shrink_to_fit();
808 
809  for(auto fe : funexp) {
810  lst fun = lst{fe.op(0).op(0), fe.op(0).op(1)};
811  lst exp = lst{fe.op(1).op(0), fe.op(1).op(1)};
812  ex rem = 1;
813  for(int i=2; i<fe.op(1).nops(); i++) {
814  auto expi = fe.op(1).op(i);
815  if(expi.info(info_flags::nonnegint)) {
816  rem *= pow(fe.op(0).op(i), expi);
817  } else {
818  fun.append(fe.op(0).op(i));
819  exp.append(fe.op(1).op(i));
820  }
821  }
822 
823  rem = collect_ex(rem, x(w));
824 
825  lst rem_lst;
826  if(is_a<add>(rem)) {
827  for(auto item : rem) rem_lst.append(item);
828  } else {
829  rem_lst.append(rem);
830  }
831 
832  for(auto item : rem_lst) {
833  auto fe2 = fe;
834  let_op_append(fe2, 0, item);
835  let_op_append(fe2, 1, 1);
836  fe2.let_op(0) = fun;
837  fe2.let_op(1) = exp;
838  FunExp.push_back(fe2);
839  }
840  }
841  }
842 
843  // Section 2.1 @ https://arxiv.org/pdf/1712.04441.pdf
844  // also refers to Feng/Thinking.pdf
846  if(IsZero) return;
847  if(Verbose > 2) cout << PRE << "\\--Scaleless: " << FunExp.size() << " :> " << flush;
848 
849  auto verb = Verbose; Verbose = 0;
850  auto sl_res =
851  GiNaC_Parallel(FunExp.size(), [this](int idx)->ex {
852  auto funexp = FunExp[idx];
853  if(funexp.nops()<3) return funexp;
854  symbol s;
855  auto fun = funexp.op(0);
856  auto exp = funexp.op(1);
857  auto deltas = funexp.op(2);
858  bool is0;
859  for(int di=0; di<deltas.nops(); di++) {
860  auto delta = ex_to<lst>(deltas.op(di));
861  is0 = false;
862  if(delta.nops()<2) continue;
863  // to make sure the integrand is projective
864  if(delta.nops()>1) {
865  lst sRepl;
866  for(int j=0; j<delta.nops(); j++) {
867  sRepl.append(delta[j]==delta[j]*s);
868  }
869 
870  bool is_s = true;
871  ex n_s = 0;
872  for(int j=0; j<fun.nops(); j++) {
873  auto tmp = expand(fun.op(j).subs(sRepl));
874  if(tmp.degree(s)!=tmp.ldegree(s)) {
875  is_s = false;
876  break;
877  }
878  n_s += tmp.degree(s) * exp.op(j);
879  }
880  if(!is_s || !normal(n_s+delta.nops()).is_zero()) continue;
881  }
882 
883  for(size_t i=1; i<ex_to<numeric>(GiNaC::pow(2,delta.nops())).to_long()-1; i++) {
884  lst sRepl;
885  auto ci = i;
886  ex n_s = 0;
887  for(int j=0; (j<delta.nops() && ci>0); j++) {
888  if((ci % 2)==1) {
889  sRepl.append(delta[j]==delta[j]*s);
890  n_s += 1;
891  }
892  ci = ci / 2;
893  }
894 
895  bool is_s = true;
896  for(int j=0; j<fun.nops(); j++) {
897  if(exp.op(j).info(info_flags::nonnegint)) continue;
898  auto tmp = collect_ex(fun.op(j).subs(sRepl),s);
899  if(tmp.degree(s)!=tmp.ldegree(s)) {
900  is_s = false;
901  break;
902  }
903  n_s += tmp.degree(s) * exp.op(j);
904  }
905  if(!is_s) continue;
906  if(!normal(n_s).is_zero()) {
907  is0 = true;
908  break;
909  }
910  }
911  if(is0) break;
912  }
913  if(!is0) return lst{fun, exp, deltas};
914  else return lst{};
915  }, "SL");
916  Verbose = verb;
917 
918  FunExp.clear();
919  FunExp.shrink_to_fit();
920  for(auto item : sl_res) {
921  lst fed = ex_to<lst>(item);
922  if(fed.nops()<1) continue;
923  FunExp.push_back(fed);
924  }
925 
926  if(Verbose > 2) cout << FunExp.size() << endl;
927  if(FunExp.size()<1) IsZero = true;
928  }
929 
930  void SecDec::RemoveDeltas() {
931  if(IsZero) return;
932 
933  while(true) {
934  bool exit = true;
935  exvector funexp = FunExp;
936  FunExp.clear();
937  FunExp.shrink_to_fit();
938  for(auto fe : funexp) {
939  if(fe.nops()<3) {
940  FunExp.push_back(fe);
941  continue;
942  }
943  auto xs = ex_to<lst>(fe.op(2).op(0));
944  lst re_deltas;
945  for(int i=1; i<fe.op(2).nops(); i++) {
946  re_deltas.append(fe.op(2).op(i));
947  }
948  if(re_deltas.nops()>0) exit = false;
949 
950  if(xs.nops()<1) {
951  FunExp.push_back(lst{fe.op(0), fe.op(1), re_deltas});
952  } else {
953  for(int i=0; i<xs.nops(); i++) {
954  auto xj = xs.op(i);
955  ex jInv = 0;
956  for(auto xi : xs) jInv += xi;
957  jInv = jInv.subs(xj==1);
958 
959  exmap repl;
960  for(int j=0; j<xs.nops(); j++) {
961  auto xxj = xs.op(j);
962  if(xxj != xj) repl[xxj] = xj*xxj;
963  }
964 
965  lst funs;
966  lst exps = ex_to<lst>(fe.op(1));
967  ex expns = 0;
968  for(int j=0; j<fe.op(0).nops(); j++) {
969  auto fun = fe.op(0).op(j);
970  fun = expand_ex(fun.subs(repl),xj);
971  if(!fun.is_polynomial(xj)) {
972  cerr << "xj: " << xj << endl;
973  cerr << "fun: " << fun << endl;
974  cerr << funexp << endl;
975  throw Error("RemoveDeltas: fun is NOT polynormial of xj.");
976  }
977  auto expn = fun.degree(xj);
978  fun = pow(xj, -expn) * fun;
979  fun = expand_ex(fun.subs(xj==1/xj),xj);
980  fun = fun.subs(xj==jInv);
981  funs.append(fun);
982  expns += expn * exps.op(j);
983  }
984 
985  funs.append(jInv);
986  exps.append(ex(0)-xs.nops()-expns);
987  if(re_deltas.nops()>0) FunExp.push_back(lst{funs, exps, re_deltas});
988  else FunExp.push_back(lst{funs, exps});
989  }
990  }
991  }
992  if(exit) break;
993  }
994  XReOrders();
995  if(use_Normalizes) Normalizes();
996  }
997 
998  void SecDec::XEnd() {
999  if(Verbose > 2) cout << PRE << "\\--BiSection: " << FunExp.size() << " :> " << flush;
1000  auto verb = Verbose; Verbose=0;
1001  GiNaC_Parallel_RM["BiSec"] = !Debug;
1002  auto funexps =
1003  GiNaC_Parallel(FunExp.size(), [this](int idx)->ex {
1004  auto fe = FunExp[idx];
1005  lst para_res_lst;
1006  if(xSign(fe.op(0).op(1))==0 && !fe.has(WRA(w))) {
1007  auto fes = AutoEnd(fe);
1008  for(auto fei : fes) para_res_lst.append(fei);
1009  } else {
1010  para_res_lst.append(fe);
1011  }
1012  return para_res_lst;
1013  }, "BiSec");
1014  Verbose = verb;
1015 
1016  FunExp.clear();
1017  FunExp.shrink_to_fit();
1018  for(auto &item : funexps) {
1019  for(auto &it : ex_to<lst>(item)) FunExp.push_back(ex_to<lst>(it));
1020  }
1021  if(Verbose > 2) cout << FunExp.size() << endl;
1022  }
1023 
1024  void SecDec::BiSection(ex xi, ex x0) {
1025  auto fes = FunExp;
1026  FunExp.clear();
1027  for(auto fe : fes) {
1028  symbol xy;
1029  auto tmp = ex_to<lst>(subs(subs(fe, lst{xi==x0*xy}), lst{xy==xi}));
1030  tmp.let_op(0).let_op(0) = tmp.op(0).op(0) * x0;
1031  FunExp.push_back(tmp);
1032 
1033  tmp = ex_to<lst>(subs(subs(fe, lst{xi==(x0-1)*xy+1}), lst{xy==xi}));
1034  tmp.let_op(0).let_op(0) = tmp.op(0).op(0) * (1-x0);
1035  FunExp.push_back(tmp);
1036  }
1037  }
1038 
1039  // after SDPrepares, Integrands can be expanded in ep safely.
1040  void SecDec::SDPrepares() {
1041  if(IsZero) return;
1042  if(FunExp.size()<1) {
1043  IsZero = true;
1044  return;
1045  }
1046  if(SecDec==NULL) SecDec = new SecDecG();
1047 
1048  MB();
1049  RemoveDeltas();
1050  if(CheckEnd) XEnd();
1051 
1052  Integrands.clear();
1053  Integrands.shrink_to_fit();
1054  auto fes = FunExp;
1055  FunExp.clear();
1056  FunExp.shrink_to_fit();
1057  for(auto &fe : fes) {
1058  bool to_add = true;
1059  for(auto item : fe.op(0)) {
1060  if(item.is_zero()) {
1061  to_add = false;
1062  break;
1063  }
1064  }
1065  if(to_add) FunExp.push_back(fe);
1066  }
1067  if(FunExp.size()<1) {
1068  IsZero = true;
1069  return;
1070  }
1071 
1072  exmap eps_map;
1073  for(auto epi : eps_lst) eps_map[epi.op(0)] = 0;
1074 
1075  SecDec->use_XMonomials = use_XMonomials;
1076  if(Verbose > 1) cout << Color_HighLight << " SDPrepares @ " << now() << RESET << endl;
1077  auto sd_res =
1078  GiNaC_Parallel(FunExp.size(), [this,&eps_map](int idx)->ex {
1079  // return a lst, element pattern: { {{x1,n1}, {x2,n2}, ...}, {{e1, n1},{e2,n2}, ...} }.
1080  auto fe = FunExp[idx];
1081  auto xns_pns = DS(fe);
1082  lst para_res_lst;
1083  for(auto const &item : xns_pns) {
1084 
1085  // take z-poles
1086  if(item.has(vz)) {
1087  auto ct = item.op(1).op(0).op(0);
1088  ct = ct.subs(lst{ CT(w)==w,FTX(w1,w2)==1 }).subs(pow(vs,vz)==1);
1089  int sNN = vsN - vsRank(ct.subs(pow(vs,vz)==1));
1090 
1091  lst zpols;
1092  // poles from Gamma(-z)
1093  for(int vn=0; vn<=sNN; vn++) zpols.append(vn);
1094 
1095  // poles from xi^{c1*z+c0} with c1<0
1096  for(auto xn : item.op(0)) {
1097  auto xn_op1 = expand(xn.op(1).normal());
1098  ex c1 = xn_op1.coeff(vz);
1099  ex c0 = xn_op1.subs(vz==0);
1100 
1101  if(!is_a<numeric>(c1)) {
1102  cerr << ErrColor << "SDPrepares: c1 is not a number: " << c1 << RESET << endl;
1103  exit(1);
1104  }
1105 
1106  if(c1<0) {
1107  int pxn = -1;
1108  while(true) {
1109  ex zp = (pxn-c0)/c1;
1110  ex zpn = zp.subs(eps_map).subs(lst{ep==0,epz==0});
1111  if(!is_a<numeric>(zpn)) {
1112  cerr << ErrColor << "SDPrepares: zpn is not a number: " << zpn << RESET << endl;
1113  exit(1);
1114  }
1115  if(zpn>sNN) break;
1116  zpols.append(zp);
1117  pxn -= 1;
1118  }
1119  }
1120  }
1121  zpols.sort();
1122  zpols.unique();
1123 
1124  symbol ss;
1125  for(auto zp : zpols) {
1126  auto item2 = item.subs(vz==ss+zp).subs(ss==vz);
1127  para_res_lst.append(item2);
1128  }
1129  } else {
1130  para_res_lst.append(item);
1131  }
1132  }
1133  return para_res_lst;
1134  }, "SD");
1135 
1136  ex min_expn = 1, min_expn2 = 10;
1137  exvector ibp_in_vec;
1138  for(auto &item : sd_res) {
1139  for(auto &it : ex_to<lst>(item)) {
1140  ex expn = 0;
1141  for(auto xn : it.op(0)) {
1142  ex nxn = xn.op(1).subs(eps_map).subs(lst{ep==0, vz==0});
1143  if(nxn<-1) expn += nxn+1;
1144  if(min_expn2>nxn) min_expn2 = nxn+1;
1145  }
1146  if(expn < min_expn) min_expn = expn;
1147 
1148  lst xns = ex_to<lst>(it.op(0));
1149  lst pns;
1150  for(int i=0; i<it.op(1).nops(); i++) {
1151  lst pn = ex_to<lst>(it.op(1).op(i));
1152  if(i<2) {
1153  pns.append(pn);
1154  continue;
1155  }
1156  if(pn.op(0).is_equal(1)) continue;
1157  if(is_zero(pn.op(1))) {
1158  if(is_zero(pn.op(0))) throw Error("SDPrepares: 0^0 found.");
1159  continue;
1160  }
1161  pns.append(pn);
1162  }
1163  ibp_in_vec.push_back(lst{xns, pns});
1164  }
1165  }
1166 
1167  if(Verbose > 1) cout << PRE << "\\--" << Color_HighLight << "Maximum x^-n: All(" << ex(0)-min_expn << "+1X), Max(" << (ex(0)-min_expn2) << "+1)" << RESET << endl;
1168 
1169  int pn = 0;
1170  exvector ibp_res_vec;
1171  while(ibp_in_vec.size()>0) {
1172  pn++;
1173  ostringstream spn;
1174  spn << "XIBP-" << (pn-1);
1175  auto ibp_res =
1176  GiNaC_Parallel(ibp_in_vec.size(), [&eps_map,&ibp_in_vec,this](int idx)->ex {
1177  // return lst
1178  // {0, element} for input with pole reached and doing nothing
1179  // {1, {element, ...}} for input whth pole NOT reached
1180  // element pattern still as { {{x1,n1}, {x2,n2}, ...}, {{e1, n1},{e2,n2}, ...} }
1181 
1182  auto xns_pns = ibp_in_vec[idx];
1183 
1184  auto xns = xns_pns.op(0);
1185  auto pns = xns_pns.op(1);
1186 
1187  exset fts;
1188  pns.op(0).find(FTX(w1,w2), fts);
1189  bool noFT = (fts.size()==1) && ( is_zero((*(fts.begin())).op(0)-1) );
1190  if(!noFT || disable_Contour) noFT = true;
1191 
1192  ex pole_requested = -1;
1193  if(noFT || PoleRequested > -1) pole_requested = PoleRequested;
1194 
1195  for(int n=0; n<xns.nops(); n++) {
1196  ex xn = xns.op(n);
1197  auto expn = xn.op(1).subs(eps_map).subs(lst{ep==0,vz==0,epz==0}).normal();
1198  if(!is_a<numeric>(expn)) throw Error("SDPrepares: expn NOT numeric: " + ex2str(expn));
1199 
1200  if(ex_to<numeric>(expn) < pole_requested) {
1201  auto xx = xn.op(0);
1202  pns.let_op(0).let_op(0) = pns.op(0).op(0) / (xn.op(1)+1);
1203 
1204  lst xns2;
1205  for(int i=0; i<xns.nops(); i++) {
1206  if(i!=n) xns2.append(xns.op(i));
1207  }
1208  lst pns2 = ex_to<lst>(pns);
1209  for(int i=0; i<pns.nops(); i++) {
1210  pns2.let_op(i).let_op(0) = pns2.op(i).op(0).subs(xx==1);
1211  }
1212 
1213  lst xns_pns_lst;
1214  xns_pns_lst.append(lst{xns2, pns2});
1215 
1216  for(int i=0; i<pns.nops(); i++) {
1217  lst xns3 = ex_to<lst>(xns);
1218  xns3.let_op(n).let_op(1) = xn.op(1)+1;
1219 
1220  ex tmp = ex(0)-pns.op(i).op(1)*diff_ex(pns.op(i).op(0),xx);
1221  if(tmp.is_zero()) continue;
1222 
1223  auto xs = get_x_from(tmp);
1224  bool tz = false;
1225  for(auto xi : xs) {
1226  if(tmp.subs(xi==0).is_zero()) {
1227  tz = true;
1228  break;
1229  }
1230  }
1231  // need collect_common_factors
1232  if(tz) tmp = collect_common_factors(expand_ex(tmp,x(w)));
1233  if(tmp.is_zero()) continue;
1234 
1235  if(tz && is_a<mul>(tmp)) {
1236  ex rem = 1;
1237  for(auto ii : tmp) {
1238  if(ii.match(pow(x(w), w2))) {
1239  bool t = true;
1240  for(int ij=0; ij<xns3.nops(); ij++) {
1241  if(xns3.op(ij).op(0)==ii.op(0)) {
1242  xns3.let_op(ij).let_op(1) += ii.op(1);
1243  t = false;
1244  break;
1245  }
1246  }
1247  if(t) xns3.append(lst{ii.op(0), ii.op(1)});
1248  } else if(ii.match(x(w))) {
1249  bool t = true;
1250  for(int ij=0; ij<xns3.nops(); ij++) {
1251  if(xns3.op(ij).op(0)==ii) {
1252  xns3.let_op(ij).let_op(1) += 1;
1253  t = false;
1254  break;
1255  }
1256  }
1257  if(t) xns3.append(lst{ii, 1});
1258  } else {
1259  rem *= ii;
1260  }
1261  }
1262  tmp = rem;
1263  }
1264 
1265  lst pns3 = ex_to<lst>(pns);
1266  if(is_zero(pns.op(i).op(1)-1)) {
1267  pns3.let_op(i).let_op(0) = tmp;
1268  } else {
1269  pns3.let_op(i).let_op(1) = pns.op(i).op(1)-1;
1270  int nn = pns.nops();
1271  if(!(pns.op(nn-1).op(1)-1).is_zero()) pns3.append(lst{ tmp, 1 });
1272  else pns3.let_op(nn-1).let_op(0) = pns.op(nn-1).op(0) * tmp;
1273  }
1274  xns_pns_lst.append(lst{xns3, pns3});
1275  }
1276  return lst{1, xns_pns_lst};
1277  }
1278  }
1279  return lst{0, xns_pns };
1280 
1281  }, spn.str().c_str());
1282 
1283  auto rets_vec =
1284  GiNaC_Parallel(ibp_res.size(), [&ibp_res](int idx)->ex {
1285  lst ret1, ret2;
1286  auto ii = ibp_res[idx];
1287  auto check = ii.op(0);
1288  if(check>0) {
1289  auto items = ii.op(1);
1290  for(auto &it : ex_to<lst>(items)) ret1.append(it);
1291  } else {
1292  exmap sub_exp;
1293  sub_exp[pow(exp(w1),w2)]=exp(w1*w2);
1294  sub_exp[sqrt(exp(w1))]=exp(w1/2);
1295  sub_exp[exp(w1)*exp(w2)*w0]=exp(w1+w2)*w0;
1296  sub_exp[exp(w1)*exp(w2)]=exp(w1+w2);
1297 
1298  auto item = ii.op(1);
1299  ex expr = 1, expr_exp = 1;
1300  for(auto pn : item.op(1)) {
1301  if(pn.op(0).has(CT(w)) || pn.op(0).has(FTX(w1,w2))) {
1302  if(pn.op(1)!=1) throw Error("SDPrepares: exponent of CT is NOT 1.");
1303  expr *= pn.op(0);
1304  } else if(xPositive(pn.op(0)) || pn.op(1).info(info_flags::integer)) {
1305  if(pn.op(0).has(exp(w))) expr_exp *= pow(pn.op(0), pn.op(1));
1306  else expr *= pow(pn.op(0), pn.op(1));
1307  } else {
1308  expr_exp *= exp(log(pn.op(0)) * pn.op(1));
1309  }
1310  }
1311  while(true) {
1312  auto expo = expr_exp.subs(sub_exp);
1313  if(is_zero(expo-expr_exp)) break;
1314  expr_exp = expo;
1315  }
1316  expr *= expr_exp;
1317  ret2.append(lst{ item.op(0), expr });
1318  }
1319  return lst{ret1, ret2};
1320  }, spn.str()+"R");
1321  ibp_in_vec.clear();
1322  ibp_in_vec.shrink_to_fit();
1323  for(auto rets : rets_vec) {
1324  for(auto item : rets.op(0)) ibp_in_vec.push_back(item);
1325  for(auto item : rets.op(1)) ibp_res_vec.push_back(item);
1326  }
1327  }
1328 
1329  auto res =
1330  GiNaC_Parallel(ibp_res_vec.size(), [&ibp_res_vec,&eps_map](int idx)->ex {
1331 
1332  // return single element in which ep/eps can be expanded safely.
1333  auto xns_expr = ibp_res_vec[idx];
1334  lst para_res_lst;
1335  auto xns = xns_expr.op(0);
1336  auto expr = xns_expr.op(1);
1337  lst exprs = { expr };
1338  symbol dx;
1339  for(auto xn : xns) {
1340  auto expn = xn.op(1).subs(eps_map).subs(lst{ep==0,vz==0,epz==0}).normal();
1341  if(!is_a<numeric>(expn)) throw Error("SDPrepares: Not a number with expn = " + ex2str(expn));
1342 
1343  lst exprs2;
1344  for(auto it : exprs) {
1345  ex rem = pow(xn.op(0), xn.op(1)) * it;
1346  if(ex_to<numeric>(expn)<=-1) {
1347  ex dit = it;
1348  ex dit0 = dit.subs(xn.op(0)==0);
1349  ex ifact = 1;
1350  if(!is_zero(dit0)) {
1351  rem -= pow(xn.op(0), xn.op(1)) * dit0 / ifact;
1352  exprs2.append(dit0/(xn.op(1)+1)/ifact);
1353  }
1354  for(int i=1; i+expn<0; i++) {
1355  dit = diff_ex(dit, xn.op(0));
1356  if(is_zero(dit)) break;
1357  dit0 = dit.subs(xn.op(0)==0);
1358  ifact *= i;
1359  if(!is_zero(dit0)) {
1360  rem -= pow(xn.op(0), xn.op(1)+i) * dit0 / ifact;
1361  exprs2.append(dit0/(xn.op(1)+i+1)/ifact);
1362  }
1363  }
1364  }
1365  exprs2.append(rem);
1366  }
1367  exprs = exprs2;
1368  }
1369 
1370  for(auto const &it : exprs) {
1371  if(it.is_zero()) continue;
1372  auto xs = get_x_from(it);
1373  lst x2y;
1374  for(int i=0; i<xs.size(); i++) x2y.append(xs[i]==y(i));
1375  para_res_lst.append(it.subs(x2y).subs(y(w)==x(w)));
1376  }
1377 
1378  return para_res_lst;
1379  }, "XTaylor");
1380 
1381  // Take z-residues
1382  bool zResides = false;
1383  exvector ints;
1384  for(auto &item : res) {
1385  for(auto it : ex_to<lst>(item)) {
1386  if(!it.is_zero()) ints.push_back(it);
1387  if(!zResides && it.has(vz)) zResides = true;
1388  }
1389  }
1390 
1391  if(zResides) {
1392  Integrands =
1393  GiNaC_Parallel(ints.size(), [&ints](int idx)->ex {
1394  auto item = ints[idx];
1395  ex it = item;
1396  if(it.has(vz)) {
1397  exset cts;
1398  it.find(CT(w),cts);
1399  lst repl;
1400  for(auto ct : cts) {
1401  ex cc = 1;
1402  ex ll = 1;
1403  lst cls;
1404  if(is_a<mul>(ct.op(0))) {
1405  for(auto ii : ct.op(0)) cls.append(ii);
1406  } else {
1407  cls.append(ct.op(0));
1408  }
1409  for(auto cl : cls) {
1410  if(cl.has(vz)) cc *= cl;
1411  else ll *= cl;
1412  }
1413  if(cc!=1) repl.append(ct==cc*CT(ll));
1414  }
1415  it = it.subs(repl);
1416  it = series_ex(it,vz,-1);
1417  it = ex(0)-it.coeff(vz, -1);
1418  }
1419  return it;
1420  }, "zResidue");
1421  } else {
1422  Integrands = ints;
1423  }
1424 
1425  }
1426 
1427  void SecDec::EpsExpands() {
1428  if(IsZero) return;
1429  if(Integrands.size()<1) {
1430  IsZero = true;
1431  return;
1432  }
1433 
1434  if(Verbose > 1) cout << Color_HighLight << " EpsExpands @ " << now() << RESET << endl;
1435 
1436  if(Verbose > 1) cout << PRE << "\\--Collecting: " << Integrands.size() << " :> " << flush;
1437  exmap int_map;
1438  for(auto &item : Integrands) {
1439  if(item.is_zero()) continue;
1440  exset cts;
1441  item.find(CT(w), cts);
1442  if(cts.size() != 1) {
1443  cerr << "cts: " << cts << endl;
1444  cerr << "item: " << item << endl;
1445  throw Error("EpsExpands: CT size is NOT 1: ");
1446  }
1447  ex ct = (*(cts.begin())).subs(CT(w)==w);
1448  auto it = item.subs(CT(w)==1);
1449  int_map[it] = int_map[it]+ct;
1450  }
1451 
1452  Integrands.clear();
1453  Integrands.shrink_to_fit();
1454  for(auto kv : int_map) {
1455  if(kv.second.is_zero()) continue;
1456  Integrands.push_back(CT(kv.second) * kv.first);
1457  }
1458  if(Verbose > 1) cout << Integrands.size() << endl;
1459 
1460  GiNaC_Parallel_RM["EpsEx"] = !Debug;
1461  auto res =
1462  GiNaC_Parallel(Integrands.size(), [this](int idx)->ex {
1463  // return { {two elements}, {two elements}, ...},
1464  // 1st: x-independent coefficient, expanded in ep/eps
1465  // 2nd: x-integrand
1466  auto item = Integrands[idx];
1467  if(item.is_zero()) return lst{ lst{0, 0} };
1468  exset cts;
1469  item.find(CT(w), cts);
1470  if(cts.size() != 1) {
1471  cerr << "cts: " << cts << endl;
1472  cerr << "item: " << item << endl;
1473  throw Error("EpsExpands: CT size is NOT 1: ");
1474  }
1475  ex ct = (*(cts.begin())).subs(CT(w)==w);
1476  auto it = item.subs(CT(w)==1);
1477 
1478  if(ct.has(epz)) {
1479  if(is_a<mul>(ct)) {
1480  ex ct0 = 1, ct1 = 1;
1481  for(auto ci : ct) {
1482  if(ci.has(epz)) ct1 *= ci;
1483  else ct0 *= ci;
1484  }
1485  ct = ct0;
1486  it *= ct1;
1487  } else {
1488  it *= ct;
1489  ct = 1;
1490  }
1491  }
1492  if(it.has(epz)) it = series_ex(it,epz,0);
1493 
1494  auto cvs = collect_lst(it, lst{epz, vs});
1495  lst para_res_lst;
1496  for(auto cv : cvs) {
1497  auto cc = cv.op(1) * ct; // multiple ct here
1498  auto vv = cv.op(0);
1499  for(auto epi : eps_lst) {
1500  auto epis = epi.op(0);
1501  if(!cc.has(epis) && !vv.has(epis)) continue;
1502  int iN = ex2int(epi.op(1));
1503  int cN = epsRank(cc, epis);
1504  vv = series_ex(vv, epis, iN-cN);
1505  int vN = vv.ldegree(epis);
1506  cc = series_ex(cc, epis, iN-vN);
1507  }
1508 
1509  lst pat_lst;
1510  for(auto item : eps_lst) pat_lst.append(item.op(0));
1511  auto ics = collect_lst(vv, pat_lst);
1512  for(auto ic : ics) {
1513  ex intg = ic.op(0);
1514  ex pref = ic.op(1)*cc;
1515  ex n1 = intg.subs(x(w)==ex(1)/7);
1516  if(is_zero(n1)) {
1517  n1 = intg.subs(x(w)==ex(1)/13);
1518  if(is_zero(n1)) intg = exnormal(intg);
1519  }
1520  for(auto epi : eps_lst) pref = series_ex(pref, epi.op(0), ex2int(epi.op(1)));
1521  para_res_lst.append(lst{pref, intg});
1522  }
1523  }
1524 
1525  return para_res_lst;
1526 
1527  }, "EpsEx");
1528 
1529  if(Verbose > 1) cout << PRE << "\\--Collecting: ";
1530  map<ex, ex, ex_is_less> int_pref;
1531  size_t ncollect = 0;
1532  ex expr_nox = 0;
1533  for(auto &item : res) {
1534  ncollect += item.nops();
1535  for(auto &kv : ex_to<lst>(item)) {
1536  if(!kv.op(1).has(x(w))) expr_nox += kv.op(0) * kv.op(1);
1537  else int_pref[kv.op(1)] += kv.op(0);
1538  }
1539  }
1540 
1541  if(Verbose > 1) cout << ncollect << " :> " << flush;
1542  expResult.clear();
1543  expResult.shrink_to_fit();
1544  for(auto kv : int_pref) {
1545  if(kv.second.is_zero()) continue;
1546  expResult.push_back(lst{kv.second, kv.first});
1547  }
1548  if(!is_zero(expr_nox)) {
1549  expr_nox = expr_nox.subs(FTX(w1,w2)==1);
1550  expResult.push_back(lst{expr_nox, 1});
1551  }
1552  if(Verbose > 1) cout << expResult.size() << endl;
1553  }
1554 
1555 }
#define RESET
Definition: BASIC.h:79
SecDec header file.
class used to wrap error message
Definition: BASIC.h:242
bool use_XMonomials
Definition: SD.h:98
virtual vector< exmap > x2y(const ex &xpol)=0
static bool VerifySD(vector< exmap > vmap, bool quick=true)
Verify the Sector Decompostion is valid.
Definition: SecDec.cpp:88
static ex XMonomials(const ex &expr)
XMonomials.
Definition: SecDec.cpp:18
SecDec by geometric method.
Definition: SD.h:106
static vector< vector< int > > RunQHull(const matrix &pts)
Definition: SecDecG.cpp:176
SecDec the main class to use Sector Decompostion method.
Definition: SD.h:413
exmap nReplacement
Definition: SD.h:425
bool IsBad(ex f, vector< exmap > vmap)
Check the x-END.
Definition: SecDec.cpp:174
SecDecBase * SecDec
Definition: SD.h:429
static bool VerifySD(vector< exmap > vmap, bool quick=true)
Definition: SecDec.cpp:123
exvector AutoEnd(ex po_ex)
Auto BiSection.
Definition: SecDec.cpp:223
bool use_XReOrders
Definition: SD.h:440
bool IsZero
Definition: SD.h:433
exvector FunExp
Definition: SD.h:426
bool disable_Contour
Definition: SD.h:424
void Scalelesses()
Definition: SecDec.cpp:845
bool use_XMonomials
Definition: SD.h:423
exvector Integrands
Definition: SD.h:427
lst BisectionPoints
Definition: SD.h:442
namespace for Numerical integration with Sector Decomposition method
Definition: AsyMB.cpp:10
exvector get_xy_from(ex pol)
Definition: Helpers.cpp:14
ex xyz_pow_simplify(const ex expr_in)
Definition: Helpers.cpp:336
exvector get_y_from(ex pol)
Definition: Helpers.cpp:44
exvector get_x_from(ex pol)
Definition: Helpers.cpp:29
int epsRank(ex expr_in, ex epi)
Definition: Helpers.cpp:195
ex expand_ex(ex const &expr_in, std::function< bool(const ex &)> has_func)
the expand like Mathematica
Definition: BASIC.cpp:1188
const char * Color_HighLight
Definition: BASIC.cpp:248
const iSymbol iEpsilon
bool xPositive(ex const expr)
check the expr is xPositive, i.e., each x-monomial item is postive
Definition: BASIC.cpp:1290
ex NN(ex expr, int digits)
the nuerical evaluation
Definition: BASIC.cpp:1278
ex w0
Definition: BASIC.h:499
map< string, bool > GiNaC_Parallel_RM
Definition: Init.cpp:148
const char * ErrColor
Definition: BASIC.cpp:246
const Symbol ep
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
ex series_ex(ex const &expr_in, ex const &s0, int sn0)
the series like Mathematica, include s^n
Definition: BASIC.cpp:968
string now(bool use_date)
date/time string
Definition: BASIC.cpp:525
bool Debug
Definition: Init.cpp:141
ex w
Definition: Init.cpp:90
const Symbol vs
ex exnormal(const ex &expr, int opt)
normalize a expression
Definition: BASIC.cpp:1906
const int o_flintf
Definition: Init.cpp:111
const Symbol epz
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
ex diff_ex(ex const expr, ex const xp, unsigned nth, bool expand)
the differential like Mathematica
Definition: BASIC.cpp:1063
int Verbose
Definition: Init.cpp:139
void let_op_append(ex &ex_in, const ex item)
append item into expression
Definition: BASIC.cpp:1324
lst vec2lst(const exvector &ev)
convert exvector to lst
Definition: BASIC.cpp:902
string ex2str(const ex &expr)
convert ex to output string, the defalut printer format will be used
Definition: BASIC.cpp:715
ex exfactor(const ex &expr, int opt)
factorize a expression
Definition: BASIC.cpp:1854
bool IsZero(const ex &e)
ex w1
Definition: BASIC.h:499
ex w3
Definition: BASIC.h:499
int xSign(ex const expr)
the always sign for expr
Definition: BASIC.cpp:1312
exvector GiNaC_Parallel(int ntotal, std::function< ex(int)> f, const string &key)
GiNaC Parallel Evaluation using fork.
Definition: BASIC.cpp:260
ex w2
Definition: BASIC.h:499
const Symbol vz
string PRE
Definition: Init.cpp:140
int ex2int(ex num)
ex to integer
Definition: BASIC.cpp:893
ex subs(const ex &e, init_list sl)
Definition: BASIC.h:51