HepLib
ChengWu.cpp
Go to the documentation of this file.
1 
6 #include "SD.h"
7 #include <cmath>
8 
9 namespace HepLib::SD {
10 
14  void SecDec::ChengWu(const ex & ft) {
15 
16  if(is_a<lst>(ft) && ft.nops()>1) {
17  for(auto item : ft) {
18  for(auto &fe : FunExp) {
19  let_op_prepend(fe, 0, item);
20  let_op_prepend(fe, 1, 0);
21  }
22  }
23 
24  int nt = ft.nops();
25  for(int i=0; i<nt; i++) {
26  auto FunExp2 = FunExp;
27  FunExp.clear();
28  FunExp.shrink_to_fit();
29  for(auto fe : FunExp2) {
30  auto ft2 = fe.op(0).op(0);
31  if(!fe.op(1).op(0).is_equal(0)) throw Error("SecDec::ChengWu exponent is NOT 0.");
32  let_op_remove_first(fe, 0);
33  let_op_remove_first(fe, 1);
34  auto fes = ChengWu::Apply(fe, ft2);
35  for(auto fe2 : fes) {
36  let_op_remove_first(fe2, 0);
37  let_op_remove_first(fe2, 1);
38  FunExp.push_back(fe2);
39  }
40  }
41  }
42  } else {
43  if(is_a<lst>(ft)) FunExp = ChengWu::Apply(FunExp, ft.op(0));
44  else FunExp = ChengWu::Apply(FunExp, ft);
45  //remove the first item in op(0) and op(1)
46  for(auto &fe : FunExp) {
47  let_op_remove_first(fe, 0);
48  let_op_remove_first(fe, 1);
49  }
50  }
51  }
52 
61  exvector ChengWu::Apply(const exvector & fe_vec, const ex & ft) {
62  exvector ret_fe_vec;
63  for(auto fe : fe_vec) {
64  ex cft = ft;
65  if(is_zero(cft)) {
66  if(fe.op(0).nops()>1) cft = fe.op(0).op(1);
67  else cft = 1;
68  }
69  if(fe.nops()<3 || xSign(cft)!=0 || cft.has(WRA(w))) {
70  let_op_prepend(fe, 0, 1);
71  let_op_prepend(fe, 1, 0);
72  ret_fe_vec.push_back(fe);
73  continue;
74  }
75  auto deltas = fe.op(2);
76  for(int di=0; di<deltas.nops(); di++) {
77  Projectivize(fe, deltas.op(di)); //make sure projective
78  }
79  let_op_prepend(fe, 0, cft);
80  let_op_prepend(fe, 1, 0);
81  auto ret = Evaluate(fe);
82  for(auto item : ret) ret_fe_vec.push_back(item);
83  }
84 
85  return exvector(std::move(ret_fe_vec));
86  }
87 
94  bool ChengWu::isProjective(const ex fe, const ex delta) {
95  static symbol s;
96  lst sRepl;
97  for(int j=0; j<delta.nops(); j++) sRepl.append(delta.op(j)==delta.op(j)*s);
98  ex over_all_sn = 0;
99  int nnn = fe.op(0).nops();
100  for(int i=0; i<nnn; i++) {
101  auto func = fe.op(0).op(i);
102  if(!func.has(x(w))) continue;
103  func = expand_ex(func.subs(sRepl),s);
104  auto sn = func.degree(s);
105  if(sn!=func.ldegree(s)) return false;
106  over_all_sn += sn*fe.op(1).op(i);
107  }
108  over_all_sn = normal(over_all_sn+delta.nops());
109  return over_all_sn.is_zero();
110  }
111 
118  void ChengWu::Projectivize(ex &fe, const ex delta, const ex xsum_in) {
119  static symbol s;
120  lst sRepl;
121  for(int j=0; j<delta.nops(); j++) sRepl.append(delta.op(j)==delta.op(j)*s);
122  ex xsum = xsum_in;
123  if(xsum.is_zero()) {
124  for(auto xi : delta) xsum += xi;
125  }
126 
127  ex over_all_sn = 0;
128  int nnn = fe.op(0).nops();
129  for(int i=0; i<nnn; i++) {
130  auto func = fe.op(0).op(i);
131  if(!func.has(x(w))) continue;
132  func = expand_ex(func.subs(sRepl),s);
133  auto sn = func.degree(s);
134  over_all_sn += sn*fe.op(1).op(i);
135  if(!is_a<add>(func)) func = lst{func};
136  ex tmp = 0;
137  for(auto ii : func) {
138  auto sni = ii.degree(s);
139  if(sni!=sn) tmp += ii.subs(s==1) * pow(xsum, sn-sni);
140  else tmp += ii.subs(s==1);
141  }
142  fe.let_op(0).let_op(i) = tmp;
143  }
144 
145  over_all_sn = normal(over_all_sn+delta.nops());
146  if(!over_all_sn.is_zero()) {
147  let_op_append(fe, 0, xsum);
148  let_op_append(fe, 1, ex(0)-over_all_sn);
149  }
150  }
151 
158  void ChengWu::Scalelize(ex &fe, const ex xi, const ex cy) {
159  if(is_a<lst>(xi)) Scalelize(fe, ex_to<lst>(xi), cy);
160  else Scalelize(fe, lst{xi}, cy);
161  }
162 
169  void ChengWu::Scalelize(ex &fe, const lst xs, const ex cy) {
170  lst x2y, y2x;
171  for(auto xi : xs) {
172  if(cy.has(xi)) {
173  throw Error("Scalelize: cy has xi: " + ex2str(cy) + "/" + ex2str(xi));
174  }
175  auto yi = xi.subs(x(w)==y(w));
176  x2y.append(xi==cy*yi);
177  y2x.append(yi==xi);
178  }
179  auto nnn = fe.op(0).nops();
180  for(int i=0; i<nnn; i++) {
181  if(!fe.op(0).op(i).has(x(w))) continue;
182  auto tmp = fe.op(0).op(i).subs(x2y).subs(y2x);
183  tmp = exfactor(tmp);
184  tmp = numer_denom(tmp);
185  if(tmp.op(1).subs(x(w)==1)<0) {
186  tmp.let_op(0) = ex(0)-tmp.op(0);
187  tmp.let_op(1) = ex(0)-tmp.op(1);
188  }
189  fe.let_op(0).let_op(i) = tmp.op(0);
190  if(tmp.op(1)!=1) {
191  let_op_append(fe, 0, tmp.op(1));
192  let_op_append(fe, 1, ex(0)-fe.op(1).op(i));
193  }
194  }
195 
196  ex det = exfactor(cy);
197  det = numer_denom(det);
198  if(det.op(1).subs(x(w)==1)<0) {
199  det.let_op(0) = ex(0)-det.op(0);
200  det.let_op(1) = ex(0)-det.op(1);
201  }
202  auto xn = xs.nops();
203  if(det.op(0)!=1) {
204  let_op_append(fe, 0, det.op(0));
205  let_op_append(fe, 1, xn);
206  }
207  if(det.op(1)!=1) {
208  let_op_append(fe, 0, det.op(1));
209  let_op_append(fe, 1, ex(0)-xn);
210  }
211  }
212 
219  exvector ChengWu::Binarize(ex const fe, ex const eqn) {
220  exvector ovec;
221  Binarize(fe, eqn, ovec);
222  return exvector(std::move(ovec));
223  }
224 
231  void ChengWu::Binarize(ex const fe, ex const eqn, exvector & ovec) {
232  auto xij = get_x_from(eqn);
233  if(xij.size()!=2) {
234  throw Error("Binarize: xij.size()!=2, " + ex2str(xij));
235  }
236  ex xi = xij[0];
237  ex xj = xij[1];
238  ex ci = eqn.coeff(xi);
239  ex cj = eqn.coeff(xj);
240  if(!((ci*xi+cj*xj-eqn).is_zero() && is_a<numeric>(ci * cj) && (ci*cj)<0)) {
241  throw Error("Binarize: ((ci*xi+cj*xj-eqn).is_zero() && is_a<numeric>(ci * cj) && (ci*cj)<0)");
242  }
243 
244  bool pos1 = (ci>0);
245 
246  ci = abs(ci);
247  cj = abs(cj);
248  symbol yi,yj;
249  // Part I: ci xi-cj xj>0, i.e., xi>cj/ci xj
250  auto f1 = ex_to<lst>(fe.op(0));
251  auto e1 = ex_to<lst>(fe.op(1));
252  ex c1 = cj/ci;
253  for(int i=0; i<f1.nops(); i++) {
254  f1.let_op(i) = f1.op(i).subs(lst{xi==c1*yj/(1+c1)+yi,xj==yj/(1+c1)}).subs(lst{yi==xi,yj==xj});
255  }
256  if(e1.op(0)==1) f1.let_op(0) = f1.op(0)/(1+c1); // Jaccobi
257  else if(e1.op(1)==1) f1.let_op(1) = f1.op(1)/(1+c1);
258  else {
259  f1.append(1/(1+c1));
260  e1.append(1);
261  }
262  auto fe1 = fe;
263  fe1.let_op(0) = f1;
264  fe1.let_op(1) = e1;
265 
266  // Part II: ci xi-cj xj<0, i.e., i.e., xj>ci/cj xi
267  auto f2 = ex_to<lst>(fe.op(0));
268  auto e2 = ex_to<lst>(fe.op(1));
269  ex c2 = ci/cj;
270  for(int i=0; i<f2.nops(); i++) {
271  f2.let_op(i) = f2.op(i).subs(lst{xj==c2*yi/(1+c2)+yj,xi==yi/(1+c2)}).subs(lst{yi==xi,yj==xj});
272  }
273  if(e2.op(0)==1) f2.let_op(0) = f2.op(0)/(1+c2); // Jaccobi
274  else if(e2.op(1)==1) f2.let_op(1) = f2.op(1)/(1+c2);
275  else {
276  f2.append(1/(1+c2));
277  e2.append(1);
278  }
279  auto fe2 = fe;
280  fe2.let_op(0) = f2;
281  fe2.let_op(1) = e2;
282 
283  // return vector
284  if(pos1) {
285  ovec.push_back(fe1);
286  ovec.push_back(fe2);
287  } else {
288  ovec.push_back(fe2);
289  ovec.push_back(fe1);
290  }
291  }
292 
300  bool ChengWu::isLinearizable(const ex ft, const ex delta, lst & xcs) {
301  for(auto xi : delta) {
302  ex f = ft;
303  if(!f.has(xi) || f.degree(xi)!=1) continue;
304  auto cxi = f.coeff(xi);
305  f = f.subs(xi==0);
306  int cxi_sgn = xSign(cxi);
307 
308  if(cxi_sgn!=0) {
309  if(cxi_sgn<0) cxi = ex(0)-cxi;
310  lst ret;
311  if(is_zero(f) || isLinearizable(f, delta, ret)) {
312  for(auto xc : ret) xcs.append(xc);
313  xcs.append(lst{xi, cxi});
314  return true;
315  }
316  }
317  }
318  return false;
319  }
320 
327  void ChengWu::Linearize(const lst xcs_in, ex & fe, ex & ft) {
328  lst xcs = xcs_in;
329  ex inv_det = 1;
330  auto nnn = xcs.nops();
331  for(int i=0; i<nnn; i++) {
332  auto xi = xcs.op(i).op(0);
333  auto s = xcs.op(i).op(1);
334  auto yi = xi.subs(x(w)==y(w));
335  inv_det *= s;
336  xcs.let_op(i).let_op(1) = yi/s;
337  for(int j=i+1; j<nnn; j++) {
338  xcs.let_op(j) = xcs.op(j).subs(xi==yi/s);
339  }
340  }
341  exmap x2y;
342  for(auto ss : xcs) x2y[ss.op(0)]=ss.op(1);
343 
344  // rescaling ft
345  if(ft.has(x(w))) ft = ft.subs(x2y).subs(y(w)==x(w));
346 
347  // rescaling each item
348  nnn = fe.op(0).nops();
349  for(int i=0; i<nnn; i++) {
350  if(!fe.op(0).op(i).has(x(w))) continue;
351  auto tmp = fe.op(0).op(i).subs(x2y);
352  tmp = tmp.subs(y(w)==x(w));
353  auto num_den = numer_denom(tmp);
354  if(num_den.op(1).subs(x(w)==1)<0) {
355  num_den.let_op(0) = ex(0)-num_den.op(0);
356  num_den.let_op(1) = ex(0)-num_den.op(1);
357  }
358  fe.let_op(0).let_op(i) = num_den.op(0);
359  if(num_den.op(1)!=1) {
360  let_op_append(fe, 0, num_den.op(1));
361  let_op_append(fe, 1, ex(0)-fe.op(1).op(i));
362  }
363  }
364 
365  // determinant
366  inv_det = inv_det.subs(y(w)==x(w));
367  auto idet_num_den = numer_denom(inv_det);
368  if(idet_num_den.op(1).subs(x(w)==1)<0) {
369  idet_num_den.let_op(0) = ex(0)-idet_num_den.op(0);
370  idet_num_den.let_op(1) = ex(0)-idet_num_den.op(1);
371  }
372  if(idet_num_den.op(0)!=1) {
373  let_op_append(fe, 0, idet_num_den.op(0));
374  let_op_append(fe, 1, -1);
375  }
376  if(idet_num_den.op(1)!=1) {
377  let_op_append(fe, 0, idet_num_den.op(1));
378  let_op_append(fe, 1, 1);
379  }
380  }
381 
396  bool ChengWu::isPartilizable(const ex ft, const ex delta, lst &xcs, int mode) {
397  for(auto xi : delta) {
398  if(!ft.has(xi) || ft.degree(xi)!=1) continue;
399 
400  auto cxi = ft.coeff(xi);
401  int cxi_sgn = xSign(cxi);
402  auto re_ft = ft.subs(xi==0);
403 
404  if(cxi_sgn!=0) {
405  if(cxi_sgn<0) cxi = ex(0)-cxi;
406 
407  lst ret;
408  ret.append(lst{xi, cxi});
409  if(is_zero(re_ft) || isPartilizable(re_ft, delta, ret, mode)) { // mode=0
410  for(int i=0; i<ret.nops(); i++) xcs.append(ret.op(i));
411  return true;
412  }
413 
414  if((mode>0) && (xSign(re_ft)!=0)) { // mode=1: x_i P_i + G_0, with P_i and G_0 same-sign
415  xcs.append(lst{xi, cxi});
416  if(re_ft.subs(x(w)==1)<0) re_ft = ex(0)-re_ft;
417  xcs.append(lst{0, re_ft});
418  return true;
419  } else if(mode>1) { // mode=2: x_i P_i + G_0, with P-i same-sign, G_0 ~ (xm-xn)^n
420  auto fflst = SecDec::XRefined_lst(re_ft);
421  if(fflst.nops()==1) {
422  symbol s;
423  auto ff = fflst.op(0).subs(x(w)==s*x(w));
424  if(get_x_from(ff).size()==2 && ff.degree(s)==1 && ff.ldegree(s)==1) {
425  xcs.append(lst{xi, cxi});
426  xcs.append(lst{0, fflst.op(0)});
427  return true;
428  }
429  }
430  }
431  } else if(mode>2) {
432  lst bilst;
433  auto cclst = SecDec::XRefined_lst(cxi);
434  if(cclst.nops()==1) {
435  symbol s;
436  auto cc = cclst.op(0).subs(x(w)==s*x(w));
437  if(get_x_from(cc).size()==2 && cc.degree(s)==1 && cc.ldegree(s)==1) {
438  bilst.append(lst{0, cclst.op(0)});
439  }
440  } else continue;
441  if(bilst.nops()!=1) continue;
442 
443  if(mode==3 && xSign(re_ft)!=0) {
444  // mode=3: x_i P_i + x_0 G_0 + Q_0, with P-i same-sign, G_0 ~ (xm-xn)^n & Q_0 same-sign
445  xcs.append(bilst.op(0));
446  return true;
447  }
448 
449  if(mode==4) {
450  // mode=4: x_i P_i + x_0 G_0 + Q_0, with P-i same-sign, G_0~(xm-xn)^n & Q_0~(xm'-xn')^n'
451  auto fflst = SecDec::XRefined_lst(re_ft);
452  if(fflst.nops()==1) {
453  symbol s;
454  auto ff = fflst.op(0).subs(x(w)==s*x(w));
455  if(get_x_from(ff).size()==2 && ff.degree(s)==1 && ff.ldegree(s)==1) {
456  bilst.append(lst{0, fflst.op(0)});
457  }
458  } else continue;
459  if(bilst.nops()!=2) continue;
460 
461  bool ok = true;
462  for(auto ii : get_x_from(bilst.op(0))) {
463  if(bilst.op(1).has(ii)) {
464  ok = false;
465  break;
466  }
467  }
468  if(ok) {
469  xcs.append(bilst.op(0));
470  xcs.append(bilst.op(1));
471  return true;
472  }
473  }
474  }
475 
476  // TODO: other modes
477  }
478  return false;
479  }
480 
488  void ChengWu::Partilize(const lst xcs, const lst delta_in, const ex fe_in, exvector & ret_lst) {
489  if(Verbose>10) {
490  cout << PRE << "\\--" << Color_HighLight << "ChengWu @xcs=" << xcs.nops() << RESET << endl;
491  }
492  lst ret = xcs;
493  auto fe = fe_in;
494  auto delta = delta_in;
495  auto ilast = ret.nops()-1;
496  ex inv_det = 1;
497  if(is_zero(get_op(ret,ilast,0))) {
498  lst rep_xs;
499  for(int i=ilast-1; i>=0; i--) rep_xs.append(get_op(ret,i,0));
500  ex xfi=0;
501  for(auto xi : delta) {
502  if(!rep_xs.has(xi) && !get_op(ret,ilast,1).has(xi)) {
503  xfi = xi;
504  break;
505  }
506  }
507  if(is_zero(xfi)) {
508  ex xs0=0;
509  for(auto xi : delta) {
510  if(!rep_xs.has(xi)) {
511  xs0 = xi;
512  break;
513  }
514  }
515  if(is_zero(xs0)) throw Error("Partilize: (xs0.is_zero())");
516  xfi = x(x_free_index(fe));
517  delta.append(xfi);
518  int dlti = -1;
519  for(int i=0; i<fe.op(2).nops(); i++) {
520  if(fe.op(2).op(i).has(xs0)) {
521  dlti = i;
522  break;
523  }
524  }
525  if(dlti<0) throw Error("Partilize: dlti<0 found.");
526  let_op_append(fe, 2, dlti, xfi);
527  let_op_append(fe, 0, xs0);
528  let_op_append(fe, 0, xfi+xs0);
529  let_op_append(fe, 1, 1);
530  let_op_append(fe, 1, -2);
531  }
532 
533  let_op(ret, ilast, 0, xfi);
534  for(int i=ilast-1; i>=0; i--) let_op(ret, i, 1, get_op(ret,i,1)*xfi);
535  }
536  lst rm_xs;
537  for(int i=ilast; i>=0; i--) {
538  auto xi = ret.op(i).op(0);
539  rm_xs.append(xi);
540  auto yi = xi.subs(x(w)==y(w));
541  auto s = ret.op(i).op(1);
542  inv_det *= s;
543  ret.let_op(i).let_op(1) = yi/s;
544  for(int j=i-1; j>=0; j--) {
545  ret.let_op(j) = ret.op(j).subs(xi==yi/s);
546  }
547  }
548  lst x2y;
549  for(auto ss : ret) x2y.append(ss.op(0)==ss.op(1));
550 
551  auto nnn = fe.op(0).nops();
552  for(int i=0; i<nnn; i++) {
553  if(!fe.op(0).op(i).has(x(w))) continue;
554  auto tmp = exfactor(fe.op(0).op(i).subs(x2y));
555  tmp = tmp.subs(y(w)==x(w));
556  auto num_den = numer_denom(tmp);
557  if(num_den.op(1).subs(x(w)==1)<0) {
558  num_den.let_op(0) = ex(0)-num_den.op(0);
559  num_den.let_op(1) = ex(0)-num_den.op(1);
560  }
561  fe.let_op(0).let_op(i) = num_den.op(0);
562  if(num_den.op(1)!=1) {
563  let_op_append(fe, 0, num_den.op(1));
564  let_op_append(fe, 1, ex(0)-fe.op(1).op(i));
565  }
566  }
567 
568  inv_det = exfactor(inv_det.subs(y(w)==x(w)));
569  auto idet_num_den = numer_denom(inv_det);
570  if(idet_num_den.op(1).subs(x(w)==1)<0) {
571  idet_num_den.let_op(0) = ex(0)-idet_num_den.op(0);
572  idet_num_den.let_op(1) = ex(0)-idet_num_den.op(1);
573  }
574  if(idet_num_den.op(0)!=1) {
575  let_op_append(fe, 0, idet_num_den.op(0));
576  let_op_append(fe, 1, -1);
577  }
578  if(idet_num_den.op(1)!=1) {
579  let_op_append(fe, 0, idet_num_den.op(1));
580  let_op_append(fe, 1, 1);
581  }
582 
583  if(!isProjective(fe,delta)) {
584  ex re_xi = 0;
585  for(auto xi : delta) {
586  if(!rm_xs.has(xi)) {
587  re_xi = xi;
588  break;
589  }
590  }
591  if(is_zero(re_xi)) {
592  cout << "xcs = " << xcs << endl;
593  cout << "fe = " << fe_in << endl;
594  throw Error("Partilize: rm_xs = " + ex2str(rm_xs) + " & delta = " + ex2str(delta));
595  }
596  Projectivize(fe, delta, re_xi);
597  }
598 
599  auto new_ft = SecDec::XRefined(get_op(fe, 0, 0));
600  symbol ss;
601  auto sft = new_ft.subs(x(w)==ss*x(w));
602  if(sft.degree(ss)!=1 || sft.ldegree(ss)!=1) throw Error("Partilize: ldegree/degree is NOT 1.");
603  auto rxs = get_x_from(new_ft);
604  lst xPos, xNeg;
605  ex sPos=0, sNeg=0;
606  for(auto xi : rxs) {
607  auto ci = new_ft.coeff(xi);
608  if(ci>0) {
609  xPos.append(xi); sPos += xi*ci;
610  } else {
611  xNeg.append(xi); sNeg -= xi*ci;
612  }
613  }
614  if(is_zero(sPos*sNeg)) throw Error("Partilize: sPos * sNeg = 0");
615  ex px = 0;
616  ex nx = 0;
617  for(int i=rm_xs.nops(); i>0; i--) {
618  auto xi = rm_xs.op(i-1);
619  if(is_zero(nx) && xNeg.has(xi)) nx = xi;
620  else if(is_zero(px) && xPos.has(xi)) px = xi;
621  if(!is_zero(px) && !is_zero(nx)) break;
622  }
623  if(sNeg<sPos) {
624  ex s = sNeg/nx;
625  if(s!=1) Scalelize(fe,xPos,s);
626  s = sPos/px;
627  if(s!=1) Scalelize(fe,nx,s);
628  } else {
629  ex s = sPos/px;
630  if(s!=1) Scalelize(fe,xNeg,s);
631  s = sNeg/nx;
632  if(s!=1) Scalelize(fe,px,s);
633  }
634  ex eqn = px-nx;
635  auto bilst = Binarize(fe, eqn);
636  for(auto item : bilst) {
637  if(!isProjective(item,delta)) throw Error("Partilize: something is wrong here.");
638  let_op(item, 0, 0, 1);
639  ret_lst.push_back(item);
640  }
641  }
642 
648  exvector ChengWu::Evaluate(const ex & in_fe) {
649  exvector fe_lst, ret_lst;
650  fe_lst.push_back(in_fe);
651  static int total_modes = 5;
652  while(true) {
653  exvector fe_lst2;
654  for(int i=0; i<fe_lst.size(); i++) {
655  auto fe = fe_lst[i];
656  auto ft = get_op(fe, 0, 0);
657 
658  // make sure, otherwise Projectivise may change things
659  if(!get_op(fe, 1, 0).is_zero()) {
660  throw Error("Evaluate: (!get_op(fe, 1, 0).is_zero())");
661  }
662  ft = SecDec::XRefined(ft);
663 
664  if(fe.nops()<3 || xSign(ft)!=0) {
665  let_op(fe, 0, 0, 1);
666  ret_lst.push_back(fe);
667  goto ok_label;
668  }
669 
670  // loop modes and deltas
671  for(int mode=0; mode<total_modes; mode++) {
672  for(int di=0; di<fe.op(2).nops(); di++) {
673  lst delta = ex_to<lst>(fe.op(2).op(di));
674  if(delta.nops()<2) continue;
675 
676  lst ret;
677  if(mode==0 && isPartilizable(ft, delta, ret, mode)) {
678  Partilize(ret, delta, fe, ret_lst);
679  goto ok_label;
680  }
681 
682  ret.remove_all();
683  if(mode==1 && isPartilizable(ft, delta, ret, mode)) {
684  Partilize(ret, delta, fe, ret_lst);
685  goto ok_label;
686  }
687 
688  ret.remove_all();
689  if(mode==2 && isPartilizable(ft, delta, ret, mode)) {
690  auto bilst = Binarize(fe, get_op(ret, ret.nops()-1, 1));
691  for(auto item : bilst) fe_lst2.push_back(item);
692  goto ok_label;
693  }
694 
695  ret.remove_all();
696  if(mode==3 && isPartilizable(ft, delta, ret, mode)) {
697  auto bilst = Binarize(fe, get_op(ret, ret.nops()-1, 1));
698  for(auto item : bilst) fe_lst2.push_back(item);
699  goto ok_label;
700  }
701 
702  ret.remove_all();
703  if(mode==4 && isPartilizable(ft, delta, ret, mode)) {
704  auto eq1 = get_op(ret, ret.nops()-1, 1);
705  auto eq2 = get_op(ret, ret.nops()-2, 1);
706  for(auto item : Binarize(fe, eq1)) {
707  for(auto ii : Binarize(item, eq2)) fe_lst2.push_back(ii);
708  }
709  goto ok_label;
710  }
711 
712  //TODO: other mode
713 
714  }}
715  let_op(fe, 0, 0, 2);
716  ret_lst.push_back(fe);
717 
718  ok_label: ;
719  }
720  if(fe_lst2.size()<1) break;
721  fe_lst = fe_lst2;
722  }
723 
724  return exvector(std::move(ret_lst));
725  }
726 
727  namespace {
728  int maxn(const ex pols, const ex xi) {
729  int max_xn = -1;
730  for(auto item : pols) {
731  auto cxn = item.degree(xi);
732  if(max_xn<cxn) max_xn = cxn;
733  }
734  return max_xn;
735  }
736  }
737 
743  exvector ChengWu::WickRotation(const exvector & fe_vec) {
744  exvector ret_vec, run_vec, run2_vec;
745  run_vec = fe_vec;
746  ReRun:
747  for(auto fe : run_vec) {
748  ex ft = fe.op(0).op(1);
749  auto xs = get_x_from(ft);
750 
751  // Case: x(cx!=0) + (c0!=0)
752  for(auto xi : xs) {
753  auto ftx = collect_ex(ft,xi);
754  auto c0 = ftx.subs(xi==0);
755 
756  auto cx = ftx-c0;
757  if(xSign(cx)!=0) {
758  if(xSign(c0)!=0) {
759  auto max_xn = maxn(fe.op(0),xi)+1;
760  auto wra = WRA(-xSign(cx) * Pi/max_xn);
761  fe.let_op(0) = fe.op(0).subs(xi==exp(I*wra)*xi);
762  if(fe.op(1).op(0)!=1) throw Error("WickRotation: fe.op(0).op(0)!=1.");
763  fe.let_op(0).let_op(0) = fe.op(0).op(0) * exp(I*wra);
764  ret_vec.push_back(fe);
765  goto next_fe;
766  } else {
767  auto item = SecDec::XRefined(c0);
768  auto xs2 = get_x_from(item);
769  for(auto xi2 : xs2) {
770  if(item.degree(xi2)!=1) continue;
771  auto xc0 = item.subs(xi2==0);
772  auto xc1 = item.coeff(xi2);
773  if(xSign(xc0)*xSign(xc1)<0) {
774  ex xx = 0;
775  for(auto xi3 : xs2) {
776  if(xi3==xi || xi3==xi2) continue;
777  xx = xi3;
778  break;
779  }
780  if(is_zero(xx)) continue;
781  Scalelize(fe,xi2,-xc0/xc1/xx);
782  auto bfe = Binarize(fe, xi2-xx);
783  for(auto item2 : bfe) run2_vec.push_back(item2);
784  goto next_fe;
785  }
786  }
787  }
788  }
789  }
790 
791  // Case: a*x^2+b*x+c with a*c<0 or b*c>0
792  for(auto xi : xs) {
793  auto ftx = collect_ex(ft,xi);
794  if(ftx.degree(xi)!=2) continue;
795 
796  auto c0 = ftx.coeff(xi,0);
797  auto c1 = ftx.coeff(xi,1);
798  auto c2 = ftx.coeff(xi,2);
799 
800  if(xSign(c2) * xSign(c0) <0) {
801  if(xSign(c1)!=0) {
802  auto max_xn = maxn(fe.op(0),xi)+1;
803  auto wra = WRA(-xSign(c2) * Pi/max_xn);
804  fe.let_op(0) = fe.op(0).subs(xi==exp(I*wra)*xi);
805  fe.let_op(0).let_op(1) = c2*pow(xi,2)*exp(I*wra)+c1*xi+c0*exp(-I*wra);
806  if(fe.op(1).op(0)!=1) throw Error("WickRotation: fe.op(0).op(0)!=1.");
807  fe.let_op(0).let_op(0) = fe.op(0).op(0) * exp(I*wra*(1+fe.op(1).op(1)));
808  ret_vec.push_back(fe);
809  goto next_fe;
810  } else {
811  auto item = SecDec::XRefined(c1);
812  auto xs2 = get_x_from(item);
813  for(auto xi2 : xs2) {
814  if(item.degree(xi2)!=1) continue;
815  auto xc0 = item.subs(xi2==0);
816  auto xc1 = item.coeff(xi2);
817  if(xSign(xc0)*xSign(xc1)<0) {
818  ex xx = 0;
819  for(auto xi3 : xs2) {
820  if(xi3==xi || xi3==xi2) continue;
821  xx = xi3;
822  break;
823  }
824  if(is_zero(xx)) continue;
825  Scalelize(fe,xi2,-xc0/xc1/xx);
826  auto bfe = Binarize(fe, xi2-xx);
827  for(auto item2 : bfe) run2_vec.push_back(item2);
828  goto next_fe;
829  }
830  }
831  }
832  }
833 
834  if(xSign(c1) * xSign(c0) >0) {
835  if(xSign(c2)!=0) {
836  auto max_xn = maxn(fe.op(0),xi)+1;
837  auto wra = WRA(xSign(c0) * Pi/max_xn);
838  fe.let_op(0) = fe.op(0).subs(xi==exp(I*wra)*xi);
839  fe.let_op(0).let_op(1) = c2*pow(xi,2)+c1*xi*exp(-I*wra)+c0*exp(-2*I*wra);
840  if(fe.op(1).op(0)!=1) throw Error("WickRotation: fe.op(0).op(0)!=1.");
841  fe.let_op(0).let_op(0) = fe.op(0).op(0) * exp(I*wra*(1+2*fe.op(1).op(1)));
842  ret_vec.push_back(fe);
843  goto next_fe;
844  } else {
845  auto item = SecDec::XRefined(c2);
846  auto xs2 = get_x_from(item);
847  for(auto xi2 : xs2) {
848  if(item.degree(xi2)!=1) continue;
849  auto xc0 = item.subs(xi2==0);
850  auto xc1 = item.coeff(xi2);
851  if(xSign(xc0)*xSign(xc1)<0) {
852  ex xx = 0;
853  for(auto xi3 : xs2) {
854  if(xi3==xi || xi3==xi2) continue;
855  xx = xi3;
856  break;
857  }
858  if(is_zero(xx)) continue;
859  Scalelize(fe,xi2,-xc0/xc1/xx);
860  auto bfe = Binarize(fe, xi2-xx);
861  for(auto item2 : bfe) run2_vec.push_back(item2);
862  goto next_fe;
863  }
864  }
865  }
866  }
867 
868  auto c0x = Factor(c0);
869  if(is_a<mul>(c0x)) {
870  ex cc = 1;
871  for(auto item : c0x) {
872  if(!is_a<numeric>(item) && !item.match(x(w))) cc *= item;
873  }
874  c0x = cc;
875  }
876  if(Factor(c0x).match(pow(w,2))) {
877  auto item = SecDec::XRefined(c0x);
878  auto xs2 = get_x_from(item);
879  for(auto xi2 : xs2) {
880  if(item.degree(xi2)!=1) continue;
881  auto xc0 = item.subs(xi2==0);
882  auto xc1 = item.coeff(xi2);
883  if(xSign(xc0)*xSign(xc1)<0) {
884  ex xx = 0;
885  for(auto xi3 : xs2) {
886  if(xi3==xi || xi3==xi2) continue;
887  xx = xi3;
888  break;
889  }
890  if(is_zero(xx)) continue;
891  Scalelize(fe,xi2,-xc0/xc1/xx);
892  auto bfe = Binarize(fe, xi2-xx);
893  for(auto item2 : bfe) run2_vec.push_back(item2);
894  goto next_fe;
895  }
896  }
897  }
898  }
899 
900  ret_vec.push_back(fe);
901  next_fe: ;
902  }
903  if(run2_vec.size()>0) {
904  run_vec = run2_vec;
905  run2_vec.clear();
906  goto ReRun;
907  }
908  return exvector(std::move(ret_vec));
909  }
910 
911 }
#define RESET
Definition: BASIC.h:79
SecDec header file.
class used to wrap error message
Definition: BASIC.h:242
static bool isPartilizable(const ex ft, const ex delta, lst &xcs, int mode=0)
isPartilize w.r.t. F-term, generized to isLinearizable
Definition: ChengWu.cpp:396
static bool isLinearizable(const ex ft, const ex delta, lst &xcs)
Linearize w.r.t. F-term.
Definition: ChengWu.cpp:300
static exvector Apply(const exvector &fe_vec, const ex &ft=0)
ChengWu-rized on the vector of fe, note that 1st element of the output, which needs to be droped,...
Definition: ChengWu.cpp:61
static exvector Binarize(ex const fe, ex const eqn)
Binarize the input fe, w.r.t. the linear eqn.
Definition: ChengWu.cpp:219
static void Projectivize(ex &fe, const ex delta, const ex xsum=0)
to Projectivize the input fe, the fe will be updated
Definition: ChengWu.cpp:118
static void Partilize(const lst xcs, const lst delta, const ex fe, exvector &ret_lst)
Partilize function, generized to Linearize.
Definition: ChengWu.cpp:488
static void Scalelize(ex &fe, const lst xs, const ex cy)
to Scalelize the input fe, the fe will be updated
Definition: ChengWu.cpp:169
static bool isProjective(const ex fe, const ex delta)
to check the input fe is Projective or NOT w.r.t. delta
Definition: ChengWu.cpp:94
static exvector Evaluate(const ex &fe)
ChengWu Internal, make sure ft in the first term, ONLY appear in ChengWu.cpp.
Definition: ChengWu.cpp:648
static void Linearize(const lst xcs, ex &fe, ex &ft)
Linearize w.r.t. xcs_in.
Definition: ChengWu.cpp:327
static exvector WickRotation(const exvector &fe_vec)
WickRotation, just check WRA exist or NOT to see successful or NOT. Still Experimental.
Definition: ChengWu.cpp:743
static lst XRefined_lst(ex const &ft)
Definition: Helpers.cpp:392
static ex XRefined(ex const &ft)
Definition: Helpers.cpp:369
void ChengWu(const ex &ft=0)
ChengWu for SecDec class.
Definition: ChengWu.cpp:14
exvector FunExp
Definition: SD.h:426
namespace for Numerical integration with Sector Decomposition method
Definition: AsyMB.cpp:10
int x_free_index(ex expr)
Definition: Helpers.cpp:165
exvector get_x_from(ex pol)
Definition: Helpers.cpp:29
ex Factor(const ex expr_in)
Definition: Helpers.cpp:244
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
void let_op_prepend(ex &ex_in, const ex item)
preppend item into expression
Definition: BASIC.cpp:1335
void let_op(ex &ex_in, int index1, int index2, const ex item)
update index1-th.index2-th of expression with item
Definition: BASIC.cpp:1560
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 w
Definition: Init.cpp:90
ex get_op(const ex ex_in, int index1, int index2)
return index1-th.index2-th of expression
Definition: BASIC.cpp:1606
int Verbose
Definition: Init.cpp:139
void let_op_append(ex &ex_in, const ex item)
append item into expression
Definition: BASIC.cpp:1324
string ex2str(const ex &expr)
convert ex to output string, the defalut printer format will be used
Definition: BASIC.cpp:715
void let_op_remove_first(ex &ex_in)
remove first from expression
Definition: BASIC.cpp:1355
ex exfactor(const ex &expr, int opt)
factorize a expression
Definition: BASIC.cpp:1854
int xSign(ex const expr)
the always sign for expr
Definition: BASIC.cpp:1312
string PRE
Definition: Init.cpp:140
ex subs(const ex &e, init_list sl)
Definition: BASIC.h:51