HepLib
CIPrepares.cpp
Go to the documentation of this file.
1 
6 #include "SD.h"
7 #include <math.h>
8 #include <cmath>
9 
10 namespace 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::IsZero
bool IsZero
Definition: SD.h:433
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:160
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:148
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:161
HepLib::now
string now(bool use_date)
date/time string
Definition: BASIC.cpp:525
HepLib::Debug
bool Debug
Definition: Init.cpp:141
HepLib::w
ex w
Definition: Init.cpp:90
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:139
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.1