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