HepLib
KillPowers.cpp
Go to the documentation of this file.
1 
6 #include "SD.h"
7 #include <cmath>
8 
9 namespace 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:90
HepLib::Verbose
int Verbose
Definition: Init.cpp:139
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:140
HepLib::subs
ex subs(const ex &e, init_list sl)
Definition: BASIC.h:51
Generated by doxygen 1.9.1