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