HepLib
Solver.cpp
Go to the documentation of this file.
1 
6 #include "BASIC.h"
7 #include "IBP.h"
8 #include <cmath>
9 
10 namespace HepLib {
11 
12  // return MR & T, such that MR = T.M0
13  static pair<matrix,matrix> RowReduce(matrix mat) {
14  Fermat &fermat = Fermat::get();
15  int &v_max = fermat.vmax;
16 
17  lst rep_vs;
18  ex tree = mat;
19  for(const_preorder_iterator i = tree.preorder_begin(); i != tree.preorder_end(); ++i) {
20  auto e = (*i);
21  if(is_a<symbol>(e) || e.match(a(w))) rep_vs.append(e);
22  }
23  rep_vs.sort();
24  rep_vs.unique();
25  //sort_lst(rep_vs);
26 
27  exmap v2f;
28  symtab st;
29  int fvi = 0;
30  for(auto vi : rep_vs) {
31  auto name = "v" + to_string(fvi);
32  v2f[vi] = Symbol(name);
33  st[name] = vi;
34  fvi++;
35  }
36 
37  stringstream ss;
38  if(fvi>111) {
39  cout << rep_vs << endl;
40  throw Error("Fermat: Too many variables.");
41  }
42  if(fvi>v_max) {
43  for(int i=v_max; i<fvi; i++) ss << "&(J=v" << i << ");" << endl;
44  fermat.Execute(ss.str());
45  ss.clear();
46  ss.str("");
47  v_max = fvi;
48  }
49 
50  int nrow = mat.rows();
51  int ncol = mat.cols();
52 
53  ss << "Array m[" << nrow << "," << ncol+1 << "];" << endl;
54  fermat.Execute(ss.str());
55  ss.clear();
56  ss.str("");
57 
58  ss << "[m]:=[(";
59  for(int c=0; c<ncol; c++) {
60  for(int r=0; r<nrow; r++) {
61  ss << mat(r,c).subs(iEpsilon==0).subs(v2f) << ",";
62  }
63  }
64  if(nrow>0) ss << "0";
65  for(int r=1; r<nrow; r++) ss << ",0";
66  ss << ")];" << endl;
67  fermat.Execute(ss.str());
68  ss.clear();
69  ss.str("");
70 
71  bool sparse = false;
72  if(nrow>1100) sparse = true;
73  if(sparse) fermat.Execute("Sparse([m]);");
74  fermat.Execute("Redrowech([m],[t]);");
75  ss << "&(U=1);" << endl; // ugly printing, the whitespace matters
76  matrix mr(nrow, ncol);
77  if(true) { // read matrix m
78  if(sparse) ss << "![m]" << endl;
79  else ss << "![m" << endl;
80  auto ostr = fermat.Execute(ss.str());
81  ss.clear();
82  ss.str("");
83 
84  // make sure last char is 0
85  if(ostr[ostr.length()-1]!='0') throw Error("RowReduce, last char is NOT 0.");
86  ostr = ostr.substr(0, ostr.length()-1);
87  string_trim(ostr);
88 
89  ostr.erase(0, ostr.find(":=")+2);
90  size_t sn = ostr.length();
91  char lc;
92  for(size_t i=0; i<sn; i++) {
93  char & c = ostr[i];
94  if(c=='[') {
95  c = '{';
96  if(sparse && i>0 && lc=='}') ostr[i-1] = ',';
97  } else if(c==']') c = '}';
98  else if(sparse && (c==' '||c=='\t'||c=='\n'||c=='\r')) continue;
99  lc = c;
100  }
101 
102  Parser fp(st);
103  auto res = fp.Read(ostr);
104  if(sparse) {
105  for(auto const & item : res) {
106  int r = -1;
107  for(auto const & it : item) {
108  if(r==-1) r = ex2int(it)-1;
109  else {
110  int c = ex2int(it.op(0))-1;
111  mr(r,c) = it.op(1);
112  }
113  }
114  }
115  } else {
116  for(int r=0; r<nrow; r++) {
117  auto cur = res.op(r);
118  for(int c=0; c<ncol; c++) mr(r,c) = cur.op(c);
119  }
120  }
121  }
122  matrix mt(nrow, nrow);
123  if(true) { // read matrix t
124  if(sparse) ss << "![t]" << endl;
125  else ss << "![t" << endl;
126  auto ostr = fermat.Execute(ss.str());
127  ss.clear();
128  ss.str("");
129 
130  // make sure last char is 0
131  if(ostr[ostr.length()-1]!='0') throw Error("Solver::Export, last char is NOT 0.");
132  ostr = ostr.substr(0, ostr.length()-1);
133  string_trim(ostr);
134 
135  ostr.erase(0, ostr.find(":=")+2);
136  size_t sn = ostr.length();
137  char lc;
138  for(size_t i=0; i<sn; i++) {
139  char & c = ostr[i];
140  if(c=='[') {
141  c = '{';
142  if(sparse && i>0 && lc=='}') ostr[i-1] = ',';
143  } else if(c==']') c = '}';
144  else if(sparse && (c==' '||c=='\t'||c=='\n'||c=='\r')) continue;
145  lc = c;
146  }
147  Parser fp(st);
148  auto res = fp.Read(ostr);
149  if(sparse) {
150  for(auto const & item : res) {
151  int r = -1;
152  for(auto const & it : item) {
153  if(r==-1) r = ex2int(it)-1;
154  else {
155  int c = ex2int(it.op(0))-1;
156  mt(r,c) = it.op(1);
157  }
158  }
159  }
160  } else {
161  for(int r=0; r<nrow; r++) {
162  auto cur = res.op(r);
163  for(int c=0; c<nrow; c++) mt(r,c) = cur.op(c);
164  }
165  }
166  }
167  // note the order, before exfactor (normal_fermat will be called again here)
168  ss << "&(U=0);" << endl; // disable ugly printing
169  ss << "@([m],[t]);" << endl;
170  ss << "&_G;" << endl;
171  fermat.Execute(ss.str());
172  ss.clear();
173  ss.str("");
174  return make_pair(mr,mt);
175  }
176 
177  // both row and column start from 0
178  typedef map<int,map<int,ex>> SparseMatrix;
179  static void RowReduce(SparseMatrix & smat, int pn) {
180  Fermat &fermat = Fermat::get();
181  int &v_max = fermat.vmax;
182 
183  lst rep_vs;
184  for(auto const & rv : smat) for(auto const & cv : rv.second) {
185  ex tree = cv.second;
186  for(const_preorder_iterator i = tree.preorder_begin(); i != tree.preorder_end(); ++i) {
187  auto e = (*i);
188  if(is_a<symbol>(e) || e.match(a(w))) rep_vs.append(e);
189  }
190  }
191  rep_vs.sort();
192  rep_vs.unique();
193  //sort_lst(rep_vs);
194 
195  exmap v2f;
196  symtab st;
197  int fvi = 0;
198  for(auto vi : rep_vs) {
199  auto name = "v" + to_string(fvi);
200  v2f[vi] = Symbol(name);
201  st[name] = vi;
202  fvi++;
203  }
204 
205  stringstream ss;
206  if(fvi>111) {
207  cout << rep_vs << endl;
208  throw Error("Fermat: Too many variables.");
209  }
210  if(fvi>v_max) {
211  for(int i=v_max; i<fvi; i++) ss << "&(J=v" << i << ");" << endl;
212  fermat.Execute(ss.str());
213  ss.clear();
214  ss.str("");
215  v_max = fvi;
216  }
217 
218  ostringstream oss;
219  int nrow=-1, ncol=-1;
220  int rcur = 0, rtot = smat.size();
221  oss << "[";
222  for(auto const & rv : smat) {
223  rcur++;
224  if(rv.second.size()<1) continue;
225  int r = rv.first;
226  if(r>nrow) nrow = r;
227  oss << "[" << r+1 << ",";
228  int ccur = 0, ctot = rv.second.size();
229  for(auto const & cv : rv.second) {
230  ccur++;
231  int c = cv.first;
232  if(c>ncol) ncol = c;
233  oss << "[" << c+1 << "," << cv.second.subs(iEpsilon==0).subs(v2f) << "]";
234  if(ccur!=ctot) oss << ",";
235  else oss << "]";
236  }
237  if(rcur!=rtot) oss << "`" << endl;
238  else oss << "];" << endl;
239  }
240  nrow = nrow+1;
241  ncol = ncol+1;
242 
243  ss << "Array m[" << nrow << "," << ncol+1 << "] Sparse;" << endl;
244  fermat.Execute(ss.str());
245  ss.clear();
246  ss.str("");
247 
248  ss << "[m]:=" << oss.str();
249  fermat.Execute(ss.str());
250  ss.clear();
251  ss.str("");
252 
253  if(pn>0) fermat.Execute("Redrowech([m],,,"+to_string(pn)+");");
254  else fermat.Execute("Redrowech([m]);");
255  ss << "&(U=1);" << endl; // ugly printing, the whitespace matters
256  smat.clear();
257  if(true) { // read matrix m
258  ss << "![m]" << endl;
259  auto ostr = fermat.Execute(ss.str());
260  ss.clear();
261  ss.str("");
262 
263  // make sure last char is 0
264  if(ostr[ostr.length()-1]!='0') throw Error("RowReduce, last char is NOT 0.");
265  ostr = ostr.substr(0, ostr.length()-1);
266  string_trim(ostr);
267 
268  ostr.erase(0, ostr.find(":=")+2);
269  size_t sn = ostr.length();
270  char lc;
271  for(size_t i=0; i<sn; i++) {
272  char & c = ostr[i];
273  if(c=='[') {
274  c = '{';
275  if(i>0 && lc=='}') ostr[i-1] = ',';
276  } else if(c==']') c = '}';
277  else if(c==' '||c=='\t'||c=='\n'||c=='\r') continue;
278  lc = c;
279  }
280 
281  Parser fp(st);
282  auto res = fp.Read(ostr);
283  for(auto const & item : res) {
284  int r = -1;
285  for(auto const & it : item) {
286  if(r==-1) r = ex2int(it)-1;
287  else {
288  int c = ex2int(it.op(0))-1;
289  smat[r][c] = it.op(1);
290  }
291  }
292  }
293  }
294  // note the order, before exfactor (normal_fermat will be called again here)
295  ss << "&(U=0);" << endl; // disable ugly printing
296  ss << "@([m]);" << endl;
297  ss << "&_G;" << endl;
298  fermat.Execute(ss.str());
299  ss.clear();
300  ss.str("");
301  }
302 
303  void Solver::SolveSparse(const ex & sector, const map<int,int> & n2n) {
304  int pdim = Propagator.nops();
305  lst ibps = IBPs;
306 
307  // add more eqn to ibps
308  if(true) {
309  if(true) {
310  int n = ibps.nops();
311  for(int i=0; i<pdim; i++) {
312  for(int j1=0; j1<n; j1++) for(int j2=0; j2<n; j2++) {
313  ibps.append(ibps.op(i).subs(lst{a(j1)==a(j1)+1}).subs(lst{a(j2)==a(j2)-1}));
314 // ibps.append(ibps.op(i).subs(lst{a(j1)==a(j1)+1}).subs(lst{a(j2)==a(j2)+1}));
315 // ibps.append(ibps.op(i).subs(lst{a(j1)==a(j1)-1}).subs(lst{a(j2)==a(j2)-1}));
316  }
317  }
318  } else {
319  int n = ibps.nops();
320  for(int i=0; i<pdim; i++) {
321  for(int j=0; j<n; j++) ibps.append(ibps.op(j).subs(a(i)==a(i)+1));
322 // for(int j=0; j<n; j++) ibps.append(ibps.op(j).subs(a(i)==a(i)+2));
323 // for(int j=0; j<n; j++) ibps.append(ibps.op(j).subs(a(i)==a(i)+3));
324 // for(int j=0; j<n; j++) ibps.append(ibps.op(j).subs(a(i)==a(i)+4));
325  }
326 
327  n = ibps.nops();
328  for(int i=0; i<pdim; i++) {
329  for(int j=0; j<n; j++) ibps.append(ibps.op(j).subs(a(i)==a(i)-1));
330 // for(int j=0; j<n; j++) ibps.append(ibps.op(j).subs(a(i)==a(i)-2));
331 // for(int j=0; j<n; j++) ibps.append(ibps.op(j).subs(a(i)==a(i)-3));
332 // for(int j=0; j<n; j++) ibps.append(ibps.op(j).subs(a(i)==a(i)-4));
333  }
334  }
335  } else {
336  exset fs;
337  find(ibps, F(w), fs);
338  int n = ibps.nops();
339  for(auto fi : fs) {
340  exmap a2a;
341  for(int i=0; i<fi.op(0).nops(); i++) a2a[a(i)] = fi.op(0).op(i);
342  for(int j=0; j<n; j++) ibps.append(ibps.op(j).subs(a2a));
343  }
344  for(auto fi : fs) {
345  exmap a2a;
346  for(int i=0; i<fi.op(0).nops(); i++) a2a[a(i)] = fi.op(0).op(i)+1;
347  for(int j=0; j<n; j++) ibps.append(ibps.op(j).subs(a2a));
348  }
349  for(auto fi : fs) {
350  exmap a2a;
351  for(int i=0; i<fi.op(0).nops(); i++) a2a[a(i)] = fi.op(0).op(i)-1;
352  for(int j=0; j<n; j++) ibps.append(ibps.op(j).subs(a2a));
353  }
354  }
355 
356  // some a's replaced by integer
357  vector<int> nfix(pdim);
358  for(int i=0; i<pdim; i++) nfix[i] = 0;
359  exmap a2n;
360  for(auto const & kv : n2n) {
361  nfix[kv.first] = 1;
362  int k = kv.first;
363  if(k<0) k = pdim + k;
364  a2n[a(k)] = kv.second;
365  }
366  if(a2n.size()>0) {
367  int nibps = ibps.nops();
368  for(int i=0; i<nibps; i++) ibps.let_op(i) = ibps.op(i).subs(a2n);
369  }
370 
371  lst cons;
372  map<ex,lst,ex_is_less> cons_vec;
373  while(ibps.nops()>0) {
374  ibps.sort();
375  ibps.unique();
376 
377  exset fset;
378  find(ibps,F(w),fset);
379  map<int,exvector> fmap_vec;
380  for(auto fi : fset) { // sort fset
381  int psum = 0, nsum = 0;
382  auto ns = fi.op(0);
383  for(int i=0; i<pdim; i++) {
384  auto ni = ns.op(i).subs(a(i)==0,nopat);
385  if(sector.op(i)==1) {
386  if(nfix[i]==1 && ni<=0) psum -= 100000;
387  else psum += ex2int(ni)-1; // -1 from FIRE
388  } else {
389  if(nfix[i]==1 && ni>0) nsum += ex2int(ni);
390  else nsum -= ex2int(ni);
391  }
392  }
393  int key = psum+nsum;
394  fmap_vec[-key].push_back(fi); // - to reverse the order
395  }
396 
397  map<ex,int,ex_is_less> f2n;
398  int next_gi = 0;
399  exvector fvec(fset.size());
400  int nr = 1, nc = 0;
401  for(auto const & ev : fmap_vec) {
402  if(next_gi<1) next_gi = ev.second.size();
403  for(auto const & v : ev.second) {
404  fvec[nc] = v;
405  f2n[v] = nc;
406  nc++;
407  }
408  nr++;
409  }
410  nr = ibps.nops();
411 
412  if(Verbose>0) {
413  cout << "+----------------------------------------" << endl;
414  cout << "| IBPs: " << nr << ", Fs: " << next_gi << "/" << nc << endl;
415  cout << "+----------------------------------------" << endl;
416  }
417 
418  SparseMatrix smat;
419  for(int r=0; r<nr; r++) {
420  auto ibp = ibps.op(r);
421  auto cvs = collect_lst(ibp,F(w));
422  for(auto & cv : cvs) {
423  int c = f2n[cv.op(1)];
424  ex cc = cv.op(0);
425  if(!cc.is_zero()) smat[r][c] = cc;
426  }
427  }
428  RowReduce(smat,next_gi);
429 
430  ibps.remove_all();
431  for(auto const & rv : smat) {
432  int next_sol_ibp = -1; // -1:init, 0:next, 1:sol, 2:ibp
433  ex sol = 0, ns, fc;
434  exvector nv(pdim);
435  exmap amap;
436  for(auto const & cv : rv.second) {
437  if(next_sol_ibp==-1) {
438  if(cv.first<next_gi) {
439  next_sol_ibp = 0;
440  if(!is_zero(cv.second-1)) throw Error("Solver::SolveSparse, NOT 1.");
441  ns = fvec[cv.first].op(0).subs(a(w)==0);
442  for(int i=0; i<pdim; i++) {
443  if(nfix[i]==0) amap[a(i)] = a(i)-ns.op(i);
444  }
445  fc = fvec[cv.first].subs(amap,nopat);
446  for(int i=0; i<pdim; i++) {
447  if(nfix[i]==1) {
448  nv[i] = ns.op(i);
449  if(sector.op(i)==1 && nv[i]<1) goto next_row; // subsector
450  else if(sector.op(i)!=1 && nv[i]>0) goto next_row; // out of this sector
451  } else nv[i] = (sector.op(i)==1 ? 1 : 0);
452  }
453  continue; // need continue
454  } else next_sol_ibp = 2; // ibp
455  }
456 
457  if(next_sol_ibp==0) {
458  if(cv.first<next_gi) goto next_row; // not a solution or ibp
459  next_sol_ibp = 1; // solution found
460  }
461 
462  if(true) {
463  ex f = fvec[cv.first];
464  ex nd = cv.second;
465  if(next_sol_ibp==1) {
466  f = f.subs(amap,nopat);
467  nd = nd.subs(amap,nopat);
468  ex num = factor_flint(nd.numer());
469  ex den = factor_flint(nd.denom());
470  ns = f.op(0).subs(a(w)==0); // reuse ns here
471  for(int i=0; i<pdim; i++) {
472  if(sector.op(i)!=1) {
473  ex nsi = ns.op(i);
474  if(nfix[i]==1) {
475  if(nsi>0) goto next_row; // out of this sector
476  else continue;
477  }
478  nsi = -nsi; // note '-' here, ni<=nsi
479  // due to possible (ai+nsi-1)*F(...,ai-nsi,...), so nsi+1 may be OK
480  while(true) {
481  ex cci = num.subs(a(i)==1+nsi);
482  if(!cci.is_zero() || nsi>=nv[i]) break;
483  nsi++;
484  }
485  if(nsi<nv[i]) nv[i] = nsi;
486  }
487  }
488 
489  if(!is_a<mul>(den)) den = lst{ den };
490  for(const auto & item : den) {
491  ex it = item;
492  if(is_a<power>(it)) it = it.op(0);
493  if(!it.has(a(w)) || it.has(d)) continue;
494  exset as;
495  find(it,a(w),as);
496  if(as.size()!=1) {
497  cout << "size is NOT 1: " << it << endl;
498  throw Error("Solver::SolveSparse, error occured.");
499  }
500  ex aw = *(as.begin());
501  if(it.degree(aw)!=1) {
502  cout << "degree is NOT 1: " << it << endl;
503  throw Error("Solver::SolveSparse, error occured.");
504  }
505  ex a0 = -it.coeff(aw,0)/it.coeff(aw,1);
506  int an = ex2int(aw.op(0));
507  if(sector.op(an)!=1) {
508  if(a0<=nv[an]) nv[an] = a0-1;
509  } else {
510  if(a0>=nv[an]) nv[an] = a0+1;
511  }
512  }
513  }
514  sol -= nd * f;
515  }
516  }
517 
518  if(next_sol_ibp==1) {
519  for(int i=0; i<pdim; i++) { // check
520  if(nfix[i]==1 && fc.op(0).op(i)!=nv[i]) {
521  cout << i << ": " << fc << endl << nv << endl;
522  throw Error("Solver::SolveSparse, error occured.");
523  } else if(nfix[i]==0 && fc.op(0).op(i)!=a(i)) {
524  cout << i << ": " << fc << endl << nv << endl;
525  throw Error("Solver::SolveSparse, error occured.");
526  }
527  }
528 
529  if(true) {
530  int psum = 0, nsum = 0;
531  for(int i=0; i<pdim; i++) {
532  if(nfix[i]==1) continue;
533  if(sector.op(i)==1 && !nv[i].is_equal(1)) psum++;
534  if(sector.op(i)!=1 && !nv[i].is_zero()) nsum++;
535  }
536  if(psum<=1 && nsum<=1) {
537  ex con = vec2lst(nv);
538  cons.append(con);
539  cons_vec[con].append(sol);
540  }
541  }
542  } else ibps.append(sol);
543 
544  next_row: ;
545  }
546 cons.sort();
547 cons.unique();
548 sort_lst(cons,false);
549 cout << cons << endl << endl;
550 for(auto & item : cons) {
551  if(item==lst{1,1,1,1,1,1,1,1,1,1,0,0}) {
552  cout << endl << cons_vec[item] << endl << endl;
553  exit(0);
554  }
555 }
556  }
557 
558  cons.sort();
559  cons.unique();
560  sort_lst(cons,false);
561 
562 // cout << endl;
563 // for(auto & k : cons) {
564 // cons_vec[k].sort();
565 // cons_vec[k].unique();
566 // cout << k << endl;
567 // cout << cons_vec[k] << endl << endl;
568 // }
569 // cout << endl << cons << endl;
570 
571 
572  }
573 
574  void Solver::IBP() {
575  int pdim = Propagator.nops();
576  lst InExternal;
577  for(auto ii : Internal) InExternal.append(ii);
578  for(auto ii : External) InExternal.append(ii);
579 
580  if(ISP.nops()<1) {
581  for(auto it : Internal) {
582  for(auto ii : InExternal) ISP.append(it*ii);
583  }
584  ISP.sort();
585  ISP.unique();
586  }
587 
588  if(ISP.nops() > pdim) {
589  cout << "ISP = " << ISP << endl;
590  cout << "Propagator = " << Propagator << endl;
591  throw Error("Solver::IBP: #(ISP) > #(Propagator).");
592  }
593 
594  lst sp2s, s2sp, ss;
595  int _pic=0;
596  for(auto item : ISP) {
597  _pic++;
598  Symbol si("P"+to_string(_pic));
599  ss.append(si);
600  sp2s.append(w*item==w*si);
601  sp2s.append(item==si);
602  s2sp.append(si==item);
603  }
604 
605  lst leqns;
606  for(int i=0; i<ISP.nops(); i++) { // note NOT pdim
607  auto eq = Propagator.op(i).expand().subs(iEpsilon==0); // drop iEpsilon
608  eq = eq.subs(sp2s);
609  eq = eq.subs(Replacement);
610  if(eq.has(iWF(w))) throw Error("Solver::IBP, iWF used in eq.");
611  leqns.append(eq == iWF(i));
612  }
613  auto s2p = lsolve(leqns, ss);
614  if(s2p.nops() != ISP.nops()) throw Error("Solver::IBP, lsolve failed.");
615 
616  // 1st version
617  if(true)
618  if(DSP.nops()<1) {
619  for(auto p1 : Internal)
620  for(auto p2 : InExternal)
621  DSP.append(lst{p1,p2});
622  }
623 
624  // Lee version
625  if(false)
626  if(DSP.nops()<1) {
627  for(int i=0; i<Internal.nops(); i++)
628  DSP.append(lst{Internal.op(i), Internal.op(i==Internal.nops()-1 ? 0 : i+1)});
629  for(auto p2 : InExternal) DSP.append(lst{Internal.op(0), p2});
630  DSP.append(lst{Internal, Internal});
631  }
632 
633  IBPs.remove_all();
634  lst nsa;
635  for(int i=0; i<pdim; i++) nsa.append(a(i));
636  for(auto sp : DSP) {
637  auto ilp = sp.op(0);
638  auto iep = sp.op(1);
639 
640  ex ibp = 0;
641  symbol ss;
642  for(int i=0; i<pdim; i++) {
643  auto ns = nsa;
644  ns.let_op(i) = nsa.op(i) + 1;
645  auto dp = Propagator.op(i).subs(ilp==ss).diff(ss).subs(ss==ilp);
646  ibp -= (a(i)+Shift[i]) * F(ns) * dp;
647  }
648 
649  ibp = ibp * iep;
650  ibp = ibp.expand();
651  ibp = ibp.subs(sp2s);
652  ibp = ibp.subs(Replacement);
653  ibp = ibp.subs(s2p);
654 
655  ex res = 0;
656  for(int i=0; i<pdim; i++) {
657  auto ci = ibp.coeff(iWF(i), 1);
658  ci = MapFunction([i](const ex &e, MapFunction &self)->ex {
659  if(!e.has(F(w))) return e;
660  else if(e.match(F(w))) {
661  lst tmp = ex_to<lst>(e.op(0));
662  tmp.let_op(i) = tmp.op(i)-1;
663  return F(tmp);
664  } else return e.map(self);
665  })(ci);
666  res += ci;
667  }
668  res += ibp.subs(lst{iWF(w)==0});
669  auto cilp = iep.coeff(ilp);
670  if(!is_zero(cilp)) res += d*cilp*F(nsa);
671  IBPs.append(res);
672  }
673  }
674 
675  void Solver::Solve(const ex & sector, const map<int,int> & n2n) {
676  int pdim = Propagator.nops();
677  lst ibps = IBPs;
678 
679  // add more eqn to ibps
680  if(false) {
681 // int n = ibps.nops();
682 // for(int i=0; i<pdim; i++) {
683 // for(int j1=0; j1<n; j1++) for(int j2=0; j2<n; j2++)
684 // ibps.append(ibps.op(i).subs(lst{a(j1)==a(j1)+1}).subs(lst{a(j2)==a(j2)-1}));
685 // }
686 
687  int n = ibps.nops();
688  for(int i=0; i<pdim; i++) {
689  for(int j=0; j<n; j++) ibps.append(ibps.op(j).subs(a(i)==a(i)+1));
690 // for(int j=0; j<n; j++) ibps.append(ibps.op(j).subs(a(i)==a(i)+2));
691 // for(int j=0; j<n; j++) ibps.append(ibps.op(j).subs(a(i)==a(i)+3));
692 // for(int j=0; j<n; j++) ibps.append(ibps.op(j).subs(a(i)==a(i)+4));
693  }
694 
695  //n = ibps.nops();
696  for(int i=0; i<pdim; i++) {
697  for(int j=0; j<n; j++) ibps.append(ibps.op(j).subs(a(i)==a(i)-1));
698 // for(int j=0; j<n; j++) ibps.append(ibps.op(j).subs(a(i)==a(i)-2));
699 // for(int j=0; j<n; j++) ibps.append(ibps.op(j).subs(a(i)==a(i)-3));
700 // for(int j=0; j<n; j++) ibps.append(ibps.op(j).subs(a(i)==a(i)-4));
701  }
702  } else {
703  exset fs;
704  find(ibps, F(w1,w2), fs);
705 cout << fs << endl;
706 exit(0);
707  }
708 
709  // some a's replaced by integer
710  vector<int> nfix(pdim);
711  for(int i=0; i<pdim; i++) nfix[i] = 0;
712  for(auto const & kv : n2n) nfix[kv.first] = 1;
713  exmap a2n;
714  for(auto const & kv : n2n) {
715  int k = kv.first;
716  if(k<0) k = pdim + k;
717  a2n[a(k)] = kv.second;
718  }
719  int nibps = ibps.nops();
720  for(int i=0; i<nibps; i++) ibps.let_op(i) = ibps.op(i).subs(a2n);
721 
722  lst cons;
723  map<ex,lst,ex_is_less> cons_vec;
724  while(ibps.nops()>0) {
725  ibps.sort();
726  ibps.unique();
727 
728  exvector fvec;
729  exset fset;
730  find(ibps,F(w),fset);
731  int nmax = -10000000, pmax = -10000000;
732  for(auto fi : fset) { // sort fset
733  int psum = 0, nsum = 0;
734  auto ns = fi.op(0);
735  for(int i=0; i<pdim; i++) {
736  auto ni = ns.op(i).subs(a(i)==0,nopat);
737  if(sector.op(i)==1) {
738  if(nfix[i]==1 && ni<=0) psum -= 100000;
739  else psum += ex2int(ni)-1; // -1 from FIRE
740  } else {
741  if(nfix[i]==1 && ni>0) nsum += ex2int(ni);
742  else nsum -= ex2int(ni);
743  }
744  }
745 
746  if(psum+nsum<pmax+nmax) continue;
747  if(psum+nsum>pmax+nmax) {
748  pmax = psum;
749  nmax = nsum;
750  fvec.clear();
751  }
752 
753  fvec.push_back(fi);
754  }
755  //sort_vec(fvec);
756 
757  int nr = ibps.nops();
758  int nc = fvec.size();
759  matrix mat(nr, nc);
760 
761  if(Verbose>0) {
762  cout << "+----------------------------------------" << endl;
763  cout << "| IBPs: " << nr << ", Fs: " << nc << endl;
764  cout << "+----------------------------------------" << endl;
765  }
766 
767  matrix bvec(nr,1);
768  for(int r=0; r<nr; r++) {
769  auto ibp = ibps.op(r);
770  for(int c=0; c<nc; c++) {
771  mat(r,c) = ibp;
772  ibp = ibp.subs(fvec[c]==0,nopat);
773  mat(r,c) = (mat(r,c)-ibp).coeff(fvec[c]);
774  }
775  bvec(r,0) = ibp;
776  }
777 
778  // note that: mr = mt.mul(mat)
779  auto rt = RowReduce(mat);
780  auto mr = rt.first;
781  auto mt = rt.second;
782  bvec = mt.mul(bvec);
783  for(int r=0; r<nr; r++) bvec(r,0) = collect_ex(bvec(r,0),F(w),o_flintfD);
784 
785  //cout << "CHECK: 1=" << ex_to<matrix>(normal(mr.sub(mt.mul(mat)))).is_zero_matrix() << endl;
786 
787  auto ibp_con_sol_vec = GiNaC_Parallel(nr, [&](const int r)->ex {
788  ex ibp=0, con=0;
789  int cc = -1;
790  for(int c=0; c<nc; c++) {
791  if(mr(r,c).is_zero()) continue;
792  if(mr(r,c)==1 && cc==-1) cc = c;
793  else goto next_row; // not a solution, skip
794  }
795 
796  if(cc!=-1) { // a solution found N -> N-1, N-2, subsectors
797  auto ns = fvec[cc].op(0).subs(a(w)==0);
798  exmap amap;
799  for(int i=0; i<pdim; i++) {
800  if(nfix[i]==0) amap[a(i)] = a(i)-ns.op(i);
801  }
802  fvec[cc] = fvec[cc].subs(amap,nopat);
803  auto sol = bvec(r,0).subs(amap,nopat);
804  sol = normal_flint(sol,o_flintfD);
805 
806  exvector nv(pdim);
807  for(int i=0; i<pdim; i++) {
808  if(nfix[i]==1) {
809  nv[i] = ns.op(i);
810  if(sector.op(i)==1 && nv[i]<1) goto next_row; // subsector
811  else if(sector.op(i)!=1 && nv[i]>0) goto next_row; // out of this sector
812  } else nv[i] = (sector.op(i)==1 ? 1 : 0);
813  }
814 
815  ex num = sol.numer();
816  exset fset;
817  find(sol,F(w),fset);
818  for(auto const & f : fset) { // make sure not to increase the sector
819  ex cc = factor_flint(num.coeff(f));
820  ns = f.op(0).subs(a(w)==0); // reuse ns here
821  for(int i=0; i<pdim; i++) {
822  if(sector.op(i)!=1) {
823  ex nsi = ns.op(i);
824  if(nfix[i]==1) {
825  if(nsi>0) goto next_row; // out of this sector
826  else continue;
827  }
828  nsi = -nsi; // note '-' here, ni<=nsi
829  // due to possible (ai+nsi-1)*F(...,ai-nsi,...), so nsi+1 may be OK
830  while(true) {
831  ex cci = cc.subs(a(i)==1+nsi);
832  if(!cci.is_zero() || nsi>=nv[i]) break;
833  nsi++;
834  }
835  if(nsi<nv[i]) nv[i] = nsi;
836  }
837  }
838  }
839 
840  ex den = sol.denom();
841  if(!is_a<mul>(den)) den = lst{den};
842  for(const auto & item : den) {
843  ex it = item;
844  if(is_a<power>(it)) it = it.op(0);
845  if(!it.has(a(w)) || it.has(d)) continue;
846  exset as;
847  find(it,a(w),as);
848  if(as.size()!=1) {
849  cout << "size is NOT 1: " << it << endl;
850  throw Error("Solver::Solve, error occured.");
851  }
852  ex aw = *(as.begin());
853  if(it.degree(aw)!=1) {
854  cout << "degree is NOT 1: " << it << endl;
855  throw Error("Solver::Solve, error occured.");
856  }
857  ex a0 = -it.coeff(aw,0)/it.coeff(aw,1);
858  int an = ex2int(aw.op(0));
859  if(sector.op(an)!=1) {
860  if(a0<=nv[an]) nv[an] = a0-1;
861  } else {
862  if(a0>=nv[an]) nv[an] = a0+1;
863  }
864  }
865 
866  for(int i=0; i<pdim; i++) { // check
867  if(nfix[i]==1 && fvec[cc].op(0).op(i)!=nv[i]) {
868  cout << fvec << endl << nv << endl;
869  throw Error("Solver::Solve, error occured.");
870  } else if(nfix[i]==0 && fvec[cc].op(0).op(i)!=a(i)) {
871  cout << fvec << endl << nv << endl;
872  throw Error("Solver::Solve, error occured.");
873  }
874  }
875 
876  //cout << "corner: " << nv << endl;
877  int psum = 0, nsum = 0;
878  for(int i=0; i<pdim; i++) {
879  if(nfix[i]==1) continue;
880  if(sector.op(i)==1 && !nv[i].is_equal(1)) psum++;
881  if(sector.op(i)!=1 && !nv[i].is_zero()) nsum++;
882  }
883  if(nsum<=1 && psum<=1) con = lst{ vec2lst(nv), sol };
884  } else if(!bvec(r,0).is_zero()) ibp = bvec(r,0).numer();
885  next_row: ;
886  return lst{ibp, con};
887  });
888 
889  ibps.remove_all();
890  for(int r=0; r<nr; r++) {
891  auto const & item = ibp_con_sol_vec[r];
892  if(!item.op(0).is_zero()) ibps.append(item.op(0));
893  if(is_a<lst>(item.op(1))) {
894  cons.append(item.op(1).op(0));
895  cons_vec[item.op(1).op(0)].append(item.op(1).op(1));
896  }
897  }
898 
899 cons.sort();
900 cons.unique();
901 cout << cons << endl << endl;
902  }
903 
904  cons.sort();
905  cons.unique();
906  sort_lst(cons,false);
907 
908  cout << endl;
909  for(auto & k : cons) {
910  cons_vec[k].sort();
911  cons_vec[k].unique();
912  cout << k << endl;
913  cout << cons_vec[k] << endl << endl;
914  }
915  cout << endl << cons << endl;
916  }
917 
921  void Solver::Export() {
922 
923  }
924 
925 
929  void Solver::Run() {
930 
931  }
932 
936  void Solver::Import() {
937 
938  }
939 
940 }
a
int * a
Definition: Functions.cpp:234
BASIC.h
Basic header file.
IBP.h
IBP header file.
HepLib::Error
class used to wrap error message
Definition: BASIC.h:242
HepLib::Fermat
interface to communicate with Fermat program
Definition: BASIC.h:797
HepLib::Fermat::get
static Fermat & get()
Definition: Fermat.cpp:7
HepLib::Fermat::Execute
string Execute(string)
Definition: Process.cpp:102
HepLib::IBP::Internal
lst Internal
Definition: IBP.h:27
HepLib::IBP::Shift
map< int, ex > Shift
Definition: IBP.h:36
HepLib::IBP::ISP
lst ISP
Definition: IBP.h:34
HepLib::IBP::Replacement
lst Replacement
Definition: IBP.h:29
HepLib::IBP::DSP
lst DSP
Definition: IBP.h:33
HepLib::IBP::External
lst External
Definition: IBP.h:28
HepLib::IBP::Propagator
lst Propagator
Definition: IBP.h:30
HepLib::MapFunction
class to wrap map_function of GiNaC
Definition: BASIC.h:632
HepLib::Parser
class to parse for string or file, helpful with Symbol class
Definition: BASIC.h:645
HepLib::Parser::Read
ex Read(const string &instr, bool s2S=true)
Definition: Functions.cpp:111
HepLib::Solver::Run
void Run() override
invoke kira program for reduction
Definition: Solver.cpp:929
HepLib::Solver::IBP
void IBP()
Definition: Solver.cpp:574
HepLib::Solver::Solve
void Solve(const ex &sector, const map< int, int > &a2n={})
Definition: Solver.cpp:675
HepLib::Solver::Import
void Import() override
import kira result
Definition: Solver.cpp:936
HepLib::Solver::Export
void Export() override
Export input data for KIRA.
Definition: Solver.cpp:921
HepLib::Solver::SolveSparse
void SolveSparse(const ex &sector, const map< int, int > &a2n={})
Definition: Solver.cpp:303
HepLib::Symbol
class extended to GiNaC symbol class, represent a positive symbol
Definition: BASIC.h:113
HepLib
HepLib namespace.
Definition: BASIC.cpp:17
HepLib::iEpsilon
const iSymbol iEpsilon
HepLib::sp
ex sp(const ex &a, const ex &b)
translated the vector dot a.b to a*b, useful in SecDec
Definition: Pair.cpp:237
HepLib::factor_flint
ex factor_flint(const ex &e, bool nd=true)
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:90
HepLib::d
const Symbol d
HepLib::SparseMatrix
map< int, map< int, ex > > SparseMatrix
Definition: Solver.cpp:178
HepLib::o_flintfD
const int o_flintfD
Definition: Init.cpp:112
HepLib::normal_flint
ex normal_flint(const ex &expr, int opt=o_flint)
HepLib::as
const Symbol as
HepLib::collect_lst
lst collect_lst(ex const &expr_in, std::function< bool(const ex &)> has_func, int opt)
the collect function like Mathematica, reture the lst { {c1,v1}, {c2,v2}, ... }
Definition: BASIC.cpp:1222
HepLib::Verbose
int Verbose
Definition: Init.cpp:139
HepLib::sort_lst
void sort_lst(lst &ilst, bool less=true)
sort the list in less order, or the reverse
Definition: Sort.cpp:79
HepLib::vec2lst
lst vec2lst(const exvector &ev)
convert exvector to lst
Definition: BASIC.cpp:902
HepLib::w1
ex w1
Definition: BASIC.h:499
HepLib::nopat
unsigned nopat
Definition: Init.cpp:88
HepLib::GiNaC_Parallel
exvector GiNaC_Parallel(int ntotal, std::function< ex(int)> f, const string &key)
GiNaC Parallel Evaluation using fork.
Definition: BASIC.cpp:260
HepLib::w2
ex w2
Definition: BASIC.h:499
HepLib::string_trim
void string_trim(string &str)
Definition: Functions.cpp:156
HepLib::ex2int
int ex2int(ex num)
ex to integer
Definition: BASIC.cpp:893
HepLib::subs
ex subs(const ex &e, init_list sl)
Definition: BASIC.h:51
Generated by doxygen 1.9.1