HepLib
Loading...
Searching...
No Matches
KillPowers.cpp
Go to the documentation of this file.
1
6#include "SD.h"
7#include <cmath>
8
9namespace HepLib::SD {
10
17 bool SecDec::KillPowerD(ex fe, int kpi) {
18 // check last element p-n is iWF(1)^1
19 auto nlast = fe.op(0).nops()-1;
20 if(fe.op(0).op(nlast)==iWF(1) && fe.op(1).op(nlast).is_zero()) {
21 FunExp.push_back(fe);
22 return false;
23 }
24 ex ft = fe.op(0).op(1);
25 lst fts = XRefined_lst(ft);
26 ex eqn;
27
28 //-----------------------------------------------------
29 bool ok2 = true;
30 for(auto ftitem : fts) {
31 auto xs = get_x_from(ftitem);
32 for(int i=0; i<xs.size(); i++) {
33 for(int j=i+1; j<xs.size(); j++) {
34 bool has_xij = false;
35 for(auto delta : ex_to<lst>(fe.op(2))) {
36 if(delta.has(xs[i]) && delta.has(xs[j])) {
37 has_xij = true;
38 break;
39 }
40 }
41 if(!has_xij) continue;
42
43 symbol xi, xj;
44 auto ftij = ftitem.subs(lst{xs[i]==xi, xs[j]==xj});
45 auto xs2 = get_x_from(ftij);
46 auto fts2 = XRefined_lst(ftij.subs(x(w)==0));
47
48 for(auto item : fts2) {
49 if(!item.has(xi) || !item.has(xj)) continue;
50 if(item.match(pow(w1,w2))) eqn = item.op(0).expand();
51 else eqn = item;
52 if(eqn.degree(xi)!=1 || eqn.degree(xj)!=1) continue;
53 ex ci = eqn.coeff(xi);
54 ex cj = eqn.coeff(xj);
55 if((ci*xi+cj*xj-eqn).is_zero() && is_a<numeric>(ci*cj) && (ci*cj)<0) {
56 eqn = eqn.subs(lst{xi==xs[i], xj==xs[j]});
57 ok2 = false;
58 goto OK2;
59 }
60 }
61 }
62 }
63 }
64
65 OK2:
66 if(!ok2) {
67 auto xij = get_x_from(eqn);
68 ex xi = xij[0];
69 ex xj = xij[1];
70 ex ci = eqn.coeff(xi);
71 ex cj = eqn.coeff(xj);
72
73 // handle eqn==ci xi - cj xj
74 if((ci*xi+cj*xj-eqn).is_zero() && is_a<numeric>(ci * cj) && (ci*cj)<0) {
75 if(Verbose>10) cout << PRE << "\\--" << Color_HighLight << "KillPowerD ["<<kpi<<"]: " << eqn << RESET << endl;
76 ci = abs(ci);
77 cj = abs(cj);
78 symbol yi,yj;
79 // Part I: ci xi-cj xj>0, i.e., xi>cj/ci xj
80 auto f1 = ex_to<lst>(fe.op(0));
81 auto e1 = ex_to<lst>(fe.op(1));
82 ex c1 = cj/ci;
83 for(int i=0; i<f1.nops(); i++) {
84 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});
85 }
86 if(e1.op(0)==1) f1.let_op(0) = f1.op(0)/(1+c1); // Jaccobi
87 else if(e1.op(1)==1) f1.let_op(1) = f1.op(1)/(1+c1);
88 else {
89 f1.append(1/(1+c1));
90 e1.append(1);
91 }
92 FunExp.push_back(lst{f1,e1,fe.op(2)});
93 // Part II: ci xi-cj xj<0, i.e., i.e., xj>ci/cj xi
94 auto f2 = ex_to<lst>(fe.op(0));
95 auto e2 = ex_to<lst>(fe.op(1));
96 ex c2 = ci/cj;
97 for(int i=0; i<f2.nops(); i++) {
98 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});
99 }
100 if(e2.op(0)==1) f2.let_op(0) = f2.op(0)/(1+c2); // Jaccobi
101 else if(e2.op(1)==1) f2.let_op(1) = f2.op(1)/(1+c2);
102 else {
103 f2.append(1/(1+c2));
104 e2.append(1);
105 }
106 FunExp.push_back(lst{f2,e2,fe.op(2)});
107 } else {
108 throw Error("KillPowerD: Not Expected Region.");
109 }
110 return true;
111 }
112
113 // when ok2 is true
114 let_op_append(fe, 0, iWF(1));
115 let_op_append(fe, 1, 0);
116 FunExp.push_back(fe);
117 return false;
118 }
119
127 bool SecDec::KillPower(ex fe, int kpi, int bits) {
128 // check last element p-n is iWF(1)^1
129 auto nlast = fe.op(0).nops()-1;
130 if(is_zero(fe.op(0).op(nlast)-iWF(1)) && fe.op(1).op(nlast).is_zero()) {
131 FunExp.push_back(fe);
132 return false;
133 }
134 ex ft = fe.op(0).op(1);
135 lst fts = XRefined_lst(ft);
136 ex eqn;
137
138 //-----------------------------------------------------
139 bool ok2 = true;
140 if(((bits/2)%2)==1) {
141 for(auto ftitem : fts) {
142 auto xs = get_x_from(ftitem);
143 for(int i=0; i<xs.size(); i++) {
144 for(int j=i+1; j<xs.size(); j++) {
145 symbol xi, xj;
146 auto ftij = ftitem.subs(lst{xs[i]==xi, xs[j]==xj});
147 auto xs2 = get_x_from(ftij);
148 for(int nn=0; nn<std::pow(2,xs2.size()); nn++) {
149 int tnn = nn;
150 lst xsubs;
151 for(int ni=0; ni<xs2.size(); ni++) {
152 if(tnn%2==1) xsubs.append(xs2[ni]==1);
153 else xsubs.append(xs2[ni]==0);
154 tnn /= 2;
155 }
156 auto fts2 = XRefined_lst(ftij.subs(xsubs));
157
158 int NN = 100;
159 for(auto item : fts2) {
160 if(!item.has(xi) || !item.has(xj)) continue;
161 if(item.match(pow(w1,w2)) || (item.degree(xi)==1 && item.degree(xj)==1)) {
162 if(item.match(pow(w1,w2))) eqn = item.op(0);
163 else eqn = item;
164 auto t1 = eqn.subs(lst{xi==1/ex(11), xj==1/ex(19)});
165 if(t1.is_zero()) t1 = eqn.subs(lst{xi==1/ex(3), xj==1/ex(23)});
166 if(t1.is_zero()) t1 = eqn.subs(lst{xi==1/ex(13), xj==1/ex(37)});
167
168 bool ook = true;
169 for(int ni=0; ni<=NN; ni++) {
170 for(int nj=0; nj<=NN; nj++) {
171 auto t2 = eqn.subs(lst{xi==ni/ex(NN), xj==nj/ex(NN)});
172 if(t1*t2 < 0) {
173 ook = false;
174 break;
175 }
176 }
177 if(!ook) break;
178 }
179 if(ook) continue;
180
181 if(item.match(pow(w1,w2)) && Verbose>0) {
182 if(eqn.degree(xi)>1 || eqn.degree(xj)>1 || eqn.coeff(xi).has(xj)) {
183 cout << WarnColor << PRE << "\\--Warning: Not handled with eqn1=" << eqn << RESET << endl;
184 continue;
185 }}
186
187 eqn = eqn.subs(lst{xi==xs[i], xj==xs[j]});
188 ok2 = false;
189 goto OK2;
190 }
191 }
192 }
193 }
194 }
195 }}
196
197 OK2:
198 if(!ok2) {
199 auto xij = get_x_from(eqn);
200 ex xi = xij[0];
201 ex xj = xij[1];
202
203 if(false)
204 if((eqn-(xi+xj-1)).is_zero() || (eqn+(xi+xj-1)).is_zero()) {
205 symbol xx;
206 auto f1 = fe.op(0);
207 auto e1 = fe.op(1);
208 for(int i=0; i<f1.nops(); i++) f1.let_op(i) = f1.op(i).subs(xi==1-xx).subs(xx==xi);
209 FunExp.push_back(lst{f1,e1});
210 return true;
211 }
212
213 ex ci = eqn.coeff(xi);
214 ex cj = eqn.coeff(xj);
215
216 // handle eqn==ci xi - cj xj
217 if((ci*xi+cj*xj-eqn).is_zero() && is_a<numeric>(ci * cj) && (ci*cj)<0) {
218 if(Verbose>10) cout << PRE << "\\--" << Color_HighLight << "KillPower ["<<kpi<<"]: " << eqn << RESET << endl;
219 ci = abs(ci);
220 cj = abs(cj);
221 if(is_zero(ci-cj)) {
222 symbol xx;
223 // Part I: xi<xj
224 auto f1 = ex_to<lst>(fe.op(0));
225 auto e1 = ex_to<lst>(fe.op(1));
226 for(int i=0; i<f1.nops(); i++) f1.let_op(i) = f1.op(i).subs(xi==xx*xj).subs(xx==xi);
227 f1.append(xj); // Jaccobi
228 e1.append(1);
229 FunExp.push_back(lst{f1,e1});
230 // Part II: xj<xi
231 auto f2 = ex_to<lst>(fe.op(0));
232 auto e2 = ex_to<lst>(fe.op(1));
233 for(int i=0; i<f2.nops(); i++) f2.let_op(i) = f2.op(i).subs(xj==xx*xi).subs(xx==xj);
234 f2.append(xi); // Jaccobi
235 e2.append(1);
236 FunExp.push_back(lst{f2,e2});
237 } else {
238 // we set c1 > c2 always
239 ex c1 = ci;
240 ex x1 = xi;
241 ex c2 = cj;
242 ex x2 = xj;
243 if(c1 < c2) {
244 c1 = cj;
245 x1 = xj;
246 c2 = ci;
247 x2 = xi;
248 }
249 symbol xx;
250 // Part I: 0<x2<1, c2/c1<x1<1
251 auto f1 = ex_to<lst>(fe.op(0));
252 auto e1 = ex_to<lst>(fe.op(1));
253 for(int i=0; i<f1.nops(); i++) f1.let_op(i) = f1.op(i).subs(x1==xx*(c1-c2)/c1+c2/c1).subs(xx==x1);
254 if(e1.op(0)==1) f1.let_op(0) = f1.op(0)*(c1-c2)/c1; // Jaccobi
255 else if(e1.op(1)==1) f1.let_op(1) = f1.op(1)*(c1-c2)/c1;
256 else {
257 f1.append((c1-c2)/c1);
258 e1.append(1);
259 }
260 FunExp.push_back(lst{f1,e1});
261 // Part II: x1>c2/c1 x2, i.e., x2<c1/c2 x1, 0<x1<c2/c1
262 auto f2 = ex_to<lst>(fe.op(0));
263 auto e2 = ex_to<lst>(fe.op(1));
264 for(int i=0; i<f2.nops(); i++) f2.let_op(i) = f2.op(i).subs(x1==xx*c2/c1).subs(xx==x1);
265 if(e2.op(0)==1) f2.let_op(0) = f2.op(0)*c2/c1; // Jaccobi
266 else if(e2.op(1)==1) f2.let_op(1) = f2.op(1)*c2/c1;
267 else {
268 f2.append(c2/c1);
269 e2.append(1);
270 }
271 // now x2<x1, 0<x1<1
272 for(int i=0; i<f2.nops(); i++) f2.let_op(i) = f2.op(i).subs(x2==xx*x1).subs(xx==x2);
273 f2.append(x1); // Jaccobi
274 e2.append(1);
275 FunExp.push_back(lst{f2,e2});
276 // Part III: x1<c2/c1 x2
277 auto f3 = ex_to<lst>(fe.op(0));
278 auto e3 = ex_to<lst>(fe.op(1));
279 for(int i=0; i<f3.nops(); i++) f3.let_op(i) = f3.op(i).subs(x1==xx*x2*c2/c1).subs(xx==x1);
280 if(e3.op(0)==1) f3.let_op(0) = f3.op(0)*c2/c1; // Jaccobi
281 else if(e3.op(1)==1) f3.let_op(1) = f3.op(1)*c2/c1;
282 else {
283 f3.append(c2/c1);
284 e3.append(1);
285 }
286 f3.append(x2);
287 e3.append(1);
288 FunExp.push_back(lst{f3,e3});
289 }
290 } else if( (eqn-(xi+xj-1)).is_zero() || (eqn+(xi+xj-1)).is_zero() ) {
291 if(Verbose>10) cout << PRE << "\\--" << Color_HighLight << "KillPower ["<<kpi<<"]: " << eqn << RESET << endl;
292 symbol xx, yy, zz;
293 // Part I: xi+xj-1>0
294 auto f1 = ex_to<lst>(fe.op(0));
295 auto e1 = ex_to<lst>(fe.op(1));
296 for(int i=0; i<f1.nops(); i++) f1.let_op(i) = f1.op(i).subs(xj==xx+1-xi).subs(xx==xj);
297 // now 0<xi<1, 0<xj<xi
298 for(int i=0; i<f1.nops(); i++) f1.let_op(i) = f1.op(i).subs(xj==xx*xi).subs(xx==xj);
299 f1.append(xi); // Jaccobi
300 e1.append(1);
301 FunExp.push_back(lst{f1,e1});
302 // Part IIa: 1-xi-xj>0, (a): xj<xi
303 auto f2 = ex_to<lst>(fe.op(0));
304 auto e2 = ex_to<lst>(fe.op(1));
305 for(int i=0; i<f2.nops(); i++) f2.let_op(i) = f2.op(i).subs(lst{xi==(1+zz)*yy/2,xj==(1-zz)*yy/2}).subs(lst{yy==xi, zz==xj});
306 f2.append(xi/2); // Jaccobi
307 e2.append(1);
308 FunExp.push_back(lst{f2,e2});
309 // Part IIb: 1-xi-xj>0, (a): xj>xi
310 auto f3 = ex_to<lst>(fe.op(0));
311 auto e3 = ex_to<lst>(fe.op(1));
312 for(int i=0; i<f3.nops(); i++) f3.let_op(i) = f3.op(i).subs(lst{xj==(1+zz)*yy/2,xi==(1-zz)*yy/2}).subs(lst{yy==xi, zz==xj});
313 f3.append(xi/2); // Jaccobi
314 e3.append(1);
315 FunExp.push_back(lst{f3,e3});
316 } else {
317 let_op_append(fe, 0, iWF(1));
318 let_op_append(fe, 1, 0);
319 FunExp.push_back(fe);
320 if(Verbose>0) cout << WarnColor << PRE << "\\--Warning: Not handled with eqn2=" << eqn << RESET << endl;
321 }
322 return true;
323 }
324
325 //-----------------------------------------------------
326 bool ok1 = true;
327 if((bits%2)==1) {
328 for(auto ftitem : fts) {
329 auto xs = get_x_from(ftitem);
330 for(int i=0; i<xs.size(); i++) {
331 symbol xi;
332 auto fti = ftitem.subs(lst{xs[i]==xi});
333 auto xs2 = get_x_from(fti);
334 for(int nn=0; nn<std::pow(2,xs2.size()); nn++) {
335 int tnn = nn;
336 lst xsubs;
337 for(int ni=0; ni<xs2.size(); ni++) {
338 if(tnn%2==1) xsubs.append(xs2[ni]==1);
339 else xsubs.append(xs2[ni]==0);
340 tnn /= 2;
341 }
342 auto fts2 = XRefined_lst(fti.subs(xsubs));
343
344 int NN = 100;
345 for(auto item : fts2) {
346 if((item.degree(xi)==1 || item.match(pow(w1,w2))) && item.has(xi)) {
347 if(item.match(pow(w1,w2))) eqn = item.op(0);
348 else if(xs2.size()<1) eqn = item;
349 else continue;
350 auto t1 = eqn.subs(lst{xi==1/ex(11)});
351 if(t1.is_zero()) t1 = eqn.subs(lst{xi==1/ex(3)});
352 if(t1.is_zero()) t1 = eqn.subs(lst{xi==1/ex(13)});
353
354 bool ook = true;
355 for(int ni=0; ni<=NN; ni++) {
356 auto t2 = eqn.subs(lst{xi==ni/ex(NN)});
357 if(t1*t2 < 0) {
358 ook = false;
359 break;
360 }
361 }
362 if(ook) continue;
363
364 if(Verbose>0 && eqn.degree(xi)>1) {
365 cout << WarnColor << PRE << "\\--Warning: Not handled with eqn3=" << eqn << RESET << endl;
366 continue;
367 }
368
369 if(!ook) {
370 eqn = eqn.subs(lst{xi==xs[i]});
371 ok1 = false;
372 goto OK1;
373 }
374 }
375 }
376 }
377 }
378 }}
379
380 OK1:
381 if(!ok1) {
382 auto xij = get_x_from(eqn);
383 ex xi = xij[0];
384 ex ci = eqn.coeff(xi);
385 ex c0 = eqn.subs(lst{xi==0});
386 // handle eqn==ci xi - c0
387 if((ci*xi+c0-eqn).is_zero() && is_a<numeric>(ci*c0) && (ci*c0)<0 && abs(c0)<abs(ci)) {
388 if(Verbose>10) cout << PRE << "\\--" << Color_HighLight << "KillPower ["<<kpi<<"]: " << eqn << RESET << endl;
389 ci = abs(ci);
390 c0 = abs(c0);
391 ex cc = c0/ci;
392 symbol xx;
393 // Part I: xi<cc
394 auto f1 = ex_to<lst>(fe.op(0));
395 auto e1 = ex_to<lst>(fe.op(1));
396 for(int i=0; i<f1.nops(); i++) f1.let_op(i) = f1.op(i).subs(xi==xx*cc).subs(xx==1-xi);
397 if(e1.op(0)==1) f1.let_op(0) = f1.op(0)*cc; // Jaccobi
398 else if(e1.op(1)==1) f1.let_op(1) = f1.op(1)*cc;
399 else {
400 f1.append(cc);
401 e1.append(1);
402 }
403 FunExp.push_back(lst{f1,e1});
404 // Part II: xi>cc
405 auto f2 = ex_to<lst>(fe.op(0));
406 auto e2 = ex_to<lst>(fe.op(1));
407 for(int i=0; i<f2.nops(); i++) f2.let_op(i) = f2.op(i).subs(xi==(1-cc)*xx+cc).subs(xx==xi);
408 if(e2.op(0)==1) f2.let_op(0) = f2.op(0)*(1-cc); // Jaccobi
409 else if(e2.op(1)==1) f2.let_op(1) = f2.op(1)*(1-cc);
410 else {
411 f2.append(1-cc);
412 e2.append(1);
413 }
414 FunExp.push_back(lst{f2,e2});
415 } else {
416 if(Verbose>0) cout << WarnColor << PRE << "\\--Warning: Not handled with eqn4=" << eqn << RESET << endl;
417 FunExp.push_back(fe);
418 }
419 return true;
420 }
421
422 // when ok2 && ok1 is true
423 let_op_append(fe, 0, iWF(1));
424 let_op_append(fe, 1, 0);
425 FunExp.push_back(fe);
426 return false;
427 }
428
433 void SecDec::KillPowers(int bits) {
434 MB(); // make sure MB first, if possible
435 int kpi = 0;
436 while(kpi<30) {
437 kpi++;
438 if(kpi>10 && (kpi % 10)==0) {
439 cout << "Warning: kip>10, (kpi=" << kpi << ") maybe a dead loop!" << endl;
440 }
441
442 exvector funexp;
443 for(auto fe : FunExp) {
444 funexp.push_back(fe);
445 }
446 FunExp.clear();
447 FunExp.shrink_to_fit();
448
449 bool repeat = false;
450 for(auto &fe : funexp) {
451 bool ret;
452 if(fe.nops()>2) ret = KillPowerD(fe, kpi);
453 else ret = KillPower(fe, kpi, bits);
454 if(ret) repeat = true;
455 }
456 if(!repeat) break;
457 }
458
459 for(auto &fe : FunExp) {
460 auto nn = fe.op(0).nops()-1;
461 if(is_zero(fe.op(0).op(nn)-iWF(1)) && fe.op(1).op(nn).is_zero()) {
462 let_op_remove_last(fe, 0);
463 let_op_remove_last(fe, 1);
464 }
465 }
466 Normalizes();
467 }
468
469
470}
RESET
#define RESET
Definition BASIC.h:79
SD.h
SecDec header file.
HepLib::SD::SecDec::XRefined_lst
static lst XRefined_lst(ex const &ft)
Definition Helpers.cpp:392
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::get_x_from
exvector get_x_from(ex pol)
Definition Helpers.cpp:29
HepLib::Color_HighLight
const char * Color_HighLight
Definition BASIC.cpp:248
HepLib::NN
ex NN(ex expr, int digits)
the nuerical evaluation
Definition BASIC.cpp:1278
HepLib::WarnColor
const char * WarnColor
Definition BASIC.cpp:247
HepLib::w
ex w
Definition Init.cpp:93
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::w1
ex w1
Definition BASIC.h:499
HepLib::let_op_remove_last
void let_op_remove_last(ex &ex_in)
remove last from expression
Definition BASIC.cpp:1345
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