HepLib
Form.cpp
Go to the documentation of this file.
1 
6 #include "HEP.h"
7 #include "ginac/parse_context.h"
8 
9 namespace HepLib {
10 
11  //-----------------------------------------------------------
12  // Extend Parser for form, copied from ginac/parser of GiNaC
13  //-----------------------------------------------------------
14  namespace {
15 
16  alignas(2) static ex SP_reader(const exvector& ev) {
17  return SP(ev[0], ev[1]).subs(SP_map);
18  }
19 
20  alignas(2) static ex LC_reader(const exvector& ev) {
21  if(ev.size()==4) return (-I)*LC(ev[0], ev[1], ev[2], ev[3]);
22  else return Eps(ev);
23  }
24 
25  alignas(2) static ex SUNT_reader(const exvector& ev) {
26  int n = ev.size();
27  if(n==3) return SUNT(ev[0], ev[1], ev[2]);
28  else if(n>3) {
29  lst as;
30  for(int i=0; i<n-2; i++) as.append(ev[i]);
31  return SUNT(as,ev[n-2],ev[n-1]);
32  } else throw Error("SUNT_reader: number of arguments less than 3.");
33  }
34 
35  alignas(2) static ex TTR_reader(const exvector& ev) {
36  if(ev.size()==1) return TTR(ev[0]);
37  lst as;
38  for(auto item : ev) as.append(item);
39  return TTR(as);
40  }
41 
42  alignas(2) static ex DGamma_reader(const exvector& ev) {
43  int n = ev.size();
44  if(n==1) return GAS(ev[0]);
45  ex ret = 1;
46  int rl = ex2int(ev[0]);
47  for(int i=1; i<n; i++) ret = ret * GAS(ev[i],rl);
48  return ret;
49  }
50 
51  alignas(2) static ex gi_reader(const exvector& ev) {
52  int n = ev.size();
53  if(n==0) return GAS(1);
54  else if(n==1) return GAS(1,ex2int(ev[0]));
55  else throw Error("DGamma_reader: number of arguments more than 2.");
56  }
57 
58  alignas(2) static ex Matrix_reader(const exvector& ev) {
59  return Matrix(ev[0],ev[1],ev[2]);
60  }
61  }
62 
63  //-----------------------------------------------------------
64  // Run Form Program
65  //-----------------------------------------------------------
66  namespace {
67 
68  struct mapGamma : public map_function {
69  public:
70  ex operator()(const ex &e) {
71  if(is_a<DGamma>(e)) return DGamma(ex_to<DGamma>(e), gline);
72  else if(!DGamma::has(e)) return e;
73  else return e.map(*this);
74  }
75  mapGamma(int _gline) : gline(_gline) { };
76  private:
77  unsigned gline = 0;
78  };
79 
80  struct mapTR : public map_function {
81  public:
82  ex operator()(const ex &e) {
83  if (e.match(TR(w))) {
84  ex trs = e.op(0);
85  gline++;
86  trs = mapGamma(gline)(trs);
87  trs = DGamma(1, gline) * trs;
88  if(glmax<gline) {
89  glmax = gline;
90  if(glmax>128) throw Error("too large index with glmax>128.");
91  }
92  return TR(trs);
93  } else if(is_a<add>(e)) {
94  ex res = 0;
95  unsigned gl = gline;
96  unsigned glmax = gline;
97  for(auto item : e) {
98  gline = gl;
99  res += (*this)(item);
100  if(glmax<gline) glmax = gline;
101  }
102  gline = glmax;
103  return res;
104  } else return e.map(*this);
105  }
106  unsigned glmax = 0;
107  private:
108  unsigned gline = 0;
109  };
110 
111  string init_script = R"EOF(
112 CFunction pow,sqrt,gamma,HF,Matrix,WF;
113 Tensor TTR(cyclic), f(antisymmetric), T, f4, colTp;
114 Symbols reX,I2R,NF,NA,d,I,Pi;
115 AutoDeclare Symbols gCF, trcN;
116 Dimension NA;
117 AutoDeclare Index colA;
118 Dimension NF;
119 AutoDeclare Index colF;
120 
121 #procedure SUNTrace
122 Dimension NF;
123 
124 repeat;
125  id,once,TTR(?a) = T(?a,colF1,colF1);
126  sum colF1;
127  repeat;
128  id,once,T(colA1?,colA2?,?a,colF1?,colF2?) = T(colA1,colF1,colF3)*T(colA2,?a,colF3,colF2);
129  sum colF3;
130  id,once,f4(colA1?,colA2?,colA3?,colA4?) = f(colA1,colA2,colA5) * f(colA5,colA3,colA4);
131  sum colA5;
132  endrepeat;
133 endrepeat;
134 
135 #do colXi = 1,1
136  if ( count(f,1) || match(T(colA1?,colF1?,colF2?)*T(colA1?,colF3?,colF4?)) ) redefine colXi "0";
137  id,once,f(colA1?,colA2?,colA3?) = 1/I2R/i_*T(colA1,colF1,colF2)*T(colA2,colF2,colF3)*T(colA3,colF3,colF1)-1/I2R/i_*T(colA3,colF1,colF2)*T(colA2,colF2,colF3)*T(colA1,colF3,colF1);
138  sum colF1,colF2,colF3;
139  id T(colA1?,colF1?,colF2?)*T(colA1?,colF3?,colF4?) = colTp(colF1,colF2,colF3,colF4);
140  #do colXj = 1,1
141  if ( count(colTp,1) ) redefine colXj "0";
142  .sort
143  id,once,colTp(colF1?,colF2?,colF3?,colF4?) = I2R*(d_(colF1,colF4)*d_(colF2,colF3)-d_(colF1,colF2)*d_(colF3,colF4)/NF);
144  #enddo
145 #enddo
146 
147 repeat;
148  id T(colA1?,?a,colF1?,colF2?)*T(colA2?,?b,colF2?,colF3?) = T(colA1,?a,colA2,?b,colF1,colF3);
149 endrepeat;
150 id T(?a,colF1?,colF1?) = TTR(?a);
151 id TTR(colA1?) = 0;
152 id TTR(colA1?,colA2?) = I2R*d_(colA1,colA2);
153 .sort
154 
155 #endprocedure
156 .global
157  )EOF";
158  }
159 
160  //-----------------------------------------------------------
161  // form
162  //-----------------------------------------------------------
163  namespace {
164  ex runform(const ex &expr_in, int verb) {
165  static map<pid_t, Form> form_map;
166  auto pid = getpid();
167  if((form_map.find(pid)==form_map.end())) { // init section
168  ostringstream ss;
169  ss << init_script << endl;
170  form_map[pid].Init("form");
171  form_map[pid].Execute(ss.str());
172  }
173  Form &fprc = form_map[pid];
174 
175  ex expr = expr_in.subs(SP_map);
176  ex all_expr = expr;
177  stringstream sss;
178  FormFormat ids(sss);
179  for(auto kv : SP_map) {
180  ids << "id " << kv.first << "=" << kv.second << ";" << endl;
181  all_expr += iWF(kv.first) * iWF(kv.second);
182  }
183  string idstr = sss.str();
184 
185  lst vec_lst, VD_lst, CF_lst, CA_lst, sym_lst;
186  for(const_preorder_iterator i = all_expr.preorder_begin(); i != all_expr.preorder_end(); ++i) {
187  if(is_a<Vector>(*i)) vec_lst.append(*i);
188  else if(is_a<Index>(*i)) {
189  if(ex_to<Index>(*i).type==Index::Type::VD) VD_lst.append(*i);
190  else if(ex_to<Index>(*i).type==Index::Type::CF) CF_lst.append(*i);
191  else if(ex_to<Index>(*i).type==Index::Type::CA) CA_lst.append(*i);
192  } else if(is_a<symbol>(*i)) sym_lst.append(*i);
193  else if(is_a<GiNaC::function>(*i)) {
194  static vector<string> fun_vec = {
195  "iWF", "WF", "TR", "sin", "cos", "HF", "TTR", "Matrix"
196  };
197  auto func = ex_to<GiNaC::function>(*i).get_name();
198  bool ok = false;
199  for(auto fi : fun_vec) {
200  if(fi==func) { ok = true; break; }
201  }
202  if(!ok) {
203  cout << (*i) << endl;
204  throw Error("runform: Functions not defined in FORM: "+func);
205  }
206  }
207  }
208  vec_lst.sort(); vec_lst.unique();
209  VD_lst.sort(); VD_lst.unique();
210  CF_lst.sort(); CF_lst.unique();
211  CA_lst.sort(); CA_lst.unique();
212  sym_lst.append(d);
213  sym_lst.sort(); sym_lst.unique();
214 
215  stringstream ss;
216  FormFormat ff(ss);
217  ff << ".store" << endl;
218  symtab st;
219  if(sym_lst.nops()>0) {
220  ff << "Symbols";
221  for(auto sx : sym_lst) {
222  auto s = ex_to<symbol>(sx);
223  ff << " " << s.get_name();
224  st[s.get_name()] = sx;
225  }
226  ff << ";" << endl;
227  }
228  if(vec_lst.nops()>0) {
229  ff << "Vectors";
230  for(auto vx : vec_lst) {
231  auto v = ex_to<Vector>(vx);
232  ff << " " << v;
233  st[v.name.get_name()] = v;
234  for(auto vvx : vec_lst) {
235  auto vv = ex_to<Vector>(vvx);
236  st[v.name.get_name()+"__"+vv.name.get_name()] = SP(v,vv).subs(SP_map);
237  }
238  }
239  ff << ";" << endl;
240  }
241  if(VD_lst.nops()>0) {
242  if(form_using_dim4) ff << "Dimension 4;" << endl;
243  else ff << "Dimension d;" << endl;
244  ff << "Indices";
245  for(auto ix : VD_lst) {
246  auto i = ex_to<Index>(ix);
247  ff << " " << i;
248  st[i.name.get_name()] = i;
249  auto itr = Index::Dimension.find(ix);
250  if(itr!=Index::Dimension.end()) ff << "=" << itr->second;
251  }
252  ff << ";" << endl;
253  }
254  if(CF_lst.nops()>0) {
255  ff << "Dimension NF;" << endl;
256  ff << "Indices";
257  for(auto ix : CF_lst) {
258  auto i = ex_to<Index>(ix);
259  ff << " " << i;
260  st[i.name.get_name()] = i;
261  auto itr = Index::Dimension.find(ix);
262  if(itr!=Index::Dimension.end()) ff << "=" << itr->second;
263  }
264  ff << ";" << endl;
265  }
266  if(CA_lst.nops()>0) {
267  ff << "Dimension NA;" << endl;
268  ff << "Indices";
269  for(auto ix : CA_lst) {
270  auto i = ex_to<Index>(ix);
271  ff << " " << i;
272  st[i.name.get_name()] = i;
273  auto itr = Index::Dimension.find(ix);
274  if(itr!=Index::Dimension.end()) ff << "=" << itr->second;
275  }
276  ff << ";" << endl;
277  }
278  ff << ".global" << endl;
279  ff << endl;
280 
281  // trace and contract
282  bool islst = is_a<lst>(expr);
283  lst expr_lst;
284  if(islst) expr_lst = ex_to<lst>(expr);
285  else expr_lst.append(expr);
286  auto total = expr_lst.nops();
287 
288  string ostr;
289  int gid = 1;
290  ostr = "{";
291  int c = 1;
292  int kid = 0;
293  map<ex,int,ex_is_less> e2i_map;
294  for(auto it : expr_lst) {
295  ff << ".store" << endl;
296  ex item = it;
297  item = MapFunction([](const ex & e, MapFunction &self)->ex{
298  if(!e.has(TR(w))) return e;
299  else if(e.match(TR(w))) {
300  ex nd = e.op(0).normal().numer_denom();
301  return TR(nd.op(0))/nd.op(1);
302  } else return e.map(self);
303  })(item);
304  item = normal(item); // TODO: normal
305  // pull out global common factor
306  item = collect_common_factors(item);
307  item = CoPat(item,[](const ex &e)->bool{return Index::has(e) || DGamma::has(e) || Eps::has(e);});
308  auto ckey = item.op(0);
309  if(!is_a<numeric>(ckey)) {
310  int ckid;
311  if(e2i_map.find(ckey)==e2i_map.end()) {
312  kid++;
313  e2i_map[ckey] = kid;
314  ckid = kid;
315  } else ckid = e2i_map[ckey];
316  string gCF = "gCF" + to_string(ckid);
317  st[gCF] = item.op(0);
318  item = Symbol(gCF) * item.op(1);
319  } else {
320  item = ckey * item.op(1);
321  }
322 
323  // pull out color factor
324  auto cv_lst = collect_lst(item, [](const ex &e)->bool{return Index::hasc(e);});
325  item=0;
326  exvector color_vec;
327  for(int i=0; i<cv_lst.nops(); i++) {
328  auto it = cv_lst.op(i);
329  auto cc = it.op(0);
330  auto vv = it.op(1);
331  color_vec.push_back(vv);
332  item += cc * Symbol("[cl"+to_string(i)+"]");
333  }
334 
335  ostringstream coss;
336  for(int i=0; i<color_vec.size(); i++) {
337  coss << " [cl" << i << "]";
338  ff << "L [cl" << i << "]=" << color_vec[i] << ";" << endl;
339  ff << ".sort" << endl;
340  ff << "#call SUNTrace" << endl;
341  ff << ".sort" << endl;
342  if(form_using_su3) {
343  ff << "id NF^reX?=3^reX;" << endl;
344  ff << "id I2R=1/2;" << endl;
345  ff << "id NA^reX?=8^reX;" << endl;
346  ff << ".sort" << endl;
347  }
348  }
349  ff << endl;
350 
351  // methods to handle TR objects
352  auto tr_mode = form_trace_mode;
353  exset trs;
354  find(item,TR(w),trs);
355  if(form_trace_mode==form_trace_auto) {
356  if(trs.size()<2) tr_mode = form_trace_all;
357  else tr_mode = form_trace_each_each;
358  }
359  if(tr_mode==form_trace_all) {
360  mapTR tr;
361  item = tr(item);
362  ff << "L [o]=" << item << ";" << endl;
363  ff << ".sort" << endl;
364  ff << "#do i = 1,1" << endl;
365  ff << "id once HF(reX?)=reX;" << endl;
366  ff << idstr;
367  ff << "if(count(HF,1)>0) redefine i \"0\";" << endl;
368  ff << ".sort" << endl;
369  ff << "#enddo" << endl;
370  ff << "contract 0;" << endl;
371  ff << ".sort" << endl;
372  ff << idstr << ".sort" << endl;
373  ff << endl;
374  for(int gl=1; gl<=tr.glmax; gl++) {
375  if(form_using_gamma5) {
376  ff << ".sort" << endl;
377  if(form_using_dim4) ff << "Dimension 4;" << endl;
378  else ff << "Dimension d;" << endl;
379  ff << "Indices [g5_i1], [g5_i2], [g5_i3], [g5_i4];" << endl;
380  ff << "id g_(" << gl << ",5_) = e_([g5_i1],[g5_i2],[g5_i3],[g5_i4])*g_(" << gl << ",[g5_i1],[g5_i2],[g5_i3],[g5_i4])/24;" << endl;
381  ff << "sum [g5_i1],[g5_i2],[g5_i3],[g5_i4];" << endl;
382  ff << ".sort" << endl;
383  }
384  if(form_using_dim4) ff << "trace4 " << gl << ";" << endl;
385  else ff << "tracen " << gl << ";" << endl;
386  ff << ".sort" << endl;
387  ff << idstr << ".sort" << endl;
388  }
389  } else if(tr_mode==form_trace_each_all) {
390  exmap tr2v;
391  int trn=0;
392  exvector trvec;
393  for(auto tr : trs) {
394  tr2v[tr] = Symbol("[tr"+to_string(trn)+"]");
395  auto tri = mapGamma(gid)(tr.op(0));
396  trvec.push_back(tri);
397  trn++;
398  }
399  item = item.subs(tr2v);
400  for(int i=0; i<trvec.size(); i++) {
401  ff << "L [tr" << i << "]=" << trvec[i] << ";" << endl;
402  if(form_using_gamma5) {
403  ff << ".sort" << endl;
404  if(form_using_dim4) ff << "Dimension 4;" << endl;
405  else ff << "Dimension d;" << endl;
406  ff << "Indices [g5_i1], [g5_i2], [g5_i3], [g5_i4];" << endl;
407  ff << "id g_(" << gid << ",5_) = e_([g5_i1],[g5_i2],[g5_i3],[g5_i4])*g_(" << gid << ",[g5_i1],[g5_i2],[g5_i3],[g5_i4])/24;" << endl;
408  ff << "sum [g5_i1],[g5_i2],[g5_i3],[g5_i4];" << endl;
409  }
410  if(form_using_dim4) ff << "trace4 " << gid << ";" << endl;
411  else ff << "tracen " << gid << ";" << endl;
412  ff << ".sort" << endl;
413  ff << idstr << ".sort" << endl;
414  }
415  ff << "L [o]=" << item << ";" << endl;
416  ff << ".sort" << endl;
417  ff << "#do i = 1,1" << endl;
418  ff << "id once HF(reX?)=reX;" << endl;
419  ff << idstr;
420  ff << "if(count(HF,1)>0) redefine i \"0\";" << endl;
421  ff << ".sort" << endl;
422  ff << "#enddo" << endl;
423  ff << idstr << ".sort" << endl;
424  ff << endl;
425  } else if(tr_mode==form_trace_each_each) {
426  exmap tr2v;
427  int trn=0;
428  exvector trvec;
429  for(auto tr : trs) {
430  tr2v[tr] = Symbol("trcN"+to_string(trn));
431  auto tri = mapGamma(gid)(tr.op(0));
432  trvec.push_back(tri);
433  trn++;
434  }
435  item = item.subs(tr2v);
436 
437  if(color_vec.size()>0) ff << "drop" << coss.str() << ";" << endl;
438  ff << "L [o]=" << item << ";" << endl;
439  ff << ".sort" << endl;
440  ff << "#do i = 1,1" << endl;
441  ff << "id once HF(reX?)=reX;" << endl;
442  ff << idstr;
443  ff << "if(count(HF,1)>0) redefine i \"0\";" << endl;
444  ff << ".sort" << endl;
445  ff << "#enddo" << endl;
446  ff << endl;
447 
448  for(int i=0; i<trvec.size(); i++) {
449  ff << "L [tr" << i << "]=" << trvec[i] << ";" << endl;
450  if(form_using_gamma5) {
451  ff << ".sort" << endl;
452  if(form_using_dim4) ff << "Dimension 4;" << endl;
453  else ff << "Dimension d;" << endl;
454  ff << "Indices [g5_i1], [g5_i2], [g5_i3], [g5_i4];" << endl;
455  ff << "id g_(" << gid << ",5_) = e_([g5_i1],[g5_i2],[g5_i3],[g5_i4])*g_(" << gid << ",[g5_i1],[g5_i2],[g5_i3],[g5_i4])/24;" << endl;
456  ff << "sum [g5_i1],[g5_i2],[g5_i3],[g5_i4];" << endl;
457  }
458  if(form_using_dim4) ff << "trace4 " << gid << ";" << endl;
459  else ff << "tracen " << gid << ";" << endl;
460  ff << ".sort" << endl;
461  ff << idstr << ".sort" << endl;
462  ff << "drop [tr" << i << "];" << endl;
463  ff << "id trcN" << i << "=[tr" << i << "];" << endl;
464  ff << ".sort" << endl;
465  ff << idstr << ".sort" << endl;
466  ff << endl;
467  }
468  } else {
469  throw Error("runform: unsupported form_trace_mode = " + to_string(form_trace_mode));
470  }
471 
472  ff << "contract 0;" << endl;
473  ff << ".sort" << endl;
474  ff << idstr << ".sort" << endl;
475 
476  if(verb==1) {
477  cout << "\r \r";
478  cout << PRE << "\\--Form Script @ " << c << " / " << total << flush;
479  } else if(verb>1) {
480  cout << "--------------------------------------" << endl;
481  cout << "Form Script @ " << c << " / " << total << endl;
482  cout << "--------------------------------------" << endl;
483  cout << ss.str() << endl;
484  }
485 
486  auto script = ss.str();
487  string_replace_all(script, "sin(", "sin_(");
488  string_replace_all(script, "cos(", "cos_(");
489  try {
490  auto otmp = fprc.Execute(script);
491  if(verb>2) {
492  cout << "--------------------------------------" << endl;
493  cout << "Form Output @ " << c << " / " << total << endl;
494  cout << "--------------------------------------" << endl;
495  cout << otmp << endl;
496  }
497 
498  ostr += otmp;
499  ss.clear();
500  ss.str("");
501 
502  if(c<total) ostr += ",";
503  c++;
504  } catch(Error& err) {
505  form_map.erase(pid);
506  throw;
507  }
508  }
509  if(verb==1) cout << endl;
510  ostr += "}";
511 
512  string_replace_all(ostr, "[", "(");
513  string_replace_all(ostr, "]", ")");
514  string_replace_all(ostr, "\\\n", "");
515  string_replace_all(ostr, " ","");
516  for(auto v : vec_lst) {
517  string pat(ex_to<Vector>(v).name.get_name());
518  string from = pat+"(";
519  string to = "d_("+pat+",";
520  string_replace_all(ostr, from, to);
521  from = pat+".";
522  to = pat+"__";
523  string_replace_all(ostr, from, to);
524  }
525 
526  string_replace_all(ostr, "d_(", "SP(");
527  string_replace_all(ostr, "e_(", "LC(");
528  string_replace_all(ostr, "sin_(", "sin(");
529  string_replace_all(ostr, "cos_(", "cos(");
530  string_replace_all(ostr, "g_(", "DG(");
531  string_replace_all(ostr, ",5_", ",5");
532  string_replace_all(ostr, ",6_", ",6");
533  string_replace_all(ostr, ",7_", ",7");
534  string_replace_all(ostr, "gi_(", "GI(");
535  string_replace_all(ostr, "TTR(", "TTRX(");
536 
537  st["I2R"] = TF;
538  st["NA"] = NA;
539  st["NF"] = NF;
540  st["I"] = I;
541  st["i_"] = I;
542 
543  Parser fp(st);
544  fp.FTable.insert({{"SP", 2}, reader_func(SP_reader)});
545  fp.FTable.insert({{"LC", 4}, reader_func(LC_reader)});
546  fp.FTable.insert({{"Matrix", 3}, reader_func(Matrix_reader)});
547  for(int i=1; i<30; i++) fp.FTable.insert({{"T", i}, reader_func(SUNT_reader)});
548  for(int i=1; i<30; i++) fp.FTable.insert({{"TTRX", i}, reader_func(TTR_reader)});
549  for(int i=1; i<30; i++) fp.FTable.insert({{"DG", i}, reader_func(DGamma_reader)});
550  fp.FTable.insert({{"GI", 0}, reader_func(gi_reader)});
551  fp.FTable.insert({{"GI", 1}, reader_func(gi_reader)});
552  ex ret = fp.Read(ostr);
553  if(!islst) ret = ret.op(0);
554  return ret;
555  }}
556 
563  ex form(const ex & iexpr, int verb) {
564  ex expr = iexpr;
565  if(!is_a<lst>(expr)) {
566  static MapFunction mf([](const ex & e, MapFunction & self)->ex{
567  if(!e.has(TR(w))) return e;
568  else if(e.match(TR(w))) {
569  if(is_a<mul>(e.op(0))) {
570  ex c = 1, v = 1;
571  for(auto const & ei : e.op(0)) {
572  if(DGamma::has(ei)) v *= ei;
573  else c *= ei;
574  }
575  return c*TR(v);
576  } else return e;
577  } else return e.map(self);
578  });
579  expr = mf(expr);
580  }
581  if(form_expand_mode==form_expand_none || is_a<lst>(expr)) return runform(expr, verb);
582  else if(form_expand_mode==form_expand_tr) {
583  auto cv_lst = collect_lst(expr.subs(SP_map), TR(w));
584  lst to_lst;
585  for(auto cv : cv_lst) to_lst.append(cv.op(0)*cv.op(1));
586  lst out_lst = ex_to<lst>(runform(to_lst, verb));
587 
588  ex ret = 0;
589  for(int i=0; i<cv_lst.nops(); i++) ret += out_lst.op(i);
590 
591  return ret.subs(SP_map);
592  } else if(form_expand_mode==form_expand_ci) {
593  auto cv_lst = collect_lst(expr.subs(SP_map), [](const ex & e)->bool {
594  return e.has(TR(w)) || Index::hasc(e);
595  });
596  lst to_lst;
597  for(auto cv : cv_lst) to_lst.append(cv.op(0)*cv.op(1));
598  lst out_lst = ex_to<lst>(runform(to_lst, verb));
599 
600  ex ret = 0;
601  for(int i=0; i<cv_lst.nops(); i++) ret += out_lst.op(i);
602 
603  return ret.subs(SP_map);
604  } else if(form_expand_mode==form_expand_li) {
605  auto cv_lst = collect_lst(expr.subs(SP_map), [](const ex & e)->bool {
606  return e.has(TR(w)) || Index::hasv(e);
607  });
608  lst to_lst;
609  for(auto cv : cv_lst) to_lst.append(cv.op(0)*cv.op(1));
610  lst out_lst = ex_to<lst>(runform(to_lst, verb));
611 
612  ex ret = 0;
613  for(int i=0; i<cv_lst.nops(); i++) ret += out_lst.op(i);
614 
615  return ret.subs(SP_map);
616  } else if(form_expand_mode==form_expand_all) {
617  auto cv_lst = collect_lst(expr.subs(SP_map), [](const ex & e)->bool {
618  return e.has(TR(w)) || SUNT::has(e) || SUNF::has(e) || SUNF4::has(e) || Index::has(e) || DGamma::has(e);
619  });
620  lst to_lst;
621  for(auto cv : cv_lst) to_lst.append(cv.op(1));
622  lst out_lst = ex_to<lst>(runform(to_lst, verb));
623 
624  ex ret = 0;
625  for(int i=0; i<cv_lst.nops(); i++) ret += cv_lst.op(i).op(0) * out_lst.op(i);
626 
627  return ret.subs(SP_map);
628  } else throw Error("form: unsupported form_expand_mode = " + to_string(form_expand_mode));
629  return 0;
630  }
631 
637  ex charge_conjugate(const ex & expr) {
638  if(expr.has(Matrix(w1,w2,w3))) throw Error("charge_conjugate: Matrix found.");
639  if(!DGamma::has(expr)) return expr;
640  if(is_a<DGamma>(expr)) {
641  DGamma g = ex_to<DGamma>(expr);
642  if(is_a<Vector>(g.pi) || is_a<Index>(g.pi)) return -expr;
643  else if(is_zero(g.pi-1) || is_zero(g.pi-5)) return expr;
644  }
645 
646  if(is_a<add>(expr)) {
647  ex ret = 0;
648  for(auto item : expr) ret += charge_conjugate(item);
649  return ret;
650  } else if(is_a<mul>(expr)) {
651  ex ret = 1;
652  for(auto item : expr) ret *= charge_conjugate(item);
653  return ret;
654  } else if(is_a<ncmul>(expr)) {
655  ex ret = 1;
656  int n = expr.nops();
657  for(int i=n-1; i>-1; i--) ret *= charge_conjugate(expr.op(i));
658  return ret;
659  }
660  cout << DGamma::has(expr) << " : " << expr << endl;
661  throw Error("charge_conjugate: unexpected region.");
662  return 0;
663  }
664 
665 }
666 
glmax
unsigned glmax
Definition: Form.cpp:106
HEP.h
HEP header file.
HepLib::DGamma
class for Dirac Gamma object
Definition: HEP.h:437
HepLib::DGamma::has
static bool has(const ex &e)
HepLib::Eps::has
static bool has(const ex &e)
HepLib::Error
class used to wrap error message
Definition: BASIC.h:242
HepLib::Index::hasc
static bool hasc(const ex &e)
Definition: Basic.cpp:228
HepLib::Index::Dimension
static GiNaC::exmap Dimension
Definition: HEP.h:109
HepLib::Index::has
static bool has(const ex &e)
HepLib::MapFunction
class to wrap map_function of GiNaC
Definition: BASIC.h:632
HepLib
HepLib namespace.
Definition: BASIC.cpp:17
HepLib::form_using_dim4
bool form_using_dim4
Definition: Init.cpp:209
HepLib::NA
const Symbol NA
HepLib::form_expand_tr
const int form_expand_tr
Definition: Init.cpp:202
HepLib::CoPat
lst CoPat(const ex &e, std::function< bool(const ex &)> f)
Definition: BASIC.h:725
HepLib::form_expand_none
const int form_expand_none
Definition: Init.cpp:201
HepLib::NF
const Symbol NF
HepLib::SP_map
exmap SP_map
Definition: Init.cpp:179
HepLib::form_expand_mode
int form_expand_mode
Definition: Init.cpp:206
HepLib::form_trace_all
const int form_trace_all
Definition: Init.cpp:196
HepLib::GAS
ex GAS(const ex &expr, unsigned rl)
function similar to GAD/GSD in FeynClac
Definition: DGamma.cpp:246
HepLib::form_expand_ci
const int form_expand_ci
Definition: Init.cpp:203
HepLib::form_trace_mode
int form_trace_mode
Definition: Init.cpp:199
HepLib::form_using_su3
bool form_using_su3
Definition: Init.cpp:208
HepLib::w
ex w
Definition: Init.cpp:90
HepLib::form_trace_each_all
const int form_trace_each_all
Definition: Init.cpp:197
HepLib::d
const Symbol d
HepLib::form_expand_all
const int form_expand_all
Definition: Init.cpp:205
HepLib::CA
const Symbol CA
HepLib::TF
const Symbol TF
HepLib::CF
const Symbol CF
HepLib::as
const Symbol as
HepLib::charge_conjugate
ex charge_conjugate(const ex &expr)
make the charge conjugate operaton, M -> C^{-1} . M^T . C w.r.t. a Matrix object
Definition: Form.cpp:637
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::form_using_gamma5
bool form_using_gamma5
Definition: Init.cpp:210
HepLib::string_replace_all
void string_replace_all(string &str, const string &from, const string &to)
Definition: Functions.cpp:148
HepLib::form_trace_each_each
const int form_trace_each_each
Definition: Init.cpp:198
HepLib::form_trace_auto
const int form_trace_auto
Definition: Init.cpp:195
HepLib::w1
ex w1
Definition: BASIC.h:499
HepLib::w3
ex w3
Definition: BASIC.h:499
HepLib::form_expand_li
const int form_expand_li
Definition: Init.cpp:204
HepLib::LC
ex LC(ex pi1, ex pi2, ex pi3, ex pi4)
function similar to LCD in FeynCalc
Definition: Eps.cpp:179
HepLib::w2
ex w2
Definition: BASIC.h:499
HepLib::form
ex form(const ex &iexpr, int verb)
evalulate expr in form program, see also the form_trace_mode and form_expand_mode
Definition: Form.cpp:563
HepLib::PRE
string PRE
Definition: Init.cpp:140
HepLib::SP
ex SP(const ex &a, bool use_map=false)
Definition: Pair.cpp:166
HepLib::ex2int
int ex2int(ex num)
ex to integer
Definition: BASIC.cpp:893
Generated by doxygen 1.9.1