16 if(is_a<lst>(ft) && ft.nops()>1) {
25 for(
int i=0; i<nt; i++) {
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.");
63 for(
auto fe : fe_vec) {
66 if(fe.op(0).nops()>1) cft = fe.op(0).op(1);
69 if(fe.nops()<3 ||
xSign(cft)!=0 || cft.has(WRA(
w))) {
72 ret_fe_vec.push_back(fe);
75 auto deltas = fe.op(2);
76 for(
int di=0; di<deltas.nops(); di++) {
82 for(
auto item : ret) ret_fe_vec.push_back(item);
85 return exvector(std::move(ret_fe_vec));
97 for(
int j=0; j<delta.nops(); j++) sRepl.append(delta.op(j)==delta.op(j)*s);
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;
104 auto sn = func.degree(s);
105 if(sn!=func.ldegree(s))
return false;
106 over_all_sn += sn*fe.op(1).op(i);
108 over_all_sn = normal(over_all_sn+delta.nops());
109 return over_all_sn.is_zero();
121 for(
int j=0; j<delta.nops(); j++) sRepl.append(delta.op(j)==delta.op(j)*s);
124 for(
auto xi : delta) xsum += xi;
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;
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};
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);
142 fe.let_op(0).let_op(i) = tmp;
145 over_all_sn = normal(over_all_sn+delta.nops());
146 if(!over_all_sn.is_zero()) {
159 if(is_a<lst>(xi))
Scalelize(fe, ex_to<lst>(xi), cy);
175 auto yi = xi.subs(x(
w)==y(
w));
176 x2y.append(xi==cy*yi);
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);
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);
189 fe.let_op(0).let_op(i) = tmp.op(0);
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);
222 return exvector(std::move(ovec));
234 throw Error(
"Binarize: xij.size()!=2, " +
ex2str(xij));
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)");
250 auto f1 = ex_to<lst>(fe.op(0));
251 auto e1 = ex_to<lst>(fe.op(1));
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});
256 if(e1.op(0)==1) f1.let_op(0) = f1.op(0)/(1+c1);
257 else if(e1.op(1)==1) f1.let_op(1) = f1.op(1)/(1+c1);
267 auto f2 = ex_to<lst>(fe.op(0));
268 auto e2 = ex_to<lst>(fe.op(1));
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});
273 if(e2.op(0)==1) f2.let_op(0) = f2.op(0)/(1+c2);
274 else if(e2.op(1)==1) f2.let_op(1) = f2.op(1)/(1+c2);
301 for(
auto xi : delta) {
303 if(!f.has(xi) || f.degree(xi)!=1)
continue;
304 auto cxi = f.coeff(xi);
306 int cxi_sgn =
xSign(cxi);
309 if(cxi_sgn<0) cxi = ex(0)-cxi;
312 for(
auto xc : ret) xcs.append(xc);
313 xcs.append(lst{xi, cxi});
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));
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);
342 for(
auto ss : xcs) x2y[ss.op(0)]=ss.op(1);
345 if(ft.has(x(
w))) ft = ft.subs(x2y).subs(y(
w)==x(
w));
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);
358 fe.let_op(0).let_op(i) = num_den.op(0);
359 if(num_den.op(1)!=1) {
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);
372 if(idet_num_den.op(0)!=1) {
376 if(idet_num_den.op(1)!=1) {
397 for(
auto xi : delta) {
398 if(!ft.has(xi) || ft.degree(xi)!=1)
continue;
400 auto cxi = ft.coeff(xi);
401 int cxi_sgn =
xSign(cxi);
402 auto re_ft = ft.subs(xi==0);
405 if(cxi_sgn<0) cxi = ex(0)-cxi;
408 ret.append(lst{xi, cxi});
410 for(
int i=0; i<ret.nops(); i++) xcs.append(ret.op(i));
414 if((mode>0) && (
xSign(re_ft)!=0)) {
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});
421 if(fflst.nops()==1) {
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)});
434 if(cclst.nops()==1) {
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)});
441 if(bilst.nops()!=1)
continue;
443 if(mode==3 &&
xSign(re_ft)!=0) {
445 xcs.append(bilst.op(0));
452 if(fflst.nops()==1) {
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)});
459 if(bilst.nops()!=2)
continue;
463 if(bilst.op(1).has(ii)) {
469 xcs.append(bilst.op(0));
470 xcs.append(bilst.op(1));
494 auto delta = delta_in;
495 auto ilast = ret.nops()-1;
497 if(is_zero(
get_op(ret,ilast,0))) {
499 for(
int i=ilast-1; i>=0; i--) rep_xs.append(
get_op(ret,i,0));
501 for(
auto xi : delta) {
502 if(!rep_xs.has(xi) && !
get_op(ret,ilast,1).has(xi)) {
509 for(
auto xi : delta) {
510 if(!rep_xs.has(xi)) {
515 if(is_zero(xs0))
throw Error(
"Partilize: (xs0.is_zero())");
519 for(
int i=0; i<fe.op(2).nops(); i++) {
520 if(fe.op(2).op(i).has(xs0)) {
525 if(dlti<0)
throw Error(
"Partilize: dlti<0 found.");
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);
537 for(
int i=ilast; i>=0; i--) {
538 auto xi = ret.op(i).op(0);
540 auto yi = xi.subs(x(
w)==y(
w));
541 auto s = ret.op(i).op(1);
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);
549 for(
auto ss : ret) x2y.append(ss.op(0)==ss.op(1));
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);
561 fe.let_op(0).let_op(i) = num_den.op(0);
562 if(num_den.op(1)!=1) {
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);
574 if(idet_num_den.op(0)!=1) {
578 if(idet_num_den.op(1)!=1) {
585 for(
auto xi : delta) {
592 cout <<
"xcs = " << xcs << endl;
593 cout <<
"fe = " << fe_in << endl;
594 throw Error(
"Partilize: rm_xs = " +
ex2str(rm_xs) +
" & delta = " +
ex2str(delta));
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.");
607 auto ci = new_ft.coeff(xi);
609 xPos.append(xi); sPos += xi*ci;
611 xNeg.append(xi); sNeg -= xi*ci;
614 if(is_zero(sPos*sNeg))
throw Error(
"Partilize: sPos * sNeg = 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;
636 for(
auto item : bilst) {
637 if(!
isProjective(item,delta))
throw Error(
"Partilize: something is wrong here.");
639 ret_lst.push_back(item);
649 exvector fe_lst, ret_lst;
650 fe_lst.push_back(in_fe);
651 static int total_modes = 5;
654 for(
int i=0; i<fe_lst.size(); i++) {
656 auto ft =
get_op(fe, 0, 0);
659 if(!
get_op(fe, 1, 0).is_zero()) {
660 throw Error(
"Evaluate: (!get_op(fe, 1, 0).is_zero())");
664 if(fe.nops()<3 ||
xSign(ft)!=0) {
666 ret_lst.push_back(fe);
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;
691 for(
auto item : bilst) fe_lst2.push_back(item);
698 for(
auto item : bilst) fe_lst2.push_back(item);
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);
716 ret_lst.push_back(fe);
720 if(fe_lst2.size()<1)
break;
724 return exvector(std::move(ret_lst));
728 int maxn(
const ex pols,
const ex xi) {
730 for(
auto item : pols) {
731 auto cxn = item.degree(xi);
732 if(max_xn<cxn) max_xn = cxn;
744 exvector ret_vec, run_vec, run2_vec;
747 for(
auto fe : run_vec) {
748 ex ft = fe.op(0).op(1);
754 auto c0 = ftx.subs(xi==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);
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);
775 for(
auto xi3 : xs2) {
776 if(xi3==xi || xi3==xi2)
continue;
780 if(is_zero(xx))
continue;
783 for(
auto item2 : bfe) run2_vec.push_back(item2);
794 if(ftx.degree(xi)!=2)
continue;
796 auto c0 = ftx.coeff(xi,0);
797 auto c1 = ftx.coeff(xi,1);
798 auto c2 = ftx.coeff(xi,2);
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);
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);
819 for(
auto xi3 : xs2) {
820 if(xi3==xi || xi3==xi2)
continue;
824 if(is_zero(xx))
continue;
827 for(
auto item2 : bfe) run2_vec.push_back(item2);
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);
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);
853 for(
auto xi3 : xs2) {
854 if(xi3==xi || xi3==xi2)
continue;
858 if(is_zero(xx))
continue;
861 for(
auto item2 : bfe) run2_vec.push_back(item2);
871 for(
auto item : c0x) {
872 if(!is_a<numeric>(item) && !item.match(x(
w))) cc *= item;
876 if(
Factor(c0x).match(pow(
w,2))) {
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);
885 for(
auto xi3 : xs2) {
886 if(xi3==xi || xi3==xi2)
continue;
890 if(is_zero(xx))
continue;
893 for(
auto item2 : bfe) run2_vec.push_back(item2);
900 ret_vec.push_back(fe);
903 if(run2_vec.size()>0) {
908 return exvector(std::move(ret_vec));
class used to wrap error message
static bool isPartilizable(const ex ft, const ex delta, lst &xcs, int mode=0)
isPartilize w.r.t. F-term, generized to isLinearizable
static bool isLinearizable(const ex ft, const ex delta, lst &xcs)
Linearize w.r.t. F-term.
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,...
static exvector Binarize(ex const fe, ex const eqn)
Binarize the input fe, w.r.t. the linear eqn.
static void Projectivize(ex &fe, const ex delta, const ex xsum=0)
to Projectivize the input fe, the fe will be updated
static void Partilize(const lst xcs, const lst delta, const ex fe, exvector &ret_lst)
Partilize function, generized to Linearize.
static void Scalelize(ex &fe, const lst xs, const ex cy)
to Scalelize the input fe, the fe will be updated
static bool isProjective(const ex fe, const ex delta)
to check the input fe is Projective or NOT w.r.t. delta
static exvector Evaluate(const ex &fe)
ChengWu Internal, make sure ft in the first term, ONLY appear in ChengWu.cpp.
static void Linearize(const lst xcs, ex &fe, ex &ft)
Linearize w.r.t. xcs_in.
static exvector WickRotation(const exvector &fe_vec)
WickRotation, just check WRA exist or NOT to see successful or NOT. Still Experimental.
static lst XRefined_lst(ex const &ft)
static ex XRefined(ex const &ft)
void ChengWu(const ex &ft=0)
ChengWu for SecDec class.
namespace for Numerical integration with Sector Decomposition method
int x_free_index(ex expr)
exvector get_x_from(ex pol)
ex Factor(const ex expr_in)
ex expand_ex(ex const &expr_in, std::function< bool(const ex &)> has_func)
the expand like Mathematica
const char * Color_HighLight
void let_op_prepend(ex &ex_in, const ex item)
preppend item into expression
void let_op(ex &ex_in, int index1, int index2, const ex item)
update index1-th.index2-th of expression with item
ex collect_ex(ex const &expr_in, std::function< bool(const ex &)> has_func, int opt)
the collect function like Mathematica
ex get_op(const ex ex_in, int index1, int index2)
return index1-th.index2-th of expression
void let_op_append(ex &ex_in, const ex item)
append item into expression
string ex2str(const ex &expr)
convert ex to output string, the defalut printer format will be used
void let_op_remove_first(ex &ex_in)
remove first from expression
ex exfactor(const ex &expr, int opt)
factorize a expression
int xSign(ex const expr)
the always sign for expr
ex subs(const ex &e, init_list sl)