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