HepLib
Loading...
Searching...
No Matches
CIPrepares.cpp
Go to the documentation of this file.
1
6#include "SD.h"
7#include <math.h>
8#include <cmath>
9
10namespace HepLib::SD {
11
17 void SecDec::CIPrepares(const string & key) {
18 if(expResult.size()<1) IsZero = true;
19 if(IsZero) return;
20 int rc;
21
22 if(Verbose > 0) cout << Color_HighLight << " CIPrepares @ " << now() << RESET << endl;
23 auto pid = getpid();
24 string spid = to_string(pid);
25 if(get_env("GiNaC_Parallel_Host").length()>0) spid = get_env("GiNaC_Parallel_Host")+"_"+spid;
26
27 GiNaC_Parallel_RM["FCI-F"] = false;
28 auto resf =
29 GiNaC_Parallel(expResult.size(), [this](int idx)->ex {
30 // return lst{ kv.op(0), kv.op(1), ft};
31 auto kv = expResult[idx];
32 auto expr = kv.op(1);
33 auto xs = get_xy_from(expr);
34 if(xs.size()<1) {
35 return lst{kv.op(0), kv.op(1), 1};
36 }
37
38 exset ftxset;
39 expr.find(FTX(w1,w2), ftxset);
40 ex ft;
41 int ftxsize = -1;
42 for(auto item : ftxset) {
43 auto xys = get_xy_from(item.op(0));
44 if((int)xys.size() > ftxsize) {
45 ft = item.op(0);
46 ftxsize = xys.size();
47 }
48 }
49
50 bool need_contour_deformation = false;
51 if(ft.has(x(w)) && !ft.has(PL(w))) {
52 exmap eps_map;
53 ex epn = ex(1)/111;
54 for(auto epi : eps_lst) {
55 eps_map[epi.op(0)] = epn;
56 epn = epn / 13;
57 }
58 auto tmp = ft.subs(nReplacement).subs(eps_map).subs(CV(w1,w2)==w2).expand();
59 if(is_a<add>(tmp)) {
60 for(auto item : tmp) {
61 if(!is_a<numeric>(item.subs(x(w)==1))) {
62 throw Error("CIPrepares: (!is_a<numeric>(item.subs(x(w)==1)))");
63 }
64 if(item.subs(x(w)==1) < 0) {
65 need_contour_deformation = true;
66 break;
67 }
68 }
69 } else {
70 if(!is_a<numeric>(tmp.subs(x(w)==1))) {
71 throw Error("CIPrepares: (!is_a<numeric>(tmp.subs(x(w)==1)))");
72 }
73 if(tmp.subs(x(w)==1) < 0) need_contour_deformation = true;
74 }
75 if(!need_contour_deformation) ft = 1; //note the difference with SDPrepare
76 } else if(!ft.has(x(w))) {
77 ft = 1;
78 }
79
80 ft = collect_common_factors(ft);
81 return lst{ kv.op(0), kv.op(1), ft};
82
83 }, "FCI-F");
84
85
86 //============================================================================================================
87 lst fts;
88 for(auto item : resf) {
89 if(item.op(2).has(x(w))) {
90 fts.append(item.op(2));
91 }
92 }
93 fts.sort();
94 fts.unique();
95
96 exvector ftnvec;
97 map<ex,int,ex_is_less> ftnmap;
98 int ft_n = 1;
99 FT_N_XN.remove_all();
100
101 for(auto item : fts) {
102 ftnvec.push_back(lst{item, ft_n});
103 ftnmap[item] = ft_n;
104 FT_N_XN.append(lst{item, ft_n, get_xy_from(item).size()});
105 ft_n++;
106 }
107 //ftnvec item: lst { ft, ft-id }
108
109 exvector res_vec;
110 exmap cf_int;
111
112 for(auto &item : resf) {
113 auto ii = ex_to<lst>(item);
114 if(ii.op(2)==1) {
115 ii.append(-1);
116 } else {
117 int ft_n = ftnmap[item.op(2)];
118 if(ft_n==0) throw Error("CIPrepares: ft_n==0, " + ex2str(item.op(2)));
119 ii.append(ft_n);
120 }
121 if(ii.op(0).is_zero() || ii.op(1).subs(FTX(w1,w2)==1).is_zero()) continue;
122 if(IBF==0) res_vec.push_back(ii);
123 else if(IBF==1) { // by coefficient and F-term
124 ex key = ii;
125 key[1] = 1;
126 cf_int[key] += ii.op(1);
127 } else { // by F-term
128 auto cvs = collect_lst(ii.op(0), [](const ex & e)->bool { return has_symbol(e); });
129 for(auto const & cv : cvs) {
130 ex key = ii;
131 key[0] = cv.op(1);
132 key[1] = 1;
133 cf_int[key] += cv.op(0)*ii.op(1);
134 }
135 }
136 }
137
138 if(IBF!=0) {
139 for(auto kv : cf_int) {
140 lst ii = ex_to<lst>(kv.first);
141 ii[1] = kv.second;
142 res_vec.push_back(ii);
143 }
144 }
145
146 //res_vec item: lst { coeff, integrand, ft, ft-id }
147
148 if(res_vec.size()<1) {
149 IsZero = true;
150 return;
151 }
152 //============================================================================================================
153
154 // Prepare FT-lambda
155 GiNaC_Parallel_RM["FCI-C"] = false;
156 GiNaC_Parallel(ftnvec.size(), [&ftnvec,spid](int idx)->ex {
157 // return nothing
158 auto kv = ftnvec[idx];
159 ex ft = kv.op(0);
160 ex ft_n = kv.op(1);
161 auto fxs = get_xy_from(ft);
162 lst las;
163
164 auto pls = get_pl_from(ft);
165 int npls = pls.size()>0 ? ex_to<numeric>(pls[pls.size()-1].subs(lst{PL(w)==w})).to_int() : -1;
166 lst plRepl;
167 for(int i=0; i<npls+1; i++) {
168 ostringstream pl;
169 pl << "pl[" << i << "]";
170 plRepl.append(PL(i) == Symbol(pl.str()));
171 }
172
173 ex DFs[fxs.size()], DDFs[fxs.size()][fxs.size()];
174 for(int i=0; i<fxs.size(); i++) {
175 auto df = diff_ex(ft, fxs[i]);
176 DFs[i] = collect_common_factors(df);
177 ostringstream ilaos;
178 ilaos << "ilas[" << i << "]";
179 symbol ila(ilaos.str());
180 for(int j=0; j<fxs.size(); j++) {
181 auto ddf = diff_ex(DFs[i], fxs[j]);
182 DDFs[i][j] = collect_common_factors(ddf);
183 }
184 }
185
186 if(!dir_exists(spid)) auto rc = system(("mkdir -p "+spid).c_str());
187 ostringstream cppfn, sofn;
188 cppfn << spid << "/F" << ft_n << ".cpp";
189 sofn << spid << "/F" << ft_n << ".o";
190 std::ofstream ofs;
191 ofs.open(cppfn.str(), ios::out);
192 if (!ofs) throw Error("failed to open *.cpp file! (1)");
193
194 lst cxRepl, czRepl;
195 for (int i=0; i<fxs.size(); i++) {
196 ostringstream sx, sz;
197 sx << "x[" << i << "]";
198 cxRepl.append(fxs[i] == Symbol(sx.str()));
199 sz << "z[" << i << "]";
200 czRepl.append(fxs[i] == Symbol(sz.str()));
201 }
202
203 /*----------------------------------------------*/
204 ofs << "#include \"NFunctions.h\"" << endl;
205 /*----------------------------------------------*/
206 auto cppL = CppFormat(ofs, "L");
207 auto cppQ = CppFormat(ofs, "Q");
208 auto cppMP = CppFormat(ofs, "MP");
209
210 // FL_fid
211 ofs << "dREAL FL_" << ft_n << "(const dREAL* x, const dREAL *pl) {" << endl;
212 ofs << "dREAL yy = ";
213 EvalL(ft.subs(plRepl).subs(cxRepl)).print(cppL);
214 ofs << ";" << endl;
215 ofs << "return yy;" << endl;
216 ofs << "}" << endl;
217 ofs << endl;
218
219 // FQ_fid
220 ofs << "qREAL FQ_" << ft_n << "(const qREAL* x, const qREAL *pl) {" << endl;
221 ofs << "qREAL yy = ";
222 EvalQ(ft.subs(plRepl).subs(cxRepl)).print(cppQ);
223 ofs << ";" << endl;
224 ofs << "return yy;" << endl;
225 ofs << "}" << endl;
226 ofs << endl;
227
228 // FMP_fid
229 ofs << "mpREAL FMP_" << ft_n << "(const mpREAL* x, const mpREAL *pl) {" << endl;
230 ofs << "mpREAL yy = ";
231 EvalMP(ft.subs(plRepl).subs(cxRepl)).print(cppMP);
232 ofs << ";" << endl;
233 ofs << "return yy;" << endl;
234 ofs << "}" << endl;
235 ofs << endl;
236
237 // D's FL_fid
238 ofs << "dREAL FL_" << ft_n << "(const int i, const dREAL* x, const dREAL *pl) {" << endl;
239 for(int i=0; i<fxs.size(); i++) {
240 ofs << "if("<<i<<"==i) return ";
241 EvalL(DFs[i].subs(plRepl).subs(cxRepl)).print(cppL);
242 ofs << ";" << endl;
243 }
244 ofs << "return 0;" << endl;
245 ofs << "}" << endl;
246 ofs << endl;
247
248 // D's FQ_fid
249 ofs << "qREAL FQ_" << ft_n << "(const int i, const qREAL* x, const qREAL *pl) {" << endl;
250 for(int i=0; i<fxs.size(); i++) {
251 ofs << "if("<<i<<"==i) return ";
252 EvalQ(DFs[i].subs(plRepl).subs(cxRepl)).print(cppQ);
253 ofs << ";" << endl;
254 }
255 ofs << "return 0;" << endl;
256 ofs << "}" << endl;
257 ofs << endl;
258
259 // D's FMP_fid
260 ofs << "mpREAL FMP_" << ft_n << "(const int i, const mpREAL* x, const mpREAL *pl) {" << endl;
261 for(int i=0; i<fxs.size(); i++) {
262 ofs << "if("<<i<<"==i) return ";
263 EvalMP(DFs[i].subs(plRepl).subs(cxRepl)).print(cppMP);
264 ofs << ";" << endl;
265 }
266 ofs << "return 0;" << endl;
267 ofs << "}" << endl;
268 ofs << endl;
269
270 // DD's FL_fid
271 ofs << "dREAL FL_" << ft_n << "(const int i, const int j, const dREAL* x, const dREAL *pl) {" << endl;
272 for(int i=0; i<fxs.size(); i++) {
273 for(int j=0; j<fxs.size(); j++) {
274 ofs << "if("<<i<<"==i && "<<j<<"==j) return ";
275 EvalL(DDFs[i][j].subs(plRepl).subs(cxRepl)).print(cppL);
276 ofs << ";" << endl;
277 }}
278 ofs << "return 0;" << endl;
279 ofs << "}" << endl;
280 ofs << endl;
281
282 // DD's FQ_fid
283 ofs << "qREAL FQ_" << ft_n << "(const int i, const int j, const qREAL* x, const qREAL *pl) {" << endl;
284 for(int i=0; i<fxs.size(); i++) {
285 for(int j=0; j<fxs.size(); j++) {
286 ofs << "if("<<i<<"==i && "<<j<<"==j) return ";
287 EvalQ(DDFs[i][j].subs(plRepl).subs(cxRepl)).print(cppQ);
288 ofs << ";" << endl;
289 }}
290 ofs << "return 0;" << endl;
291 ofs << "}" << endl;
292 ofs << endl;
293
294 // DD's FMP_fid
295 ofs << "mpREAL FMP_" << ft_n << "(const int i, const int j, const mpREAL* x, const mpREAL *pl) {" << endl;
296 for(int i=0; i<fxs.size(); i++) {
297 for(int j=0; j<fxs.size(); j++) {
298 ofs << "if("<<i<<"==i && "<<j<<"==j) return ";
299 EvalMP(DDFs[i][j].subs(plRepl).subs(cxRepl)).print(cppMP);
300 ofs << ";" << endl;
301 }}
302 ofs << "return 0;" << endl;
303 ofs << "}" << endl;
304 ofs << endl;
305
306 // X2ZL_fid
307 ofs << "void X2ZL_" << ft_n << "(const dREAL* x, dCOMPLEX* z, dCOMPLEX* r, dREAL* dff, const dREAL* pl, const dREAL* las) {" << endl;
308 ofs << "X2Z("<<fxs.size()<<",FL_"<<ft_n<<",FL_"<<ft_n<<",x,z,r,dff,pl,las);" << endl;
309 ofs << "}" << endl;
310 ofs << endl;
311
312 // X2ZQ_fid
313 ofs << "void X2ZQ_" << ft_n << "(const qREAL* x, qCOMPLEX* z, qCOMPLEX* r, qREAL* dff, const qREAL* pl, const qREAL* las) {" << endl;
314 ofs << "X2Z("<<fxs.size()<<",FQ_"<<ft_n<<",FQ_"<<ft_n<<",x,z,r,dff,pl,las);" << endl;
315 ofs << "}" << endl;
316 ofs << endl;
317
318 // X2ZMP_fid
319 ofs << "void X2ZMP_" << ft_n << "(const mpREAL* x, mpCOMPLEX* z, mpCOMPLEX* r, mpREAL* dff, const mpREAL* pl, const mpREAL* las) {" << endl;
320 ofs << "X2Z("<<fxs.size()<<",FMP_"<<ft_n<<",FMP_"<<ft_n<<",x,z,r,dff,pl,las);" << endl;
321 ofs << "}" << endl;
322 ofs << endl;
323
324 // MatL_id
325 ofs << "void MatL_"<<ft_n<<"(dCOMPLEX* mat, const dREAL* x, const dREAL* dff, const dREAL* pl, const dREAL* las) {" << endl;
326 ofs << "Mat("<<fxs.size()<<",FL_"<<ft_n<<",mat,x,dff,pl,las);" << endl;
327 ofs << "}" << endl;
328 ofs << endl;
329
330 // MatQ_fid
331 ofs << "void MatQ_"<<ft_n<<"(qCOMPLEX *mat, const qREAL* x, const qREAL* dff, const qREAL *pl, const qREAL *las) {" << endl;
332 ofs << "Mat("<<fxs.size()<<",FQ_"<<ft_n<<",mat,x,dff,pl,las);" << endl;
333 ofs << "}" << endl;
334 ofs << endl;
335
336 // MatMP_fid
337 ofs << "void MatMP_"<<ft_n<<"(mpCOMPLEX *mat, const mpREAL* x, const mpREAL* dff, const mpREAL *pl, const mpREAL *las) {" << endl;
338 ofs << "Mat("<<fxs.size()<<",FMP_"<<ft_n<<",mat,x,dff,pl,las);" << endl;
339 ofs << "}" << endl;
340 ofs << endl;
341
342 // for Minimization of F(z)-image, x[xn-1] is the lambda
343 ofs << "extern \"C\" " << endl;
344 ofs << "dREAL imgF_"<<ft_n<<"(const int xn, const dREAL* x, const dREAL *pl, const dREAL *las_in) {" << endl;
345 ofs << "dREAL las[xn-1];" <<endl;
346 ofs << "for(int i=0; i<xn-1; i++) las[i] = las_in[i]*x[xn-1];" <<endl;
347 ofs << "dCOMPLEX z[xn], r[xn];" << endl;
348 ofs << "dREAL dff[xn+1];" << endl;
349 ofs << "X2ZL_"<<ft_n<<"(x,z,r,dff,pl,las);" << endl;
350 ofs << "dCOMPLEX zf = ";
351 EvalL(ft.subs(plRepl).subs(czRepl)).print(cppL);
352 ofs << ";" << endl;
353 ofs << "return -zf.imag()/x[xn-1];" << endl; // find max image part, check with 0
354 ofs << "}" << endl;
355 ofs << endl;
356
357 // for Minimization of F
358 ofs << "extern \"C\" " << endl;
359 ofs << "dREAL minF_"<<ft_n<<"(const int xn, const dREAL* x, const dREAL *pl, const dREAL *las_in) {" << endl;
360 ofs << "return FL_"<<ft_n<<"(x,pl);" << endl;
361 ofs << "}" << endl;
362 ofs << endl;
363
364 // for Minimization of -F
365 ofs << "extern \"C\" " << endl;
366 ofs << "dREAL minFM_"<<ft_n<<"(const int xn, const dREAL* x, const dREAL *pl, const dREAL *las_in) {" << endl;
367 ofs << "return 0.L-FL_"<<ft_n<<"(x,pl);" << endl;
368 ofs << "}" << endl;
369 ofs << endl;
370
371 // for Minimization of DF-i
372 for(int i=0; i<fxs.size(); i++) {
373 ofs << "extern \"C\" " << endl;
374 ofs << "dREAL dirC_"<<ft_n<<"_"<<i<<"(const int xn, const dREAL* x, const dREAL *pl, const dREAL *las) {" << endl;
375 ofs << "int i = " << i << ";" << endl;
376 ofs << "dREAL yy = x[i]*(1-x[i])*FL_"<<ft_n<<"(i, x, pl);" << endl;
377 ofs << "return -fabs(yy);" << endl;
378 ofs << "}" << endl;
379 ofs << endl;
380 }
381
382 ostringstream cmd;
383 cmd << cpp << " -fPIC -c " << INC_FLAGS << " -include NFunctions.h -o " << sofn.str() << " " << cppfn.str();
384 auto rc = system(cmd.str().c_str());
385
386 if(!file_exists(sofn.str().c_str())) {
387 cmd.clear();
388 cmd.str("");
389 cmd << "cp " << sofn.str() << " . ";
390 rc = system(cmd.str().c_str());
391 }
392
393 return 0;
394
395 }, "FCI-C");
396
397 bool hasF = (ftnvec.size()>0);
398 if(hasF) {
399 ostringstream sofn, cmd;
400 if(key != "") {
401 sofn << key << "F.so";
402 } else {
403 sofn << spid << "F.so";
404 }
405 cmd << cpp << " " << LIB_FLAGS << " -Wl,-rpath,. -rdynamic -fPIC -shared -lHepLib -lquadmath -lmpfr -lgmp " << " -o " << sofn.str() << " " << spid << "/F*.o";
406 cmd << " -lHepLib -lquadmath -lmpfr -lgmp";
407 cmd << " 1> /dev/null 2> /dev/null";
408 rc = system(cmd.str().c_str());
409
410 cmd.clear();
411 cmd.str("");
412 cmd << "rm -rf " << spid;
413 if(!Debug) rc = system(cmd.str().c_str());
414 }
415
416
417 //============================================================================================================
418 // Compile the null.o
419 if(true) {
420 ostringstream cmd;
421 if(!dir_exists(spid)) cmd << "mkdir -p " << spid;
422 rc = system(cmd.str().c_str());
423 cmd.clear();
424 cmd.str("");
425 cmd << "echo ''>" << spid << "/null.cpp;";
426 cmd << cpp << " -fPIC -c -o " << spid << "/null.o " << spid << "/null.cpp";
427 rc = system(cmd.str().c_str());
428 }
429
430 // Prepare Integrand
431 GiNaC_Parallel_RM["FCI-I"] = false;
432 auto res =
433 GiNaC_Parallel(res_vec.size(), [&res_vec,spid](int idx)->ex {
434 // return lst{ no-x-result, xn, x-indepent prefactor, ft_n }
435 // or lst{ id(SD(D|Q)_id in .so), xn, x-indepent prefactor, ft_n }
436
437 static symbol xwra("xwra");
438 auto kvf = res_vec[idx];
439 auto expr = kvf.op(1);
440 auto xs = get_xy_from(expr);
441 auto ft_n = kvf.op(3);
442 bool hasF = (ft_n>0);
443
444 if(xs.size()<1) {
445 ostringstream cmd;
446 cmd << "cp " << spid << "/null.o " << spid << "/" << idx << ".o";
447 auto rc = system(cmd.str().c_str());
448 return lst{ expr.subs(lst{FTX(w1,w2)==1}), xs.size(), kvf.op(0), -1};
449 }
450
451 if(xs.size()<1) xs.push_back(x(0));
452
453 auto ft = kvf.op(2);
454 auto fxs = get_xy_from(ft);
455
456 exset ftxset;
457 expr.find(FTX(w1,w2), ftxset);
458 lst ftxlst;
459 for(auto it : ftxset) ftxlst.append(it);
460 expr = collect_ex(expr, FTX(w1,w2));
461 vector<pair<ex,ex>> ft_expr;
462 for(auto item : ftxlst) {
463 ft_expr.push_back(make_pair(item.op(1), expr.coeff(item)));
464 }
465
466 lst cxRepl, czRepl, czzRepl, plRepl;
467 for (int i=0; i<fxs.size(); i++) {
468 ostringstream xs, zs, zzs;
469 xs << "x[" << i << "]";
470 zs << "z[" << i << "]";
471 zzs << "zz[" << i << "]";
472 cxRepl.append(fxs[i] == Symbol(xs.str()));
473 czRepl.append(fxs[i] == Symbol(zs.str()));
474 czzRepl.append(fxs[i] == Symbol(zzs.str()));
475 }
476 int count = fxs.size();
477 for(auto xi : xs) {
478 auto xii = xi.subs(czRepl);
479 if(is_zero(xii-xi)) {
480 ostringstream xs, zs;
481 xs << "x[" << count << "]";
482 cxRepl.append(xi == Symbol(xs.str()));
483 czRepl.append(xi == Symbol(xs.str()));
484 czzRepl.append(xi == Symbol(xs.str()));
485 count++;
486 }
487 }
488 if(count!=xs.size()) throw Error("CIPrepares: (count!=xs.size())");
489 auto pls = get_pl_from(expr);
490 int npls = pls.size()>0 ? ex_to<numeric>(pls[pls.size()-1].subs(lst{PL(w)==w})).to_int() : -1;
491 for(int i=0; i<npls+1; i++) {
492 ostringstream pl;
493 pl << "pl[" << i << "]";
494 plRepl.append(PL(i) == Symbol(pl.str()));
495 }
496
497 if(!dir_exists(spid)) auto rc = system(("mkdir -p "+spid).c_str());
498 ostringstream cppfn;
499 cppfn << spid << "/" << idx << ".cpp";
500 std::ofstream ofs;
501 ofs.open(cppfn.str(), ios::out);
502 if (!ofs) throw Error("failed to open *.cpp file! (2)");
503
504 /*----------------------------------------------*/
505 ofs << "#include \"NFunctions.h\"" << endl;
506 /*----------------------------------------------*/
507 if(hasF) {
508 ofs << "qREAL FQ_"<<ft_n<<"(const qREAL*, const qREAL*);" << endl; // for FT only
509 ofs << "dREAL FL_"<<ft_n<<"(const int, const dREAL*, const dREAL*);" << endl;
510 ofs << "qREAL FQ_"<<ft_n<<"(const int, const qREAL*, const qREAL*);" << endl;
511 ofs << "qREAL FMP_"<<ft_n<<"(const int, const mpREAL*, const mpREAL*);" << endl;
512 ofs << "dREAL FL_"<<ft_n<<"(const int, const int, const dREAL*, const dREAL*);" << endl;
513 ofs << "qREAL FQ_"<<ft_n<<"(const int, const int, const qREAL*, const qREAL*);" << endl;
514 ofs << "qREAL FMP_"<<ft_n<<"(const int, const int, const mpREAL*, const mpREAL*);" << endl;
515 ofs << "void X2ZL_"<<ft_n<<"(const dREAL*, dCOMPLEX*, dCOMPLEX*, dREAL*, const dREAL*, const dREAL*);" << endl;
516 ofs << "void X2ZQ_"<<ft_n<<"(const qREAL*, qCOMPLEX*, qCOMPLEX*, qREAL*, const qREAL*, const qREAL*);" << endl;
517 ofs << "void X2ZMP_"<<ft_n<<"(const mpREAL*, mpCOMPLEX*, mpCOMPLEX*, mpREAL*, const mpREAL*, const mpREAL*);" << endl;
518 ofs << "void MatL_"<<ft_n<<"(dCOMPLEX*, const dREAL*, const dREAL*, const dREAL*, const dREAL*);" << endl;
519 ofs << "void MatQ_"<<ft_n<<"(qCOMPLEX*, const qREAL*, const qREAL*, const qREAL*, const qREAL*);" << endl;
520 ofs << "void MatMP_"<<ft_n<<"(mpCOMPLEX*, const mpREAL*, const mpREAL*, const mpREAL*, const mpREAL*);" << endl;
521 ofs << endl << endl;
522 }
523
524 /*----------------------------------------------*/
525 // long double
526 /*----------------------------------------------*/
527 auto cppL = CppFormat(ofs, "L");
528 // alwasy export non-complex function
529 if(true) {
530 ofs << "extern \"C\" " << endl;
531 ofs << "int SDD_"<<idx<<"(const unsigned int xn, const dREAL x[], const unsigned int yn, dREAL y[], const dREAL pl[], const dREAL las[]) {" << endl;
532 if(Debug) {
533 auto tmp = expr.subs(FTX(w1,w2)==1).subs(cxRepl).subs(plRepl);
534 ofs << "//for debug, intg: " << endl << "//" << tmp << endl;
535 }
536
537 auto intg = expr.subs(FTX(w1,w2)==1);
538 intg = xyz_pow_simplify(intg);
539 bool hasF2 = intg.has(iEpsilon) || intg.has(I);
540
541 // WickRotation
542 hasF2 = hasF2 || intg.has(WRA(w));
543 if(intg.has(WRA(w))) {
544 exset wras;
545 find(intg, WRA(w), wras);
546 if(wras.size()!=1) throw Error("CIPrepares: Too many WRA(w).");
547 auto wra = (*(wras.begin())).op(0);
548 intg = intg.subs(WRA(w)==xwra);
549 ofs << "dREAL wra = "; EvalL(wra).print(cppL); ofs << ";" << endl;
550
551 exset pows_set;
552 find(intg, pow(w1,w2), pows_set);
553 find(intg, sqrt(w), pows_set);
554 lst pow_subs;
555 for(auto item : pows_set) {
556 if(item.has(xwra)) {
557 if(item.match(pow(w1,w2)) && !item.op(1).info(info_flags::integer)) {
558 pow_subs.append(item == exp(item.op(1)*log(item.op(0))));
559 } else if(item.match(sqrt(w))) {
560 pow_subs.append(item == exp(log(item.op(0))/2));
561 }
562 }
563 }
564 if(pow_subs.nops()>0) intg = exp_simplify(intg);
565
566 exset logs_set;
567 find(intg, log(w), logs_set);
568 lst logs;
569 for(auto item : logs_set) {
570 if(item.has(xwra)) logs.append(item.op(0));
571 }
572
573 if(logs.nops()>0) {
574 ofs << "int nlog = "<<logs.nops()<<";" << endl;
575 ofs << "dCOMPLEX CLog[nlog], LogZ[nlog][NRCLog+1];" << endl;
576 cseParser cse;
577 lst clogs = ex_to<lst>(cse.Parse(logs));
578
579 ofs << "for(int ti=0; ti<=NRCLog; ti++) {" << endl;
580 ofs << "dREAL xwra;" << endl;
581 ofs << "if(ti==0) xwra = wra/(25*NRCLog);" << endl;
582 ofs << "else xwra = ti*wra/NRCLog;" << endl;
583 ofs << "dCOMPLEX "<<cse.oc<<"[" << cse.on()+1 << "];" << endl;
584 for(auto kv : cse.os()) {
585 ofs <<cse.oc<< "["<<kv.first<<"] = ";
586 EvalL(kv.second.subs(cxRepl).subs(plRepl)).print(cppL);
587 ofs << ";" << endl;
588 }
589
590 exmap log_subs;
591 for(int i=0; i<clogs.nops(); i++) {
592 ofs << "LogZ["<<i<<"][ti] = ";
593 EvalL(clogs.op(i).subs(cxRepl).subs(plRepl)).print(cppL);
594 ofs << ";" << endl;
595 log_subs[log(logs.op(i))] = Symbol("CLog["+to_string(i)+"]");
596 }
597 ofs << "}" << endl;
598 ofs << "for(int li=0; li<nlog; li++) CLog[li] = RCLog(LogZ[li],NRCLog);" << endl;
599 intg = intg.subs(log_subs);
600 }
601 ofs << "dREAL xwra = wra;" << endl;
602 }
603
604 cseParser cse;
605 intg = cse.Parse(intg);
606 if(hasF2) ofs << "dCOMPLEX "<<cse.oc<<"[" << cse.on()+1 << "];" << endl;
607 else ofs << "dREAL "<<cse.oc<<"[" << cse.on()+1 << "];" << endl;
608 for(auto kv : cse.os()) {
609 ofs <<cse.oc<< "["<<kv.first<<"] = ";
610 EvalL(kv.second.subs(cxRepl).subs(plRepl)).print(cppL);
611 ofs << ";" << endl;
612 }
613
614 ofs << "dCOMPLEX yy = ";
615 EvalL(intg.subs(cxRepl).subs(plRepl)).print(cppL);
616 ofs << ";" << endl;
617
618 ofs << "y[0] = yy.real();" << endl;
619 ofs << "y[1] = yy.imag();" << endl;
620 ofs << "return 0;" << endl;
621 ofs << "}" << endl;
622 ofs << endl;
623 }
624
625 if(hasF) {
626 ofs << "extern \"C\" " << endl;
627 ofs << "int CSDD_"<<idx<<"(const unsigned int xn, const dREAL x[], const unsigned int yn, dREAL y[], const dREAL pl[], const dREAL las[]) {" << endl;
628 ofs << "dREAL x0[xn];" << endl;
629 ofs << "dCOMPLEX z[xn],zz[xn],r[xn];" << endl;
630 ofs << "dREAL dff[xn+1];" << endl;
631 ofs << "dCOMPLEX yy=0, ytmp, det;" << endl;
632 ofs << "int ii, nfxs="<<fxs.size()<<";" << endl;
633 ofs << "dCOMPLEX mat[nfxs*nfxs];" << endl;
634 for(auto &kv : ft_expr) {
635 ofs << "{" << endl;
636 lst xs0;
637 for(int ii=0; ii<fxs.size(); ii++) {
638 if(!kv.first.has(fxs[ii])) xs0.append(ii);
639 }
640 ofs << "for(int i=0; i<xn; i++) z[i] = x0[i] = x[i];" << endl;
641 for(auto x0i : xs0) ofs << "x0["<<x0i<<"]=0;" << endl;
642 ofs << "X2ZL_"<<ft_n<<"(x0,z,r,dff,pl,las);" << endl;
643 ofs << "MatL_"<<ft_n<<"(mat,x0,dff,pl,las);" << endl;
644 for(auto x0i : xs0) {
645 ofs << "ii = " << x0i << ";" << endl;
646 ofs << "z[ii] = x[ii]-x[ii]*(1.L-x[ii])*r[ii];" << endl;
647 ofs << "for(int j=0; j<nfxs;j++) mat[nfxs*ii+j] = 0;" << endl;
648 ofs << "for(int i=0; i<nfxs;i++) mat[nfxs*i+ii] = 0;" << endl;
649 ofs << "mat[ii*nfxs+ii] = 1.L-(1.L-2.L*x[ii])*r[ii];" << endl;
650 }
651 ofs << "det = MatDet(mat, nfxs);" << endl;
652
653 ex intg = kv.second;
654
655 if(true) { // RCLog
656 exset pows_set;
657 find(intg, pow(w1,w2), pows_set);
658 find(intg, sqrt(w), pows_set);
659 lst pow_subs;
660 for(auto item : pows_set) {
661 if(!item.op(0).match(x(w))) {
662 if(item.match(pow(w1,w2)) && !item.op(1).info(info_flags::integer)) {
663 pow_subs.append(item == exp(item.op(1)*log(item.op(0))));
664 } else if(item.match(sqrt(w))) {
665 pow_subs.append(item == exp(log(item.op(0))/2));
666 }
667 }
668 }
669 if(pow_subs.nops()>0) intg = exp_simplify(intg);
670
671 exset logs_set;
672 find(intg, log(w), logs_set);
673 if(logs_set.size()>0) {
674 lst logs;
675 for(auto item : logs_set) {
676 if(!item.op(0).match(x(w))) logs.append(item.op(0));
677 }
678 if(logs.nops()>0) {
679 ofs << "int nlog = "<<logs.nops()<<";" << endl;
680 ofs << "dCOMPLEX CLog[nlog], LogZ[nlog][NRCLog+1];" << endl;
681 cseParser cse;
682 lst clogs = ex_to<lst>(cse.Parse(logs));
683
684 ofs << "for(int ti=0; ti<=NRCLog; ti++) {" << endl;
685 ofs << "if(ti==0) { for(int i=0; i<xn; i++) zz[i] = complex<dREAL>(z[i].real(), z[i].imag()/(25*NRCLog)); }" << endl;
686 ofs << "else { for(int i=0; i<xn; i++) zz[i] = complex<dREAL>(z[i].real(), ti*z[i].imag()/NRCLog); }" << endl;
687 ofs << "dCOMPLEX "<<cse.oc<<"[" << cse.on()+1 << "];" << endl;
688 for(auto kv : cse.os()) {
689 ofs <<cse.oc<< "["<<kv.first<<"] = ";
690 EvalL(kv.second.subs(czzRepl).subs(plRepl)).print(cppL);
691 ofs << ";" << endl;
692 }
693
694 exmap log_subs;
695 for(int i=0; i<clogs.nops(); i++) {
696 ofs << "LogZ["<<i<<"][ti] = ";
697 EvalL(clogs.op(i).subs(czzRepl).subs(plRepl)).print(cppL);
698 ofs << ";" << endl;
699 log_subs[log(logs.op(i))] = Symbol("CLog["+to_string(i)+"]");
700 }
701 ofs << "}" << endl;
702 ofs << "for(int li=0; li<nlog; li++) CLog[li] = RCLog(LogZ[li],NRCLog);" << endl;
703 intg = intg.subs(log_subs);
704 }
705 }
706 }
707
708 cseParser cse;
709 intg = cse.Parse(intg);
710 ofs << "dCOMPLEX "<<cse.oc<<"[" << cse.on()+1 << "];" << endl;
711 for(auto kv : cse.os()) {
712 ofs <<cse.oc<< "["<<kv.first<<"] = ";
713 EvalL(kv.second.subs(czRepl).subs(plRepl)).print(cppL);
714 ofs << ";" << endl;
715 }
716
717 ofs << "ytmp = ";
718 EvalL(intg.subs(czRepl).subs(plRepl)).print(cppL);
719 ofs << ";" << endl;
720 ofs << "yy += det * ytmp;" << endl << endl;
721 ofs << "}" << endl;
722 }
723
724 ofs << "y[0] = yy.real();" << endl;
725 ofs << "y[1] = yy.imag();" << endl;
726 ofs << "return 0;" << endl;
727 ofs << "}" << endl;
728 ofs << endl;
729 }
730
731
732
733 /*----------------------------------------------*/
734 // Quadruple
735 /*----------------------------------------------*/
736 auto cppQ = CppFormat(ofs, "Q");
737
738 // always export non-complex function
739 if(true) {
740 ofs << "extern \"C\" " << endl;
741 ofs << "int SDQ_"<<idx<<"(const unsigned int xn, const qREAL x[], const int unsigned yn, qREAL y[], const qREAL pl[], const qREAL las[]) {" << endl;
742
743 auto intg = expr.subs(FTX(w1,w2)==1);
744 intg = xyz_pow_simplify(intg);
745 bool hasF2 = intg.has(iEpsilon) || intg.has(I);
746
747 // WickRotation
748 hasF2 = hasF2 || intg.has(WRA(w));
749 if(intg.has(WRA(w))) {
750 exset wras;
751 find(intg, WRA(w), wras);
752 if(wras.size()!=1) throw Error("CIPrepares: Too many WRA(w).");
753 auto wra = (*(wras.begin())).op(0);
754 intg = intg.subs(WRA(w)==xwra);
755 ofs << "qREAL wra = "; EvalQ(wra).print(cppQ); ofs << ";" << endl;
756
757 exset pows_set;
758 find(intg, pow(w1,w2), pows_set);
759 find(intg, sqrt(w), pows_set);
760 lst pow_subs;
761 for(auto item : pows_set) {
762 if(item.has(xwra)) {
763 if(item.match(pow(w1,w2)) && !item.op(1).info(info_flags::integer)) {
764 pow_subs.append(item == exp(item.op(1)*log(item.op(0))));
765 } else if(item.match(sqrt(w))) {
766 pow_subs.append(item == exp(log(item.op(0))/2));
767 }
768 }
769 }
770 if(pow_subs.nops()>0) intg = exp_simplify(intg);
771
772 exset logs_set;
773 find(intg, log(w), logs_set);
774 lst logs;
775 for(auto item : logs_set) {
776 if(item.has(xwra)) logs.append(item.op(0));
777 }
778
779 if(logs.nops()>0) {
780 ofs << "int nlog = "<<logs.nops()<<";" << endl;
781 ofs << "qCOMPLEX CLog[nlog], LogZ[nlog][NRCLog+1];" << endl;
782 cseParser cse;
783 lst clogs = ex_to<lst>(cse.Parse(logs));
784
785 ofs << "for(int ti=0; ti<=NRCLog; ti++) {" << endl;
786 ofs << "qREAL xwra;" << endl;
787 ofs << "if(ti==0) xwra = wra/(25*NRCLog);" << endl;
788 ofs << "else xwra = ti*wra/NRCLog;" << endl;
789 ofs << "qCOMPLEX "<<cse.oc<<"[" << cse.on()+1 << "];" << endl;
790 for(auto kv : cse.os()) {
791 ofs <<cse.oc<< "["<<kv.first<<"] = ";
792 EvalQ(kv.second.subs(cxRepl).subs(plRepl)).print(cppQ);
793 ofs << ";" << endl;
794 }
795
796 exmap log_subs;
797 for(int i=0; i<clogs.nops(); i++) {
798 ofs << "LogZ["<<i<<"][ti] = ";
799 EvalQ(clogs.op(i).subs(cxRepl).subs(plRepl)).print(cppQ);
800 ofs << ";" << endl;
801 log_subs[log(logs.op(i))] = Symbol("CLog["+to_string(i)+"]");
802 }
803 ofs << "}" << endl;
804 ofs << "for(int li=0; li<nlog; li++) CLog[li] = RCLog(LogZ[li],NRCLog);" << endl;
805 intg = intg.subs(log_subs);
806 }
807 ofs << "qREAL xwra = wra;" << endl;
808 }
809
810 cseParser cse;
811 intg = cse.Parse(intg);
812 if(hasF2) ofs << "qCOMPLEX "<<cse.oc<<"[" << cse.on()+1 << "];" << endl;
813 else ofs << "qREAL "<<cse.oc<<"[" << cse.on()+1 << "];" << endl;
814 for(auto kv : cse.os()) {
815 ofs <<cse.oc<< "["<<kv.first<<"] = ";
816 EvalQ(kv.second.subs(cxRepl).subs(plRepl)).print(cppQ);
817 ofs << ";" << endl;
818 }
819
820 ofs << "qCOMPLEX yy = ";
821 EvalQ(intg.subs(cxRepl).subs(plRepl)).print(cppQ);
822 ofs << ";" << endl;
823 ofs << "y[0] = crealq(yy);" << endl;
824 ofs << "y[1] = cimagq(yy);" << endl;
825 ofs << "return 0;" << endl;
826 ofs << "}" << endl;
827 ofs << endl;
828 }
829
830 if(hasF) {
831 ofs << "extern \"C\" " << endl;
832 ofs << "int CSDQ_"<<idx<<"(const unsigned int xn, const qREAL x[], const int unsigned yn, qREAL y[], const qREAL pl[], const qREAL las[]) {" << endl;
833 ofs << "qREAL x0[xn];" << endl;
834 ofs << "qCOMPLEX z[xn],zz[xn],r[xn];" << endl;
835 ofs << "qREAL dff[xn+1];" << endl;
836 ofs << "qCOMPLEX yy=0, ytmp, det;" << endl;
837 ofs << "int ii, nfxs="<<fxs.size()<<";" << endl;
838 ofs << "qCOMPLEX mat[nfxs*nfxs];" << endl;
839 for(auto &kv : ft_expr) {
840 ofs << "{" << endl;
841 lst xs0;
842 for(int ii=0; ii<fxs.size(); ii++) {
843 if(!kv.first.has(fxs[ii])) xs0.append(ii);
844 }
845 ofs << "for(int i=0; i<xn; i++) z[i] = x0[i] = x[i];" << endl;
846 for(auto x0i : xs0) ofs << "x0["<<x0i<<"]=0;" << endl;
847 ofs << "X2ZQ_"<<ft_n<<"(x0,z,r,dff,pl,las);" << endl;
848 ofs << "MatQ_"<<ft_n<<"(mat,x0,dff,pl,las);" << endl;
849 for(auto x0i : xs0) {
850 ofs << "ii = " << x0i << ";" << endl;
851 ofs << "z[ii] = x[ii]-x[ii]*(1.Q-x[ii])*r[ii];" << endl;
852 ofs << "for(int j=0; j<nfxs;j++) mat[nfxs*ii+j] = 0;" << endl;
853 ofs << "for(int i=0; i<nfxs;i++) mat[nfxs*i+ii] = 0;" << endl;
854 ofs << "mat[ii*nfxs+ii] = 1.Q-(1.Q-2.Q*x[ii])*r[ii];" << endl;
855 }
856 ofs << "det = MatDet(mat, nfxs);" << endl;
857
858 ex intg = kv.second;
859
860 if(true) { // RCLog
861 exset pows_set;
862 find(intg, pow(w1,w2), pows_set);
863 find(intg, sqrt(w), pows_set);
864 lst pow_subs;
865 for(auto item : pows_set) {
866 if(!item.op(0).match(x(w))) {
867 if(item.match(pow(w1,w2)) && !item.op(1).info(info_flags::integer)) {
868 pow_subs.append(item == exp(item.op(1)*log(item.op(0))));
869 } else if(item.match(sqrt(w))) {
870 pow_subs.append(item == exp(log(item.op(0))/2));
871 }
872 }
873 }
874 if(pow_subs.nops()>0) intg = exp_simplify(intg);
875
876 exset logs_set;
877 find(intg, log(w), logs_set);
878 if(logs_set.size()>0) {
879 lst logs;
880 for(auto item : logs_set) {
881 if(!item.op(0).match(x(w))) logs.append(item.op(0));
882 }
883 if(logs.nops()>0) {
884 ofs << "int nlog = "<<logs.nops()<<";" << endl;
885 ofs << "qCOMPLEX CLog[nlog], LogZ[nlog][NRCLog+1];" << endl;
886 cseParser cse;
887 lst clogs = ex_to<lst>(cse.Parse(logs));
888
889 ofs << "for(int ti=0; ti<=NRCLog; ti++) {" << endl;
890 ofs << "if(ti==0) { for(int i=0; i<xn; i++) zz[i] = crealq(z[i]) + 1.Qi*cimagq(z[i])/(25.Q*NRCLog); }" << endl;
891 ofs << "else { for(int i=0; i<xn; i++) zz[i] = crealq(z[i]) + 1.Qi*ti*cimagq(z[i])/NRCLog; }" << endl;
892 ofs << "qCOMPLEX "<<cse.oc<<"[" << cse.on()+1 << "];" << endl;
893 for(auto kv : cse.os()) {
894 ofs <<cse.oc<< "["<<kv.first<<"] = ";
895 EvalQ(kv.second.subs(czzRepl).subs(plRepl)).print(cppQ);
896 ofs << ";" << endl;
897 }
898
899 exmap log_subs;
900 for(int i=0; i<clogs.nops(); i++) {
901 ofs << "LogZ["<<i<<"][ti] = ";
902 EvalQ(clogs.op(i).subs(czzRepl).subs(plRepl)).print(cppQ);
903 ofs << ";" << endl;
904 log_subs[log(logs.op(i))] = Symbol("CLog["+to_string(i)+"]");
905 }
906 ofs << "}" << endl;
907 ofs << "for(int li=0; li<nlog; li++) CLog[li] = RCLog(LogZ[li],NRCLog);" << endl;
908 intg = intg.subs(log_subs);
909 }
910 }
911 }
912
913 cseParser cse;
914 intg = cse.Parse(intg);
915 ofs << "qCOMPLEX "<<cse.oc<<"[" << cse.on()+1 << "];" << endl;
916 for(auto kv : cse.os()) {
917 ofs <<cse.oc<< "["<<kv.first<<"] = ";
918 EvalQ(kv.second.subs(czRepl).subs(plRepl)).print(cppQ);
919 ofs << ";" << endl;
920 }
921
922 ofs << "ytmp = ";
923 EvalQ(intg.subs(czRepl).subs(plRepl)).print(cppQ);
924 ofs << ";" << endl;
925 ofs << "yy += det * ytmp;" << endl << endl;
926 ofs << "}" << endl;
927 }
928 ofs << "y[0] = crealq(yy);" << endl;
929 ofs << "y[1] = cimagq(yy);" << endl;
930 ofs << "return 0;" << endl;
931 ofs << "}" << endl;
932 ofs << endl;
933
934 // Export the F-term, only Quadruple-type
935 ofs << "extern \"C\" " << endl;
936 ofs << "qREAL FT_"<<idx<<"(const qREAL x[], const qREAL pl[]) {" << endl;
937 ofs << "qREAL yy = FQ_" << ft_n << "(x, pl);" << endl;
938 ofs << "return yy;" << endl;
939 ofs << "}" << endl;
940 ofs << endl;
941 }
942
943
944 /*----------------------------------------------*/
945 // Multiple Precision
946 /*----------------------------------------------*/
947 auto cppMP = CppFormat(ofs, "MP");
948
949 // always export non-complex function
950 if(true) {
951 ofs << "extern \"C\" " << endl;
952 ofs << "int SDMP_"<<idx<<"(const unsigned int xn, const mpREAL x[], const unsigned int yn, mpREAL y[], const mpREAL pl[], const mpREAL las[]) {" << endl;
953 auto intg = expr.subs(FTX(w1,w2)==1);
954 intg = xyz_pow_simplify(intg);
955 bool hasF2 = intg.has(iEpsilon) || intg.has(I);
956
957 // WickRotation
958 hasF2 = hasF2 || intg.has(WRA(w));
959 if(intg.has(WRA(w))) {
960 exset wras;
961 find(intg, WRA(w), wras);
962 if(wras.size()!=1) throw Error("CIPrepares: Too many WRA(w).");
963 auto wra = (*(wras.begin())).op(0);
964 intg = intg.subs(WRA(w)==xwra);
965 ofs << "mpREAL wra = "; EvalMP(wra).print(cppMP); ofs << ";" << endl;
966
967 exset pows_set;
968 find(intg, pow(w1,w2), pows_set);
969 find(intg, sqrt(w), pows_set);
970 lst pow_subs;
971 for(auto item : pows_set) {
972 if(item.has(xwra)) {
973 if(item.match(pow(w1,w2)) && !item.op(1).info(info_flags::integer)) {
974 pow_subs.append(item == exp(item.op(1)*log(item.op(0))));
975 } else if(item.match(sqrt(w))) {
976 pow_subs.append(item == exp(log(item.op(0))/2));
977 }
978 }
979 }
980 if(pow_subs.nops()>0) intg = exp_simplify(intg);
981
982 exset logs_set;
983 find(intg, log(w), logs_set);
984 lst logs;
985 for(auto item : logs_set) {
986 if(item.has(xwra)) logs.append(item.op(0));
987 }
988
989 if(logs.nops()>0) {
990 ofs << "int nlog = "<<logs.nops()<<";" << endl;
991 ofs << "mpCOMPLEX CLog[nlog], LogZ[nlog][NRCLog+1];" << endl;
992 cseParser cse;
993 lst clogs = ex_to<lst>(cse.Parse(logs));
994
995 ofs << "for(int ti=0; ti<=NRCLog; ti++) {" << endl;
996 ofs << "mpREAL xwra;" << endl;
997 ofs << "if(ti==0) xwra = wra/(25*NRCLog);" << endl;
998 ofs << "else xwra = ti*wra/NRCLog;" << endl;
999 ofs << "mpCOMPLEX "<<cse.oc<<"[" << cse.on()+1 << "];" << endl;
1000 for(auto kv : cse.os()) {
1001 ofs <<cse.oc<< "["<<kv.first<<"] = ";
1002 EvalMP(kv.second.subs(cxRepl).subs(plRepl)).print(cppMP);
1003 ofs << ";" << endl;
1004 }
1005
1006 exmap log_subs;
1007 for(int i=0; i<clogs.nops(); i++) {
1008 ofs << "LogZ["<<i<<"][ti] = ";
1009 EvalMP(clogs.op(i).subs(cxRepl).subs(plRepl)).print(cppMP);
1010 ofs << ";" << endl;
1011 log_subs[log(logs.op(i))] = Symbol("CLog["+to_string(i)+"]");
1012 }
1013 ofs << "}" << endl;
1014 ofs << "for(int li=0; li<nlog; li++) CLog[li] = RCLog(LogZ[li],NRCLog);" << endl;
1015 intg = intg.subs(log_subs);
1016 }
1017 ofs << "mpREAL xwra = wra;" << endl;
1018 }
1019
1020 cseParser cse;
1021 intg = cse.Parse(intg);
1022 if(hasF2) ofs << "mpCOMPLEX "<<cse.oc<<"[" << cse.on()+1 << "];" << endl;
1023 else ofs << "mpREAL "<<cse.oc<<"[" << cse.on()+1 << "];" << endl;
1024 for(auto kv : cse.os()) {
1025 ofs <<cse.oc<< "["<<kv.first<<"] = ";
1026 EvalMP(kv.second.subs(cxRepl).subs(plRepl)).print(cppMP);
1027 ofs << ";" << endl;
1028 }
1029
1030 ofs << "mpCOMPLEX yy = ";
1031 EvalMP(intg.subs(cxRepl).subs(plRepl)).print(cppMP);
1032 ofs << ";" << endl;
1033 ofs << "y[0] = yy.real();" << endl;
1034 ofs << "y[1] = yy.imag();" << endl;
1035 ofs << "return 0;" << endl;
1036 ofs << "}" << endl;
1037 ofs << endl;
1038 }
1039
1040 if(hasF) {
1041 ofs << "extern \"C\" " << endl;
1042 ofs << "int CSDMP_"<<idx<<"(const unsigned int xn, const mpREAL x[], const unsigned int yn, mpREAL y[], const mpREAL pl[], const mpREAL las[]) {" << endl;
1043 ofs << "mpREAL x0[xn];" << endl;
1044 ofs << "mpCOMPLEX z[xn],zz[xn],r[xn];" << endl;
1045 ofs << "mpREAL dff[xn+1];" << endl;
1046 ofs << "mpCOMPLEX yy=mpREAL(0), ytmp, det;" << endl;
1047 ofs << "int ii, nfxs="<<fxs.size()<<";" << endl;
1048 ofs << "mpCOMPLEX mat[nfxs*nfxs];" << endl;
1049 for(auto &kv : ft_expr) {
1050 ofs << "{" << endl;
1051 lst xs0;
1052 for(int ii=0; ii<fxs.size(); ii++) {
1053 if(!kv.first.has(fxs[ii])) xs0.append(ii);
1054 }
1055 ofs << "for(int i=0; i<xn; i++) z[i] = x0[i] = x[i];" << endl;
1056 for(auto x0i : xs0) ofs << "x0["<<x0i<<"]=0;" << endl;
1057 ofs << "X2ZMP_"<<ft_n<<"(x0,z,r,dff,pl,las);" << endl;
1058 ofs << "MatMP_"<<ft_n<<"(mat,x0,dff,pl,las);" << endl;
1059 for(auto x0i : xs0) {
1060 ofs << "ii = " << x0i << ";" << endl;
1061 ofs << "z[ii] = x[ii]-x[ii]*(1-x[ii])*r[ii];" << endl;
1062 ofs << "for(int j=0; j<nfxs;j++) mat[nfxs*ii+j] = mpREAL(0);" << endl;
1063 ofs << "for(int i=0; i<nfxs;i++) mat[nfxs*i+ii] = mpREAL(0);" << endl;
1064 ofs << "mat[ii*nfxs+ii] = mpREAL(1)-(1-2*x[ii])*r[ii];" << endl;
1065 }
1066 ofs << "det = MatDet(mat, nfxs);" << endl;
1067
1068 ex intg = kv.second;
1069
1070 if(true) { // RCLog
1071 exset pows_set;
1072 find(intg, pow(w1,w2), pows_set);
1073 find(intg, sqrt(w), pows_set);
1074 lst pow_subs;
1075 for(auto item : pows_set) {
1076 if(!item.op(0).match(x(w))) {
1077 if(item.match(pow(w1,w2)) && !item.op(1).info(info_flags::integer)) {
1078 pow_subs.append(item == exp(item.op(1)*log(item.op(0))));
1079 } else if(item.match(sqrt(w))) {
1080 pow_subs.append(item == exp(log(item.op(0))/2));
1081 }
1082 }
1083 }
1084 if(pow_subs.nops()>0) intg = exp_simplify(intg);
1085
1086 exset logs_set;
1087 find(intg, log(w), logs_set);
1088 if(logs_set.size()>0) {
1089 lst logs;
1090 for(auto item : logs_set) {
1091 if(!item.op(0).match(x(w))) logs.append(item.op(0));
1092 }
1093 if(logs.nops()>0) {
1094 ofs << "int nlog = "<<logs.nops()<<";" << endl;
1095 ofs << "mpCOMPLEX CLog[nlog], LogZ[nlog][NRCLog+1];" << endl;
1096 cseParser cse;
1097 lst clogs = ex_to<lst>(cse.Parse(logs));
1098
1099 ofs << "for(int ti=0; ti<=NRCLog; ti++) {" << endl;
1100 ofs << "if(ti==0) { for(int i=0; i<xn; i++) zz[i] = complex<mpREAL>(z[i].real(), z[i].imag()/(25*NRCLog)); }" << endl;
1101 ofs << "else { for(int i=0; i<xn; i++) zz[i] = complex<mpREAL>(z[i].real(), ti*z[i].imag()/NRCLog); }" << endl;
1102 ofs << "mpCOMPLEX "<<cse.oc<<"[" << cse.on()+1 << "];" << endl;
1103 for(auto kv : cse.os()) {
1104 ofs <<cse.oc<< "["<<kv.first<<"] = ";
1105 EvalMP(kv.second.subs(czzRepl).subs(plRepl)).print(cppMP);
1106 ofs << ";" << endl;
1107 }
1108
1109 exmap log_subs;
1110 for(int i=0; i<clogs.nops(); i++) {
1111 ofs << "LogZ["<<i<<"][ti] = ";
1112 EvalMP(clogs.op(i).subs(czzRepl).subs(plRepl)).print(cppMP);
1113 ofs << ";" << endl;
1114 log_subs[log(logs.op(i))] = Symbol("CLog["+to_string(i)+"]");
1115 }
1116 ofs << "}" << endl;
1117 ofs << "for(int li=0; li<nlog; li++) CLog[li] = RCLog(LogZ[li],NRCLog);" << endl;
1118 intg = intg.subs(log_subs);
1119 }
1120 }
1121 }
1122
1123 cseParser cse;
1124 intg = cse.Parse(intg);
1125 ofs << "mpCOMPLEX "<<cse.oc<<"[" << cse.on()+1 << "];" << endl;
1126 for(auto kv : cse.os()) {
1127 ofs <<cse.oc<< "["<<kv.first<<"] = ";
1128 EvalMP(kv.second.subs(czRepl).subs(plRepl)).print(cppMP);
1129 ofs << ";" << endl;
1130 }
1131
1132 ofs << "ytmp = ";
1133 EvalMP(intg.subs(czRepl).subs(plRepl)).print(cppMP);
1134 ofs << ";" << endl;
1135 ofs << "yy += det * ytmp;" << endl << endl;
1136 ofs << "}" << endl;
1137 }
1138
1139 ofs << "y[0] = yy.real();" << endl;
1140 ofs << "y[1] = yy.imag();" << endl;
1141 ofs << "return 0;" << endl;
1142 ofs << "}" << endl;
1143 ofs << endl;
1144 }
1145
1146 ofs.close();
1147 ostringstream ofn, cmd;
1148 ofn << spid << "/" << idx << ".o";
1149 cmd << cpp << " -pipe -fPIC " << INC_FLAGS << " -include NFunctions.h -c -o " << ofn.str() << " " << cppfn.str();
1150 auto rc = system(cmd.str().c_str());
1151 if(!Debug) remove(cppfn.str().c_str());
1152 return lst{ idx, xs.size(), kvf.op(0), ft_n };
1153 }, "FCI-I");
1154
1155// if(true) { // try other complilation method
1156// ostringstream cmd;
1157// cmd << "cd " << spid << ";find . -name '*.cpp' | xargs -n 100 -P " << CpuCores() << " " << cpp << " -pipe -fPIC " << INC_FLAGS << " -include NFunctions.h -c";
1158// auto rc = system(cmd.str().c_str());
1159// //if(!Debug) remove(cppfn.str().c_str());
1160// }
1161
1162 //============================================================================================================
1163
1164 ostringstream fsofn, sofn, garfn, cmd;
1165 if(key != "") {
1166 fsofn << key << "F.so";
1167 sofn << key << ".so";
1168 garfn << key << ".ci.gar";
1169 lst gar_res;
1170 for(auto &item : res) gar_res.append(item);
1171 archive ar;
1172 ar.archive_ex(gar_res, "res");
1173 ar.archive_ex(19790923, "c");
1174 ar.archive_ex(FT_N_XN, "ftnxn");
1175 ar.archive_ex(soLimit,"soLimit");
1176 ar.archive_ex(eps_lst,"eps_lst");
1177 ofstream out(garfn.str());
1178 out << ar;
1179 out.close();
1180 cmd.clear();
1181 cmd.str("");
1182 cmd << "rm -f " << key << "X*.so";
1183 rc = system(cmd.str().c_str());
1184 } else {
1185 fsofn << spid << "F.so";
1186 sofn << spid << ".so";
1187 for(auto &item : res) ciResult.push_back(ex_to<lst>(item));
1188 cmd.clear();
1189 cmd.str("");
1190 cmd << "rm -f " << spid << "X*.so";
1191 rc = system(cmd.str().c_str());
1192 }
1193
1194 int res_size = res.size();
1195 if(soLimit<100) soLimit = 100;
1196 if(res_size>soLimit) {
1197 cmd.clear();
1198 cmd.str("");
1199 cmd << cpp << " " << LIB_FLAGS << " -Wl,-rpath,. -rdynamic -fPIC -shared -lHepLib -lquadmath -lmpfr -lgmp";
1200 if(hasF) cmd << " " << fsofn.str();
1201 cmd << " -o " << sofn.str() << " $(seq -f '" << spid << "/%g.o' 0 " << (soLimit-1) << ")";
1202 cmd << " -lHepLib -lquadmath -lmpfr -lgmp";
1203 cmd << " 1> /dev/null 2> /dev/null";
1204 rc = system(cmd.str().c_str());
1205
1206 for(int n=1; true; n++) {
1207 int start = n*soLimit;
1208 int end = (n+1)*soLimit-1;
1209 if(end>res_size-1) end = res_size-1;
1210 sofn.clear();
1211 sofn.str("");
1212 if(key != "") sofn << key << "X" << n << ".so";
1213 else sofn << spid << "X" << n << ".so";
1214 cmd.clear();
1215 cmd.str("");
1216 cmd << cpp << " " << LIB_FLAGS << " -Wl,-rpath,. -rdynamic -fPIC -shared -lHepLib -lquadmath -lmpfr -lgmp";
1217 if(hasF) cmd << " " << fsofn.str();
1218 cmd << " -o " << sofn.str() << " $(seq -f '" << spid << "/%g.o' " << start << " " << end << ")";
1219 if(hasF) cmd << " " << fsofn.str();
1220 cmd << " -lHepLib -lquadmath -lmpfr -lgmp";
1221 cmd << " 1> /dev/null 2> /dev/null";
1222 rc = system(cmd.str().c_str());
1223 if(end>=res_size-1) break;
1224 }
1225 } else {
1226 cmd.clear();
1227 cmd.str("");
1228 cmd << cpp << " " << LIB_FLAGS << " -Wl,-rpath,. -rdynamic -fPIC -shared -lHepLib -lquadmath -lmpfr -lgmp";
1229 if(hasF) cmd << " " << fsofn.str();
1230 cmd << " -o " << sofn.str() << " " << spid << "/*.o";
1231 if(hasF) cmd << " " << fsofn.str();
1232 cmd << " -lHepLib -lquadmath -lmpfr -lgmp";
1233 cmd << " 1> /dev/null 2> /dev/null";
1234 rc = system(cmd.str().c_str());
1235 }
1236 cmd.clear();
1237 cmd.str("");
1238 if(Debug) cmd << "rm -rf " << key << "_debug;mv -f " << spid << " " << key << "_debug;rm -f " << key << "_debug/*.o";
1239 else cmd << "rm -rf " << spid;
1240 rc = system(cmd.str().c_str());
1241
1242 }
1243
1244
1245}
RESET
#define RESET
Definition BASIC.h:79
intg
lst intg
Definition Others.cpp:26
SD.h
SecDec header file.
HepLib::Error
class used to wrap error message
Definition BASIC.h:242
HepLib::SD::SecDec::expResult
exvector expResult
Definition SD.h:427
HepLib::SD::SecDec::nReplacement
exmap nReplacement
Definition SD.h:424
HepLib::SD::SecDec::CIPrepares
void CIPrepares(const string &key="")
Prepare for the Contours and Integrates calls .so will be generated, for more detailed exported funct...
Definition CIPrepares.cpp:17
HepLib::SD
namespace for Numerical integration with Sector Decomposition method
Definition AsyMB.cpp:10
HepLib::SD::get_xy_from
exvector get_xy_from(ex pol)
Definition Helpers.cpp:14
HepLib::SD::xyz_pow_simplify
ex xyz_pow_simplify(const ex expr_in)
Definition Helpers.cpp:349
HepLib::SD::exp_simplify
ex exp_simplify(const ex expr_in)
Definition Helpers.cpp:315
HepLib::SD::get_pl_from
exvector get_pl_from(ex pol)
Definition Helpers.cpp:74
HepLib::Color_HighLight
const char * Color_HighLight
Definition BASIC.cpp:248
HepLib::INC_FLAGS
string INC_FLAGS
Definition Init.cpp:165
HepLib::iEpsilon
const iSymbol iEpsilon
HepLib::EvalL
ex EvalL(ex expr)
Definition BASIC.cpp:1258
HepLib::GiNaC_Parallel_RM
map< string, bool > GiNaC_Parallel_RM
Definition Init.cpp:151
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:1204
HepLib::EvalQ
ex EvalQ(ex expr)
Definition BASIC.cpp:1263
HepLib::file_exists
bool file_exists(string fn)
Definition BASIC.h:289
HepLib::LIB_FLAGS
string LIB_FLAGS
Definition Init.cpp:166
HepLib::now
string now(bool use_date)
date/time string
Definition BASIC.cpp:527
HepLib::Debug
bool Debug
Definition Init.cpp:144
HepLib::w
ex w
Definition Init.cpp:93
HepLib::EvalMP
ex EvalMP(ex expr)
Definition BASIC.cpp:1268
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:1224
HepLib::diff_ex
ex diff_ex(ex const expr, ex const xp, unsigned nth, bool expand)
the differential like Mathematica
Definition BASIC.cpp:1065
HepLib::Verbose
int Verbose
Definition Init.cpp:142
HepLib::ex2str
string ex2str(const ex &expr)
convert ex to output string, the defalut printer format will be used
Definition BASIC.cpp:717
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::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::get_env
string get_env(const string &name)
Definition BASIC.cpp:2371
HepLib::subs
ex subs(const ex &e, init_list sl)
Definition BASIC.h:51
Generated by doxygen 1.9.8