HepLib
Loading...
Searching...
No Matches
Form.cpp
Go to the documentation of this file.
1
6#include "HEP.h"
7#include "ginac/parse_context.h"
8
9namespace 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 GMat_reader(const exvector& ev) {
59 return GMat(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(
112CFunction pow,sqrt,gamma,HF,GMat,WF;
113Tensor TTR(cyclic), f(antisymmetric), T, f4, colTp;
114Symbols reX,I2R,NF,NA,d,I,Pi;
115AutoDeclare Symbols gCF, trcN;
116Dimension NA;
117AutoDeclare Index colA;
118Dimension NF;
119AutoDeclare Index colF;
120
121#procedure SUNTrace
122Dimension NF;
123
124repeat;
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;
133endrepeat;
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
147repeat;
148 id T(colA1?,?a,colF1?,colF2?)*T(colA2?,?b,colF2?,colF3?) = T(colA1,?a,colA2,?b,colF1,colF3);
149endrepeat;
150id T(?a,colF1?,colF1?) = TTR(?a);
151id TTR(colA1?) = 0;
152id 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", "GMat"
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, "\n", "");
516 string_replace_all(ostr, " ","");
517 for(auto v : vec_lst) {
518 string pat(ex_to<Vector>(v).name.get_name());
519 string from = pat+"(";
520 string to = "d_("+pat+",";
521 string_replace_all(ostr, from, to);
522 from = pat+".";
523 to = pat+"__";
524 string_replace_all(ostr, from, to);
525 }
526
527 string_replace_all(ostr, "d_(", "SP(");
528 string_replace_all(ostr, "e_(", "LC(");
529 string_replace_all(ostr, "sin_(", "sin(");
530 string_replace_all(ostr, "cos_(", "cos(");
531 string_replace_all(ostr, "g_(", "DG(");
532 string_replace_all(ostr, ",5_", ",5");
533 string_replace_all(ostr, ",6_", ",6");
534 string_replace_all(ostr, ",7_", ",7");
535 string_replace_all(ostr, "gi_(", "GI(");
536 string_replace_all(ostr, "TTR(", "TTRX(");
537
538 st["I2R"] = TF;
539 st["NA"] = NA;
540 st["NF"] = NF;
541 st["I"] = I;
542 st["i_"] = I;
543
544 Parser fp(st);
545 fp.FTable.insert({{"SP", 2}, reader_func(SP_reader)});
546 //fp.FTable.insert({{"LC", 4}, reader_func(LC_reader)});
547 for(int i=3; i<10; i++) fp.FTable.insert({{"LC", i}, reader_func(LC_reader)});
548 fp.FTable.insert({{"GMat", 3}, reader_func(GMat_reader)});
549 for(int i=1; i<30; i++) fp.FTable.insert({{"T", i}, reader_func(SUNT_reader)});
550 for(int i=1; i<30; i++) fp.FTable.insert({{"TTRX", i}, reader_func(TTR_reader)});
551 for(int i=1; i<30; i++) fp.FTable.insert({{"DG", i}, reader_func(DGamma_reader)});
552 fp.FTable.insert({{"GI", 0}, reader_func(gi_reader)});
553 fp.FTable.insert({{"GI", 1}, reader_func(gi_reader)});
554 ex ret = fp.Read(ostr);
555 if(!islst) ret = ret.op(0);
556 return ret;
557 }}
558
565 ex form(const ex & iexpr, int verb) {
566 ex expr = iexpr;
567 if(!is_a<lst>(expr)) {
568 static MapFunction mf([](const ex & e, MapFunction & self)->ex{
569 if(!e.has(TR(w))) return e;
570 else if(e.match(TR(w))) {
571 if(is_a<mul>(e.op(0))) {
572 ex c = 1, v = 1;
573 for(auto const & ei : e.op(0)) {
574 if(DGamma::has(ei)) v *= ei;
575 else c *= ei;
576 }
577 return c*TR(v);
578 } else return e;
579 } else return e.map(self);
580 });
581 expr = mf(expr);
582 }
583 if(form_expand_mode==form_expand_none || is_a<lst>(expr)) return runform(expr, verb);
584 else if(form_expand_mode==form_expand_tr) {
585 auto cv_lst = collect_lst(expr.subs(SP_map), TR(w));
586 lst to_lst;
587 for(auto cv : cv_lst) to_lst.append(cv.op(0)*cv.op(1));
588 lst out_lst = ex_to<lst>(runform(to_lst, verb));
589
590 ex ret = 0;
591 for(int i=0; i<cv_lst.nops(); i++) ret += out_lst.op(i);
592
593 return ret.subs(SP_map);
594 } else if(form_expand_mode==form_expand_ci) {
595 auto cv_lst = collect_lst(expr.subs(SP_map), [](const ex & e)->bool {
596 return e.has(TR(w)) || Index::hasc(e);
597 });
598 lst to_lst;
599 for(auto cv : cv_lst) to_lst.append(cv.op(0)*cv.op(1));
600 lst out_lst = ex_to<lst>(runform(to_lst, verb));
601
602 ex ret = 0;
603 for(int i=0; i<cv_lst.nops(); i++) ret += out_lst.op(i);
604
605 return ret.subs(SP_map);
606 } else if(form_expand_mode==form_expand_li) {
607 auto cv_lst = collect_lst(expr.subs(SP_map), [](const ex & e)->bool {
608 return e.has(TR(w)) || Index::hasv(e);
609 });
610 lst to_lst;
611 for(auto cv : cv_lst) to_lst.append(cv.op(0)*cv.op(1));
612 lst out_lst = ex_to<lst>(runform(to_lst, verb));
613
614 ex ret = 0;
615 for(int i=0; i<cv_lst.nops(); i++) ret += out_lst.op(i);
616
617 return ret.subs(SP_map);
618 } else if(form_expand_mode==form_expand_all) {
619 auto cv_lst = collect_lst(expr.subs(SP_map), [](const ex & e)->bool {
620 return e.has(TR(w)) || SUNT::has(e) || SUNF::has(e) || SUNF4::has(e) || Index::has(e) || DGamma::has(e);
621 });
622 lst to_lst;
623 for(auto cv : cv_lst) to_lst.append(cv.op(1));
624 lst out_lst = ex_to<lst>(runform(to_lst, verb));
625
626 ex ret = 0;
627 for(int i=0; i<cv_lst.nops(); i++) ret += cv_lst.op(i).op(0) * out_lst.op(i);
628
629 return ret.subs(SP_map);
630 } else throw Error("form: unsupported form_expand_mode = " + to_string(form_expand_mode));
631 return 0;
632 }
633
639 ex charge_conjugate(const ex & expr) {
640 if(expr.has(GMat(w1,w2,w3))) throw Error("charge_conjugate: GMat found.");
641
642 if(!DGamma::has(expr) && !AsGamma::has(expr)) return expr;
643 if(is_a<DGamma>(expr)) {
644 DGamma g = ex_to<DGamma>(expr);
645 if(is_a<Vector>(g.pi) || is_a<Index>(g.pi)) return -expr;
646 else if(is_zero(g.pi-1) || is_zero(g.pi-5)) return expr;
647 } else if(is_a<AsGamma>(expr)) {
648 AsGamma ag = ex_to<AsGamma>(expr);
649 ex sign = 1;
650 lst pis_lst;
651 for(int i=ag.pis.size()-1; i>=0; --i) {
652 ex pi = ag.pis[i];
653 if(is_a<Vector>(pi) || is_a<Index>(pi)) sign *= -1;
654 else if(is_zero(pi-1) || is_zero(pi-5)) sign *= 1;
655 else throw Error("charge_conjugate: only p/i/1/5 supported.");
656 pis_lst.append(pi);
657 }
658 return sign * AsGamma::from(pis_lst);
659 }
660
661 if(is_a<add>(expr)) {
662 ex ret = 0;
663 for(auto item : expr) ret += charge_conjugate(item);
664 return ret;
665 } else if(is_a<mul>(expr)) {
666 ex ret = 1;
667 for(auto item : expr) ret *= charge_conjugate(item);
668 return ret;
669 } else if(is_a<ncmul>(expr)) {
670 ex ret = 1;
671 int n = expr.nops();
672 for(int i=n-1; i>-1; i--) ret *= charge_conjugate(expr.op(i));
673 return ret;
674 }
675 cout << expr << endl;
676 throw Error("charge_conjugate: unexpected region.");
677 return 0;
678 }
679
685 ex gamma_transpose(const ex & expr) {
686 if(expr.has(GMat(w1,w2,w3))) throw Error("charge_conjugate: GMat found.");
687 if(AsGamma::has(expr)) throw Error("gamma_transpose: Not supported for AsGamma.");
688
689 if(!DGamma::has(expr)) return expr;
690 if(is_a<DGamma>(expr)) {
691 if(expr.is_equal(DGamma::C)) return (-1)*DGamma::C;
692 DGamma g = ex_to<DGamma>(expr);
693 g.isTr = !g.isTr;
694 return g;
695 }
696
697 if(is_a<add>(expr)) {
698 ex ret = 0;
699 for(auto item : expr) ret += gamma_transpose(item);
700 return ret;
701 } else if(is_a<mul>(expr)) {
702 ex ret = 1;
703 for(auto item : expr) ret *= gamma_transpose(item);
704 return ret;
705 } else if(is_a<ncmul>(expr)) {
706 ex ret = 1;
707 int n = expr.nops();
708 for(int i=n-1; i>-1; i--) ret *= gamma_transpose(expr.op(i));
709 return ret;
710 }
711 cout << "gamma_transpose input:" << expr << endl;
712 throw Error("gamma_transpose: unexpected region.");
713 return 0;
714 }
715
716}
717
glmax
unsigned glmax
Definition Form.cpp:106
HEP.h
HEP header file.
HepLib::DGamma::has
static bool has(const ex &e)
HepLib::Eps::has
static bool has(const ex &e)
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
HepLib namespace.
Definition BASIC.cpp:17
HepLib::form_using_dim4
bool form_using_dim4
Definition Init.cpp:212
HepLib::NA
const Symbol NA
HepLib::CoPat
lst CoPat(const ex &e, std::function< bool(const ex &)> f)
Definition BASIC.h:725
HepLib::NF
const Symbol NF
HepLib::SP_map
exmap SP_map
Definition Init.cpp:182
HepLib::form_trace_all
const int form_trace_all
Definition Init.cpp:199
HepLib::GAS
ex GAS(const ex &expr, unsigned rl)
function similar to GAD/GSD in FeynClac
Definition DGamma.cpp:281
HepLib::form_trace_mode
int form_trace_mode
Definition Init.cpp:202
HepLib::form_using_su3
bool form_using_su3
Definition Init.cpp:211
HepLib::w
ex w
Definition Init.cpp:93
HepLib::form_trace_each_all
const int form_trace_each_all
Definition Init.cpp:200
HepLib::d
const Symbol d
HepLib::TF
const Symbol TF
HepLib::gamma_transpose
ex gamma_transpose(const ex &expr)
make the transpose operaton M --> M^T
Definition Form.cpp:683
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 GMat 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:213
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:201
HepLib::form_trace_auto
const int form_trace_auto
Definition Init.cpp:198
HepLib::LC
ex LC(ex pi1, ex pi2, ex pi3, ex pi4)
function similar to LCD in FeynCalc
Definition Eps.cpp:179
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:143
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.8