HepLib
Loading...
Searching...
No Matches
FIRE.cpp
Go to the documentation of this file.
1
6#include "IBP.h"
7#include <cmath>
8
9#include "exFlint.h"
10
11namespace HepLib {
12
13 // FROM FIRE6.m
14 namespace {
15 unsigned long long HardCodedPrimes[128] = {
16 2017llu, 18446744073709551557llu, 18446744073709551533llu,
17 18446744073709551521llu, 18446744073709551437llu, 18446744073709551427llu,
18 18446744073709551359llu, 18446744073709551337llu, 18446744073709551293llu,
19 18446744073709551263llu, 18446744073709551253llu, 18446744073709551191llu,
20 18446744073709551163llu, 18446744073709551113llu, 18446744073709550873llu,
21 18446744073709550791llu, 18446744073709550773llu, 18446744073709550771llu,
22 18446744073709550719llu, 18446744073709550717llu, 18446744073709550681llu,
23 18446744073709550671llu, 18446744073709550593llu, 18446744073709550591llu,
24 18446744073709550539llu, 18446744073709550537llu, 18446744073709550381llu,
25 18446744073709550341llu, 18446744073709550293llu, 18446744073709550237llu,
26 18446744073709550147llu, 18446744073709550141llu, 18446744073709550129llu,
27 18446744073709550111llu, 18446744073709550099llu, 18446744073709550047llu,
28 18446744073709550033llu, 18446744073709550009llu, 18446744073709549951llu,
29 18446744073709549861llu, 18446744073709549817llu, 18446744073709549811llu,
30 18446744073709549777llu, 18446744073709549757llu, 18446744073709549733llu,
31 18446744073709549667llu, 18446744073709549621llu, 18446744073709549613llu,
32 18446744073709549583llu, 18446744073709549571llu, 18446744073709549519llu,
33 18446744073709549483llu, 18446744073709549441llu, 18446744073709549363llu,
34 18446744073709549331llu, 18446744073709549327llu, 18446744073709549307llu,
35 18446744073709549237llu, 18446744073709549153llu, 18446744073709549123llu,
36 18446744073709549067llu, 18446744073709549061llu, 18446744073709549019llu,
37 18446744073709548983llu, 18446744073709548899llu, 18446744073709548887llu,
38 18446744073709548859llu, 18446744073709548847llu, 18446744073709548809llu,
39 18446744073709548703llu, 18446744073709548599llu, 18446744073709548587llu,
40 18446744073709548557llu, 18446744073709548511llu, 18446744073709548503llu,
41 18446744073709548497llu, 18446744073709548481llu, 18446744073709548397llu,
42 18446744073709548391llu, 18446744073709548379llu, 18446744073709548353llu,
43 18446744073709548349llu, 18446744073709548287llu, 18446744073709548271llu,
44 18446744073709548239llu, 18446744073709548193llu, 18446744073709548119llu,
45 18446744073709548073llu, 18446744073709548053llu, 18446744073709547821llu,
46 18446744073709547797llu, 18446744073709547777llu, 18446744073709547731llu,
47 18446744073709547707llu, 18446744073709547669llu, 18446744073709547657llu,
48 18446744073709547537llu, 18446744073709547521llu, 18446744073709547489llu,
49 18446744073709547473llu, 18446744073709547471llu, 18446744073709547371llu,
50 18446744073709547357llu, 18446744073709547317llu, 18446744073709547303llu,
51 18446744073709547117llu, 18446744073709547087llu, 18446744073709547003llu,
52 18446744073709546897llu, 18446744073709546879llu, 18446744073709546873llu,
53 18446744073709546841llu, 18446744073709546739llu, 18446744073709546729llu,
54 18446744073709546657llu, 18446744073709546643llu, 18446744073709546601llu,
55 18446744073709546561llu, 18446744073709546541llu, 18446744073709546493llu,
56 18446744073709546429llu, 18446744073709546409llu, 18446744073709546391llu,
57 18446744073709546363llu, 18446744073709546337llu, 18446744073709546333llu,
58 18446744073709546289llu, 18446744073709546271llu};
59 }
60
61 inline lst syms(const exvector & ev) {
62 exset ss;
63 for(auto e : ev) {
64 for(const_preorder_iterator i=e.preorder_begin(); i!=e.preorder_end(); ++i)
65 if(is_a<symbol>(*i)) ss.insert(*i);
66 }
67 lst ls;
68 for(auto item : ss) ls.append(item);
69 return ls;
70 }
71
72 void FIRE::RRTables(const string & filename, int pnum) {
73 if(pnum<1) return;
74 exvector _tables(pnum);
75 unsigned int nmin = 0;
76 for(int i=0; i<pnum; i++) {
77 string fn = filename;
78 string_replace_all(fn, ".tables", "-"+to_string(i+1)+".tables");
79 ifstream ifs(fn);
80 string ostr((istreambuf_iterator<char>(ifs)), (istreambuf_iterator<char>()));
81 ifs.close();
82 string_replace_all(ostr, "\"", "");
83 Parser tp;
84 _tables[i] = tp.Read(ostr);
85 auto nn = _tables[i].op(1).nops();
86 if(i==0 || nmin>nn) nmin = nn;
87 }
88 exvector tables;
89 vector<numeric> primes;
90 ex convs = 0;
91 for(int i=0; i<_tables.size(); i++) {
92 auto nn = _tables[i].op(1).nops();
93 if(nmin<nn) {
94 cout << "RRTables: current table will be dropped." << endl;
95 continue;
96 }
97 if(!is_a<lst>(convs)) convs = _tables[i].op(1); // second part is equal
98 else if(!convs.is_equal(_tables[i].op(1))) {
99 throw Error("RRT: convs are different.");
100 }
101 tables.push_back(_tables[i].op(0));
102 primes.push_back(HardCodedPrimes[i+1]);
103 }
104
105 vector<pair<ex,vector<pair<ex,exvector>>>> int_mi_cs_vec;
106 int ci = 0;
107 for(int i=0; i<tables.size(); i++) {
108 if(i==0) {
109 for(auto ti : tables[0]) {
110 pair<ex,vector<pair<ex,exvector>>> int_mi_cs;
111 int_mi_cs.first = ti.op(0);
112 for(auto tii : ti.op(1)) {
113 pair<ex,exvector> mi_cs;
114 mi_cs.first = tii.op(0);
115 mi_cs.second.push_back(tii.op(1));
116 int_mi_cs.second.push_back(mi_cs);
117 }
118 int_mi_cs_vec.push_back(int_mi_cs);
119 }
120 } else {
121 for(int i2=0; i2<int_mi_cs_vec.size(); i2++) {
122 auto ti = tables[i].op(i2);
123 auto & int_mi_cs = int_mi_cs_vec[i2];
124 if(int_mi_cs.first != ti.op(0)) throw Error("E1");
125 for(int j=0; j<int_mi_cs.second.size(); j++) {
126 auto & mi_cs = int_mi_cs.second[j];
127 if(mi_cs.first != ti.op(1).op(j).op(0)) throw Error("E2");
128 mi_cs.second.push_back(ti.op(1).op(j).op(1));
129 }
130 }
131 }
132 }
133
134 std::ofstream ofs;
135 ofs.open(filename, ios::out);
136 ofs << "{" << endl;
137 ofs << " {" << endl;
138 for(int i=0; i<int_mi_cs_vec.size(); i++) {
139 auto imc = int_mi_cs_vec[i];
140 ofs << " {" << imc.first << "," << endl;
141 ofs << " {" << endl;
142 for(int j=0; j<imc.second.size(); j++) {
143 auto item = imc.second[j];
144 ofs << " {" << item.first << ", ";
145 vector<numeric> nvec;
146 for(auto e : item.second) nvec.push_back(ex_to<numeric>(e));
147 ofs << "\"" << RationalReconstruct(nvec,primes) << "\"}";
148 if(j+1<imc.second.size()) ofs << ",";
149 ofs << endl;
150 }
151 ofs << " }" << endl;
152 ofs << " }";
153 if(i+1<int_mi_cs_vec.size()) ofs << ",";
154 ofs << endl;
155 }
156 ofs << " }," << endl << " {" << endl;
157 for(int i=0; i<convs.nops(); i++) {
158 auto item = convs.op(i);
159 ofs << " {" << item.op(0) << ", " << item.op(1) << "}";
160 if(i+1<convs.nops()) ofs << ",";
161 ofs << endl;
162 }
163 ofs << " }" << endl;
164 ofs << "}" << endl;
165 ofs.close();
166 }
167
168 void FIRE::ThieleTables(const string & filename, int si, int ei) {
169 int pnum = 40;
170 map<int,ex> _tables;
171 int nd;
172 unsigned int nmin = 0;
173 for(int i=si; i<=ei; i++) {
174 cout << "\r \r" << flush;
175 cout << "Reading tables: " << ei << "|" << i << flush;
176 string fn = filename;
177 string_replace_all(fn, ".tables", "-"+to_string(i)+".tables");
178 ifstream ifs(fn);
179 string ostr((istreambuf_iterator<char>(ifs)), (istreambuf_iterator<char>()));
180 ifs.close();
181 for(auto & c : ostr) if(c=='\"') c = ' ';
182 Parser tp;
183 _tables[i] = tp.Read(ostr);
184 auto nn = _tables[i].op(1).nops();
185 if(i==si || nmin>nn) nmin = nn;
186 }
187 cout << endl;
188
189 exvector tables;
190 exvector keys;
191 ex convs = 0;
192 for(int i=si; i<=ei; i++) {
193 auto nn = _tables[i].op(1).nops();
194 if(nmin<nn) {
195 cout << "ThieleTables: current table will be dropped." << endl;
196 continue;
197 }
198 if(!is_a<lst>(convs)) convs = _tables[i].op(1); // second part is equal
199 else if(!convs.is_equal(_tables[i].op(1))) {
200 throw Error("RRT: convs are different.");
201 }
202 tables.push_back(_tables[i].op(0));
203 keys.push_back(i);
204 }
205
206 vector<pair<ex,vector<pair<ex,exvector>>>> int_mi_cs_vec;
207 int ci = 0;
208 for(int i=0; i<tables.size(); i++) {
209 if(i==0) {
210 for(auto ti : tables[0]) {
211 pair<ex,vector<pair<ex,exvector>>> int_mi_cs;
212 int_mi_cs.first = ti.op(0);
213 for(auto tii : ti.op(1)) {
214 pair<ex,exvector> mi_cs;
215 mi_cs.first = tii.op(0);
216 mi_cs.second.push_back(tii.op(1));
217 int_mi_cs.second.push_back(mi_cs);
218 }
219 int_mi_cs_vec.push_back(int_mi_cs);
220 }
221 } else {
222 for(int i2=0; i2<int_mi_cs_vec.size(); i2++) {
223 auto ti = tables[i].op(i2);
224 auto & int_mi_cs = int_mi_cs_vec[i2];
225 if(int_mi_cs.first != ti.op(0)) throw Error("E1");
226 for(int j=0; j<int_mi_cs.second.size(); j++) {
227 auto & mi_cs = int_mi_cs.second[j];
228 if(mi_cs.first != ti.op(1).op(j).op(0)) throw Error("E2");
229 mi_cs.second.push_back(ti.op(1).op(j).op(1));
230 }
231 }
232 }
233 }
234
235 std::ofstream ofs;
236 ofs.open(filename, ios::out);
237 ofs << "{" << endl;
238 ofs << " {" << endl;
239 for(int i=0; i<int_mi_cs_vec.size(); i++) {
240 cout << "\r \r" << flush;
241 cout << "Thiele: " << int_mi_cs_vec.size() << "|" << i+1 << flush;
242 auto imc = int_mi_cs_vec[i];
243 ofs << " {" << imc.first << "," << endl;
244 ofs << " {" << endl;
245
246 auto nnn = imc.second.size();
247 vector<fmpz_mpoly_q_struct*> res_vec(nnn);
248 vector<lst> xs_vec(nnn);
249 #pragma omp parallel for schedule(dynamic, 1)
250 for(int j=0; j<nnn; j++) {
251 auto item = imc.second[j];
252 xs_vec[j] = syms(item.second);
253 fmpz_mpoly_ctx_t ctx;
254 fmpz_mpoly_ctx_init(ctx, xs_vec[j].nops(), ORD_LEX);
255 res_vec[j] = (fmpz_mpoly_q_struct*) flint_malloc(sizeof(fmpz_mpoly_q_struct));
256 fmpz_mpoly_q_init(res_vec[j], ctx);
257 int n = keys.size();
258 vector<fmpz_mpoly_q_struct*> key_vec(n);
259 vector<fmpz_mpoly_q_struct*> val_vec(n);
260 vector<fmpz_mpoly_q_struct*> coeff_vec(n);
261 for(int k=0; k<n; k++) {
262 key_vec[j] = (fmpz_mpoly_q_struct*) flint_malloc(sizeof(fmpz_mpoly_q_struct));
263 fmpz_mpoly_q_init(key_vec[k], ctx);
264 _to_(xs_vec[j], key_vec[k], ctx, item.first.op(k));
265 val_vec[j] = (fmpz_mpoly_q_struct*) flint_malloc(sizeof(fmpz_mpoly_q_struct));
266 fmpz_mpoly_q_init(val_vec[k], ctx);
267 _to_(xs_vec[j], val_vec[j], ctx, item.second[k]);
268 coeff_vec[j] = (fmpz_mpoly_q_struct*) flint_malloc(sizeof(fmpz_mpoly_q_struct));
269 fmpz_mpoly_q_init(coeff_vec[k], ctx);
270 fmpz_mpoly_q_zero(coeff_vec[k], ctx);
271 }
272
273 fmpz_mpoly_q_t t1, t2;
274 fmpz_mpoly_q_init(t1, ctx);
275 fmpz_mpoly_q_init(t2, ctx);
276 for(int i=0; i<n; i++) {
277 fmpz_mpoly_q_set(t1, val_vec[i], ctx);
278 for(int j=0; j<i; j++) {
279 if(fmpz_mpoly_q_equal(t1, coeff_vec[j], ctx)) {
280 n = i-1;
281 goto out;
282 }
283 fmpz_mpoly_q_sub(t2, key_vec[i], key_vec[j], ctx);
284 fmpz_mpoly_q_sub(t1, t1, coeff_vec[j], ctx);
285 fmpz_mpoly_q_div(t1, t2, t1, ctx);
286 }
287 fmpz_mpoly_q_set(coeff_vec[i], t1, ctx);
288 }
289 throw Error("Thiele unstable!");
290 out: ;
291 fmpz_mpoly_q_set(t1, coeff_vec[n], ctx);
292 fmpz_mpoly_q_t dx;
293 fmpz_mpoly_q_init(dx, ctx);
294 _to_(xs_vec[j], dx, ctx, d);
295 for(int i=n-1; i>=0; i--) {
296 fmpz_mpoly_q_sub(t2, dx, key_vec[i], ctx);
297 fmpz_mpoly_q_div(t2, t2, t1, ctx);
298 fmpz_mpoly_q_add(t1, coeff_vec[i], t2, ctx);
299 }
300 fmpz_mpoly_q_clear(dx, ctx);
301 fmpz_mpoly_q_set(res_vec[i], t1, ctx);
302 for(int k=0; k<n; k++) {
303 fmpz_mpoly_q_clear(key_vec[k], ctx);
304 fmpz_mpoly_q_clear(val_vec[k], ctx);
305 fmpz_mpoly_q_clear(coeff_vec[k], ctx);
306 }
307 fmpz_mpoly_q_clear(t1, ctx);
308 fmpz_mpoly_q_clear(t2, ctx);
309 fmpz_mpoly_ctx_clear(ctx);
310 }
311 for(int j=0; j<imc.second.size(); j++) {
312 fmpz_mpoly_ctx_t ctx;
313 fmpz_mpoly_ctx_init(ctx, xs_vec[j].nops(), ORD_LEX);
314 auto item = imc.second[j];
315 ex res;
316 res = _to_(xs_vec[j], res_vec[j], ctx);
317 fmpz_mpoly_q_clear(res_vec[j], ctx);
318 ofs << " {" << item.first << ", ";
319 ofs << "\"" << res << "\"}";
320 if(j+1<imc.second.size()) ofs << ",";
321 ofs << endl;
322 }
323 ofs << " }" << endl;
324 ofs << " }";
325 if(i+1<int_mi_cs_vec.size()) ofs << ",";
326 ofs << endl;
327 }
328 cout << endl;
329 ofs << " }," << endl << " {" << endl;
330 for(int i=0; i<convs.nops(); i++) {
331 auto item = convs.op(i);
332 ofs << " {" << item.op(0) << ", " << item.op(1) << "}";
333 if(i+1<convs.nops()) ofs << ",";
334 ofs << endl;
335 }
336 ofs << " }" << endl;
337 ofs << "}" << endl;
338 ofs.close();
339 }
340
344 void FIRE::Export() {
345 if(WorkingDir=="") WorkingDir = to_string(getpid())+"IBP";
346 int pdim = Propagator.nops();
347 string spn = to_string(ProblemNumber);
348 int pn = 0; // to avoid unsigned short overflow in FIRE
349
350 if(!file_exists(WorkingDir+"/"+spn+".start") || !file_exists(WorkingDir+"/"+spn+".config")) {
351
352 if(Integral.nops()<1) return;
353 lst InExternal;
354 for(auto ii : Internal) InExternal.append(ii);
355 for(auto ii : External) InExternal.append(ii);
356
357 if(ISP.nops()<1) {
358 for(auto it : Internal) {
359 for(auto ii : InExternal) ISP.append(it*ii);
360 }
361 ISP.sort();
362 ISP.unique();
363 }
364
365 if(ISP.nops() > pdim) {
366 cout << "ISP = " << ISP << endl;
367 cout << "Propagator = " << Propagator << endl;
368 throw Error("FIRE::Export: #(ISP) > #(Propagator).");
369 }
370
371 lst sp2s, s2sp, ss;
372 int _pic=0;
373 for(auto item : ISP) {
374 _pic++;
375 Symbol si("P"+to_string(_pic));
376 ss.append(si);
377 sp2s.append(w*item==w*si);
378 sp2s.append(item==si);
379 s2sp.append(si==item);
380 }
381
382 if(!has_w(Replacement)) {
383 lst repl;
384 for(auto item : Replacement) {
385 repl.append(item);
386 repl.append(w*item.op(0) == w*item.op(1));
387 }
388 Replacement = repl;
389 }
390
391 lst eqns;
392 for(int i=0; i<ISP.nops(); i++) { // note NOT pdim
393 auto eq = expand(Propagator.op(i)).subs(iEpsilon==0); // drop iEpsilon
394 eq = eq.subs(sp2s);
395 eq = eq.subs(Replacement);
396 if(eq.has(iWF(w))) throw Error("FIRE::Export, iWF used in eq.");
397 eqns.append(eq == iWF(i));
398 }
399 auto s2p = lsolve(eqns, ss);
400 if(s2p.nops() != ISP.nops()) {
401 cout << "ISP=" << ISP << endl;
402 cout << "Propagator=" << Propagator << endl;
403 cout << "eqns=" << eqns << endl;
404 throw Error("FIRE::Export: lsolve failed.");
405 }
406
407 // FIRE version
408 if(true)
409 if(DSP.nops()<1) {
410 for(auto p1 : Internal)
411 for(auto p2 : InExternal)
412 DSP.append(lst{p1,p2});
413 }
414
415 // Lee version
416 if(false)
417 if(DSP.nops()<1) {
418 for(int i=0; i<Internal.nops(); i++)
419 DSP.append(lst{Internal.op(i), Internal.op(i==Internal.nops()-1 ? 0 : i+1)});
420 for(auto p2 : InExternal) DSP.append(lst{Internal.op(0), p2});
421 DSP.append(lst{Internal, Internal});
422 allIBP = true;
423 }
424
425 vector<exmap> ibps;
426 exvector IBPvec;
427 lst ns0;
428 for(int i=0; i<pdim; i++) ns0.append(0);
429
430 for(auto sp : DSP) {
431 auto ilp_lst = sp.op(0);
432 auto iep_lst = sp.op(1);
433 if(!is_a<lst>(ilp_lst)) {
434 ilp_lst = lst{ ilp_lst };
435 iep_lst = lst{ iep_lst };
436 }
437
438 exmap nc_map;
439 for(int ilst=0; ilst<ilp_lst.nops(); ilst++) {
440 symbol ss;
441 auto ilp = ilp_lst.op(ilst);
442 auto iep = iep_lst.op(ilst);
443 lst dp_lst;
444 for(int i=0; i<pdim; i++) {
445 dp_lst.append(Propagator.op(i).subs(ilp==ss).diff(ss).subs(ss==ilp));
446 }
447
448 for(int i=0; i<pdim; i++) { // diff on each propagator
449 auto ns = ns0;
450 ns.let_op(i) = ns.op(i)+1; // note the covention
451 auto tmp = dp_lst.op(i) * iep;
452 tmp = expand(tmp);
453 tmp = tmp.subs(Replacement);
454 tmp = tmp.subs(sp2s);
455 tmp = tmp.subs(s2p);
456 tmp = tmp.subs(Replacement);
457
458 tmp = ex(0) - (Shift[i+1]+a(i+1))*tmp; // note Shift here
459
460 for(int j=0; j<pdim; j++) {
461 auto cj = tmp.coeff(iWF(j));
462 if(is_zero(cj)) continue;
463 lst cns = ns;
464 cns.let_op(j) = cns.op(j)-1; // note the covention
465 nc_map[cns] = nc_map[cns] + cj;
466 }
467 tmp = tmp.subs(iWF(w)==0); // constant term
468 if(!is_zero(tmp)) nc_map[ns] = nc_map[ns] + tmp;
469 }
470
471 auto cilp = iep.coeff(ilp);
472 if(!is_zero(cilp)) nc_map[ns0] = nc_map[ns0] + d*cilp;
473 }
474
475 bool ok = false;
476 for(auto nc : nc_map) {
477 if(!is_zero(nc.second)) {
478 ok = true;
479 IBPvec.push_back(nc.second);
480 }
481 }
482 if(ok) ibps.push_back(nc_map);
483 }
484
485 auto Variables = gather_symbols(IBPvec);
486 for(auto e : External) {
487 if(Variables.has(e)) {
488 cout << "IBPvec: " << IBPvec << endl;
489 cout << "Replacement: " << Replacement << endl;
490 cout << "Variables: " << Variables << endl;
491 throw Error("FIRE: Replacement NOT work!");
492 }
493 }
494
495 ostringstream start;
496
497 start << "ExampleDimension[" << pn << "]=" << pdim << endl << endl;
498 start << "ProblemNumber=" << pn << endl << endl;
499
500 // .start - SBasisL
501 if(Version==5) {
502 PermutationsR(2, pdim, [pdim,pn,&start](const int *ns) {
503 start << "SBasisL[" << pn << ",{";
504 for(int i=0; i<pdim; i++) start << (ns[i]<1 ? -1 : 1) << (i<pdim-1 ? "," : "");
505 start << "}]=0" << endl << endl;
506 });
507 }
508
509 // .start - SBasis0L, SBasis0D && SBasis0C
510 start << "SBasis0L[" << pn << "]=" << ibps.size() << endl << endl;
511 ostringstream oss;
512 for(int i=0; i<ibps.size(); i++) {
513 start << "SBasis0D[" << pn << "," << (i+1) << "]=";
514 lst items;
515 for(auto kv : ibps[i]) {
516 items.append(kv.first);
517 if(Version==5) {
518 oss << "SBasis0C[" << pn << "," << (i+1) << "," << kv.first << "]=" <<
519 collect_common_factors(kv.second.normal()) << endl << endl;
520 } else {
521 lst olst;
522 auto cv_lst = collect_lst(kv.second, a(w));
523 for(auto item : cv_lst) {
524 auto cc = item.op(0);
525 auto cv = item.op(1);
526 cc = collect_common_factors(cc.normal());
527 if(is_zero(cc)) continue;
528 if(is_zero(cv-1)) cv=0;
529 else cv = cv.subs(a(w)==w);
530 olst.append(lst{cc, cv});
531 }
532 oss << "SBasis0C[" << pn << "," << (i+1) << "," << kv.first << "]=" << olst << endl << endl;
533 }
534 }
535 start << items << endl << endl;
536 }
537 start << oss.str();
538
539 // .start - SBasisS
540 start << "SBasisS[" << pn << "]={{{";
541 // 1,2,3,...
542 for(int i=0; i<pdim; i++) start << (i+1) << (i<pdim-1 ? "," : "");
543 start << "},{";
544 // 1,1,1,...
545 for(int i=0; i<pdim; i++) start << 1 << (i<pdim-1 ? "," : "");
546 start << "},{";
547 // 0,0,0,...
548 for(int i=0; i<pdim; i++) start << 0 << (i<pdim-1 ? "," : "");
549 start << "}}}" << endl << endl;
550
551 // .start - SBasisR
552 lst Rlst;
553 Rlst.append(lst{});
554 for(int i=0; i<pdim; i++) {
555 let_op_append(Rlst, 0, -1);
556 }
557 for(auto lpi : Internal) {
558 vector<int> ns_vec;
559 lst ns0;
560 for(int i=0; i<pdim; i++) ns0.append(1);
561 for(int i=0; i<pdim; i++) {
562 if(Propagator.op(i).has(lpi)) ns0.let_op(i) = -1;
563 else ns_vec.push_back(i);
564 }
565 size_t tot = std::pow(2LL,ns_vec.size());
566 for(size_t n=0; n<tot; n++) {
567 int cn = n;
568 lst ns1 = ns0;
569 for(int j=0; j<ns_vec.size(); j++) {
570 if((cn%2)==1) ns1.let_op(ns_vec[j]) = -1;
571 cn /= 2;
572 }
573 Rlst.append(ns1);
574 }
575 }
576
577 // Lee Zero Sector
578 if(true) {
579 exset sectors;
580 IsAlwaysZero = true;
581 lst ns0;
582 for(int i=0; i<pdim; i++) ns0.append(1);
583 size_t tot = std::pow(2LL,pdim);
584
585 for(size_t n=0; n<tot; n++) {
586 int cn = n;
587 lst sector = ns0;
588 for(int j=0; j<pdim; j++) {
589 if((cn%2)==1) sector.let_op(j) = 0;
590 cn /= 2;
591 }
592 if(SECTOR.nops()==pdim) {
593 for(int j=0; j<pdim; j++) if(sector.op(j)>SECTOR.op(j)) goto add_sector_done;
594 }
595 sectors.insert(sector);
596 add_sector_done: ;
597 }
598
599 while(!sectors.empty()) {
600 auto first = *(sectors.begin());
601 sectors.erase(first);
602
603 auto sector = ex_to<lst>(first);
604 if(IsZero(sector)) { // from LiteRed: all subsector is zero
605 lst ns1 = sector;
606 for(int j=0; j<pdim; j++) {
607 if(sector.op(j).is_zero()) ns1.let_op(j) = -1;
608 }
609 Rlst.append(ns1);
610 exset subsets;
611 for(auto item : sectors) {
612 bool ok = true;
613 for(int j=0; j<pdim; j++) {
614 if(sector.op(j)<item.op(j)) { ok = false; break; }
615 }
616 if(ok) subsets.insert(item);
617 }
618 for(auto item : subsets) {
619 sectors.erase(item);
620 lst ns1 = ex_to<lst>(item);
621 for(int j=0; j<pdim; j++) {
622 if(item.op(j).is_zero()) ns1.let_op(j) = -1;
623 }
624 Rlst.append(ns1);
625 }
626 } else { // from LiteRed: all supersector is non-zero
627 if(IsAlwaysZero) IsAlwaysZero = false;
628 exset supsets;
629 for(auto item : sectors) {
630 bool ok = true;
631 for(int j=0; j<pdim; j++) {
632 if(sector.op(j)>item.op(j)) { ok = false; break; }
633 }
634 if(ok) supsets.insert(item);
635 }
636 for(auto item : supsets) sectors.erase(item);
637 }
638 }
639
640 if(IsAlwaysZero) {
641 Rules.remove_all();
642 for(auto ii : Integral) Rules.append(F(ProblemNumber, ii)==0);
643 return;
644 }
645 }
646
647 // handle Cut Propagator
648 if(Cut.nops()>0) {
649 for(auto cx : Cut) {
650 int ci = ex_to<numeric>(cx-1).to_int(); // start from 1 in Cut
651 lst ns0;
652 for(int i=0; i<pdim; i++) ns0.append(1);
653 ns0.let_op(ci) = -1;
654 for(int n=0; n<std::pow(2,pdim-1); n++) {
655 int cn = n;
656 lst ns1 = ns0;
657 for(int j=0; j<pdim; j++) {
658 if(ci==j) continue;
659 if((cn%2)==1) ns1.let_op(j) = -1;
660 cn /= 2;
661 }
662 Rlst.append(ns1);
663 }
664 }
665 }
666
667 Rlst.sort();
668 Rlst.unique();
669 //sort_lst(Rlst);
670
671 for(auto iR : Rlst) {
672 start << "SBasisR[" << pn << ",{";
673 for(int i=0; i<pdim; i++) start << iR.op(i) << (i<pdim-1 ? "," : "");
674 start << "}]=True" << endl << endl;
675 }
676
677 // .start - Others
678 start << "SBasisRL[" << pn << "]=0" << endl << endl;
679 start << "HPI[" << pn << "]={}" << endl << endl;
680
681 string sss = start.str();
682 string_replace_all(sss, "=", " = ");
683 string_replace_all(sss, ",", ", ");
684
685 if(!dir_exists(WorkingDir)) auto rc = system(("mkdir -p " + WorkingDir).c_str());
686
687 ofstream start_out(WorkingDir+"/"+spn+".start");
688 start_out << sss << endl;
689 start_out.close();
690
691 // .config
692 {
693 ostringstream config;
694
695 config << "#variables ";
696 bool first = true;
697 exvector ev_sort;
698 for(auto v : Variables) {
699 auto fw = fermat_weight.find(v);
700 if(fw!=fermat_weight.end()) ev_sort.push_back(lst{ fw->second, v });
701 else ev_sort.push_back(lst{ 0, v });
702 }
703 sort_vec(ev_sort);
704 for(auto nv : ev_sort) {
705 const symbol & s = ex_to<symbol>(nv.op(1));
706 if(false && !islower(s.get_name()[0])) {
707 cout << "Replacement: " << Replacement << endl;
708 cout << "IBPvec: " << IBPvec << endl;
709 throw Error("FIRE: Fermat requires a name must begin with a lower case letter: "+s.get_name());
710 }
711 config << (first ? "" : ",") << s;
712 auto itr = NVariables.find(nv.op(1));
713 if(itr!=NVariables.end()) config << "->" << itr->second;
714 first=false;
715 }
716 config << endl;
717
718 if(PosPref!=1) config << "#pos_pref "<< PosPref << endl;
719 config << "#database db" << ProblemNumber << endl;
720 config << "##prime 111" << endl; // ## for comment
721 if(allIBP) config << "#allIBP" << endl;
722 config << "#start" << endl;
723
724 if(SECTOR.nops()==pdim) {
725 int pmax = pdim;
726 for(int i=pdim-1; i>=0; i--) {
727 if(SECTOR.op(i)!=0) {
728 pmax = i+1;
729 break;
730 }
731 }
732 int pmin = -1;
733 for(int i=0; i<pdim; i++) {
734 if(SECTOR.op(i)!=0) {
735 pmin = i+1;
736 break;
737 }
738 }
739 config << "#problem " << pn << " ";
740 if(pmin>1) config << "|" << pmin << "," << pmax << "|";
741 else if(pmax!=pdim) config << "|" << pmax << "|";
742 config << ProblemNumber << ".start" << endl;
743 } else config << "#problem " << pn << " " << ProblemNumber << ".start" << endl;
744
745 config << "#output " << ProblemNumber << ".tables" << endl;
746
747 if(PIntegral.nops()>0) {
748 ostringstream oss;
749 oss << "{";
750 int nn = PIntegral.nops();
751 for(int i=0; i<nn; i++) {
752 if(PIntegral.op(i).nops()!=pdim) throw Error("FIRE::Export@1, Index dimension NOT match Propagator.");
753 oss << "{" << pn << "," << PIntegral.op(i) << (i<nn-1 ? "}," : "}");
754 }
755 oss << "}";
756 ofstream pref_out(WorkingDir+"/"+spn+".pref");
757 pref_out << oss.str() << endl;
758 pref_out.close();
759 config << "#preferred " << ProblemNumber << ".pref" << endl;
760 }
761
762 config << "#integrals " << ProblemNumber << ".intg" << endl;
763
764 ofstream config_out(WorkingDir+"/"+spn+".config");
765 config_out << config.str() << endl;
766 config_out.close();
767 }
768
769 }
770
771 // *.intg
772 if(!file_exists(WorkingDir+"/"+spn+".intg")) {
773 ostringstream intg;
774 intg << "{";
775 for(int i=0; i<Integral.nops(); i++) {
776 if(Integral[i].nops()!=pdim) throw Error("FIRE::Export@2, Index dimension NOT match Propagator.");
777 intg << "{" << pn << "," << Integral[i] << (i<Integral.nops()-1 ? "}," : "}");
778 }
779 intg << "}" << endl;
780
781 ofstream intg_out(WorkingDir+"/"+spn+".intg");
782 intg_out << intg.str() << endl;
783 intg_out.close();
784 }
785
786 // export .gar here
787 if(!file_exists(WorkingDir+"/"+spn+".gar")) {
788 garWrite(TO(), WorkingDir+"/"+spn+".gar");
789 }
790 }
791
795 void FIRE::Run() {
796 if(IsAlwaysZero) return;
797 if(Execute=="") Execute = InstallPrefix + "/FIRE/M/FIRE";
798 if(file_exists(WorkingDir + "/" + to_string(ProblemNumber) + ".tables")) return;
799 if(Integral.nops()<1) return;
800 int tried = 0;
801 while(tried<3) {
802 tried++;
803 ostringstream cmd;
804 cmd << "cd " << WorkingDir << " && " << Execute;
805 if(Version>5) {
806 if(opt=="") cmd << " -silent -t " << T1 << "," << T2 << " -lt " << LT1 << "," << LT2;
807 else cmd << " " << opt;
808 }
809 cmd << " -c " << ProblemNumber << " 1>/dev/null 2>&1";
810 auto rc = system(cmd.str().c_str());
811 rc = system(("rm -rf "+WorkingDir+"/db"+to_string(ProblemNumber)).c_str());
812 if(file_exists(WorkingDir + "/" + to_string(ProblemNumber) + ".tables")) break;
813 sleep(3);
814 }
815 if(!file_exists(WorkingDir + "/" + to_string(ProblemNumber) + ".tables")) {
816 cout << "Propagator: " << Propagator << endl;
817 throw Error("FIRE::Run failed!");
818 }
819 }
820
824 void FIRE::Import() {
825 if(IsAlwaysZero) return;
826 if(Integral.nops()<1) return;
827 string spn = to_string(ProblemNumber);
828 ifstream ifs(WorkingDir+"/"+spn+".tables");
829 string ostr((istreambuf_iterator<char>(ifs)), (istreambuf_iterator<char>()));
830 ifs.close();
831
832 //string_replace_all(ostr, "\"", "");
833 for(auto & c : ostr) if(c=='"') c = ' ';
834
835 Parser tp;
836 auto tp_lst = tp.Read(ostr);
837 exmap id2F;
838 for(auto item : tp_lst.op(1)) {
839 if(!is_zero(item.op(1).op(0))) throw Error("FIRE::Reduce: pn is NOT 0.");
840 id2F[item.op(0)] = F(ProblemNumber, item.op(1).op(1));
841 }
842
843 for(auto item : tp_lst.op(0)) {
844 ex left = item.op(0).subs(id2F);
845 ex right = 0;
846 for(auto it : item.op(1)) {
847 right += it.op(0).subs(id2F) * it.op(1);
848 }
849 if(left.is_equal(right)) MIntegral.append(left);
850 else Rules.append(left==right);
851 }
852 MIntegral.sort();
853 MIntegral.unique();
854 sort_lst(MIntegral);
855
856 // handle Cut not equal 1, using #preferred
857 if(reCut && Cut.nops()>0) {
858 lst ois, vis;
859 auto mis = MIntegral;
860 MIntegral.remove_all();
861 int nrun = -1;
862 while(true) {
863 nrun++;
864 if(nrun>100) break;
865
866 lst iis, pis;
867 bool allOK = true;
868 for(auto item : mis) {
869 lst mi = ex_to<lst>(item.op(1));
870 bool isOK = true;
871 for(auto cx : Cut) {
872 if(!is_zero(mi.op(ex_to<numeric>(cx).to_int()-1)-1)) {
873 isOK = false;
874 allOK = false;
875 break;
876 }
877 }
878 if(nrun==0 && isOK) MIntegral.append(item);
879 else if(!isOK) {
880 if(nrun==0) {
881 ois.append(item);
882 vis.append(item);
883 }
884 iis.append(mi);
885 auto pi = mi;
886 vector<int> ipos, ineg;
887 for(int i=0; i<mi.nops(); i++) {
888 bool isCut = false;
889 for(auto cx : Cut) {
890 if(is_zero(cx-1-i)) {
891 isCut = true;
892 break;
893 }
894 }
895 if(isCut) pi.let_op(i)=1;
896 else if(mi.op(i)<=0) ineg.push_back(i);
897 else ipos.push_back(i);
898 }
899
900 lst tis;
901 int max = 2;
902 ex total = pow(numeric(max), ineg.size());
903 for(numeric in=0; in<total; in++) {
904 auto cin = in;
905 auto pi2 = pi;
906 for(int i=0; i<ineg.size(); i++) {
907 int re = mod(cin,max).to_int();
908 pi2.let_op(ineg[i]) = ex(0)-re;
909 cin = (cin-re)/max;
910 }
911 tis.append(pi2);
912 }
913
914 int max2 = 2*max+1;
915 total = pow(numeric(max2), ipos.size());
916 for(auto pi : tis) {
917 for(numeric in=0; in<total; in++) {
918 auto cin = in;
919 auto pi2 = pi;
920 for(int i=0; i<ipos.size(); i++) {
921 int re = mod(cin,max2).to_int();
922 pi2.let_op(ipos[i]) = re-max;
923 cin = (cin-re)/max2;
924 }
925 pis.append(pi2);
926 }
927 }
928 }
929 }
930 if(allOK) break;
931
932 // Reduce again
933 pis.sort();
934 pis.unique();
935 iis.sort();
936 iis.unique();
937 if(pis.nops()>0) {
938 FIRE fire;
939 fire.Propagator = Propagator;
940 fire.Internal = Internal;
941 fire.External = External;
942 fire.Replacement = Replacement;
943 fire.ProblemNumber = ProblemNumber;
944 fire.ISP = ISP;
945 fire.DSP = DSP;
946 fire.Cut = Cut;
947 fire.reCut = false;
948 fire.Shift = Shift;
949 fire.Integral = iis;
950 fire.PIntegral = pis;
951 fire.WorkingDir = WorkingDir + "_C"+to_string(ProblemNumber);
952 fire.Reduce();
953 auto rc = system(("rm -rf " + fire.WorkingDir).c_str());
954
955 auto vis_chk = vis;
956 vis_chk = ex_to<lst>(subs(vis_chk, fire.Rules));
957 if(vis_chk.is_equal(vis)) break;
958 vis = vis_chk;
959 exset fs;
960 find(vis,F(w1,w2),fs);
961 mis.remove_all();
962 for(auto item : fs) mis.append(item);
963 } else break;
964 }
965
966 exmap c2m;
967 for(int i=0; i<ois.nops(); i++) c2m[ois.op(i)] = vis.op(i);
968 for(auto item : mis) MIntegral.append(item);
969 MIntegral.sort();
970 MIntegral.unique();
971 sort_lst(MIntegral);
972
973 auto rules = Rules;
974 Rules.remove_all();
975 lst rIntegral;
976 for(auto r : rules) {
977 auto ri = r.op(0);
978 bool isMI = false;
979 for(auto item : MIntegral) {
980 if(item.is_equal(ri)) {
981 isMI = true;
982 rIntegral.append(ri);
983 break;
984 }
985 }
986 if(isMI) continue;
987 auto rv = r.op(1);
988 rIntegral.append(rv);
989 Rules.append(ri==rv.subs(c2m));
990 }
991
992 rIntegral.sort();
993 rIntegral.unique();
994 exset fs;
995 find(rIntegral,F(w1,w2),fs);
996 MIntegral.remove_all();
997 for(auto fi : fs) MIntegral.append(fi);
998 sort_lst(MIntegral);
999 }
1000 }
1001
1002}
a
int * a
Definition Functions.cpp:234
IBP.h
IBP header file.
file_exists
bool file_exists(const char *fn)
Definition Process.cpp:9
HepLib::Error
class used to wrap error message
Definition BASIC.h:242
HepLib::FIRE
IBP reduction using FIRE program.
Definition IBP.h:68
HepLib::FIRE::RRTables
static void RRTables(const string &filename, int pnum)
Definition FIRE.cpp:72
HepLib::FIRE::ThieleTables
static void ThieleTables(const string &filename, int si, int ei)
Definition FIRE.cpp:168
HepLib::FIRE::Export
void Export() override
Export start config intgral etc. files.
Definition FIRE.cpp:344
HepLib::IBP::Internal
lst Internal
Definition IBP.h:27
HepLib::IBP::Shift
map< int, ex > Shift
Definition IBP.h:36
HepLib::IBP::PIntegral
lst PIntegral
Definition IBP.h:40
HepLib::IBP::ISP
lst ISP
Definition IBP.h:34
HepLib::IBP::Replacement
lst Replacement
Definition IBP.h:29
HepLib::IBP::Reduce
void Reduce()
HepLib::IBP::DSP
lst DSP
Definition IBP.h:33
HepLib::IBP::reCut
bool reCut
Definition IBP.h:37
HepLib::IBP::Cut
lst Cut
Definition IBP.h:32
HepLib::IBP::Integral
lst Integral
Definition IBP.h:31
HepLib::IBP::External
lst External
Definition IBP.h:28
HepLib::IBP::Propagator
lst Propagator
Definition IBP.h:30
HepLib::IBP::ProblemNumber
int ProblemNumber
Definition IBP.h:39
HepLib::IBP::WorkingDir
string WorkingDir
Definition IBP.h:38
HepLib::IBP::Rules
lst Rules
Definition IBP.h:43
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::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::gather_symbols
lst gather_symbols(const ex &e)
get all symbols from input expression
Definition BASIC.cpp:539
HepLib::Version
string Version
Definition Init.cpp:85
HepLib::fermat_weight
map< ex, long long, ex_is_less > fermat_weight
Definition Init.cpp:155
HepLib::RationalReconstruct
numeric RationalReconstruct(numeric a, numeric p)
Definition Rational.cpp:8
HepLib::file_exists
bool file_exists(string fn)
Definition BASIC.h:289
HepLib::has_w
bool has_w(const ex &e)
Definition BASIC.cpp:2233
HepLib::syms
lst syms(const exvector &ev)
Definition FIRE.cpp:61
HepLib::w
ex w
Definition Init.cpp:93
HepLib::d
const Symbol d
HepLib::InstallPrefix
string InstallPrefix
Definition Init.cpp:162
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::sort_vec
void sort_vec(exvector &ivec, bool less=true)
sort the list in less order, or the reverse
Definition Sort.cpp:54
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::string_replace_all
void string_replace_all(string &str, const string &from, const string &to)
Definition Functions.cpp:148
HepLib::garWrite
void garWrite(const string &garfn, const map< string, ex > &resMap)
garWrite to write the string-key map to the archive
Definition BASIC.cpp:639
HepLib::let_op_append
void let_op_append(ex &ex_in, const ex item)
append item into expression
Definition BASIC.cpp:1324
HepLib::dir_exists
bool dir_exists(string dir)
Definition BASIC.h:297
HepLib::IsZero
bool IsZero(const ex &e)
HepLib::w1
ex w1
Definition BASIC.h:499
HepLib::w2
ex w2
Definition BASIC.h:499
HepLib::PermutationsR
void PermutationsR(int n, int m, std::function< void(const int *)> f)
Definition Functions.cpp:268
HepLib::subs
ex subs(const ex &e, init_list sl)
Definition BASIC.h:51
Generated by doxygen 1.9.8