HepLib
Loading...
Searching...
No Matches
Solver.cpp
Go to the documentation of this file.
1
6#include "BASIC.h"
7#include "IBP.h"
8#include <cmath>
9
10namespace 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 fermat.Execute("&(u=3);");
254
255 if(pn>0) fermat.Execute("Redrowech([m],,,"+to_string(pn)+");");
256 else fermat.Execute("Redrowech([m]);");
257 ss << "&(U=1);" << endl; // ugly printing, the whitespace matters
258 smat.clear();
259 if(true) { // read matrix m
260 ss << "![m]" << endl;
261 auto ostr = fermat.Execute(ss.str());
262 ss.clear();
263 ss.str("");
264
265 // make sure last char is 0
266 if(ostr[ostr.length()-1]!='0') throw Error("RowReduce, last char is NOT 0.");
267 ostr = ostr.substr(0, ostr.length()-1);
268 string_trim(ostr);
269
270 ostr.erase(0, ostr.find(":=")+2);
271 size_t sn = ostr.length();
272 char lc;
273 for(size_t i=0; i<sn; i++) {
274 char & c = ostr[i];
275 if(c=='[') {
276 c = '{';
277 if(i>0 && lc=='}') ostr[i-1] = ',';
278 } else if(c==']') c = '}';
279 else if(c==' '||c=='\t'||c=='\n'||c=='\r') continue;
280 lc = c;
281 }
282
283 Parser fp(st);
284 auto res = fp.Read(ostr);
285 for(auto const & item : res) {
286 int r = -1;
287 for(auto const & it : item) {
288 if(r==-1) r = ex2int(it)-1;
289 else {
290 int c = ex2int(it.op(0))-1;
291 smat[r][c] = it.op(1);
292 }
293 }
294 }
295 }
296 // note the order, before exfactor (normal_fermat will be called again here)
297 ss << "&(U=0);" << endl; // disable ugly printing
298 ss << "@([m]);" << endl;
299 ss << "&_G;" << endl;
300 fermat.Execute(ss.str());
301 ss.clear();
302 ss.str("");
303 }
304
305 void Solver::SolveSparse(const ex & sector, const map<int,int> & n2n) {
306 int pdim = Propagator.nops();
307 lst ibps = IBPs;
308
309 // add more eqn to ibps
310 if(false) {
311 if(true) {
312 int n = ibps.nops();
313 for(int i=0; i<pdim; i++) {
314 for(int j1=0; j1<n; j1++) for(int j2=0; j2<n; j2++) {
315 ibps.append(ibps.op(i).subs(lst{a(j1)==a(j1)+1}).subs(lst{a(j2)==a(j2)-1}));
316 //ibps.append(ibps.op(i).subs(lst{a(j1)==a(j1)+1}).subs(lst{a(j2)==a(j2)+1}));
317 //ibps.append(ibps.op(i).subs(lst{a(j1)==a(j1)-1}).subs(lst{a(j2)==a(j2)-1}));
318 }
319 }
320 } else {
321 int n = ibps.nops();
322 for(int i=0; i<pdim; i++) {
323 for(int j=0; j<n; j++) ibps.append(ibps.op(j).subs(a(i)==a(i)+1));
324// for(int j=0; j<n; j++) ibps.append(ibps.op(j).subs(a(i)==a(i)+2));
325// for(int j=0; j<n; j++) ibps.append(ibps.op(j).subs(a(i)==a(i)+3));
326// for(int j=0; j<n; j++) ibps.append(ibps.op(j).subs(a(i)==a(i)+4));
327 }
328
329 n = ibps.nops();
330 for(int i=0; i<pdim; i++) {
331 for(int j=0; j<n; j++) ibps.append(ibps.op(j).subs(a(i)==a(i)-1));
332// for(int j=0; j<n; j++) ibps.append(ibps.op(j).subs(a(i)==a(i)-2));
333// for(int j=0; j<n; j++) ibps.append(ibps.op(j).subs(a(i)==a(i)-3));
334// for(int j=0; j<n; j++) ibps.append(ibps.op(j).subs(a(i)==a(i)-4));
335 }
336 }
337 } else {
338 {
339 exset fs;
340 find(ibps, F(w), fs);
341 int n = ibps.nops();
342 for(auto fi : fs) {
343 exmap a2a;
344 for(int i=0; i<fi.op(0).nops(); i++) a2a[a(i)] = a(i)-fi.op(0).op(i).subs(a(i)==0);
345 for(int j=0; j<n; j++) ibps.append(ibps.op(j).subs(a2a));
346 }
347 }
348// for(auto fi : fs) {
349// exmap a2a;
350// for(int i=0; i<fi.op(0).nops(); i++) a2a[a(i)] = fi.op(0).op(i)+1;
351// for(int j=0; j<n; j++) ibps.append(ibps.op(j).subs(a2a));
352// }
353// for(auto fi : fs) {
354// exmap a2a;
355// for(int i=0; i<fi.op(0).nops(); i++) a2a[a(i)] = fi.op(0).op(i)-1;
356// for(int j=0; j<n; j++) ibps.append(ibps.op(j).subs(a2a));
357// }
358 }
359
360 // some a's replaced by integer
361 vector<int> nfix(pdim);
362 for(int i=0; i<pdim; i++) nfix[i] = 0;
363 exmap a2n;
364 for(auto const & kv : n2n) {
365 nfix[kv.first] = 1;
366 int k = kv.first;
367 if(k<0) k = pdim + k;
368 a2n[a(k)] = kv.second;
369 }
370 if(a2n.size()>0) {
371 int nibps = ibps.nops();
372 for(int i=0; i<nibps; i++) ibps.let_op(i) = ibps.op(i).subs(a2n);
373 }
374
375 lst cons;
376 map<ex,lst,ex_is_less> cons_vec;
377 while(ibps.nops()>0) {
378 ibps.sort();
379 ibps.unique();
380
381 exset fset;
382 find(ibps,F(w),fset);
383 map<int,exvector> fmap_vec;
384 for(auto fi : fset) { // sort fset
385 int psum = 0, nsum = 0;
386 auto ns = fi.op(0);
387 for(int i=0; i<pdim; i++) {
388 auto ni = ns.op(i).subs(a(i)==0,nopat);
389 if(sector.op(i)==1) {
390 if(nfix[i]==1 && ni<=0) psum -= 100000;
391 else psum += ex2int(ni)-1; // -1 from FIRE
392 } else {
393 if(nfix[i]==1 && ni>0) nsum += ex2int(ni);
394 else nsum -= ex2int(ni);
395 }
396 }
397 int key = psum+nsum;
398 fmap_vec[-key].push_back(fi); // - to reverse the order
399 }
400
401 map<ex,int,ex_is_less> f2n;
402 int next_gi = 0;
403 exvector fvec(fset.size());
404 int nr = 1, nc = 0;
405 for(auto const & ev : fmap_vec) {
406 if(next_gi<1) next_gi = ev.second.size();
407 for(auto const & v : ev.second) {
408 fvec[nc] = v;
409 f2n[v] = nc;
410 nc++;
411 }
412 nr++;
413 }
414 nr = ibps.nops();
415
416 if(Verbose>0) {
417 cout << "+----------------------------------------" << endl;
418 cout << "| IBPs: " << nr << ", Fs: " << next_gi << "/" << nc << endl;
419 cout << "+----------------------------------------" << endl;
420 }
421
422 SparseMatrix smat;
423 for(int r=0; r<nr; r++) {
424 auto ibp = ibps.op(r);
425 auto cvs = collect_lst(ibp,F(w));
426 for(auto & cv : cvs) {
427 int c = f2n[cv.op(1)];
428 ex cc = cv.op(0);
429 if(!cc.is_zero()) smat[r][c] = cc;
430 }
431 }
432 RowReduce(smat,next_gi);
433
434 ibps.remove_all();
435 for(auto const & rv : smat) {
436 int next_sol_ibp = -1; // -1:init, 0:next, 1:sol, 2:ibp
437 ex sol = 0, ns, fc;
438 exvector nv(pdim);
439 exmap amap;
440 for(auto const & cv : rv.second) {
441 if(next_sol_ibp==-1) {
442 if(cv.first<next_gi) {
443 next_sol_ibp = 0;
444 if(!is_zero(cv.second-1)) throw Error("Solver::SolveSparse, NOT 1.");
445 ns = fvec[cv.first].op(0).subs(a(w)==0);
446 for(int i=0; i<pdim; i++) {
447 if(nfix[i]==0) amap[a(i)] = a(i)-ns.op(i);
448 }
449 fc = fvec[cv.first].subs(amap,nopat);
450 for(int i=0; i<pdim; i++) {
451 if(nfix[i]==1) {
452 nv[i] = ns.op(i);
453 if(sector.op(i)==1 && nv[i]<1) goto next_row; // subsector
454 else if(sector.op(i)!=1 && nv[i]>0) goto next_row; // out of this sector
455 } else nv[i] = (sector.op(i)==1 ? 1 : 0);
456 }
457 continue; // need continue
458 } else next_sol_ibp = 2; // ibp
459 }
460
461 if(next_sol_ibp==0) {
462 if(cv.first<next_gi) goto next_row; // not a solution or ibp
463 next_sol_ibp = 1; // solution found
464 }
465
466 if(true) {
467 ex f = fvec[cv.first];
468 ex nd = cv.second;
469 if(next_sol_ibp==1) {
470 f = f.subs(amap,nopat);
471 nd = nd.subs(amap,nopat).numer_denom();
472 ex num = exfactor(nd.op(0));
473 ex den = exfactor(nd.op(1));
474 ns = f.op(0).subs(a(w)==0); // reuse ns here
475 for(int i=0; i<pdim; i++) {
476 if(sector.op(i)!=1) {
477 ex nsi = ns.op(i);
478 if(nfix[i]==1) {
479 if(nsi>0) goto next_row; // out of this sector
480 else continue;
481 }
482 nsi = -nsi; // note '-' here, ni<=nsi
483 // due to possible (ai+nsi-1)*F(...,ai-nsi,...), so nsi+1 may be OK
484 while(true) {
485 ex cci = num.subs(a(i)==1+nsi);
486 if(!cci.is_zero() || nsi>=nv[i]) break;
487 nsi++;
488 }
489 if(nsi<nv[i]) nv[i] = nsi;
490 }
491 }
492
493 if(!is_a<mul>(den)) den = lst{ den };
494 for(const auto & item : den) {
495 ex it = item;
496 if(is_a<power>(it)) it = it.op(0);
497 if(!it.has(a(w)) || it.has(d)) continue;
498 exset as;
499 find(it,a(w),as);
500 if(as.size()!=1) {
501 cout << "size is NOT 1: " << it << endl;
502 throw Error("Solver::SolveSparse, error occured.");
503 }
504 ex aw = *(as.begin());
505 if(it.degree(aw)!=1) {
506 cout << "degree is NOT 1: " << it << endl;
507 throw Error("Solver::SolveSparse, error occured.");
508 }
509 ex a0 = -it.coeff(aw,0)/it.coeff(aw,1);
510 int an = ex2int(aw.op(0));
511 if(sector.op(an)!=1) {
512 if(a0<=nv[an]) nv[an] = a0-1;
513 } else {
514 if(a0>=nv[an]) nv[an] = a0+1;
515 }
516 }
517 }
518 sol -= nd * f;
519 }
520 }
521
522 if(next_sol_ibp==1) {
523 for(int i=0; i<pdim; i++) { // check
524 if(nfix[i]==1 && fc.op(0).op(i)!=nv[i]) {
525 cout << i << ": " << fc << endl << nv << endl;
526 throw Error("Solver::SolveSparse, error occured.");
527 } else if(nfix[i]==0 && fc.op(0).op(i)!=a(i)) {
528 cout << i << ": " << fc << endl << nv << endl;
529 throw Error("Solver::SolveSparse, error occured.");
530 }
531 }
532
533 if(true) {
534 int psum = 0, nsum = 0;
535 for(int i=0; i<pdim; i++) {
536 if(nfix[i]==1) continue;
537 if(sector.op(i)==1 && !nv[i].is_equal(1)) psum++;
538 if(sector.op(i)!=1 && !nv[i].is_zero()) nsum++;
539 }
540 if(psum<=1 && nsum<=1) {
541 ex con = vec2lst(nv);
542 cons.append(con);
543 cons_vec[con].append(sol);
544 }
545 }
546 } else ibps.append(sol);
547
548 next_row: ;
549 }
550cons.sort();
551cons.unique();
552sort_lst(cons,false);
553cout << cons << endl << endl;
554if(n2n.size()>0) {
555 for(auto & item : cons) {
556 for(auto kv : n2n) {
557 if(kv.second==0) continue;
558 if(item.op(kv.first)!=kv.second) goto done;
559 }
560 cout << endl << item << endl << cons_vec[item] << endl << endl;
561 done: ;
562 }
563}
564 }
565
566 cons.sort();
567 cons.unique();
568 sort_lst(cons,false);
569
570 cout << endl;
571 for(auto & k : cons) {
572 cons_vec[k].sort();
573 cons_vec[k].unique();
574 cout << k << endl;
575 cout << cons_vec[k] << endl << endl;
576 }
577 cout << endl << cons << endl;
578
579
580 }
581
582 void Solver::IBP() {
583 int pdim = Propagator.nops();
584 lst InExternal;
585 for(auto ii : Internal) InExternal.append(ii);
586 for(auto ii : External) InExternal.append(ii);
587
588 if(ISP.nops()<1) {
589 for(auto it : Internal) {
590 for(auto ii : InExternal) ISP.append(it*ii);
591 }
592 ISP.sort();
593 ISP.unique();
594 }
595
596 if(ISP.nops() > pdim) {
597 cout << "ISP = " << ISP << endl;
598 cout << "Propagator = " << Propagator << endl;
599 throw Error("Solver::IBP: #(ISP) > #(Propagator).");
600 }
601
602 lst sp2s, s2sp, ss;
603 int _pic=0;
604 for(auto item : ISP) {
605 _pic++;
606 Symbol si("P"+to_string(_pic));
607 ss.append(si);
608 sp2s.append(w*item==w*si);
609 sp2s.append(item==si);
610 s2sp.append(si==item);
611 }
612
613 lst leqns;
614 for(int i=0; i<ISP.nops(); i++) { // note NOT pdim
615 auto eq = Propagator.op(i).expand().subs(iEpsilon==0); // drop iEpsilon
616 eq = eq.subs(sp2s);
617 eq = eq.subs(Replacement);
618 if(eq.has(iWF(w))) throw Error("Solver::IBP, iWF used in eq.");
619 leqns.append(eq == iWF(i));
620 }
621 auto s2p = lsolve(leqns, ss);
622 if(s2p.nops() != ISP.nops()) throw Error("Solver::IBP, lsolve failed.");
623
624 // 1st version
625 if(true)
626 if(DSP.nops()<1) {
627 for(auto p1 : Internal)
628 for(auto p2 : InExternal)
629 DSP.append(lst{p1,p2});
630 }
631
632 // Lee version
633 if(false)
634 if(DSP.nops()<1) {
635 for(int i=0; i<Internal.nops(); i++)
636 DSP.append(lst{Internal.op(i), Internal.op(i==Internal.nops()-1 ? 0 : i+1)});
637 for(auto p2 : InExternal) DSP.append(lst{Internal.op(0), p2});
638 DSP.append(lst{Internal, Internal});
639 }
640
641 IBPs.remove_all();
642 lst nsa;
643 for(int i=0; i<pdim; i++) nsa.append(a(i));
644 for(auto sp : DSP) {
645 auto ilp = sp.op(0);
646 auto iep = sp.op(1);
647
648 ex ibp = 0;
649 symbol ss;
650 for(int i=0; i<pdim; i++) {
651 auto ns = nsa;
652 ns.let_op(i) = nsa.op(i) + 1;
653 auto dp = Propagator.op(i).subs(ilp==ss).diff(ss).subs(ss==ilp);
654 ibp -= (a(i)+Shift[i]) * F(ns) * dp;
655 }
656
657 ibp = ibp * iep;
658 ibp = ibp.expand();
659 ibp = ibp.subs(sp2s);
660 ibp = ibp.subs(Replacement);
661 ibp = ibp.subs(s2p);
662
663 ex res = 0;
664 for(int i=0; i<pdim; i++) {
665 auto ci = ibp.coeff(iWF(i), 1);
666 ci = MapFunction([i](const ex &e, MapFunction &self)->ex {
667 if(!e.has(F(w))) return e;
668 else if(e.match(F(w))) {
669 lst tmp = ex_to<lst>(e.op(0));
670 tmp.let_op(i) = tmp.op(i)-1;
671 return F(tmp);
672 } else return e.map(self);
673 })(ci);
674 res += ci;
675 }
676 res += ibp.subs(lst{iWF(w)==0});
677 auto cilp = iep.coeff(ilp);
678 if(!is_zero(cilp)) res += d*cilp*F(nsa);
679 IBPs.append(res);
680 }
681 }
682
683 void Solver::Solve(const ex & sector, const map<int,int> & n2n) {
684 int pdim = Propagator.nops();
685 lst ibps = IBPs;
686
687 // add more eqn to ibps
688 if(false) {
689// int n = ibps.nops();
690// for(int i=0; i<pdim; i++) {
691// for(int j1=0; j1<n; j1++) for(int j2=0; j2<n; j2++)
692// ibps.append(ibps.op(i).subs(lst{a(j1)==a(j1)+1}).subs(lst{a(j2)==a(j2)-1}));
693// }
694
695 int 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
703 //n = ibps.nops();
704 for(int i=0; i<pdim; i++) {
705 for(int j=0; j<n; j++) ibps.append(ibps.op(j).subs(a(i)==a(i)-1));
706// for(int j=0; j<n; j++) ibps.append(ibps.op(j).subs(a(i)==a(i)-2));
707// for(int j=0; j<n; j++) ibps.append(ibps.op(j).subs(a(i)==a(i)-3));
708// for(int j=0; j<n; j++) ibps.append(ibps.op(j).subs(a(i)==a(i)-4));
709 }
710 } else {
711 exset fs;
712 find(ibps, F(w1,w2), fs);
713cout << fs << endl;
714exit(0);
715 }
716
717 // some a's replaced by integer
718 vector<int> nfix(pdim);
719 for(int i=0; i<pdim; i++) nfix[i] = 0;
720 for(auto const & kv : n2n) nfix[kv.first] = 1;
721 exmap a2n;
722 for(auto const & kv : n2n) {
723 int k = kv.first;
724 if(k<0) k = pdim + k;
725 a2n[a(k)] = kv.second;
726 }
727 int nibps = ibps.nops();
728 for(int i=0; i<nibps; i++) ibps.let_op(i) = ibps.op(i).subs(a2n);
729
730 lst cons;
731 map<ex,lst,ex_is_less> cons_vec;
732 while(ibps.nops()>0) {
733 ibps.sort();
734 ibps.unique();
735
736 exvector fvec;
737 exset fset;
738 find(ibps,F(w),fset);
739 int nmax = -10000000, pmax = -10000000;
740 for(auto fi : fset) { // sort fset
741 int psum = 0, nsum = 0;
742 auto ns = fi.op(0);
743 for(int i=0; i<pdim; i++) {
744 auto ni = ns.op(i).subs(a(i)==0,nopat);
745 if(sector.op(i)==1) {
746 if(nfix[i]==1 && ni<=0) psum -= 100000;
747 else psum += ex2int(ni)-1; // -1 from FIRE
748 } else {
749 if(nfix[i]==1 && ni>0) nsum += ex2int(ni);
750 else nsum -= ex2int(ni);
751 }
752 }
753
754 if(psum+nsum<pmax+nmax) continue;
755 if(psum+nsum>pmax+nmax) {
756 pmax = psum;
757 nmax = nsum;
758 fvec.clear();
759 }
760
761 fvec.push_back(fi);
762 }
763 //sort_vec(fvec);
764
765 int nr = ibps.nops();
766 int nc = fvec.size();
767 matrix mat(nr, nc);
768
769 if(Verbose>0) {
770 cout << "+----------------------------------------" << endl;
771 cout << "| IBPs: " << nr << ", Fs: " << nc << endl;
772 cout << "+----------------------------------------" << endl;
773 }
774
775 matrix bvec(nr,1);
776 for(int r=0; r<nr; r++) {
777 auto ibp = ibps.op(r);
778 for(int c=0; c<nc; c++) {
779 mat(r,c) = ibp;
780 ibp = ibp.subs(fvec[c]==0,nopat);
781 mat(r,c) = (mat(r,c)-ibp).coeff(fvec[c]);
782 }
783 bvec(r,0) = ibp;
784 }
785
786 // note that: mr = mt.mul(mat)
787 auto rt = RowReduce(mat);
788 auto mr = rt.first;
789 auto mt = rt.second;
790 bvec = mt.mul(bvec);
791 for(int r=0; r<nr; r++) bvec(r,0) = collect_ex(bvec(r,0),F(w),o_flintfD);
792
793 //cout << "CHECK: 1=" << ex_to<matrix>(normal(mr.sub(mt.mul(mat)))).is_zero_matrix() << endl;
794
795 auto ibp_con_sol_vec = GiNaC_Parallel(nr, [&](const int r)->ex {
796 ex ibp=0, con=0;
797 int cc = -1;
798 for(int c=0; c<nc; c++) {
799 if(mr(r,c).is_zero()) continue;
800 if(mr(r,c)==1 && cc==-1) cc = c;
801 else goto next_row; // not a solution, skip
802 }
803
804 if(cc!=-1) { // a solution found N -> N-1, N-2, subsectors
805 auto ns = fvec[cc].op(0).subs(a(w)==0);
806 exmap amap;
807 for(int i=0; i<pdim; i++) {
808 if(nfix[i]==0) amap[a(i)] = a(i)-ns.op(i);
809 }
810 fvec[cc] = fvec[cc].subs(amap,nopat);
811 auto sol = bvec(r,0).subs(amap,nopat);
812 sol = normal_flint(sol,o_flintfD);
813
814 exvector nv(pdim);
815 for(int i=0; i<pdim; i++) {
816 if(nfix[i]==1) {
817 nv[i] = ns.op(i);
818 if(sector.op(i)==1 && nv[i]<1) goto next_row; // subsector
819 else if(sector.op(i)!=1 && nv[i]>0) goto next_row; // out of this sector
820 } else nv[i] = (sector.op(i)==1 ? 1 : 0);
821 }
822
823 ex num = sol.numer();
824 exset fset;
825 find(sol,F(w),fset);
826 for(auto const & f : fset) { // make sure not to increase the sector
827 ex cc = factor_flint(num.coeff(f));
828 ns = f.op(0).subs(a(w)==0); // reuse ns here
829 for(int i=0; i<pdim; i++) {
830 if(sector.op(i)!=1) {
831 ex nsi = ns.op(i);
832 if(nfix[i]==1) {
833 if(nsi>0) goto next_row; // out of this sector
834 else continue;
835 }
836 nsi = -nsi; // note '-' here, ni<=nsi
837 // due to possible (ai+nsi-1)*F(...,ai-nsi,...), so nsi+1 may be OK
838 while(true) {
839 ex cci = cc.subs(a(i)==1+nsi);
840 if(!cci.is_zero() || nsi>=nv[i]) break;
841 nsi++;
842 }
843 if(nsi<nv[i]) nv[i] = nsi;
844 }
845 }
846 }
847
848 ex den = sol.denom();
849 if(!is_a<mul>(den)) den = lst{den};
850 for(const auto & item : den) {
851 ex it = item;
852 if(is_a<power>(it)) it = it.op(0);
853 if(!it.has(a(w)) || it.has(d)) continue;
854 exset as;
855 find(it,a(w),as);
856 if(as.size()!=1) {
857 cout << "size is NOT 1: " << it << endl;
858 throw Error("Solver::Solve, error occured.");
859 }
860 ex aw = *(as.begin());
861 if(it.degree(aw)!=1) {
862 cout << "degree is NOT 1: " << it << endl;
863 throw Error("Solver::Solve, error occured.");
864 }
865 ex a0 = -it.coeff(aw,0)/it.coeff(aw,1);
866 int an = ex2int(aw.op(0));
867 if(sector.op(an)!=1) {
868 if(a0<=nv[an]) nv[an] = a0-1;
869 } else {
870 if(a0>=nv[an]) nv[an] = a0+1;
871 }
872 }
873
874 for(int i=0; i<pdim; i++) { // check
875 if(nfix[i]==1 && fvec[cc].op(0).op(i)!=nv[i]) {
876 cout << fvec << endl << nv << endl;
877 throw Error("Solver::Solve, error occured.");
878 } else if(nfix[i]==0 && fvec[cc].op(0).op(i)!=a(i)) {
879 cout << fvec << endl << nv << endl;
880 throw Error("Solver::Solve, error occured.");
881 }
882 }
883
884 //cout << "corner: " << nv << endl;
885 int psum = 0, nsum = 0;
886 for(int i=0; i<pdim; i++) {
887 if(nfix[i]==1) continue;
888 if(sector.op(i)==1 && !nv[i].is_equal(1)) psum++;
889 if(sector.op(i)!=1 && !nv[i].is_zero()) nsum++;
890 }
891 if(nsum<=1 && psum<=1) con = lst{ vec2lst(nv), sol };
892 } else if(!bvec(r,0).is_zero()) ibp = bvec(r,0).numer();
893 next_row: ;
894 return lst{ibp, con};
895 });
896
897 ibps.remove_all();
898 for(int r=0; r<nr; r++) {
899 auto const & item = ibp_con_sol_vec[r];
900 if(!item.op(0).is_zero()) ibps.append(item.op(0));
901 if(is_a<lst>(item.op(1))) {
902 cons.append(item.op(1).op(0));
903 cons_vec[item.op(1).op(0)].append(item.op(1).op(1));
904 }
905 }
906
907cons.sort();
908cons.unique();
909cout << cons << endl << endl;
910 }
911
912 cons.sort();
913 cons.unique();
914 sort_lst(cons,false);
915
916 cout << endl;
917 for(auto & k : cons) {
918 cons_vec[k].sort();
919 cons_vec[k].unique();
920 cout << k << endl;
921 cout << cons_vec[k] << endl << endl;
922 }
923 cout << endl << cons << endl;
924 }
925
929 void Solver::Export() {
930
931 }
932
933
937 void Solver::Run() {
938
939 }
940
944 void Solver::Import() {
945
946 }
947
948}
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::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::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::SolveSparse
void SolveSparse(const ex &sector, const map< int, int > &a2n={})
Definition Solver.cpp:305
HepLib::Symbol
class extended to GiNaC symbol class, represent a positive symbol
Definition BASIC.h:113
HepLib
HepLib namespace.
Definition BASIC.cpp:17
HepLib::exfactor
ex exfactor(const ex &expr_in, int opt)
factorize a expression
Definition BASIC.cpp:1854
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:93
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:115
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:142
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:91
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.8